1 /*-------------------------------------------------------------------
2 Copyright 2011 Ravishankar Sundararaman, Kendra Letchworth-Weaver
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 <fluid/ErfFMTweight.h>
21 #include <fluid/Molecule.h>
22 #include <core/Operators.h>
23 #include <electronic/ColumnBundle.h>
24
Site(string name,int atomicNumber)25 Molecule::Site::Site(string name, int atomicNumber) : name(name), Rhs(0), atomicNumber(atomicNumber), Znuc(0), sigmaNuc(0), Zelec(0), aElec(0), Zsite(0), deltaS(0), sigmaElec(0), rcElec(0), alpha(0), aPol(0), initialized(false)
26 {
27 }
28
~Site()29 Molecule::Site::~Site()
30 { free();
31 }
32
free()33 void Molecule::Site::free()
34 { if(initialized)
35 { elecKernel.free();
36 chargeKernel.free();
37 polKernel.free();
38 if(Rhs)
39 { w0.free();
40 w1.free();
41 w2.free();
42 w3.free();
43 w1v.free();
44 w2m.free();
45 }
46 }
47 }
48
49 //Functions for use in calculating site charge kernels for fluid-electron coupling
50 //--- normalized real space exponential function
Exponential(double r,double norm,double a)51 inline double Exponential(double r, double norm, double a)
52 { return norm/(8.*M_PI*pow(a,3))*exp(-r/a);
53 }
54 //--- un-normalized real space peaked exponential function
PeakedExponential(double r,double a,double rc,double sigma)55 inline double PeakedExponential(double r, double a, double rc, double sigma)
56 { return 0.5*exp(-r/a)*erfc((rc-r)/sigma);
57 }
58
59 //Create the radial kernel and initialize its weights (expensive) only once and then reuse
getLogGridKernel()60 RadialFunctionR getLogGridKernel()
61 { static RadialFunctionR KernelR;
62 if(!KernelR.r.size())
63 { int NRpts = 10000; //number of realspace log grid points used to represent electron density kernel
64 double rVecMin = 1e-7; //minimum value in realspace log grid
65 double rLogScale = 1.005; //rLogScale=r[i+1]/r[i]
66 KernelR.r.resize(NRpts);
67 for(int i=0; i<NRpts; i++)
68 KernelR.r[i] = i ? (KernelR.r[i-1] * rLogScale) : rVecMin;
69 KernelR.dr.resize(NRpts);
70 KernelR.f.resize(NRpts);
71 KernelR.initWeights();
72 }
73 return KernelR;
74 }
75
setup(const GridInfo & gInfo)76 void Molecule::Site::setup(const GridInfo& gInfo)
77 { if(initialized) free();
78 double dG = gInfo.dGradial;
79 int nGridLoc = int(ceil(gInfo.GmaxGrid/dG))+5;
80 logPrintf(" Initializing site '%s'\n", name.c_str());
81
82 //Initialize electron density kernel:
83 if(elecFilename.length() || elecFilenameG.length() || Zelec)
84 { logPrintf(" Electron density: ");
85 if(!Zelec || !sigmaElec)
86 { logPrintf("cuspless exponential with width %lg and norm %lg\n", aElec, Zelec);
87 elecKernel.init(0, dG, gInfo.GmaxGrid, RadialFunctionG::cusplessExpTilde, Zelec, aElec);
88 deltaS -= 12.*M_PI*Zelec*pow(aElec,2);
89 }
90 else
91 { logPrintf("proportional to exp(-r/%lg)*erfc((r-%lg)/%lg) with norm %lg\n", aElec, rcElec, sigmaElec, Zelec);
92 //--- construct the radial function in real space:
93 RadialFunctionR KernelR(getLogGridKernel());
94 for(unsigned i=0; i<KernelR.r.size(); i++)
95 KernelR.f[i] = PeakedExponential(KernelR.r[i], aElec, rcElec, sigmaElec);
96 //--- normalize the function to Zelec and store in reciprocal space:
97 double scale_factor = Zelec / KernelR.transform(0, 0.);
98 for(double& f: KernelR.f)
99 f *= scale_factor;
100 KernelR.transform(0, dG, nGridLoc, elecKernel);
101 //--- calculate deltaS that accounts for mismatched charge kernels:
102 //where deltaS = -2*pi/3*int(r^2*n(r))dV (realspace formula)
103 for(unsigned i=0; i<KernelR.r.size(); i++)
104 { double r = KernelR.r[i];
105 KernelR.f[i] = (-2*M_PI/3)*(r*r) * (KernelR.f[i] - Exponential(r,Zelec,aElec));
106 }
107 deltaS -= 8*M_PI*Zelec*pow(aElec,2); //apply the potential correction for the analytical exponential function
108 deltaS += KernelR.transform(0, 0.); //apply the potential correction for the difference between radial and analytical functions
109 }
110
111 if(elecFilename.length())
112 { //--- Read from file to RadialFunctionR:
113 ifstream ifs(elecFilename.c_str());
114 if(!ifs.is_open()) die("\tCan't open radial electron density file '%s' for reading.\n", elecFilename.c_str());
115 logPrintf(" Adding realspace radial model from '%s':\n", elecFilename.c_str());
116 std::vector<double> rVec, nVec; //to store values of radius and density
117 rVec.reserve(1000); nVec.reserve(1000); //to prevent repeated reallocations
118 while(!ifs.eof())
119 { string line; getline(ifs, line);
120 if(ifs.eof()) break;
121 double r,n; istringstream(line) >> r >> n;
122 if(rVec.size() && r<rVec.back())
123 { die("\tERROR reading electron density model for %s site from %s:\n"
124 "\tRadial gridpoints must be in ascending order\n", name.c_str(), elecFilename.c_str());
125 }
126 if(ifs.fail() || ifs.bad())
127 die("\tERROR reading electron density model for %s site from %s\n", name.c_str(), elecFilename.c_str());
128 rVec.push_back(r);
129 nVec.push_back(n);
130 }
131 RadialFunctionR KernelR(rVec.size());
132 KernelR.r = rVec;
133 KernelR.f = nVec;
134 KernelR.initWeights();
135 //--- Put kernel in G space and add to existing elecKernel:
136 RadialFunctionG KernelG;
137 KernelR.transform(0, dG, nGridLoc, KernelG);
138 std::vector<double> newKernel;
139 for (int i=0; i<nGridLoc; i++)
140 { double G = i*dG;
141 newKernel.push_back(elecKernel(G)+KernelG(G));
142 }
143 elecKernel.init(0,newKernel,dG); //reinitialize elecKernel with new radial function added on.
144 Znuc = elecKernel(0)-Zsite;
145 logPrintf(" Adjusting Znuc to %lg to ensure correct site charge.\n",Znuc);
146 //--- Update deltaS correction for new piece
147 double norm = KernelG(0.0);
148 for(unsigned i=0; i<rVec.size(); i++)
149 { double r = KernelR.r[i];
150 KernelR.f[i] = (-2*M_PI/3)*(r*r) * (KernelR.f[i] - Exponential(r,norm,aElec));
151 }
152 deltaS -= 8*M_PI*norm*pow(aElec,2); //apply the potential correction for the analytical exponential function
153 deltaS += KernelR.transform(0, 0.); //apply the potential correction for the difference between radial and analytical functions
154 }
155
156 if(elecFilenameG.length())
157 { //--- Read from file to RadialFunctionG:
158 ifstream ifs(elecFilenameG.c_str());
159 if(!ifs.is_open())
160 die("\tCan't open radial Gspace electron density file '%s' for reading.\n", elecFilenameG.c_str());
161 logPrintf(" Adding Gspace radial model from '%s':\n", elecFilenameG.c_str());
162 std::vector<double> GVec,nVec; //to store values of G and density
163 GVec.reserve(1000); nVec.reserve(1000); //to prevent repeated reallocations
164 while (!ifs.eof())
165 { string line; getline(ifs, line);
166 if(ifs.eof()) break;
167 double G,n; istringstream(line) >> G >> n;
168 if(GVec.size()==1)
169 { const double dGmin = 0.1;
170 double dGfile = G - GVec.back();
171 if(dGfile <= 0.0)
172 die("\tERROR reading electron density model for %s site from %s:\n"
173 "\tGspace gridpoints must be in ascending order\n", name.c_str(), elecFilenameG.c_str())
174 if(dGfile > dGmin)
175 die("\tERROR reading electron density model for %s site from %s:\n"
176 "\tGspace grid spacing must be less than %lg\n", name.c_str(), elecFilenameG.c_str(), dGmin)
177 }
178 else if(GVec.size()>1)
179 { if(fabs( (G-GVec.back()) - (GVec[1]-GVec[0]) )>1e-8)
180 die("\tERROR reading electron density model for %s site from %s:\n"
181 "\tNon-uniform G spacing is currently unsupported\n", name.c_str(), elecFilenameG.c_str())
182 }
183 if( ifs.fail() || ifs.bad())
184 die("\tERROR reading electron density model for %s site from %s\n", name.c_str(), elecFilenameG.c_str())
185 GVec.push_back(G);
186 nVec.push_back(n);
187 }
188 RadialFunctionG KernelG;
189 KernelG.init(0, nVec, GVec[1]-GVec[0]);
190 //--- add to existing elecKernel:
191 std::vector<double> newKernel;
192 for(int i=0; i<nGridLoc; i++)
193 { double G = i*dG;
194 newKernel.push_back(elecKernel(G)+KernelG(G));
195 }
196 elecKernel.init(0,newKernel,dG); //reinitialize elecKernel with new radial function added on.
197 logPrintf("\tWARNING: Uncompensated charge kernel mismatch for fluid-fluid and electron-fluid interactions\n");
198 }
199 }
200
201 //Initialize charge kernel:
202 if(elecKernel || Znuc)
203 { logPrintf(" Charge density: gaussian nuclear width %lg", sigmaNuc);
204 std::vector<double> samples(nGridLoc);
205 for(unsigned iG=0; iG<samples.size(); iG++)
206 { double G = iG * dG;
207 if(elecKernel) samples[iG] += elecKernel(G);
208 if(Znuc) samples[iG] -= RadialFunctionG::gaussTilde(G, Znuc, sigmaNuc);
209 }
210 chargeKernel.init(0, samples, dG);
211 logPrintf(" with net site charge %lg\n", chargeKernel(0));
212 deltaS += 2.*M_PI*Znuc*pow(sigmaNuc,2);
213 }
214
215 //Initialize polarizability kernel:
216 if(alpha)
217 { logPrintf(" Polarizability: cuspless exponential with width %lg and norm %lg\n", aPol, alpha);
218 polKernel.init(0, dG, gInfo.GmaxGrid, RadialFunctionG::cusplessExpTilde, 1., aPol);
219 }
220
221 if(Rhs)
222 { logPrintf(" Hard sphere radius: %lg bohrs\n", Rhs);
223 //Initialize hard sphere weights:
224 ErfFMTweight erfFMTweight(Rhs, 0.);
225 std::vector<double> w0(nGridLoc), w1(nGridLoc), w2(nGridLoc), w3(nGridLoc), w1v(nGridLoc), w2m(nGridLoc);
226 for(int iG=0; iG<nGridLoc; iG++)
227 erfFMTweight(iG*dG, w0[iG], w1[iG], w2[iG], w3[iG], w1v[iG], w2m[iG]);
228 this->w0.init(0, w0, dG);
229 this->w1.init(0, w1, dG);
230 this->w2.init(0, w2, dG);
231 this->w3.init(0, w3, dG);
232 this->w1v.init(0, w1v, dG);
233 this->w2m.init(0, w2m, dG);
234 }
235
236 logPrintf(" Positions in reference frame:\n");
237 for(vector3<> r: positions) logPrintf(" [ %+.6lf %+.6lf %+.6lf ]\n", r[0], r[1], r[2]);
238 initialized = true;
239 }
240
Molecule(string name)241 Molecule::Molecule(string name) : name(name), initialized(false)
242 {
243 }
244
~Molecule()245 Molecule::~Molecule()
246 { if(initialized)
247 { mfKernel.free();
248 }
249 }
250
sphericalShellTilde(double G,double R)251 inline double sphericalShellTilde(double G, double R)
252 { return bessel_jl(0, G*R);
253 }
254
255 //Fourier transform of gaussian
gaussTilde(double G,double norm,double sigma)256 inline double gaussTilde(double G, double norm, double sigma)
257 { double sigmaG = sigma*G;
258 return norm * exp(-0.5*sigmaG*sigmaG);
259 }
260
261 //Lischner10 Coulomb Kernel
setCoulombCutoffKernel(double G)262 inline double setCoulombCutoffKernel(double G)
263 {
264 const double Gc = 0.33;
265 return 1/sqrt(1 + pow(G/Gc,4));
266 }
267
setup(const GridInfo & gInfo,double Rmf)268 void Molecule::setup(const GridInfo& gInfo, double Rmf)
269 { logPrintf(" Initializing fluid molecule '%s'\n", name.c_str());
270 for(auto& site: sites) site->setup(gInfo);
271 logPrintf(" Net charge: %lg dipole magnitude: %lg\n", checkCharge(), getDipole().length());
272
273 if(getCharge())
274 { //Ions use gaussian mfKernel while neutral solvent molecules use spherical shell mfKernel
275 double Res = Rmf ? Rmf/sqrt(2) : pow(3*getVhs()/(4.*M_PI), 1./3)/sqrt(2);
276 mfKernel.init(0, gInfo.dGradial, gInfo.GmaxGrid, gaussTilde, 1.0, Res);
277 logPrintf(" Initializing gaussian mfKernel with width: %lg Bohr\n",Res);
278 for(auto& site: sites) site->deltaS += 2.*M_PI*site->Zsite*pow(Res,2);
279 }
280 else
281 { //Sundararaman style charge kernel (from JCP 2014)
282 double Res = Rmf ? Rmf : pow(3*getVhs()/(4.*M_PI), 1./3);
283 mfKernel.init(0, gInfo.dGradial, gInfo.GmaxGrid, sphericalShellTilde, Res);
284 logPrintf(" Initializing spherical shell mfKernel with radius %lg Bohr\n",Res);
285 for(auto& site: sites) site->deltaS += 2.*M_PI/3.*site->Zsite*pow(Res,2);
286 //Debugging option:
287 //Lischner style charge kernel
288 // mfKernel.init(0, gInfo.dGradial, gInfo.GmaxGrid, setCoulombCutoffKernel);
289 // logPrintf("\tInitializing Lischner style mfKernel.\n");
290 //no contribution to deltaS
291 }
292 logPrintf(" deltaS corrections:\n");
293 for(auto& site: sites)
294 logPrintf(" site '%s': %lg\n", site->name.c_str(), site->deltaS);
295 initialized = true;
296 }
297
298
isMonoatomic() const299 bool Molecule::isMonoatomic() const
300 { return (sites.size()==1) && (sites[0]->positions.size()==1);
301 }
302
getCharge() const303 double Molecule::getCharge() const
304 { double Q = 0.0;
305 for(const auto& site: sites)
306 if(site->chargeKernel)
307 Q += site->chargeKernel(0) * site->positions.size();
308 if(fabs(Q) < 1e-12) return 0.; //simplify neutrality testing
309 else return Q;
310 }
311
312 //At first initialization of molecule absorb very small net charges
313 //into first site chargeKernel and elecKernel
checkCharge()314 double Molecule::checkCharge()
315 {
316 double Q = getCharge();
317 if(Q==0.)
318 return 0.;
319 if(fabs(Q) < 1e-2)
320 {
321 int ChargedSite = 0;
322 while (!sites[ChargedSite]->chargeKernel(0))
323 {
324 ChargedSite += 1;
325 }
326 std::shared_ptr<Site> s0 = sites[ChargedSite];
327 int nG = s0->chargeKernel.nCoeff;
328 double dG = 1./s0->chargeKernel.dGinv;
329 std::vector<double> chargeKernel, elecKernel;
330 for (int i=0; i < nG; i++)
331 {
332 double G = double(i)*dG;
333 chargeKernel.push_back(s0->chargeKernel(G));
334 elecKernel.push_back(s0->elecKernel(G));
335 }
336 chargeKernel[0] -= Q / s0->positions.size();
337 s0->chargeKernel.init(0,chargeKernel,dG);
338 elecKernel[0] -= Q / s0->positions.size();
339 s0->elecKernel.init(0,elecKernel,dG);
340 logPrintf(" WARNING: Molecule had net charge %lg, adjusting %s site charge by %lg to compensate.\n",Q,s0->name.c_str(),-Q/s0->positions.size());
341 return 0.;
342 }
343 else return Q;
344 }
345
346
getDipole() const347 vector3<> Molecule::getDipole() const
348 { vector3<> P;
349 for(const auto& site: sites)
350 if(site->chargeKernel)
351 for(const vector3<>& r: site->positions)
352 P += site->chargeKernel(0) * r;
353 if(P.length() < 1e-12) return vector3<>(); //simplify polarity testing
354 else return P;
355 }
356
getVhs() const357 double Molecule::getVhs() const
358 { double Vhs = 0;
359 for(const auto& site: sites)
360 if(site->w3)
361 Vhs += site->w3(0) * site->positions.size();
362 return Vhs;
363 }
364
getAlphaTot() const365 double Molecule::getAlphaTot() const
366 { double alphaTot = 0.;
367 for(const auto& site: sites)
368 if(site->polKernel)
369 alphaTot += site->alpha * site->positions.size();
370 return alphaTot;
371 }
372
getBonds() const373 std::map<double,int> Molecule::getBonds() const
374 { std::map<double,int> bond;
375 for(const auto& site1: sites)
376 { double R1 = site1->Rhs;
377 if(R1) for(vector3<> pos1: site1->positions)
378 { for(const auto& site2: sites)
379 { double R2 = site2->Rhs;
380 if(R2) for(vector3<> pos2: site2->positions)
381 { if(fabs(R1+R2-(pos1-pos2).length()) < 1e-6*(R1+R2))
382 bond[R1*R2/(R1+R2)]++;
383 }
384 }
385 }
386 }
387 //correct for double counting:
388 for(auto& bondEntry: bond)
389 { assert(bondEntry.second % 2 == 0);
390 bondEntry.second /= 2;
391 }
392 return bond;
393 }
394
setModelMonoatomic(string name,double Q,double Rhs)395 void Molecule::setModelMonoatomic(string name, double Q, double Rhs)
396 { sites.clear();
397 this->name = name;
398 auto site = std::make_shared<Molecule::Site>(name);
399 site->Znuc = Q; site->sigmaNuc = (1./6)*Rhs;
400 site->Rhs = Rhs;
401 sites.push_back(site);
402 site->positions.push_back(vector3<>());
403 }
404