Preconditioner.C 4.64 KB
Newer Older
1 2
////////////////////////////////////////////////////////////////////////////////
//
Francois Gygi committed
3 4 5 6
// Copyright (c) 2008 The Regents of the University of California
//
// This file is part of Qbox
//
Francois Gygi committed
7 8
// Qbox is distributed under the terms of the GNU General Public License
// as published by the Free Software Foundation, either version 2 of
Francois Gygi committed
9 10 11 12 13 14
// the License, or (at your option) any later version.
// See the file COPYING in the root directory of this distribution
// or <http://www.gnu.org/licenses/>.
//
////////////////////////////////////////////////////////////////////////////////
//
15 16 17 18 19 20
// Preconditioner.C
//
////////////////////////////////////////////////////////////////////////////////

#include "Preconditioner.h"
#include "Basis.h"
Francois Gygi committed
21
#include "Wavefunction.h"
22 23
#include "EnergyFunctional.h"
#include "ConfinementPotential.h"
24
#include "SlaterDet.h"
25
using namespace std;
26 27

////////////////////////////////////////////////////////////////////////////////
28 29
Preconditioner::Preconditioner(const Wavefunction& wf, EnergyFunctional& ef,
  double ecutprec) : ef_(ef), ecutprec_(ecutprec)
30
{
Francois Gygi committed
31 32 33 34 35 36 37 38 39 40 41 42 43 44 45
  kpg2_.resize(wf.nspin());
  ekin_.resize(wf.nspin());
  for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
  {
    kpg2_[ispin].resize(wf.nkp());
    ekin_[ispin].resize(wf.nkp());
    for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
    {
      const SlaterDet& sd = *(wf.sd(ispin,ikp));
      const Basis& wfbasis = sd.basis();
      kpg2_[ispin][ikp] = wfbasis.kpg2_ptr();
      ekin_[ispin][ikp].resize(sd.nstloc());
    }
  }
  update(wf);
46 47 48
}

////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
49
double Preconditioner::diag(int ispin, int ikp, int n, int ig) const
50
{
51
  const valarray<double>& fstress = ef_.confpot(ikp)->fstress();
Francois Gygi committed
52 53
  if ( ecutprec_ == 0.0 )
  {
54 55 56
    double ekin_n = ekin_[ispin][ikp][n];
    // if ekin_n == 0 (occurs for first wf, G=0, when starting without
    // randomizing wfs) replace ekin_n by fixed value 1.0
Francois Gygi committed
57
    if ( ekin_n == 0.0 )
58
      ekin_n = 1.0;
59
    const double q2 = kpg2_[ispin][ikp][ig] + fstress[ig];
60
#if 0
Francois Gygi committed
61 62 63 64 65 66 67 68 69 70 71 72
#if 0
    // Payne Teter Allan adaptive preconditioner
    const double x = 0.5 * q2 / ekin_n;
    const double num = 27.0 + x * ( 18.0 + x * ( 12 + x * 8 ));
    const double den = num + 16.0 * x*x*x*x;
    return (num/den);
#else
    // modified Payne Teter Allan adaptive preconditioner (Kresse)
    const double x = 0.5 * q2 / ( 1.5 * ekin_n );
    const double num = 27.0 + x * ( 18.0 + x * ( 12 + x * 8 ));
    const double den = num + 16.0 * x*x*x*x;
    return ( 2.0 / (1.5 *  ekin_n )) * (num/den);
73 74 75 76 77
#endif
#else
    // basic adaptive preconditioner: use ekin_n for the value of ecutprec
    double e = 0.5 * ( kpg2_[ispin][ikp][ig] + fstress[ig] );
    return ( e < ekin_n ) ? 0.5 / ekin_n : 0.5 / e;
Francois Gygi committed
78 79 80 81
#endif
  }
  else
  {
82
    double e = 0.5 * ( kpg2_[ispin][ikp][ig] + fstress[ig] );
Francois Gygi committed
83 84 85
    return ( e < ecutprec_ ) ? 0.5 / ecutprec_ : 0.5 / e;
  }
}
86

Francois Gygi committed
87 88 89 90
////////////////////////////////////////////////////////////////////////////////
void Preconditioner::update(const Wavefunction& wf)
{
  // update the kinetic energy ekin_[ispin][ikp][n] of states in wf
91 92 93 94
  for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
  {
    for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
    {
Francois Gygi committed
95 96
      const SlaterDet& sd = *(wf.sd(ispin,ikp));
      const Basis& wfbasis = sd.basis();
97 98
      // factor fac in next lines: 2.0 for G and -G (if basis is real)
      const double fac = wfbasis.real() ? 2.0 : 1.0;
Francois Gygi committed
99 100 101

      const ComplexMatrix& c = sd.c();
      const Context& sdctxt = sd.context();
102

Francois Gygi committed
103 104 105 106 107 108 109 110
      const int ngwloc = wfbasis.localsize();
      const complex<double>* p = c.cvalptr();
      const int mloc = c.mloc();
      const int nloc = c.nloc();

      valarray<double> buf(2*nloc);
      const double* pkpg2 = kpg2_[ispin][ikp];
      for ( int n = 0; n < nloc; n++ )
Francois Gygi committed
111
      {
Francois Gygi committed
112 113
        double sum_norm = 0.0;
        double sum_ekin = 0.0;
Francois Gygi committed
114
        for ( int ig = 0; ig < ngwloc; ig++ )
115
        {
Francois Gygi committed
116 117 118
          const double psi2 = norm(p[ig+n*mloc]);
          sum_norm += psi2;
          sum_ekin += psi2 * pkpg2[ig];
119
        }
Francois Gygi committed
120 121 122 123 124 125
        // correct for double counting of G=0 in norm
        if ( sdctxt.myrow() == 0 )
          sum_norm -= 0.5 * norm(p[n*mloc]);
        // store norm in buf[n] and ekin in buf[n+nloc]
        buf[n] = fac * sum_norm;
        buf[n+nloc] = fac * sum_ekin;
Francois Gygi committed
126
      }
Francois Gygi committed
127
      sdctxt.dsum('C',2*nloc,1,&buf[0],2*nloc);
128
      // factor 0.5 in next line: 1/(2m)
Francois Gygi committed
129
      for ( int n = 0; n < nloc; n++ )
130
        ekin_[ispin][ikp][n] = 0.5*buf[n] > 0.0 ? 0.5*buf[n+nloc]/buf[n] : 0;
Francois Gygi committed
131 132 133

#ifdef DEBUG
      if ( sdctxt.onpe0() )
Francois Gygi committed
134
      {
Francois Gygi committed
135 136 137
        for ( int n = 0; n < nloc; n++ )
          cout << "Preconditioner::update ekin[" << n << "] = "
               << ekin_[ispin][ikp][n] << endl;
138
      }
Francois Gygi committed
139
#endif
140 141 142
    }
  }
}