1 /*-------------------------------------------------------------------
2 Copyright 2014 Deniz Gunceler, Ravishankar Sundararaman
3
4 This file is part of JDFTx.
5
6 JDFTx is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 3 of the License, or
9 (at your option) any later version.
10
11 JDFTx is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with JDFTx. If not, see <http://www.gnu.org/licenses/>.
18 -------------------------------------------------------------------*/
19
20 #include <electronic/Dump.h>
21 #include <electronic/Everything.h>
22 #include <electronic/ColumnBundle.h>
23 #include <electronic/ColumnBundleTransform.h>
24 #include <electronic/Dump_internal.h>
25 #include <core/WignerSeitz.h>
26 #include <core/Operators.h>
27 #include <core/Units.h>
28 #include <core/LatticeUtils.h>
29
30 //---------------------- Excitations -----------------------------------------------
31
dumpExcitations(const Everything & e,const char * filename)32 void dumpExcitations(const Everything& e, const char* filename)
33 {
34 const GridInfo& g = e.gInfo;
35
36 struct excitation
37 { int q,o,u;
38 double dE;
39 double dreal, dimag, dnorm;
40
41 excitation(int q, int o, int u, double dE, double dreal, double dimag, double dnorm): q(q), o(o), u(u), dE(dE), dreal(dreal), dimag(dimag), dnorm(dnorm){};
42
43 inline bool operator<(const excitation& other) const {return dE<other.dE;}
44 void print(FILE* fp) const { fprintf(fp, "%5i %3i %3i %12.5e %12.5e %12.5e %12.5e\n", q, o, u, dE, dreal, dimag, dnorm); }
45 };
46 std::vector<excitation> excitations;
47
48 double maxHOMO=-DBL_MAX, minLUMO=DBL_MAX; // maximum (minimum) of all HOMOs (LUMOs) in all qnums
49 int maxHOMOq=0, minLUMOq=0, maxHOMOn=0, minLUMOn=0; //Indices and energies for the indirect gap
50
51 //Select relevant eigenvals:
52 std::vector<diagMatrix> eigsQP;
53 if(e.exCorr.orbitalDep && e.dump.count(std::make_pair(DumpFreq_End, DumpOrbitalDep)))
54 { //Search for an eigenvalsQP file:
55 string fname = e.dump.getFilename("eigenvalsQP");
56 FILE* fp = fopen(fname.c_str(), "r");
57 if(fp)
58 { fclose(fp);
59 eigsQP.resize(e.eInfo.nStates);
60 e.eInfo.read(eigsQP, fname.c_str());
61 }
62 }
63 const std::vector<diagMatrix>& eigs = eigsQP.size() ? eigsQP : e.eVars.Hsub_eigs;
64
65 // Integral kernel's for Fermi's golden rule
66 VectorField r;
67 { nullToZero(r, g);
68 logSuspend();
69 WignerSeitz ws(e.gInfo.R);
70 logResume();
71 const vector3<>& x0 = e.coulombParams.embedCenter; //origin
72
73 vector3<double*> rData(r.data());
74 matrix3<> invS = inv(Diag(vector3<>(e.gInfo.S)));
75 vector3<int> iv;
76 size_t i = 0;
77 for(iv[0]=0; iv[0]<e.gInfo.S[0]; iv[0]++)
78 for(iv[1]=0; iv[1]<e.gInfo.S[1]; iv[1]++)
79 for(iv[2]=0; iv[2]<e.gInfo.S[2]; iv[2]++)
80 { vector3<> x = x0 + ws.restrict(invS*iv - x0); //lattice coordinates wrapped to WS centered on x0
81 storeVector(e.gInfo.R*x, rData, i++);
82 }
83 }
84
85 //Find and cache all excitations in system (between same qnums)
86 bool insufficientBands = false;
87 for(int q=e.eInfo.qStart; q<e.eInfo.qStop; q++)
88 { //Find local HOMO and check band sufficiency:
89 int HOMO = e.eInfo.findHOMO(q);
90 if(HOMO+1>=e.eInfo.nBands) { insufficientBands=true; break; }
91
92 //Update global HOMO and LUMO of current process:
93 if(eigs[q][HOMO] > maxHOMO) { maxHOMOq = q; maxHOMOn = HOMO; maxHOMO = eigs[q][HOMO]; }
94 if(eigs[q][HOMO+1] < minLUMO) { minLUMOq = q; minLUMOn = HOMO+1; minLUMO = eigs[q][HOMO+1]; }
95
96 for(int o=HOMO; o>=0; o--)
97 { for(int u=(HOMO+1); u<e.eInfo.nBands; u++)
98 { vector3<> dreal, dimag, dnorm;
99 for(int iDir=0; iDir<3; iDir++)
100 { complex xi = integral(I(e.eVars.C[q].getColumn(u,0)) * r[iDir] * I(e.eVars.C[q].getColumn(o,0)));
101 dreal[iDir] = xi.real();
102 dimag[iDir] = xi.imag();
103 dnorm[iDir] = xi.abs();
104 }
105 double dE = eigs[q][u]-eigs[q][o]; //Excitation energy
106 excitations.push_back(excitation(q, o, u, dE, dreal.length_squared(), dimag.length_squared(), dnorm.length_squared()));
107 }
108 }
109 }
110 mpiWorld->allReduce(insufficientBands, MPIUtil::ReduceLOr);
111 if(insufficientBands)
112 { logPrintf("Insufficient bands to calculate excited states!\n");
113 logPrintf("Increase the number of bands (elec-n-bands) and try again!\n");
114 return;
115 }
116
117 //Transmit results to head process:
118 if(mpiWorld->isHead())
119 { excitations.reserve(excitations.size() * mpiWorld->nProcesses());
120 for(int jProcess=1; jProcess<mpiWorld->nProcesses(); jProcess++)
121 { //Receive data:
122 size_t nExcitations; mpiWorld->recv(nExcitations, jProcess, 0);
123 std::vector<int> msgInt(4 + nExcitations*3);
124 std::vector<double> msgDbl(2 + nExcitations*4);
125 mpiWorld->recvData(msgInt, jProcess, 1);
126 mpiWorld->recvData(msgDbl, jProcess, 2);
127 //Unpack:
128 std::vector<int>::const_iterator intPtr = msgInt.begin();
129 std::vector<double>::const_iterator dblPtr = msgDbl.begin();
130 //--- globals:
131 int j_maxHOMOq = *(intPtr++); int j_maxHOMOn = *(intPtr++); double j_maxHOMO = *(dblPtr++);
132 int j_minLUMOq = *(intPtr++); int j_minLUMOn = *(intPtr++); double j_minLUMO = *(dblPtr++);
133 if(j_maxHOMO > maxHOMO) { maxHOMOq=j_maxHOMOq; maxHOMOn=j_maxHOMOn; maxHOMO=j_maxHOMO; }
134 if(j_minLUMO < minLUMO) { minLUMOq=j_minLUMOq; minLUMOn=j_minLUMOn; minLUMO=j_minLUMO; }
135 //--- excitation array:
136 for(size_t iExcitation=0; iExcitation<nExcitations; iExcitation++)
137 { int q = *(intPtr++); int o = *(intPtr++); int u = *(intPtr++);
138 double dE = *(dblPtr++);
139 double dreal = *(dblPtr++); double dimag = *(dblPtr++); double dnorm = *(dblPtr++);
140 excitations.push_back(excitation(q, o, u, dE, dreal, dimag, dnorm));
141 }
142 }
143 }
144 else
145 { //Pack data:
146 std::vector<int> msgInt; std::vector<double> msgDbl;
147 size_t nExcitations = excitations.size();
148 msgInt.reserve(4 + nExcitations*3);
149 msgDbl.reserve(2 + nExcitations*4);
150 msgInt.push_back(maxHOMOq); msgInt.push_back(maxHOMOn); msgDbl.push_back(maxHOMO);
151 msgInt.push_back(minLUMOq); msgInt.push_back(minLUMOn); msgDbl.push_back(minLUMO);
152 for(const excitation& e: excitations)
153 { msgInt.push_back(e.q); msgInt.push_back(e.o); msgInt.push_back(e.u);
154 msgDbl.push_back(e.dE);
155 msgDbl.push_back(e.dreal); msgDbl.push_back(e.dimag); msgDbl.push_back(e.dnorm);
156 }
157 //Send data:
158 mpiWorld->send(nExcitations, 0, 0);
159 mpiWorld->sendData(msgInt, 0, 1);
160 mpiWorld->sendData(msgDbl, 0, 2);
161 }
162
163 //Process and print excitations:
164 if(!mpiWorld->isHead()) return;
165
166 FILE* fp = fopen(filename, "w");
167 if(!fp) die_alone("Error opening %s for writing.\n", filename);
168
169 std::sort(excitations.begin(), excitations.end());
170 const excitation& opt = excitations.front();
171 fprintf(fp, "Using %s eigenvalues. HOMO: %.5f LUMO: %.5f \n", eigsQP.size() ? "discontinuity-corrected QP" : "KS", maxHOMO, minLUMO);
172 fprintf(fp, "Optical (direct) gap: %.5e (from n = %i to %i in qnum = %i)\n", opt.dE, opt.o, opt.u, opt.q);
173 fprintf(fp, "Indirect gap: %.5e (from (%i, %i) to (%i, %i))\n\n", minLUMO-maxHOMO, maxHOMOq, maxHOMOn, minLUMOq, minLUMOn);
174
175 fprintf(fp, "Optical excitation energies and corresponding electric dipole transition strengths\n");
176 fprintf(fp, "qnum i f dE |<psi1|r|psi2>|^2 (real, imag, norm)\n");
177 for(const excitation& e: excitations) e.print(fp);
178 fclose(fp);
179 }
180
181 //---------------------------- FCI ------------------------------------
182
183 //Split doublet index i12 that ends on N^2 or N(N+1)/2 depending on sym = false or true, into i1 and i2
splitDoubletIndex(size_t i12,size_t N,bool sym,size_t & i1,size_t & i2)184 inline void splitDoubletIndex(size_t i12, size_t N, bool sym, size_t& i1, size_t& i2)
185 { if(sym) //N(N+1)/2 case
186 { i1 = 0; while((i1+1)*(i1+2)<=2*i12) i1++;
187 i2 = i12 - (i1*(i1+1))/2;
188 }
189 else //N^2 case
190 { i1 = i12 / N;
191 i2 = i12 - i1 * N;
192 }
193 }
194
dumpFCI(const Everything & e,const char * filename)195 void dumpFCI(const Everything& e, const char* filename)
196 { static StopWatch watch("dumpFCI"); watch.start();
197
198 const size_t nBands = e.eInfo.nBands;
199 const int nStates = e.eInfo.nStates;
200 const int nOrbitals = nBands * nStates;
201 const int nSpinor = e.eInfo.spinorLength();
202 const double Kcut = 1e-8; //threshold for non-zero matrix element
203
204 FILE* fp = NULL;
205 if(mpiWorld->isHead())
206 { fp = fopen(filename, "w");
207 if(!fp) die_alone("Error opening %s for writing.\n", filename);
208
209 //Write header:
210 fprintf(fp, " &FCI\n NORB=%d,\n NELEC=%d,\n MS2=%d,\n UHF=TRUE,\n ORBSYM= ",
211 nOrbitals, int(round(e.eInfo.nElectrons)), int(round(fabs(e.eInfo.Minitial))) );
212 for(int i=0; i<nOrbitals; i++)
213 fprintf(fp, "0, ");
214 fprintf(fp, "\n ISYM=1,\n NPROP= 1 1 1,\n PROPBITLEN= 15\n &END\n");
215 }
216
217 //Coulomb integrals:
218 #define GETwfns(a) \
219 ColumnBundle Ctmp##a; \
220 if(not e.eInfo.isMine(q##a)) \
221 Ctmp##a.init(nBands, e.basis[q##a].nbasis*nSpinor, &(e.basis[q##a]), &(e.eInfo.qnums[q##a]), isGpuEnabled()); \
222 const ColumnBundle& C##a = e.eInfo.isMine(q##a) ? e.eVars.C[q##a] : Ctmp##a; \
223 mpiWorld->bcastData((ColumnBundle&)C##a, e.eInfo.whose(q##a));
224 #define UPDATE_conjIpsi(a) \
225 for(int s=0; s<nSpinor; s++) \
226 conjIpsi##a[s] = conj(I(C##a.getColumn(i##a,s)));
227 #define UPDATE_nPair(a,b, dk) \
228 { complexScalarField nPair; \
229 for(int s=0; s<nSpinor; s++) \
230 nPair += conjIpsi##a[s] * I(C##b.getColumn(i##b,s)); \
231 if((dk).length_squared() > symmThresholdSq) \
232 multiplyBlochPhase(nPair, dk); /*correct for k-wrapping upto reciprocal lattice vector*/ \
233 n##a##b = J(nPair); \
234 }
235 #define UPDATE_12 \
236 { complexScalarFieldTilde n12; \
237 UPDATE_nPair(1,2, vector3<>()) \
238 Kn12 = O((*e.coulombWfns)(n12, dk21, 0.)); \
239 }
240 #define UPDATE_34 \
241 UPDATE_nPair(3,4, dk43-dk21)
242
243 //--- Quadruple loop over k-spin indices such that 12 and 34 have same spin indices
244 for(int q1=0; q1<nStates; q1++)
245 { GETwfns(1)
246 for(int q2=0; q2<=q1; q2++) if(e.eInfo.qnums[q1].spin == e.eInfo.qnums[q2].spin)
247 { GETwfns(2)
248 vector3<> dk21 = e.eInfo.qnums[q2].k - e.eInfo.qnums[q1].k;
249 for(int q3=0; q3<nStates; q3++)
250 { GETwfns(3)
251 for(int q4=0; q4<=q3; q4++) if(e.eInfo.qnums[q3].spin == e.eInfo.qnums[q4].spin)
252 { vector3<> dk43 = e.eInfo.qnums[q4].k - e.eInfo.qnums[q3].k;
253 if(circDistanceSquared(dk21, dk43) > symmThresholdSq) continue; //momentum conservation
254 GETwfns(4)
255
256 //Split band quadruplets over MPI:
257 size_t N12 = (q1==q2 ? (nBands*(nBands+1))/2 : nBands*nBands);
258 size_t N34 = (q3==q4 ? (nBands*(nBands+1))/2 : nBands*nBands);
259 size_t N1234 = N12 * N34;
260 size_t i1234start=0, i1234stop=0;
261 TaskDivision(N1234, mpiWorld).myRange(i1234start, i1234stop);
262
263 //Initial indices on this process:
264 size_t i1234 = i1234start; //quaruplet index
265 size_t i12 = i1234 / N34; //12 doublet index
266 size_t i34 = i1234 - i12*N34; //34 doublet index
267 size_t i1, i2, i3, i4; //individual band indices
268 splitDoubletIndex(i12, nBands, q1==q2, i1, i2);
269 splitDoubletIndex(i34, nBands, q3==q4, i3, i4);
270
271 //Initialize wavefunctions, pair densities for initial indices:
272 std::vector<complexScalarField> conjIpsi1(nSpinor), conjIpsi3(nSpinor);
273 complexScalarFieldTilde Kn12, n34;
274 //--- 12 pair
275 UPDATE_conjIpsi(1)
276 UPDATE_12
277 //--- 34 pair
278 UPDATE_conjIpsi(3)
279 UPDATE_34
280
281 //Loop over band quadruplets for this process:
282 ostringstream oss; //buffer containing output from this process
283 while(i1234 < i1234stop)
284 {
285 complex K1234 = dot(Kn12, n34);
286 if(K1234.abs() > Kcut)
287 { if(fabs(K1234.imag()) > Kcut)
288 oss << "( " << K1234.real() << ", " << K1234.imag() << ") "; //complex element
289 else
290 oss << K1234.real() << ' '; //real element
291 oss << i1*nStates+q1+1 << ' '
292 << i2*nStates+q2+1 << ' '
293 << i3*nStates+q3+1 << ' '
294 << i4*nStates+q4+1 << '\n';
295 }
296
297 i1234++; if(i1234==i1234stop) break;
298 i4++;
299 if(i4==(q3==q4 ? i3+1 : nBands))
300 { i4=0;
301 i3++;
302 if(i3==nBands)
303 { i3=0;
304 i2++;
305 if(i2==(q1==q2 ? i1+1 : nBands))
306 { i2=0;
307 i1++;
308 UPDATE_conjIpsi(1)
309 }
310 UPDATE_12
311 }
312 UPDATE_conjIpsi(3)
313 }
314 UPDATE_34
315 }
316
317 //Write matrix elements from HEAD:
318 if(mpiWorld->isHead())
319 { fputs(oss.str().c_str(), fp); //local data
320 oss.clear();
321 for(int jProcess=1; jProcess<mpiWorld->nProcesses(); jProcess++)
322 { string buf; mpiWorld->recv(buf, jProcess, 0);
323 fputs(buf.c_str(), fp); //data from other processes
324 }
325 fflush(fp);
326 }
327 else mpiWorld->send(oss.str(), 0, 0);
328 }
329 }
330 }
331 }
332 #undef GETwfns
333 #undef UPDATE_conjIpsi
334 #undef UPDATE_nPair
335 #undef UPDATE_12
336 #undef UPDATE_34
337
338 //Compute and collect one-particle matrix elements / eigenvalues:
339 std::vector<matrix> H0(nStates);
340 std::vector<diagMatrix> E(nStates);
341 for(int q=e.eInfo.qStart; q<e.eInfo.qStop; q++)
342 {
343 const ColumnBundle& Cq = e.eVars.C[q];
344 ScalarFieldArray Vloc(e.eInfo.nDensities); nullToZero(Vloc, e.gInfo);
345 for(unsigned s=0; s<Vloc.size(); s++)
346 if(s<2) Vloc[s] = Jdag(O(e.iInfo.Vlocps));
347 ColumnBundle H0Cq = Idag_DiagV_I(Cq, Vloc) - 0.5*L(Cq); //local PS and KE
348 { //Add nonlocal PS contribution:
349 std::vector<matrix> HVdagCq(e.iInfo.species.size());
350 e.iInfo.EnlAndGrad(e.eInfo.qnums[q], e.eVars.F[q], e.eVars.VdagC[q], HVdagCq);
351 e.iInfo.projectGrad(HVdagCq, Cq, H0Cq); //HCq += (nonlocalPS projectors) * Cq
352 }
353 H0[q] = Cq ^ H0Cq; //matrix elements of KE + net pseudopotential
354 E[q] = e.eVars.Hsub_eigs[q];
355
356 if(not mpiWorld->isHead())
357 { mpiWorld->sendData(E[q], 0, q);
358 mpiWorld->sendData(H0[q], 0, q);
359 }
360 }
361 if(mpiWorld->isHead())
362 { for(int q=0; q<nStates; q++)
363 if(!e.eInfo.isMine(q))
364 { H0[q].init(nBands, nBands);
365 E[q].resize(nBands);
366 mpiWorld->recvData(E[q], e.eInfo.whose(q), q);
367 mpiWorld->recvData(H0[q], e.eInfo.whose(q), q);
368 }
369 }
370
371 //Output one-particle matrix elements / eigenvalues:
372 if(mpiWorld->isHead())
373 { //One-particle
374 for(int q=0; q<nStates; q++)
375 for(size_t b1=0; b1<nBands; b1++)
376 for(size_t b2=0; b2<=b1; b2++)
377 { complex M = H0[q](b1,b2);
378 if(M.abs() > Kcut)
379 { if(fabs(M.imag()) > Kcut)
380 fprintf(fp, "( %lg, %lg )", M.real(), M.imag()); //complex element
381 else
382 fprintf(fp, "%lg", M.real()); //real element
383 fprintf(fp, " %lu %lu 0 0\n", b1*nStates+q+1, b2*nStates+q+1);
384 }
385 }
386 //Eigenvalues
387 for(int q=0; q<nStates; q++)
388 for(size_t b=0; b<nBands; b++)
389 fprintf(fp, "%lg %lu 0 0 0\n", E[q][b], b*nStates+q+1);
390 //Ewald energy
391 double Eewald = e.ener.E["Eewald"];
392 if(fabs(Eewald) > Kcut)
393 fprintf(fp, "%lg 0 0 0 0\n", Eewald);
394 fclose(fp);
395 }
396 watch.stop();
397 }
398
399 //---------------------------- Moments ------------------------------------
400
dumpMoment(const Everything & e,const char * filename)401 void dumpMoment(const Everything& e, const char* filename)
402 { logSuspend();
403 WignerSeitz ws(e.gInfo.R);
404 logResume();
405 const vector3<>& x0 = e.coulombParams.embedCenter; //origin
406
407 //Compute electronic moment in lattice coordinates:
408 vector3<> elecMoment;
409 ScalarField n = e.eVars.get_nTot();
410 double* nData = n->data();
411 matrix3<> invS = inv(Diag(vector3<>(e.gInfo.S)));
412 vector3<int> iv;
413 for(iv[0]=0; iv[0]<e.gInfo.S[0]; iv[0]++)
414 for(iv[1]=0; iv[1]<e.gInfo.S[1]; iv[1]++)
415 for(iv[2]=0; iv[2]<e.gInfo.S[2]; iv[2]++)
416 { vector3<> x = x0 + ws.restrict(invS*iv - x0); //lattice coordinates wrapped to WS centered on x0
417 elecMoment += x * (*(nData)++); //collect in lattice coordinates
418 }
419 elecMoment *= e.gInfo.dV; //integration weight
420
421 //Compute ionic moment in lattice coordinates:
422 vector3<> ionMoment;
423 for(auto sp: e.iInfo.species)
424 for(vector3<> pos: sp->atpos)
425 { vector3<> x = x0 + ws.restrict(pos - x0); //wrap atom position to WS centered on x0
426 ionMoment -= sp->Z * x;
427 }
428
429 //Zero contributions in periodic directions:
430 for(int iDir=0; iDir<3; iDir++)
431 if(not e.coulombParams.isTruncated()[iDir])
432 { elecMoment[iDir] = 0.;
433 ionMoment[iDir] = 0.;
434 }
435
436 //Convert to Cartesian and compute total:
437 elecMoment = e.gInfo.R * elecMoment;
438 ionMoment = e.gInfo.R * ionMoment;
439 vector3<> totalMoment = elecMoment + ionMoment;
440
441 //Output:
442 FILE* fp = fopen(filename, "w");
443 if(!fp) die("Error opening %s for writing.\n", filename);
444 fprintf(fp, "#Contribution %4s %12s %12s\n", "x", "y", "z");
445 fprintf(fp, "%10s %12.6lf %12.6lf %12.6lf\n", "Electronic", elecMoment[0], elecMoment[1], elecMoment[2]);
446 fprintf(fp, "%10s %12.6lf %12.6lf %12.6lf\n", "Ionic", ionMoment[0], ionMoment[1], ionMoment[2]);
447 fprintf(fp, "%10s %12.6lf %12.6lf %12.6lf\n", "Total", totalMoment[0], totalMoment[1], totalMoment[2]);
448 fclose(fp);
449 }
450
451
452 //---------------------------- XC analysis --------------------------------
453
454 namespace XC_Analysis
455 {
456 static const double cutoff = 1e-8;
457
regularize(int i,vector3<> r,double * tau)458 void regularize(int i, vector3<> r, double* tau)
459 { tau[i] = (tau[i]<cutoff ? 0. : tau[i]);
460 }
spness_kernel(int i,vector3<> r,double * tauW,double * tau,double * spness)461 void spness_kernel(int i, vector3<> r, double* tauW, double* tau, double* spness)
462 { spness[i] = (tau[i]>cutoff ? tauW[i]/tau[i] : 1.);
463 }
464
tauWeizsacker(const Everything & e)465 ScalarFieldArray tauWeizsacker(const Everything& e)
466 { ScalarFieldArray tauW(e.eVars.n.size());
467 for(size_t j=0; j<e.eVars.n.size(); j++)
468 tauW[j] = (1./8.)*lengthSquared(gradient(e.eVars.n[j]))*pow(e.eVars.n[j], -1);
469 return tauW;
470 }
471
spness(const Everything & e)472 ScalarFieldArray spness(const Everything& e)
473 {
474 const auto& tau = (e.exCorr.needsKEdensity() ? e.eVars.tau : e.eVars.KEdensity());
475
476 ScalarFieldArray tauW = tauWeizsacker(e);
477 ScalarFieldArray spness(e.eVars.n.size()); nullToZero(spness, e.gInfo);
478 for(size_t j=0; j<e.eVars.n.size(); j++)
479 { applyFunc_r(e.gInfo, regularize, tau[j]->data());
480 applyFunc_r(e.gInfo, regularize, tauW[j]->data());
481 //spness[j] = tauW[j]*pow(tau[j], -1);
482 applyFunc_r(e.gInfo, spness_kernel, tauW[j]->data(), tau[j]->data(), spness[j]->data());
483 }
484
485 return spness;
486 }
487
sHartree(const Everything & e)488 ScalarFieldArray sHartree(const Everything& e)
489 {
490 ScalarFieldArray sVh(e.eVars.n.size());
491 ScalarFieldTildeArray nTilde = J(e.eVars.n);
492 for(size_t s=0; s<e.eVars.n.size(); s++)
493 sVh[s] = I((*(e.coulomb))(nTilde[s]));
494
495 return sVh;
496 }
497 }
498
499
500 //-------------------------- Solvation radii ----------------------------------
501
set_rInv(size_t iStart,size_t iStop,const vector3<int> & S,const matrix3<> & RTR,const WignerSeitz * ws,double * rInv)502 inline void set_rInv(size_t iStart, size_t iStop, const vector3<int>& S, const matrix3<>& RTR, const WignerSeitz* ws, double* rInv)
503 { vector3<> invS; for(int k=0; k<3; k++) invS[k] = 1./S[k];
504 matrix3<> meshMetric = Diag(invS) * RTR * Diag(invS);
505 const double rWidth = 1.;
506 THREAD_rLoop(
507 double r = sqrt(meshMetric.metric_length_squared(ws->restrict(iv, S, invS)));
508 if(r < rWidth) //Polynomial with f',f" zero at origin and f,f',f" matched at rWidth
509 { double t = r/rWidth;
510 rInv[i] = (1./rWidth) * (1. + 0.5*(1.+t*t*t*(-2.+t)) + (1./12)*(1.+t*t*t*(-4.+t*3.))*2. );
511 }
512 else rInv[i] = 1./r;
513 )
514 }
515
dumpRsol(ScalarField nbound,string fname)516 void Dump::dumpRsol(ScalarField nbound, string fname)
517 {
518 //Compute normalization factor for the partition:
519 int nAtomsTot = 0; for(const auto& sp: e->iInfo.species) nAtomsTot += sp->atpos.size();
520 const double nFloor = 1e-5/nAtomsTot; //lower cap on densities to prevent Nyquist noise in low density regions
521 ScalarField nAtomicTot;
522 for(const auto& sp: e->iInfo.species)
523 { RadialFunctionG nRadial;
524 logSuspend(); sp->getAtom_nRadial(0,0, nRadial, true); logResume();
525 for(unsigned atom=0; atom<sp->atpos.size(); atom++)
526 { ScalarField nAtomic = radialFunction(e->gInfo, nRadial, sp->atpos[atom]);
527 double nMin, nMax; callPref(eblas_capMinMax)(e->gInfo.nr, nAtomic->dataPref(), nMin, nMax, nFloor);
528 nAtomicTot += nAtomic;
529 }
530 }
531 ScalarField nboundByAtomic = (nbound*nbound) * inv(nAtomicTot);
532
533 ScalarField rInv0(ScalarFieldData::alloc(e->gInfo));
534 { logSuspend(); WignerSeitz ws(e->gInfo.R); logResume();
535 threadLaunch(set_rInv, e->gInfo.nr, e->gInfo.S, e->gInfo.RTR, &ws, rInv0->data());
536 }
537
538 //Compute bound charge 1/r and 1/r^2 expectation values weighted by atom-density partition:
539 FILE* fp = fopen(fname.c_str(), "w");
540 fprintf(fp, "#Species rMean +/- rSigma [bohrs] (rMean +/- rSigma [Angstrom]) sqrt(Int|nbound^2|) in partition\n");
541 for(const auto& sp: e->iInfo.species)
542 { RadialFunctionG nRadial;
543 logSuspend(); sp->getAtom_nRadial(0,0, nRadial, true); logResume();
544 for(unsigned atom=0; atom<sp->atpos.size(); atom++)
545 { ScalarField w = radialFunction(e->gInfo, nRadial, sp->atpos[atom]) * nboundByAtomic;
546 //Get r centered at current atom:
547 ScalarFieldTilde trans; nullToZero(trans, e->gInfo); initTranslation(trans, e->gInfo.R*sp->atpos[atom]);
548 ScalarField rInv = I(trans * J(rInv0));
549 //Compute moments:
550 double wNorm = integral(w);
551 double rInvMean = integral(w * rInv) / wNorm;
552 double rInvSqMean = integral(w * rInv * rInv) / wNorm;
553 double rInvSigma = sqrt(rInvSqMean - rInvMean*rInvMean);
554 double rMean = 1./rInvMean;
555 double rSigma = rInvSigma / (rInvMean*rInvMean);
556 //Print stats:
557 fprintf(fp, "Rsol %s %.2lf +/- %.2lf ( %.2lf +/- %.2lf A ) Qrms: %.1le\n", sp->name.c_str(),
558 rMean, rSigma, rMean/Angstrom, rSigma/Angstrom, sqrt(wNorm));
559 }
560 }
561 fclose(fp);
562 }
563
dumpUnfold()564 void Dump::dumpUnfold()
565 { const GridInfo& gInfo = e->gInfo;
566 const ElecInfo& eInfo = e->eInfo;
567 //Header
568 //--- unit cell grid:
569 const matrix3<int>& M = Munfold;
570 matrix3<> invM = inv(matrix3<>(M));
571 int nUnits = abs(det(M));
572 GridInfo gInfoUnit;
573 gInfoUnit.R = gInfo.R * invM;
574 gInfoUnit.Gmax = sqrt(2*e->cntrl.Ecut); //Ecut = 0.5 Gmax^2
575 logSuspend();
576 gInfoUnit.initialize();
577 logResume();
578 //--- list of unit cell k-points for each supercell k:
579 std::vector<std::vector<vector3<>>> kUnit(eInfo.nStates);
580 for(int q=0; q<eInfo.nStates; q++)
581 { PeriodicLookup<vector3<>> plook(kUnit[q], gInfoUnit.GGT, nUnits);
582 for(const vector3<int>& iG: e->basis[q].iGarr)
583 { vector3<> k = (eInfo.qnums[q].k + iG) * invM; //corresponding unit cell k
584 for(int iDir=0; iDir<3; iDir++) k[iDir] -= ceil(-0.5+k[iDir]); //wrap to (-0.5,0.5]
585 if(plook.find(k) == string::npos) //not encountered yet
586 { plook.addPoint(kUnit[q].size(), k); //update periodic look-up table
587 kUnit[q].push_back(k);
588 if(int(kUnit[q].size())==nUnits) break; //all unit cell k for this supercell k found
589 }
590 }
591 std::sort(kUnit[q].begin(), kUnit[q].end());
592 }
593 //--- write header:
594 string fname = getFilename("bandUnfoldHeader");
595 logPrintf("Dumping '%s' ... ", fname.c_str()); logFlush();
596 if(mpiWorld->isHead())
597 { FILE* fp = fopen(fname.c_str(), "w");
598 for(int q=0; q<eInfo.nStates; q++)
599 { fprintf(fp, "Supercell kpoint ");
600 eInfo.kpointPrint(fp, q, true);
601 fprintf(fp, "\n");
602 for(vector3<>& k: kUnit[q])
603 fprintf(fp, "%+.7f %+.7f %+.7f\n", k[0], k[1], k[2]);
604 fprintf(fp, "\n");
605 }
606 fclose(fp);
607 }
608 logPrintf("done\n"); logFlush();
609 //Weights of each unit cell k:
610 fname = getFilename("bandUnfold");
611 logPrintf("Dumping '%s' ... ", fname.c_str()); logFlush();
612 std::vector<diagMatrix> unitWeights(eInfo.nStates);
613 for(int q=eInfo.qStart; q<eInfo.qStop; q++)
614 { unitWeights[q].reserve(nUnits * eInfo.nBands);
615 for(const vector3<>& k: kUnit[q])
616 { //Create basis and qnum for unit cell k:
617 Basis basisUnit;
618 logSuspend();
619 basisUnit.setup(gInfoUnit, e->iInfo, e->cntrl.Ecut, k);
620 logResume();
621 QuantumNumber qnumUnit; qnumUnit.k = k;
622 //Setup wave function transformation:
623 int nSpinor = eInfo.spinorLength();
624 ColumnBundleTransform cbt(k, basisUnit, eInfo.qnums[q].k, e->basis[q], nSpinor, SpaceGroupOp(), +1, M);
625 //Extract unit cell wavefunctions:
626 const ColumnBundle& C = e->eVars.C[q];
627 ColumnBundle OC = O(C);
628 ColumnBundle Cunit(1, basisUnit.nbasis*nSpinor, &basisUnit, &qnumUnit, isGpuEnabled());
629 ColumnBundle OCunit = Cunit.similar();
630 for(int b=0; b<eInfo.nBands; b++)
631 { Cunit.zero(); cbt.gatherAxpy(1., C,b, Cunit,0);
632 OCunit.zero(); cbt.gatherAxpy(1., OC,b, OCunit,0);
633 unitWeights[q].push_back(trace(Cunit^OCunit).real());
634 }
635 }
636 }
637 eInfo.write(unitWeights, fname.c_str(), nUnits*eInfo.nBands);
638 logPrintf("done\n"); logFlush();
639 }
640