1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (c) 2016 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 // ResponseCmd.cpp
16 //
17 ////////////////////////////////////////////////////////////////////////////////
18 
19 #include "MPIdata.h"
20 #include "BOSampleStepper.h"
21 #include "ExternalPotential.h"
22 #include "FourierTransform.h"
23 #include "ResponseCmd.h"
24 #include "release.h"
25 #include "isodate.h"
26 #include "Species.h"
27 #include "Function3d.h"
28 #include<cassert>
29 #include<vector>
30 #include<string>
31 #include<fstream>
32 using namespace std;
33 
34 ////////////////////////////////////////////////////////////////////////////////
action(int argc,char ** argv)35 int ResponseCmd::action(int argc, char **argv)
36 {
37   if ( s->wf.nst() == 0 )
38   {
39     if ( ui->onpe0() )
40       cout << "  ResponseCmd: no states, cannot run" << endl;
41     return 1;
42   }
43   if ( s->wf.ecut() == 0.0 )
44   {
45     if ( ui->onpe0() )
46       cout << "  ResponseCmd: ecut = 0.0, cannot run" << endl;
47     return 1;
48   }
49   if ( s->wf.nel() != 2 * s->wf.nst() )
50   {
51     if ( ui->onpe0() )
52       cout << "  ResponseCmd: cannot run with fractionally\n"
53               "  occupied or empty states" << endl;
54     return 1;
55   }
56 
57   bool rpa = false;
58   bool ipa = false;
59   double amplitude = 0.0;
60 
61   int iarg = 1;
62   if ( !strcmp(argv[iarg],"-vext") )
63   {
64     // response to vext
65     if ( s->vext )
66     {
67       if ( ui->onpe0() )
68         cout << "  ResponseCmd: cannot run when vext is already set" << endl;
69       return 1;
70     }
71     iarg++;
72     if ( iarg >= argc )
73     {
74       if ( ui->onpe0() )
75         cout << help_msg();
76       return 1;
77     }
78     string filename = argv[iarg];
79     iarg++;
80     string fmt = "xml";
81     if ( iarg >= argc )
82     {
83       if ( ui->onpe0() )
84         cout << help_msg();
85       return 1;
86     }
87     if ( !strcmp(argv[iarg],"-cube") )
88     {
89       fmt = "cube";
90       iarg++;
91     }
92     if ( iarg >= argc )
93     {
94       if ( ui->onpe0() )
95         cout << help_msg();
96       return 1;
97     }
98     if ( !strcmp(argv[iarg],"-RPA") )
99     {
100       rpa = true;
101       iarg++;
102     }
103     if ( iarg >= argc )
104     {
105       if ( ui->onpe0() )
106         cout << help_msg();
107       return 1;
108     }
109     if ( !strcmp(argv[iarg],"-IPA") )
110     {
111       ipa = true;
112       iarg++;
113     }
114     if ( rpa && ipa )
115     {
116       if ( ui->onpe0() )
117         cout << " Only one of -RPA or -IPA can be specified" << endl;
118       return 1;
119     }
120     if ( iarg >= argc )
121     {
122       if ( ui->onpe0() )
123         cout << help_msg();
124       return 1;
125     }
126     if ( !strcmp(argv[iarg],"-amplitude") )
127     {
128       iarg++;
129       if ( iarg >= argc )
130       {
131         if ( ui->onpe0() )
132           cout << help_msg();
133         return 1;
134       }
135       amplitude = atof(argv[iarg]);
136       iarg++;
137     }
138 
139     int nitscf = 0;
140     if ( iarg >= argc )
141     {
142       if ( ui->onpe0() )
143         cout << help_msg();
144       return 1;
145     }
146     nitscf = atoi(argv[iarg]);
147     iarg++;
148     int nite = 0;
149     if ( iarg < argc )
150       nite = atoi(argv[iarg]);
151 
152     if ( nitscf == 0 )
153     {
154       if ( ui->onpe0() )
155         cout << "  ResponseCmd: nitscf = 0, cannot run" << endl;
156       return 1;
157     }
158 
159     s->vext = new ExternalPotential(*s, filename, fmt);
160     s->vext->set_amplitude(amplitude);
161     responseVext(rpa, ipa, nitscf, nite, fmt);
162     delete s->vext;
163     s->vext = 0;
164   }
165   else
166   {
167     // polarizability calculation in constant field
168     // response amplitude [-RPA|-IPA] nitscf [nite]
169     if ( argc < 3 || argc > 5 )
170     {
171       if ( ui->onpe0() )
172         cout << help_msg();
173       return 1;
174     }
175     amplitude = atof(argv[iarg]);
176     iarg++;
177     if ( !strcmp(argv[iarg],"-RPA") )
178     {
179       rpa = true;
180       iarg++;
181     }
182     if ( iarg >= argc )
183     {
184       if ( ui->onpe0() )
185         cout << help_msg();
186       return 1;
187     }
188     if ( !strcmp(argv[iarg],"-IPA") )
189     {
190       ipa = true;
191       iarg++;
192     }
193     if ( rpa && ipa )
194     {
195       if ( ui->onpe0() )
196         cout << " Only one of -RPA or -IPA can be specified" << endl;
197       return 1;
198     }
199     if ( iarg >= argc )
200     {
201       if ( ui->onpe0() )
202         cout << help_msg();
203       return 1;
204     }
205     int nitscf = atoi(argv[iarg]);
206     iarg++;
207     int nite = 0;
208     if ( iarg < argc )
209       nite = atoi(argv[iarg]);
210 
211     if ( ui->onpe0() )
212       cout << " ResponseCmd: polarizability with "
213            << " amplitude " << argv[1] << endl;
214     responseEfield(amplitude,rpa,ipa,nitscf,nite);
215   }
216 
217   return 0;
218 }
219 
220 ////////////////////////////////////////////////////////////////////////////////
responseEfield(double amplitude,bool rpa,bool ipa,int nitscf,int nite)221 void ResponseCmd::responseEfield(double amplitude, bool rpa, bool ipa,
222   int nitscf, int nite)
223 {
224   // compute dipole change
225   BOSampleStepper* stepper = new BOSampleStepper(*s,nitscf,nite);
226   assert(stepper!=0);
227 
228   if ( rpa )
229     stepper->set_update_vxc(false);
230 
231   if ( ipa )
232   {
233     stepper->set_update_vh(false);
234     stepper->set_update_vxc(false);
235   }
236 
237   ElectricEnthalpy* el_enth = stepper->ef().el_enth();
238   if ( el_enth == 0 )
239   {
240     if ( ui->onpe0() )
241       cout << "ResponseCmd:responseEfield: ElectricEnthalpy not defined\n"
242            << "The polarization variable may be set to OFF" << endl;
243     return;
244   }
245 
246   D3vector e_field[3] = { D3vector(amplitude,0.0,0.0),
247                           D3vector(0.0,amplitude,0.0),
248                           D3vector(0.0,0.0,amplitude) };
249   D3vector dipole_p[3], dipole_m[3];
250 
251   // save wave functions
252   Wavefunction wf0(s->wf);
253   wf0 = s->wf;
254 
255   // compute change in dipole in 3 directions by finite difference
256   D3vector e_field_base = el_enth->e_field();
257   for ( int idir = 0; idir < 3; idir++ )
258   {
259     el_enth->set_e_field(e_field_base+e_field[idir]);
260     stepper->step(0);
261     dipole_p[idir] = el_enth->dipole_total();
262     s->wf = wf0;
263 
264     el_enth->set_e_field(e_field_base-e_field[idir]);
265     stepper->step(0);
266     dipole_m[idir] = el_enth->dipole_total();
267     s->wf = wf0;
268   }
269 
270   D3vector ddx = dipole_p[0] - dipole_m[0];
271   D3vector ddy = dipole_p[1] - dipole_m[1];
272   D3vector ddz = dipole_p[2] - dipole_m[2];
273 
274   const UnitCell& cell = s->wf.cell();
275   cell.fold_in_ws(ddx);
276   cell.fold_in_ws(ddy);
277   cell.fold_in_ws(ddz);
278 
279   ddx *= 0.5 / amplitude;
280   ddy *= 0.5 / amplitude;
281   ddz *= 0.5 / amplitude;
282 
283   if ( ui->onpe0() )
284   {
285     cout << "<polarizability>" << endl;
286     cout << setprecision(6);
287     cout << " <a_xx> " << setw(10) << ddx.x << " </a_xx>"
288             " <a_yx> " << setw(10) << ddx.y << " </a_yx>"
289             " <a_zx> " << setw(10) << ddx.z << " </a_zx>" << endl;
290     cout << " <a_xy> " << setw(10) << ddy.x << " </a_xy>"
291             " <a_yy> " << setw(10) << ddy.y << " </a_yy>"
292             " <a_zy> " << setw(10) << ddy.z << " </a_zy>" << endl;
293     cout << " <a_xz> " << setw(10) << ddz.x << " </a_xz>"
294             " <a_yz> " << setw(10) << ddz.y << " </a_yz>"
295             " <a_zz> " << setw(10) << ddz.z << " </a_zz>" << endl;
296     double a_iso = (ddx.x+ddy.y+ddz.z)/3.0;
297     cout << " <a_iso> " << setw(10) << a_iso << " </a_iso>" << endl;
298     double beta[9] = { ddx.x-a_iso, ddx.y, ddx.z,
299                        ddy.x, ddy.y-a_iso, ddy.z,
300                        ddz.x, ddz.y, ddz.z-a_iso };
301     double a_aniso = 0.0;
302     for ( int i = 0; i < 9; i++ )
303       a_aniso += beta[i] * beta[i];
304     a_aniso *= 2.0/15.0;
305     cout << " <a_aniso> " << setw(10) << a_aniso << " </a_aniso>" << endl;
306     cout << "</polarizability>" << endl;
307   }
308 
309   delete stepper;
310 }
311 
312 ////////////////////////////////////////////////////////////////////////////////
responseVext(bool rpa,bool ipa,int nitscf,int nite,string fmt)313 void ResponseCmd::responseVext(bool rpa, bool ipa, int nitscf, int nite,
314    string fmt)
315 {
316   s->wf.info(cout, "wavefunction");
317 
318   const int nspin = s->wf.nspin();
319   BOSampleStepper *stepper = new BOSampleStepper(*s, nitscf, nite);
320 
321   if ( rpa )
322     stepper->set_update_vxc(false);
323 
324   if ( ipa )
325   {
326     stepper->set_update_vh(false);
327     stepper->set_update_vxc(false);
328   }
329 
330   assert(fmt == "xml" || fmt == "cube");
331 
332   // save a copy of initial wave functions
333   Wavefunction wf0(s->wf);
334   wf0 = s->wf;
335 
336   ChargeDensity &cd = stepper->cd();
337   MPI_Comm vcomm = MPIdata::g_comm();
338 
339   stepper->step(0);
340   vector<vector<double> > rhor1; // density with +Vext
341   rhor1 = cd.rhor;
342 
343   s->wf = wf0;
344   s->vext->reverse();
345 
346   stepper->step(0);
347   vector<vector<double> > rhor2; // density with -Vext
348   rhor2 = cd.rhor;
349   s->vext->reverse();
350 
351   // restore initial wave functions
352   s->wf = wf0;
353 
354   // compute drho_r = rhor1 - rhor2
355   const int np012loc = cd.vft()->np012loc();
356   const double omega = cd.vbasis()->cell().volume();
357   const double omega_inv = 1.0 / omega;
358   vector<vector<double> > drho_r;
359   vector<vector<complex<double> > > drho_r_tmp;
360   drho_r.resize(nspin);
361   drho_r_tmp.resize(nspin);
362   for (int ispin = 0; ispin < nspin; ispin++)
363   {
364     drho_r[ispin].resize(np012loc);
365     drho_r_tmp[ispin].resize(np012loc);
366     for (int ir = 0; ir < np012loc; ir++)
367     {
368       drho_r[ispin][ir] = rhor1[ispin][ir] - rhor2[ispin][ir];
369       drho_r_tmp[ispin][ir] = complex<double>(drho_r[ispin][ir] * omega, 0.0);
370     }
371   }
372 
373   // Fourier (forward) transform drho_r to the basis compatible
374   // with the vext grid size
375   const UnitCell &cell = cd.vbasis()->cell();
376   const double ecut = s->vext->ecut();
377   const int np0 = cd.vft()->np0();
378   const int np1 = cd.vft()->np1();
379   const int np2 = cd.vft()->np2();
380   Basis basis(vcomm, D3vector(0, 0, 0));
381   basis.resize(cell, cell, ecut);
382   FourierTransform ft1(basis, np0, np1, np2);
383 
384   vector<vector<complex<double> > > drho_g;
385   drho_g.resize(nspin);
386   for (int ispin = 0; ispin < nspin; ispin++)
387   {
388     drho_g[ispin].resize(basis.localsize());
389     ft1.forward(&drho_r_tmp[ispin][0], &drho_g[ispin][0]);
390   }
391 
392   // Fourier (backward) transform drho_g to drho_r on the
393   // vext grid (n0,n1,n2)
394   const int n0 = s->vext->n(0);
395   const int n1 = s->vext->n(1);
396   const int n2 = s->vext->n(2);
397   FourierTransform ft2(basis, n0, n1, n2);
398   const int n012loc = ft2.np012loc();
399 
400   for (int ispin = 0; ispin < nspin; ispin++)
401   {
402     drho_r[ispin].resize(n012loc);
403     drho_r_tmp[ispin].resize(n012loc);
404     ft2.backward(&drho_g[ispin][0], &drho_r_tmp[ispin][0]);
405     for (int ir = 0; ir < drho_r[ispin].size(); ir++)
406     {
407       drho_r[ispin][ir] = 0.5 * real(drho_r_tmp[ispin][ir] * omega_inv) /
408                           s->vext->amplitude();
409     }
410   }
411   // at this point, drho_r has been computed and distributed on column context
412 
413   // now write drho_r to disk
414   const Context& ctxt = s->wf.sd_context();
415   int nprow = ctxt.nprow();
416   int myrow = ctxt.myrow();
417   int mycol = ctxt.mycol();
418   Timer tm_write_drho;
419 
420   // serial write xml file or cube file
421   // processors on row 1 collect drho from other rows
422   // processor at row 1, column 1 write drho to disk
423   vector<double> drho_r_gathered;
424   if (myrow == 0)
425     drho_r_gathered.resize(ft2.np012());
426 
427   vector<int> rcounts(nprow, 0);
428   vector<int> displs(nprow, 0);
429   int displ = 0;
430   for (int iproc = 0; iproc < nprow; iproc++)
431   {
432     displs[iproc] = displ;
433     displ += ft2.np012loc(iproc);
434     rcounts[iproc] = ft2.np012loc(iproc);
435   }
436 
437   for (int ispin = 0; ispin < nspin; ispin++)
438   {
439     MPI_Gatherv(&drho_r[ispin][0], ft2.np012loc(), MPI_DOUBLE,
440                 &drho_r_gathered[0], &rcounts[0], &displs[0],
441                 MPI_DOUBLE, 0, vcomm);
442 
443     if ( myrow == 0 && mycol == 0 )
444     {
445       string filename;
446       if (nspin == 1)
447         filename = s->vext->filename() + ".response";
448       else
449         filename = s->vext->filename() + ".response."
450                    + ((ispin == 0) ? "spin0" : "spin1");
451       if (fmt == "xml")
452       {
453         tm_write_drho.start();
454         Function3d f;
455         f.name = "delta_rho";
456         f.nx = n0;
457         f.ny = n1;
458         f.nz = n2;
459         f.a = s->atoms.cell().a(0);
460         f.b = s->atoms.cell().a(1);
461         f.c = s->atoms.cell().a(2);
462         f.val = drho_r_gathered;
463         f.write(filename);
464         tm_write_drho.stop();
465       }
466       else
467       {
468         // write cube file
469         tm_write_drho.start();
470 
471         ofstream os(filename.c_str(), ios::out);
472         // comment lines (first 2 lines of cube file)
473         os << "Created " << isodate() << " by qbox-" << release() << "\n";
474         os << "Charge density response under external potential "
475            << s->vext->filename() << "\n";
476         // atoms and unit cell
477         int natoms = s->atoms.size();
478         D3vector a0 = s->atoms.cell().a(0);
479         D3vector a1 = s->atoms.cell().a(1);
480         D3vector a2 = s->atoms.cell().a(2);
481         os << natoms << " " << -0.5 * (a0 + a1 + a2) << "\n";
482         os << n0 << " " << a0 / n0 << endl;
483         os << n1 << " " << a1 / n1 << endl;
484         os << n2 << " " << a2 / n2 << endl;
485         const int nsp = s->atoms.nsp();
486         for (int is = 0; is < nsp; is++)
487         {
488           Species *sp = s->atoms.species_list[is];
489           const int z = sp->atomic_number();
490           const int na = s->atoms.na(is);
491           for (int ia = 0; ia < na; ia++)
492           {
493             Atom *ap = s->atoms.atom_list[is][ia];
494             os << setprecision(5);
495             os << z << " " << ((double) z) << " " << ap->position() << "\n";
496           }
497         }
498         // charge density response
499         os << setprecision(6) << std::scientific;
500         for (int nx = 0; nx < n0; nx++)
501           for (int ny = 0; ny < n1; ny++)
502             for (int nz = 0; nz < n2; nz++)
503             {
504               const int ir = nx + ny * n0 + nz * n0 * n1;
505               os << drho_r_gathered[ir] << "\n";
506             }
507         os.close();
508 
509         tm_write_drho.stop();
510       } // if fmt
511     } //if ( myrow == 0 && mycol == 0 )
512   } // for ispin
513 
514   delete stepper;
515 }
516