1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (c) 2008 The Regents of the University of California
4 //
5 // This file is part of Qbox
6 //
7 // Qbox is distributed under the terms of the GNU General Public License
8 // as published by the Free Software Foundation, either version 2 of
9 // the License, or (at your option) any later version.
10 // See the file COPYING in the root directory of this distribution
11 // or <http://www.gnu.org/licenses/>.
12 //
13 ////////////////////////////////////////////////////////////////////////////////
14 //
15 // MDWavefunctionStepper.cpp
16 //
17 ////////////////////////////////////////////////////////////////////////////////
18 
19 #include "MDWavefunctionStepper.h"
20 #include "Wavefunction.h"
21 #include "SlaterDet.h"
22 #include "Sample.h"
23 #include <iostream>
24 using namespace std;
25 
26 ////////////////////////////////////////////////////////////////////////////////
MDWavefunctionStepper(Wavefunction & wf,Wavefunction * wfv,double dt,double dt2bye,TimerMap & tmap)27 MDWavefunctionStepper::MDWavefunctionStepper(Wavefunction& wf,
28   Wavefunction *wfv, double dt, double dt2bye, TimerMap& tmap) :
29   wfv_(wfv), dt_(dt), dt2bye_(dt2bye), WavefunctionStepper(wf,tmap)
30 {
31   assert(wfv!=0);
32 
33   tmap_["md_update_wf"].reset();
34   tmap_["riccati"].reset();
35   tmap_["riccati"].reset();
36   tmap_["ekin_e"].reset();
37 }
38 
39 ////////////////////////////////////////////////////////////////////////////////
update(Wavefunction & dwf)40 void MDWavefunctionStepper::update(Wavefunction& dwf)
41 {
42   // Verlet update of wf using force dwf and wfm stored in *wfv
43   for ( int isp_loc = 0; isp_loc < wf_.nsp_loc(); ++isp_loc )
44   {
45     for ( int ikp_loc = 0; ikp_loc < wf_.nkp_loc(); ++ikp_loc )
46     {
47       tmap_["md_update_wf"].start();
48       // Verlet update of wf
49       // cp = c + (c - cm) - dt2/m * hpsi
50       // This is implemented (for each coefficient) as:
51       // cp = 2*c - cm - dt2bye * hpsi
52       // cm = c
53       // c = cp
54       SlaterDet* sd = wf_.sd(isp_loc,ikp_loc);
55       SlaterDet* sdv = wfv_->sd(isp_loc,ikp_loc);
56       double* cptr = (double*) sd->c().valptr();
57       double* cptrm = (double*) sdv->c().valptr();
58       const double* dcptr =
59         (const double*) dwf.sd(isp_loc,ikp_loc)->c().cvalptr();
60       const int mloc = sd->c().mloc();
61       const int nloc = sd->c().nloc();
62       for ( int n = 0; n < nloc; n++ )
63       {
64         // note: double mloc length for complex<double> indices
65         double* c = &cptr[2*mloc*n];
66         double* cm = &cptrm[2*mloc*n];
67         const double* dc = &dcptr[2*mloc*n];
68         for ( int i = 0; i < mloc; i++ )
69         {
70           const double ctmp = c[2*i];
71           const double ctmp1 = c[2*i+1];
72           const double cmtmp = cm[2*i];
73           const double cmtmp1 = cm[2*i+1];
74           const double dctmp = dc[2*i];
75           const double dctmp1 = dc[2*i+1];
76           const double cptmp = 2.0*ctmp -  cmtmp -  dt2bye_ * dctmp;
77           const double cptmp1= 2.0*ctmp1 - cmtmp1 - dt2bye_ * dctmp1;
78 
79           c[2*i]    = cptmp;
80           c[2*i+1]  = cptmp1;
81           cm[2*i]   = ctmp;
82           cm[2*i+1] = ctmp1;
83         }
84       }
85       tmap_["md_update_wf"].stop();
86 
87       tmap_["riccati"].start();
88       sd->riccati(*sdv);
89       tmap_["riccati"].stop();
90     }
91   }
92   ekin_em_ = ekin_ep_;
93   ekin_ep_ = ekin_eh();
94 }
95 
96 ////////////////////////////////////////////////////////////////////////////////
compute_wfm(Wavefunction & dwf)97 void MDWavefunctionStepper::compute_wfm(Wavefunction& dwf)
98 {
99   // Compute wfm for first MD step using wf, wfv and dwf (= Hpsi)
100   // Replace then wfv by wfm
101 
102   // Compute cm using c and wavefunction velocity
103   // cm = c - dt * v - 0.5 * dt2/m * hpsi
104   // replace wfv by wfm
105   const double half_dt2bye = 0.5 * dt2bye_;
106   for ( int isp_loc = 0; isp_loc < wf_.nsp_loc(); ++isp_loc )
107   {
108     for ( int ikp_loc = 0; ikp_loc < wf_.nkp_loc(); ++ikp_loc )
109     {
110       SlaterDet* sd = wf_.sd(isp_loc,ikp_loc);
111       SlaterDet* sdv = wfv_->sd(isp_loc,ikp_loc);
112 
113       double* cptr = (double*) sd->c().valptr();
114       double* cptrv = (double*) sdv->c().valptr();
115       const double* dcptr =
116         (const double*) dwf.sd(isp_loc,ikp_loc)->c().cvalptr();
117       const int mloc = sd->c().mloc();
118       const int nloc = sd->c().nloc();
119       for ( int n = 0; n < nloc; n++ )
120       {
121         // note: double mloc length for complex<double> indices
122         double* c = &cptr[2*mloc*n];
123         double* cv = &cptrv[2*mloc*n];
124         const double* dc = &dcptr[2*mloc*n];
125         for ( int i = 0; i < mloc; i++ )
126         {
127           const double ctmp = c[2*i];
128           const double ctmp1 = c[2*i+1];
129           const double cvtmp = cv[2*i];
130           const double cvtmp1 = cv[2*i+1];
131           const double dctmp = dc[2*i];
132           const double dctmp1 = dc[2*i+1];
133           cv[2*i]    = ctmp  - dt_ * cvtmp  - half_dt2bye * dctmp;
134           cv[2*i+1]  = ctmp1 - dt_ * cvtmp1 - half_dt2bye * dctmp1;
135         }
136       }
137       tmap_["riccati"].start();
138       sdv->riccati(*sd);
139       tmap_["riccati"].stop();
140     }
141   }
142   ekin_em_ = 0.0;
143   ekin_ep_ = ekin_eh();
144   // Note: ekin_ep is a first-order approximation of ekin_e using wf and wfm
145   // only. This makes ekin_e consistent with the following update() call
146   // Note: *wfv_ now contains wf(t-dt)
147 }
148 
149 ////////////////////////////////////////////////////////////////////////////////
compute_wfv(Wavefunction & dwf)150 void MDWavefunctionStepper::compute_wfv(Wavefunction& dwf)
151 {
152   // Compute wfv = (wf - wfm)/dt - 0.5*dtbye*dwf
153 
154   assert(dt_!=0.0);
155   for ( int isp_loc = 0; isp_loc < wf_.nsp_loc(); ++isp_loc )
156   {
157     for ( int ikp_loc = 0; ikp_loc < wf_.nkp_loc(); ++ikp_loc )
158     {
159       // compute final velocity wfv
160       // v = ( c - cm ) / dt - 0.5 * dt/m * hpsi
161 
162       // Note: At this point, *wfv_ contains wf(t-dt)
163 
164       // hpsi must be orthogonal to the subspace spanned by c
165       // compute descent direction H psi - psi (psi^T H psi)
166 
167       SlaterDet* sd = wf_.sd(isp_loc,ikp_loc);
168       if ( sd->basis().real() )
169       {
170         // proxy real matrices c, cp
171         DoubleMatrix c(sd->c());
172         DoubleMatrix cp(dwf.sd(isp_loc,ikp_loc)->c());
173 
174         DoubleMatrix a(c.context(),c.n(),c.n(),c.nb(),c.nb());
175 
176         // factor 2.0 in next line: G and -G
177         a.gemm('t','n',2.0,c,cp,0.0);
178         // rank-1 update correction
179         a.ger(-1.0,c,0,cp,0);
180 
181         // cp = cp - c * a
182         cp.gemm('n','n',-1.0,c,a,1.0);
183       }
184       else
185       {
186         ComplexMatrix& c = sd->c();
187         ComplexMatrix& cp = dwf.sd(isp_loc,ikp_loc)->c();
188         ComplexMatrix a(c.context(),c.n(),c.n(),c.nb(),c.nb());
189         a.gemm('c','n',1.0,c,cp,0.0);
190         // cp = cp - c * a
191         cp.gemm('n','n',-1.0,c,a,1.0);
192       }
193 
194       const double dt_inv = 1.0/dt_;
195       const double half_dtbye = 0.5 * dt2bye_ / dt_;
196       double* cptr = (double*) sd->c().valptr();
197       double* cptrv = (double*) wfv_->sd(isp_loc,ikp_loc)->c().valptr();
198       const double* dcptr =
199         (const double*) dwf.sd(isp_loc,ikp_loc)->c().cvalptr();
200       const int mloc = sd->c().mloc();
201       const int nloc = sd->c().nloc();
202       for ( int n = 0; n < nloc; n++ )
203       {
204         // note: double mloc length for complex<double> indices
205         double* c = &cptr[2*mloc*n];
206         double* cv = &cptrv[2*mloc*n];
207         const double* dc = &dcptr[2*mloc*n];
208         for ( int i = 0; i < mloc; i++ )
209         {
210           const double ctmp = c[2*i];
211           const double ctmp1 = c[2*i+1];
212           const double cmtmp = cv[2*i];
213           const double cmtmp1 = cv[2*i+1];
214           const double dctmp = dc[2*i];
215           const double dctmp1 = dc[2*i+1];
216 
217           cv[2*i]   = ( ctmp  - cmtmp  ) * dt_inv
218                       - half_dtbye * dctmp;
219           cv[2*i+1] = ( ctmp1 - cmtmp1 ) * dt_inv
220                       - half_dtbye * dctmp1;
221         }
222       }
223       // Note: *wfv_ now contains the wavefunction velocity
224     }
225   }
226 }
227 
228 ////////////////////////////////////////////////////////////////////////////////
ekin_eh(void)229 double MDWavefunctionStepper::ekin_eh(void)
230 {
231   // compute ekin at time t - 0.5*dt using wf and wfm
232   tmap_["ekin_e"].start();
233   double ekin_e = 0.0;
234   // assume that wfv contains wfm
235   for ( int isp_loc = 0; isp_loc < wf_.nsp_loc(); ++isp_loc )
236   {
237     for ( int ikp_loc = 0; ikp_loc < wf_.nkp_loc(); ++ikp_loc )
238     {
239       const int ikpg = wf_.ikp_global(ikp_loc);
240       const double weight = wf_.weight(ikpg);
241       SlaterDet* sd = wf_.sd(isp_loc,ikp_loc);
242       double* cptr = (double*) sd->c().valptr();
243       double* cptrm = (double*) wfv_->sd(isp_loc,ikp_loc)->c().valptr();
244       const int mloc = sd->c().mloc();
245       const int nloc = sd->c().nloc();
246       // compute electronic kinetic energy at time t-1/2
247       const bool onrow0 = ( wf_.sd_context().myrow() == 0 );
248       const vector<double>& occ = sd->occ();
249       for ( int n = 0; n < nloc; n++ )
250       {
251         const int nglobal = sd->c().j(0,n);
252         const double occn = occ[nglobal];
253         // note: double mloc length for complex<double> indices
254         double* c = &cptr[2*mloc*n];
255         double* cm = &cptrm[2*mloc*n];
256         double tmpsum = 0.0;
257         if ( sd->basis().real() && onrow0 )
258         {
259           // correct for double counting of G=0 element
260           // i=0 coefficient is real, count only real part
261           const double ctmp = c[0];
262           const double cmtmp = cm[0];
263           tmpsum -= 0.5 * (ctmp - cmtmp)*(ctmp - cmtmp);
264         }
265         for ( int i = 0; i < mloc; i++ )
266         {
267           const double ctmp = c[2*i];
268           const double ctmp1 = c[2*i+1];
269           const double cmtmp = cm[2*i];
270           const double cmtmp1 = cm[2*i+1];
271 
272           tmpsum += (ctmp -cmtmp )*(ctmp -cmtmp ) +
273                     (ctmp1-cmtmp1)*(ctmp1-cmtmp1);
274         }
275         if ( sd->basis().real() )
276           // Note: 2 in next line: from (G,-G)
277           ekin_e += weight * ( 2.0 * occn / dt2bye_ ) * tmpsum;
278         else
279           ekin_e += weight * ( occn / dt2bye_ ) * tmpsum;
280       }
281     }
282   }
283   double tsum;
284   MPI_Allreduce(&ekin_e,&tsum,1,MPI_DOUBLE,MPI_SUM,MPIdata::comm());
285   ekin_e = tsum;
286   tmap_["ekin_e"].stop();
287   return ekin_e;
288 }
289