1 # include "wpmd.h"
2
3
4
5 // Calculates derivative overlap matrix IDD
calc_der_overlap(bool self,cdouble cc1,cdouble c2)6 void OverlapDeriv::calc_der_overlap(bool self, cdouble cc1, cdouble c2){
7 cVector_3 I3 = I1 * ((bb_4a + 2.5) / w12.a);
8 cdouble I4 = I0 * ( bb_4a *(bb_4a + 5.) + 3.75 ) / w12.a / w12.a;
9
10 // calculate derivatives <(phi_k)'_q_k | (phi_l)'_q_l>:
11 IDD.set(0, 0, I4 - (d1.l + d2.l)*I2 + d1.l*d2.l*I0 ); // over a_k_re and a_l_re
12 IDD.set(0, 1, i_unit*( I4 - (d1.l + d2.m)*I2 + d1.l*d2.m*I0 ) ); // over a_k_re and a_l_im
13 if(!self)
14 IDD.set(1, 0, i_unit1*( I4 + (d1.m - d2.l)*I2 - d1.m*d2.l*I0 ) ); // over a_k_im and a_l_re
15 else
16 IDD.set(1,0, conj(IDD(0,1)));
17 IDD.set(1, 1, I4 + (d1.m - d2.m)*I2 - d1.m*d2.m*I0 ); // over a_k_im and a_l_im
18
19 for(int i=0;i<3;i++){
20 IDD.set(0, (i+1)*2, -I3[i] + d1.l*I1[i] + d2.u[i]*(d1.l*I0 - I2) ); // over a_k_re and b_l_re
21 IDD.set(0, (i+1)*2+1, i_unit1*( I3[i] - d1.l*I1[i] + d2.v[i]*(I2 - d1.l*I0) ) ); // over a_k_re and b_l_im
22 IDD.set(1, (i+1)*2, i_unit *( I3[i] + d1.m*I1[i] + d2.u[i]*(I2 + d1.m*I0) ) ); // over a_k_im and b_l_re
23 IDD.set(1, (i+1)*2+1, -I3[i] - d1.m*I1[i] - d2.v[i]*(d1.m*I0 + I2) ); // over a_k_im and b_l_im
24 if(!self) {
25 IDD.set((i+1)*2, 0, -I3[i] + d2.l*I1[i] + d1.u[i]*(d2.l*I0 - I2) ); // over b_k_re and a_l_re
26 IDD.set((i+1)*2+1, 0, i_unit *( I3[i] - d2.l*I1[i] - d1.v[i]*(I2 - d2.l*I0) ) ); // over b_k_im and a_l_re
27 IDD.set((i+1)*2, 1, i_unit1*( I3[i] - d2.m*I1[i] + d1.u[i]*(I2 - d2.m*I0) ) ); // over b_k_re and a_l_im
28 IDD.set((i+1)*2+1, 1, -I3[i] + d2.m*I1[i] - d1.v[i]*(d2.m*I0 - I2) ); // over b_k_im and a_l_im
29 }
30 else{
31 IDD.set((i+1)*2, 0, conj(IDD(0,(i+1)*2)) ); // over b_k_re and a_l_re
32 IDD.set((i+1)*2+1, 0, conj(IDD(0,(i+1)*2+1)) ); // over b_k_im and a_l_re
33 IDD.set((i+1)*2, 1, conj(IDD(1,(i+1)*2)) ); // over b_k_re and a_l_im
34 IDD.set((i+1)*2+1, 1, conj(IDD(1,(i+1)*2+1)) ); // over b_k_im and a_l_im
35 }
36
37 for(int j=0;j<3;j++){
38 if(!self || j>=i){
39 cdouble I2ij = I0 / w12.a *
40 (i==j ? w12.b[i]*w12.b[i] / w12.a / 4 + 0.5
41 : w12.b[i]*w12.b[j] / w12.a / 4);
42 // over b_k_re and b_l_re
43 IDD.set((j+1)*2, (i+1)*2, I2ij + d1.u[i]*I1[j] + d2.u[j]*(I1[i] + d1.u[i]*I0) );
44 // over b_k_re and b_l_im
45 IDD.set((j+1)*2, (i+1)*2+1, i_unit *( I2ij + d1.u[i]*I1[j] + d2.v[j]*(I1[i] + d1.u[i]*I0) ) );
46 // over b_k_im and b_l_re
47 if(!self || i!=j)
48 IDD.set((j+1)*2+1, (i+1)*2, i_unit1*( I2ij - d1.v[i]*I1[j] + d2.u[j]*(I1[i] - d1.v[i]*I0) ) );
49 else
50 IDD.set((j+1)*2+1, (i+1)*2, conj(IDD((i+1)*2,(j+1)*2+1)));
51 // over b_k_im and b_l_im
52 IDD.set((j+1)*2+1,(i+1)*2+1, I2ij - d1.v[i]*I1[j] + d2.v[j]*(I1[i] - d1.v[i]*I0) );
53 }
54 else{ // self && j<i
55 // over b_k_re and b_l_re
56 IDD.set((j+1)*2, (i+1)*2, conj(IDD((i+1)*2, (j+1)*2)) );
57 // over b_k_re and b_l_im
58 IDD.set((j+1)*2, (i+1)*2+1, conj(IDD((i+1)*2+1,(j+1)*2)) );
59 // over b_k_im and b_l_re
60 IDD.set((j+1)*2+1, (i+1)*2, conj(IDD((i+1)*2,(j+1)*2+1)) );
61 // over b_k_im and b_l_im
62 IDD.set((j+1)*2+1,(i+1)*2+1, conj(IDD((i+1)*2+1,(j+1)*2+1 )) );
63 }
64 } // j
65 } // i
66
67 if(real(cc1)){ // adding terms for split-packet
68
69 IDD.set(8, 0, c2*da2_re() ); // over c_1_re and a_2_re
70 IDD.set(8, 1, c2*da2_im() ); // over c_1_re and a_2_im
71 IDD.set(9, 0, -i_unit*c2*da2_re() ); // over c_1_im and a_2_re
72 IDD.set(9, 1, -i_unit*c2*da2_im() ); // over c_1_im and a_2_im
73
74 IDD.set(0, 8, cc1*da1_re() ); // over c_2_re and a_1_re
75 IDD.set(1, 8, cc1*da1_im() ); // over c_2_re and a_1_im
76 IDD.set(0, 9, i_unit*cc1*da1_re() ); // over c_2_im and a_1_re
77 IDD.set(1, 9, i_unit*cc1*da1_im() ); // over c_2_im and a_1_im
78
79 for(int i=0;i<3;i++){
80 IDD.set(8, 2+2*i, c2*db2_re(i) ); // over c_1_re and b_2_re
81 IDD.set(8, 2+2*i+1, c2*db2_im(i) ); // over c_1_re and b_2_im
82 IDD.set(9, 2+2*i, -i_unit*c2*db2_re(i) ); // over c_1_im and b_2_re
83 IDD.set(9, 2+2*i+1, -i_unit*c2*db2_im(i) ); // over c_1_im and b_2_im
84
85 IDD.set(2+2*i, 8, cc1*db1_re(i) ); // over c_2_re and b_1_re
86 IDD.set(2+2*i+1, 8, cc1*db1_im(i) ); // over c_2_re and b_1_im
87 IDD.set(2+2*i, 9, i_unit*cc1*db1_re(i) ); // over c_2_im and i_1_re
88 IDD.set(2+2*i+1, 9, i_unit*cc1*db1_im(i) ); // over c_2_im and a_1_im
89 }
90
91 IDD.set(8, 8, I0 ); // over c_1_re and c_2_re
92 IDD.set(8, 9, i_unit*I0 ); // over c_1_re and c_2_im
93 IDD.set(9, 8, -i_unit*I0 ); // over c_1_im and c_2_re
94 IDD.set(9, 9, I0 ); // over c_1_im and c_2_im
95 }
96 }
97
98
99
create_wp(Vector_3 & x,Vector_3 & v,double & w,double & pw,double mass)100 WavePacket AWPMD::create_wp(Vector_3 &x, Vector_3 &v, double &w, double &pw, double mass){
101 if(mass<0)
102 mass=me;
103 if(constraint==FIX){
104 if(w0>0)
105 w=w0;
106 pw=0.;
107 }
108
109 double rw;
110 if(Lextra>0){ // width PBC, keeping the width are within [0,Lextra]
111 w=fmod(w,Lextra);
112 if(w<0) w+=Lextra;
113 rw=w; // WP width for energy evaluation is within [0, L/2]
114 if(rw > Lextra/2) rw = Lextra - rw;
115 }
116 else
117 rw=w;
118
119 WavePacket wp;
120 wp.init(rw,x,v*mass*one_h,pw*one_h);
121 return wp;
122 }
123
124
resize(int flag)125 void AWPMD::resize(int flag){
126 for(int s=0;s<2;s++){
127 //0. resizing matrices
128 Y[s].init(ne[s],1);
129 O[s].init(ne[s],1);
130 Oflg[s].init(ne[s],1);
131 //Te[s].init(nel,1);
132 //Tei[s].init(nel,1);
133 Eep[s].assign((size_t)nwp[s],0);
134 Eeip[s].assign((size_t)nwp[s],0);
135 Eeep[s].assign((size_t)nwp[s],0);
136 Ewp[s].assign((size_t)nwp[s],0);
137
138 if(flag&(0x8|0x4) && approx!=HARTREE){ //electron forces, L and M are needed
139 M[s].init(ne[s],nvar[s]);
140 L[s].init(ne[s],nvar[s]);
141 }
142 }
143 Eiep.assign((size_t)ni,0);
144 Eiip.assign((size_t)ni,0);
145
146
147 }
148
149
150 //e sets Periodic Boundary Conditions
151 //e using bit flags: 0x1 -- PBC along X
152 //e 0x2 -- PBC along Y
153 //e 0x4 -- PBC along Z
154 //e cell specifies the lengths of the simulation box in all directions
155 //e if PBCs are used, the corresponding coordinates of electrons and ions
156 //e in periodic directions must be within a range [0, cell[per_dir])
157 //e @returns 1 if OK
set_pbc(const Vector_3P pcell,int pbc_)158 int AWPMD::set_pbc(const Vector_3P pcell, int pbc_){
159 if(!pcell)
160 pbc=0;
161 else{
162 pbc=pbc_;
163 cell=*pcell;
164 }
165 return 1;
166 }
167
168 //e setup elctrons: forms internal wave packet representations
169 //e if PBCs are used the coords must be within a range [0, cell)
set_electrons(int s,int n,Vector_3P x,Vector_3P v,double * w,double * pw,double mass,double * q)170 int AWPMD::set_electrons(int s, int n, Vector_3P x, Vector_3P v, double* w, double* pw, double mass, double *q)
171 {
172 if(s < 0 || s > 1)
173 return LOGERR(-1,logfmt("AWPMD.set_electrons: invaid s setting (%d)!",s),LINFO);
174
175 norm_matrix_state[s] = NORM_UNDEFINED;
176 nwp[s]=ne[s]=n;
177 nvar[s]=8*n;
178 wp[s].resize(n);
179
180 partition1[s].clear();
181 for(int i=0;i<n;i++){
182 wp[s][i]=create_wp(x[i],v[i],w[i],pw[i], mass);
183 // assign default partition
184 partition1[s].push_back(i+1);
185 }
186
187 // assign electronic charge
188 if(q)
189 qe[s].assign(q,q+nwp[s]);
190 else
191 qe[s].assign(nwp[s],-1);
192
193
194 return 1;
195 }
196
197 //e setup ion charges and coordinates
198 //e if PBCs are used the coords must be within a range [0, cell)
set_ions(int n,double * q,Vector_3P x)199 int AWPMD::set_ions(int n, double* q, Vector_3P x)
200 {
201 ni = n;
202 qi.resize(n);
203 xi.resize(n);
204 partition1[2].clear();
205 for(int i=0;i<n;i++){
206 qi[i] = q[i], xi[i] = x[i];
207 // assign default partition for ions
208 partition1[2].push_back(i+1);
209 }
210
211 return 1;
212 }
213
214 //e same as interaction, but using Hartee factorization (no antisymmetrization)
interaction_hartree(int flag,Vector_3P fi,Vector_3P fe_x,Vector_3P fe_p,double * fe_w,double * fe_pw,Vector_2P fe_c)215 int AWPMD::interaction_hartree(int flag, Vector_3P fi, Vector_3P fe_x,
216 Vector_3P fe_p, double *fe_w, double *fe_pw, Vector_2P fe_c){
217
218 // 0. resizing the arrays if needed
219 enum APPROX tmp=HARTREE;
220 swap(tmp,approx); // do not need large matrices
221 resize(flag);
222 swap(tmp,approx);
223 //1. clearing forces
224 clear_forces(flag,fi,fe_x,fe_p,fe_w,fe_pw,fe_c);
225
226 Eee = Ew = 0.;
227 for(int s1=0;s1<2;s1++){
228 Ee[s1]=0.;
229 Eei[s1]=0.;
230 for(int c1=0;c1<ne[s1];c1++){
231 // width part
232 double w1=wp[s1][c1].get_width();
233 /*double sgn1=1;
234 if(Lextra>0){ // width PBC
235 if(w1>Lextra-w1){
236 w1=-(Lextra-w1); // '-' is to change derivative sign
237 sgn1=-1;
238 }
239 }*/
240 Vector_3 r1=wp[s1][c1].get_r();
241 Vector_3 p=wp[s1][c1].get_p()*h_plank;
242 Vector_3 pw=wp[s1][c1].get_pwidth()*h_plank;
243 // energy contribution
244 Ee[s1] += (p.norm2()+pw.norm2())/(2*me);
245 Ew += h2_me*9./(8.*w1*w1);
246 if(constraint == HARM) Ew += harm_w0_4 * w1*w1;
247 // width force contribution
248 //double dE=2*Epot/w;
249 //if(d->fw1)d->fw1[c1]+=dE;
250 //if(fw2 && fw2!=fw1)fw2[c1]+=dE;
251
252 // e-e interaction
253 for(int s2=s1;s2<2;s2++){
254 for(int c2=(s1==s2 ? c1+1 : 0) ;c2<ne[s2];c2++){
255 double w2=wp[s2][c2].get_width();
256 Vector_3 v12=wp[s2][c2].get_r()-r1;
257 // position PBC
258 v12=v12.rcell1(cell,pbc);
259 double r12=v12.normalize();
260 /*double sgn2=1; // signs
261 if(Lextra>0){ // width PBC
262 if(w2>Lextra-w2){
263 w2=-(Lextra-w2); // '-' is to change derivative sign
264 sgn2=-1;
265 }
266 }*/
267 double wsq=w1*w1+w2*w2;
268 double argw=sqrt((2./3.)*wsq);
269
270 //double arg=r12/argw;
271 //double erfa=erf(arg);
272 double Epot=coul_pref*erf_div(r12,1./argw); //erfa/r12;
273 Eee+=Epot;
274
275 // force contribution
276 /*double dEw=coul_pref*two_over_sqr_pi*exp(-arg*arg)/argw;
277 double dEpot=(Epot-dEw)/r12;
278 if(!d->fixw){
279 dEw/=wsq;
280 if(d->fw1 && c1>=0){
281 d->fw1[c1]+=sgn1*dEw*w1;
282 }
283 if(d->fw2){
284 d->fw2[c2]+=sgn2*dEw*w2;
285 }
286 }*/
287 }
288 }
289 // e-i interaction
290 double wsq1=w1*w1;
291 double argw=sqr_2_over_3*w1;
292 for(int i=0;i<ni;i++){
293 Vector_3 v12=xi[i]-r1;
294 // position PBC
295 v12=v12.rcell1(cell,pbc);
296 double r12=v12.normalize();
297
298 //double arg=r12/argw;
299 //double erfa=erf(arg);
300 double cel=-coul_pref*qi[i]; // electron charge is always -1
301 double Epot=cel*erf_div(r12,1./argw); //erfa/r12;
302 Eei[s1]+=Epot;
303 //printf("(%g %g %g)- (%g %g %g)\n",r1[0],r1[1],r1[2],xi[i][0],xi[i][1],xi[i][2]);
304 //printf("awp(%d,%d:%d)[%g]: %g\n",i,s1,c1,r12,Epot);
305 // force contribution
306 if(flag&0x3){
307 double arg=r12/argw;
308 double dEw=cel*two_over_sqr_pi*exp(-arg*arg)/argw;
309 double dEpot=(Epot-dEw)/r12;
310 fi[i]+=v12*dEpot; // ionic force
311 }
312 // electron force
313 /*if(!d->fixw){
314 dEw/=wsq;
315 if(d->fw1 && c1>=0){
316 d->fw1[c1]+=sgn1*dEw*w1;
317 }
318 }*/
319 }
320 }
321 }
322 if(calc_ii)
323 interaction_ii(flag,fi);
324 return 1;
325 }
326
327
328 //e initializes internal buffers for calculations (set_electrons must be called first)
329 //int init(){}
330
331 //e calculates interaction in the system of ni ions + electrons
332 //e the electonic subsystem must be previously setup by set_electrons, ionic by set_ions
333 //e the iterators are describing ionic system only
334 // 0x1 -- give back ion forces
335 // 0x2 -- add ion forces to the existing set
336 // 0x4 -- calculate derivatives for electronic time step (NOT IMPLEMENTED)
337 //e if PBCs are used the coords must be within a range [0, cell)
interaction(int flag,Vector_3P fi,Vector_3P fe_x,Vector_3P fe_p,double * fe_w,double * fe_pw,Vector_2P fe_c)338 int AWPMD::interaction(int flag, Vector_3P fi, Vector_3P fe_x,
339 Vector_3P fe_p, double *fe_w, double *fe_pw, Vector_2P fe_c){
340 if(approx==HARTREE)
341 return interaction_hartree(flag,fi,fe_x,fe_p,fe_w,fe_pw,fe_c);
342 // 0. resizing the arrays if needed
343 resize(flag);
344 // 1. clearing forces
345 clear_forces(flag,fi,fe_x,fe_p,fe_w,fe_pw,fe_c);
346
347 //2. calculating overlap matrix
348 for(int s=0;s<2;s++){
349 int nes = ne[s];
350 if(nes == 0) continue;
351
352 for(int k=0;k<nes;k++){
353 Y[s].set(k,k,1.); // Diagonal elements (=1)
354 Oflg[s](k,k) = 1;
355 for(int l=k+1;l<nes;l++){
356 cdouble I0kl = pbc_mul(wp[s][l],wp[s][k]).integral(); // Non-diagonal elements
357 Y[s].set(k,l,I0kl);
358 Oflg[s](k,l) = (norm(I0kl) > ovl_tolerance);
359 }
360 }
361 O[s] = Y[s]; // save overlap matrix
362
363 //3. inverting the overlap matrix
364 int info=0;
365 if(nes){
366 /*FILE *f1=fopen(logfmt("matrO_%d.d",s),"wt");
367 fileout(f1,Y[s],"%15g");
368 fclose(f1);8*/
369
370 ZPPTRF("L",&nes,Y[s].arr,&info);
371 // analyze return code here
372 if(info<0)
373 return LOGERR(info,logfmt("AWPMD.interacton: call to ZPTRF failed (exitcode %d)!",info),LINFO);
374 ZPPTRI("L",&nes,Y[s].arr,&info);
375 if(info<0)
376 return LOGERR(info,logfmt("AWPMD.interacton: call to ZPTRI failed (exitcode %d)!",info),LINFO);
377
378
379 /*f1=fopen(logfmt("matrY_%d.d",s),"wt");
380 fileout(f1,Y[s],"%15g");
381 fclose(f1);*/
382 }
383
384 // Clearing matrices for electronic forces
385 if(flag&0x4){
386 Te[s].Set(0.);
387 Tei[s].Set(0.);
388 }
389 }
390
391 Vector_3 ndr;
392 // calculating single particle contribution
393 for(int s=0;s<2;s++){
394 Ee[s]=Eei[s]=0.;
395 for(int k=0;k<ne[s];k++){
396 for(int l=k;l<ne[s];l++){
397
398 if( !Oflg[s](k,l) ) continue; // non-overlapping WPs
399
400 // electrons kinetic energy
401 WavePacket wk=wp[s][k];
402 WavePacket& wl=wp[s][l];
403 if(pbc)
404 ndr=move_to_image(wl,wk);
405
406 WavePacket wkl=wl*conj(wk);
407 //Vector_3 rrkl=wkl.get_r();
408 cVector_3 v1=wl.b*conj(wk.a)-conj(wk.b)*wl.a;
409 cdouble v=(v1*v1)/wkl.a;
410 v-=6.*conj(wk.a)*wl.a;
411 v/=wkl.a;
412 cdouble I0kl = O[s](k,l);
413 cdouble dE=-h2_me*I0kl*v/2;
414 if(flag&0x4) // matrix needed only for electronic forces
415 Te[s].set(k,l,dE);
416 // energy component (trace)
417 dE*=Y[s](l,k);
418 Ee[s]+=(l==k ? 1. : 2.)*real(dE);
419
420 cVector_3 dkl=wkl.b/(2.*wkl.a);
421
422 // e-i energy
423 cdouble sum(0.,0.);
424 for(int i=0;i<ni;i++){ // ions loop
425 cVector_3 gkli=dkl - cVector_3(xi[i]);
426
427 if(pbc) // correcting the real part (distance) according to PBC
428 gkli=rcell1(gkli,cell,pbc);
429 //-Igor- gkli=cVector_3(real(gkli).rcell1(cell,pbc),imag(gkli));
430
431 cdouble ngkli=gkli.norm();
432 cdouble c=sqrt(wkl.a);
433 //cdouble ttt = cerf_div(ngkli,c);
434 dE=-coul_pref*(qi[i])*I0kl*cerf_div(ngkli,c);
435
436 sum+=dE;
437 if(flag&0x3){// calculate forces on ions
438 if(fabs(real(ngkli))+fabs(imag(ngkli))>1e-10){
439 cdouble arg=ngkli*c;
440 cdouble dEw=-coul_pref*qi[i]*I0kl*two_over_sqr_pi*exp(-arg*arg)*c;
441 dE=(dE-dEw)/ngkli;
442 dE*=Y[s](l,k);
443 Vector_3 dir=-real(gkli);
444 dir.normalize();
445 fi[i]+=(l==k ? 1. : 2.)*real(dE)*dir;
446 }
447 }
448 }
449 dE=sum;
450 if(flag&0x4) // matrix needed only for electronic forces
451 Tei[s].set(k,l,dE);
452 // energy component (trace)
453 dE*=Y[s](l,k);
454 Eei[s]+=(l==k ? 1. : 2.)*real(dE);
455 }
456 }
457 }
458
459 // calculating e-e interaction
460 Eee = Ew = 0.;
461 // same spin
462 for(int s=0;s<2;s++){ // spin
463 for(int k=0;k<ne[s];k++){ //c1
464 for(int l=k+1;l<ne[s];l++){ //c3
465 for(int m=k;m<ne[s];m++){ //c2
466
467 if( Oflg[s](k,m) ) {
468 WavePacket wkm=pbc_mul(wp[s][m],wp[s][k]);
469 cVector_3 dkm=wkm.b/(2*wkm.a);
470 cdouble I0km=O[s](k,m);
471
472 // kl-mn
473 for(int n=l;n<ne[s];n++){ //c4
474 if(n<=m || !Oflg[s](l,n)) continue;
475
476 WavePacket wln=pbc_mul(wp[s][n],wp[s][l]);
477 if(pbc) // reducing the pair to elementary cell
478 ndr=move_to_image(wkm,wln); // mind the derivative: wln.b+=wln.a*ndr, wln.lz+=-wln.a*ndr^2-i*wln.old_p*ndr;
479 //Vector_3 rln=wln.get_r();
480
481 cVector_3 dln=wln.b/(2*wln.a);
482 cdouble dd=(dkm-dln).norm();
483 cdouble c=1./sqrt(1./wln.a+1./wkm.a);
484 cdouble Vklmn=coul_pref*I0km*O[s](l,n)*cerf_div(dd,c);
485
486 //cdouble arge=dkm*dkm*wkm.a+dln*dln*wln.a+wkm.lz+wln.lz;
487 //cdouble Vklmn=0.5*coul_pref*M_PI*M_PI*M_PI*exp(arge)*cerf_div(dd,c)/pow(wln.a*wkm.a,3./2.);
488 cdouble dE=Vklmn*(Y[s](m,k)*Y[s](n,l)-Y[s](m,l)*Y[s](n,k));
489 double rdE=real(dE);
490 if(m!=k || n!=l) // not the same pair
491 rdE*=2;
492 Eee+=rdE;
493 }//n
494 }
495
496 if( Oflg[s](l,m) ) {
497 WavePacket wlm=pbc_mul(wp[s][m],wp[s][l]);
498 cVector_3 dlm=wlm.b/(2*wlm.a);
499 cdouble I0lm=O[s](l,m);
500
501 // kl-nm
502 for(int n=l;n<ne[s];n++){
503 if(n<=m || !Oflg[s](k,n)) continue;
504
505 WavePacket wkn=pbc_mul(wp[s][n],wp[s][k]);
506 if(pbc) // reducing the pair to elementary cell
507 ndr=move_to_image(wlm,wkn); // mind the derivative: wln.b+=wln.a*ndr, wln.lz+=-wln.a*ndr^2-i*wln.old_p*ndr;
508
509 cVector_3 dkn=wkn.b/(2*wkn.a);
510 cdouble dd=(dkn-dlm).norm();
511 cdouble c=1./sqrt(1./wkn.a+1./wlm.a);
512 cdouble Vklnm=coul_pref*I0lm*O[s](k,n)*cerf_div(dd,c);
513
514 cdouble dE=Vklnm*(Y[s](n,k)*Y[s](m,l)-Y[s](n,l)*Y[s](m,k));
515 double rdE=real(dE);
516 if(m!=k || n!=l) // not the same pair
517 rdE*=2;
518 Eee+=rdE;
519 }//n
520 }
521 }// m
522 }// l
523 }// k
524 }// s
525
526 // different spin
527 for(int k=0;k<ne[0];k++){ // skm=0 //c1
528 for(int l=0;l<ne[1];l++){ // sln=1 //c3
529 for(int m=k;m<ne[0];m++){ //c2
530 if( Oflg[0](k,m) ) {
531 WavePacket wkm=pbc_mul(wp[0][m],wp[0][k]);
532 cVector_3 dkm=wkm.b/(2*wkm.a);
533 cdouble I0km=O[0](k,m);
534
535 for(int n=l;n<ne[1];n++){ // km-ln //c4
536 if( Oflg[1](n,l) ) {
537 WavePacket wln=pbc_mul(wp[1][l],wp[1][n]);
538 if(pbc) // reducing the pair to elementary cell
539 ndr=move_to_image(wkm,wln); // mind the derivative: wln.b+=wln.a*ndr, wln.lz+=-wln.a*ndr^2-i*wln.old_p*ndr;
540 //Vector_3 rln=wln.get_r();
541
542 cVector_3 dln=wln.b/(2*wln.a);
543 cdouble dd=(dkm-dln).norm();
544 cdouble c=1./sqrt(1./wln.a+1./wkm.a);
545 cdouble Vklmn=coul_pref*I0km*wln.integral()*cerf_div(dd,c);
546
547 cdouble dE=Vklmn*Y[0](m,k)*Y[1](n,l);
548 int Mkm=(m==k ? 1: 2);
549 int Mln=(n==l ? 1: 2);
550 double rdE=Mkm*Mln*real(dE); //!!!
551 Eee+=rdE;
552 } //if
553 } // n
554 } //if
555 } // m
556 }// l
557 }// k
558 if(calc_ii)
559 interaction_ii(flag,fi);
560 return 1;
561 }
562
563 //e Calculates Norm matrix and performs LU-factorization
564 //e The result is saved in AWPMD::Norm[s]
norm_matrix(int s)565 void AWPMD::norm_matrix(int s){
566 // Internal variables
567 int k, l, i, j, qi, qj;
568 int nes = ne[s], nes8 = nes*8, nnes8 = nes*nes8;
569
570 if(!nes) return;
571
572 // References for frequently used arrays
573 sqmatrix<double>& Norms = Norm[s];
574 chmatrix& Ys = Y[s];
575 smatrix<unsigned char>& Oflgs = Oflg[s];
576
577 // Allocate of vectors and matrices
578 Norms.init(nes8,1);
579 IDD.init(nes8,1);
580 if(ID.size() != nnes8)
581 ID.resize(nnes8), IDYs.resize(nnes8), ipiv.resize(nes8);
582
583 // Calculate first and second derivatives
584 for(k=0;k<nes;k++){
585 int k8 = k*8;
586 WavePacket& wk = wp[s][k];
587 NormDeriv dk(wk);
588 dk = conj(dk); // conjugate: mu -> -mu, v -> -v !!!
589
590 for(l=0;l<nes;l++){
591 if( !Oflgs(k,l) ) continue; // non-overlapping WPs
592
593 int l8 = l*8;
594 WavePacket wl = wp[s][l];
595 if(pbc) move_to_image(wk,wl);
596 WavePacket wkl=conj(wk)*wl;
597 NormDeriv dl(wl);
598
599 cdouble I0 = O[s](k,l);
600 cVector_3 I1 = wkl.b * (I0 / wkl.a / 2);
601 cdouble bb_4a = wkl.b.norm2() / wkl.a / 4;
602 cdouble I2 = I0 * (bb_4a + 1.5) / wkl.a;
603
604 // calculate derivatives <phi_k | (phi_l)'_q_l>:
605 int idx = k + l*nes8;
606 if(k != l) {
607 ID[idx] = dl.l*I0 - I2; // over a_l_re
608 ID[idx+nes] = i_unit*(dl.m*I0 - I2); // over a_l_im
609 for(i=0;i<3;i++){
610 ID[idx+((i+1)*2)*nes] = dl.u[i]*I0 + I1[i]; // over b_l_re
611 ID[idx+((i+1)*2+1)*nes] = i_unit*(dl.v[i]*I0 + I1[i]); // over b_l_im
612 }
613 } else { // k == l
614 ID[idx] = i_unit*imag(dl.l); // over a_l_re
615 ID[idx+nes] = i_unit*(dl.m - I2); // over a_l_im
616 for(i=0;i<3;i++){
617 ID[idx+((i+1)*2)*nes] = dl.u[i] + I1[i]; // over b_l_re
618 ID[idx+((i+1)*2+1)*nes] = 0.; // over b_l_im
619 }
620 }
621
622 if(k <= l) {
623 cVector_3 I3 = I1 * ((bb_4a + 2.5) / wkl.a);
624 cdouble I4 = I0 * ( bb_4a *(bb_4a + 5.) + 3.75 ) / wkl.a / wkl.a;
625
626 // calculate derivatives <(phi_k)'_q_k | (phi_l)'_q_l>:
627 IDD.set(k8, l8, I4 - (dk.l + dl.l)*I2 + dk.l*dl.l*I0 ); // over a_k_re and a_l_re
628 IDD.set(k8, l8+1, i_unit*( I4 - (dk.l + dl.m)*I2 + dk.l*dl.m*I0 ) ); // over a_k_re and a_l_im
629 if(k != l) IDD.set(k8+1, l8, i_unit1*( I4 + (dk.m - dl.l)*I2 - dk.m*dl.l*I0 ) ); // over a_k_im and a_l_re
630 IDD.set(k8+1, l8+1, I4 + (dk.m - dl.m)*I2 - dk.m*dl.m*I0 ); // over a_k_im and a_l_im
631
632 for(i=0;i<3;i++){
633 IDD.set(k8, l8+(i+1)*2, -I3[i] + dk.l*I1[i] + dl.u[i]*(dk.l*I0 - I2) ); // over a_k_re and b_l_re
634 IDD.set(k8, l8+(i+1)*2+1, i_unit1*( I3[i] - dk.l*I1[i] + dl.v[i]*(I2 - dk.l*I0) ) ); // over a_k_re and b_l_im
635 IDD.set(k8+1, l8+(i+1)*2, i_unit *( I3[i] + dk.m*I1[i] + dl.u[i]*(I2 + dk.m*I0) ) ); // over a_k_im and b_l_re
636 IDD.set(k8+1, l8+(i+1)*2+1, -I3[i] - dk.m*I1[i] - dl.v[i]*(dk.m*I0 + I2) ); // over a_k_im and b_l_im
637 if(k != l) {
638 IDD.set(k8+(i+1)*2, l8, -I3[i] + dl.l*I1[i] + dk.u[i]*(dl.l*I0 - I2) ); // over b_k_re and a_l_re
639 IDD.set(k8+(i+1)*2+1, l8, i_unit *( I3[i] - dl.l*I1[i] - dk.v[i]*(I2 - dl.l*I0) ) ); // over b_k_im and a_l_re
640 IDD.set(k8+(i+1)*2, l8+1, i_unit1*( I3[i] - dl.m*I1[i] + dk.u[i]*(I2 - dl.m*I0) ) ); // over b_k_re and a_l_im
641 IDD.set(k8+(i+1)*2+1, l8+1, -I3[i] + dl.m*I1[i] - dk.v[i]*(dl.m*I0 - I2) ); // over b_k_im and a_l_im
642 }
643
644 for(j=0;j<3;j++){
645 cdouble I2ij = I0 / wkl.a *
646 (i==j ? wkl.b[i]*wkl.b[i] / wkl.a / 4 + 0.5
647 : wkl.b[i]*wkl.b[j] / wkl.a / 4);
648 // over b_k_re and b_l_re
649 IDD.set(k8+(j+1)*2, l8+(i+1)*2, I2ij + dk.u[i]*I1[j] + dl.u[j]*(I1[i] + dk.u[i]*I0) );
650 // over b_k_re and b_l_im
651 IDD.set(k8+(j+1)*2, l8+(i+1)*2+1, i_unit *( I2ij + dk.u[i]*I1[j] + dl.v[j]*(I1[i] + dk.u[i]*I0) ) );
652 // over b_k_im and b_l_re
653 if(k != l) IDD.set(k8+(j+1)*2+1, l8+(i+1)*2, i_unit1*( I2ij - dk.v[i]*I1[j] + dl.u[j]*(I1[i] - dk.v[i]*I0) ) );
654 // over b_k_im and b_l_im
655 IDD.set(k8+(j+1)*2+1, l8+(i+1)*2+1, I2ij - dk.v[i]*I1[j] + dl.v[j]*(I1[i] - dk.v[i]*I0) );
656 } // j
657 } // i
658 } // if(k <= l)
659 } // k
660 } // l
661
662 // Calculate matrix product IDYs_(k,q_j) = Ys_(k,l) * <phi_l | (phi_j)'_q_j>
663 for(qj=0; qj<nes8; qj++){
664 j = qj / 8;
665 int idx = qj*nes;
666 for(k=0;k<nes;k++) {
667 cdouble sum = 0.;
668 for(l=0;l<nes;l++)
669 if( Oflgs(l,j) ) sum += ID[idx+l] * Ys(k,l);
670 IDYs[idx+k] = sum;
671 }
672 }
673
674 // Calculate Norm-matrix
675 for(qi=0; qi<nes8; qi++){
676 i = qi / 8;
677 int idxqi = qi*nes;
678
679 Norms(qi,qi) = 0.; // zero diagonal elements
680
681 for(qj=qi+1; qj<nes8; qj++){
682 j = qj / 8;
683 int idxqj = qj*nes;
684
685 // Calculate matrix product sum = <(phi_i)'_q_i | phi_k> * IDYs_(k,q_j)
686 cdouble sum = 0.;
687 for(k=0;k<nes;k++)
688 if( Oflgs(i,k) )
689 sum += IDYs[idxqj+k] * conj(ID[idxqi+k]);
690
691 // Update norm-matrix taking into account its anti-symmetry
692 double a = Oflgs(i,j) ? // IDD = 0 for non-overlapping WPs
693 h_plank2 * imag( (sum - IDD(qi,qj))*Ys(j,i) ) :
694 h_plank2 * imag( sum*Ys(j,i) );
695 Norms(qi,qj) = a;
696 Norms(qj,qi) = -a;
697 } // qj
698 } // qi
699
700 # if 1
701 // transform norm matrix to the physical variables
702 for(i=0;i<nes;i++){
703 WavePacket wi=wp[s][i];
704 for(k=0;k<8;k++){
705 // iterator to list all N(8*i+k,*) with fixed 8*i+k
706 sqmatrix<double>::iterator mi=Norms.fix_first(8*i+k,0);
707 for(j=i ;j<nes;j++){ // TO DO: run this loop from i+1 and take simplectic form for (i,i) block
708 WavePacket wj=wp[s][j];
709 wj.int2phys_der< eq_second >(mi+8*j,mi+8*j,mi+8*j+3,mi+8*j+6,mi+8*j+7);
710 }
711 }// finished line of blocks from right
712 for(k= 8*i;k<nes8;k++){ // TO DO: run this loop from 8*i+8 and take simplectic form for (i,i) block
713 // iterator to list all N(8i+*,k) by fixed k
714 sqmatrix<double>::iterator mi=Norms.fix_second(8*i,k);
715 wi.int2phys_der< eq_second >(mi,mi,mi+3,mi+6,mi+7);
716 }// finished line of blocks from left
717
718 for(k=0;k<8;k++){ // filling the lower triangle according to antisymmetry
719 for(j=8*i+8;j<nes8;j++)
720 Norms(j,8*i+k)=-Norms(8*i+k,j);
721 }
722 }
723 # endif
724 # if 0
725
726 // transform norm matrix to the physical variables
727 for(i=0;i<nes;i++){
728 WavePacket wi=wp[s][i];
729 for(j=i;j<nes;j++){
730 WavePacket wj=wp[s][j];
731 for(k=0;k<8;k++){
732 // iterator to list all N(8*i+k,*) with fixed 8*i+k
733 sqmatrix<double>::iterator mi=Norms.fix_first(8*i+k,8*j);
734 wj.int2phys_der< eq_second >(mi,mi,mi+3,mi+6,mi+7);
735 }
736 for(k=0;k<8;k++){
737 // iterator to list all N(8*i+k,*) with fixed 8*i+k
738 sqmatrix<double>::iterator mi=Norms.fix_second(8*i,8*j+k);
739 wi.int2phys_der< eq_second >(mi,mi,mi+3,mi+6,mi+7);
740 }
741 if(i!=j){
742 for(int k1=0;k1<8;k1++){
743 for(int k2=0;k2<8;k2++)
744 Norms(8*j+k1,8*i+k2)=-Norms(8*i+k2,8*j+k1);
745 }
746 }
747 }
748 }
749 # endif
750
751 norm_matrix_state[s] = NORM_CALCULATED;
752 }
753
754 //e Norm matrix LU-factorization
norm_factorize(int s)755 void AWPMD::norm_factorize(int s) {
756 if( norm_matrix_state[s] != NORM_CALCULATED) norm_matrix(s);
757
758 int nes8 = ne[s]*8, info;
759 DGETRF(&nes8, &nes8, Norm[s].arr, &nes8, &ipiv[0], &info);
760 if(info < 0)
761 LOGERR(info,logfmt("AWPMD.norm_factorize: call to DGETRF failed (exitcode %d)!",info),LINFO);
762
763 norm_matrix_state[s] = NORM_FACTORIZED;
764 }
765
766
767 //e Norm matrix inversion
norm_invert(int s)768 void AWPMD::norm_invert(int s) {
769 if( norm_matrix_state[s] != NORM_FACTORIZED) norm_factorize(s);
770
771 int nes8 = ne[s]*8, info;
772 int IDD_size = (int)IDD.get_datasize(nes8);
773
774 DGETRI(&nes8, Norm[s].arr, &nes8, &ipiv[0], (double*)IDD.arr, &IDD_size, &info); // use IDD for work storage
775 if(info < 0)
776 LOGERR(info,logfmt("AWPMD.norm_invert: call to DGETRI failed (exitcode %d)!",info),LINFO);
777
778 norm_matrix_state[s] = NORM_INVERTED;
779 }
780
781
782 //e Get the determinant of the norm-matrix for the particles with spin s
norm_matrix_det(int s)783 double AWPMD::norm_matrix_det(int s) {
784 double det = 1.;
785 int nes8 = ne[s]*8;
786
787 if(!nes8) return det;
788 if(norm_matrix_state[s] != NORM_FACTORIZED) norm_factorize(s);
789
790 sqmatrix<double>& Norms = Norm[s];
791 for(int i=0; i<nes8; i++)
792 det *= Norms(i, i); // product of the diagonal elements
793
794 return det;
795 }
796
797
798 //e Get the determinant logarithm of the norm-matrix for the particles with spin s
norm_matrix_detl(int s)799 double AWPMD::norm_matrix_detl(int s) {
800 double detl = 0.;
801 int nes8 = ne[s]*8;
802
803 if(!nes8) return detl;
804 if(norm_matrix_state[s] != NORM_FACTORIZED) norm_factorize(s);
805
806 sqmatrix<double>& Norms = Norm[s];
807 for(int i=0; i<nes8; i++)
808 detl += log(fabs( Norms(i, i) )); // product of the diagonal elements
809
810 return detl;
811 }
812
813
get_energy()814 double AWPMD::get_energy(){
815 double res=Eee + Ew;
816 for(int s=0;s<2;s++)
817 res+=Eei[s]+Ee[s];
818 if(calc_ii)
819 res+=Eii;
820 return res;
821 }
822
823 //e makes timestep of electronic component (NOT IMPLEMENTED)
step(double dt)824 int AWPMD::step(double dt){
825 return -1;
826 }
827
828
829 //e gets current electronic coordinates
get_electrons(int spin,Vector_3P x,Vector_3P v,double * w,double * pw,double mass)830 int AWPMD::get_electrons(int spin, Vector_3P x, Vector_3P v, double* w, double* pw, double mass){
831 if(spin<0 || spin >1)
832 return -1; // invalid spin: return LOGERR(-1,logfmt("AWPMD.get_electrons: invaid spin setting (%d)!",spin),LINFO);
833 if(mass<0)
834 mass=me;
835 for(int i=0;i<ni;i++){
836 w[i]=sqrt(3./(4*real(wp[spin][i].a)));
837 pw[i]=-2*w[i]*imag(wp[spin][i].a)/one_h; //pw[i]=-h_plank2*w[i]*imag(wp[spin][i].a);
838 x[i]=real(wp[spin][i].b)/(2*real(wp[spin][i].a));
839 v[i]=(pw[i]*x[i]/w[i] + imag(wp[spin][i].b)/one_h)/mass; //v[i]=(pw[i]*x[i]/w[i] + h_plank*imag(wp[spin][i].b))/m_electron;
840 }
841 return 1;
842 }
843
844
845
clear_forces(int flag,Vector_3P fi,Vector_3P fe_x,Vector_3P fe_p,double * fe_w,double * fe_pw,Vector_2P fe_c)846 void AWPMD::clear_forces(int flag,Vector_3P fi, Vector_3P fe_x,
847 Vector_3P fe_p, double *fe_w, double *fe_pw, Vector_2P fe_c){
848 if(flag&0x1){
849 for(int i=0;i<ni;i++)
850 fi[i]=Vector_3(0.);
851 }
852 if(flag&0x4 && !(flag&0x10)){ // electron forces requested in physical representation
853 for(int s1=0;s1<2;s1++){ // clearing forces
854 for(int c1=0;c1<ne[s1];c1++){
855 fe_x[c1]=Vector_3(0,0,0);
856 fe_p[c1]=Vector_3(0,0,0);
857 fe_w[c1]=0;
858 fe_pw[c1]=0;
859 }
860 }
861 }
862 }
863
864
865
interaction_ii(int flag,Vector_3P fi)866 int AWPMD::interaction_ii(int flag,Vector_3P fi){
867 Eii=0.;
868 for(int i=0;i<ni;i++){
869 for(int j=i+1;j<ni;j++){
870 double M12e, M12f;
871 _mytie(M12e,M12f)=check_part1ii(i,j);
872 if(M12f){
873 Vector_3 rij=xi[i]-xi[j];
874 double r=rij.norm();
875 double dE=coul_pref*qi[i]*qi[j]/r;
876 Eii+=M12e*dE;
877
878 Eiip[i]+=0.5*M12e*dE;
879 Eiip[j]+=0.5*M12e*dE;
880
881 if(flag&0x3){ // ion forces needed
882 Vector_3 df=-M12f*dE*rij/(r*r);
883 fi[i]+=df;
884 fi[j]-=df;
885 }
886 }
887 }
888 }
889 return 1;
890 }
891