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