JDFTx  1.0.0
PCM_internal.h
1 /*-------------------------------------------------------------------
2 Copyright 2011 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 #ifndef JDFTX_ELECTRONIC_PCM_INTERNAL_H
21 #define JDFTX_ELECTRONIC_PCM_INTERNAL_H
22 
23 #include <core/vector3.h>
24 #include <electronic/RadialFunction.h>
25 
26 //----------- Common PCM functions (top level interface not seen by .cu files) ------------
27 #ifndef __in_a_cu_file__
28 
29 #include <core/Operators.h>
30 #include <core/EnergyComponents.h>
31 #include <fluid/FluidSolverParams.h>
32 
33 namespace ShapeFunction
34 {
36  void compute(const ScalarField& n, ScalarField& shape, double nc, double sigma);
37 
39  void propagateGradient(const ScalarField& n, const ScalarField& E_shape, ScalarField& E_n, double nc, double sigma);
40 }
41 
42 namespace ShapeFunctionCANDLE
43 {
45  void compute(const ScalarField& n, const ScalarFieldTilde& phi,
46  ScalarField& shape, double nc, double sigma, double pCavity);
47 
49  void propagateGradient(const ScalarField& n, const ScalarFieldTilde& phi, const ScalarField& E_shape,
50  ScalarField& E_n, ScalarFieldTilde& E_phi, double& E_pCavity, double nc, double sigma, double pCavity);
51 }
52 
53 namespace ShapeFunctionSGA13
54 {
56  void expandDensity(const RadialFunctionG& w, double R, const ScalarField& n, ScalarField& nEx, const ScalarField* A_nEx=0, ScalarField* A_n=0);
57 }
58 
59 namespace ShapeFunctionSCCS
60 {
62  void compute(const ScalarField& n, ScalarField& shape, double rhoMin, double rhoMax, double epsBulk);
63 
65  void propagateGradient(const ScalarField& n, const ScalarField& E_shape, ScalarField& E_n, double rhoMin, double rhoMax, double epsBulk);
66 }
67 
68 #endif
69 
70 
71 //--------- Compute kernels (shared by CPU and GPU implementations) --------
72 
73 //Cavity shape function and gradient
74 namespace ShapeFunction
75 {
76  __hostanddev__ void compute_calc(int i, const double* nCavity, double* shape, const double nc, const double sigma)
77  { shape[i] = erfc(sqrt(0.5)*log(fabs(nCavity[i])/nc)/sigma)*0.5;
78  }
79  __hostanddev__ void propagateGradient_calc(int i, const double* nCavity, const double* grad_shape, double* grad_nCavity, const double nc, const double sigma)
80  { grad_nCavity[i] += (-1.0/(nc*sigma*sqrt(2*M_PI))) * grad_shape[i]
81  * exp(0.5*(pow(sigma,2) - pow(log(fabs(nCavity[i])/nc)/sigma + sigma, 2)));
82  }
83 }
84 
85 namespace ShapeFunctionCANDLE
86 {
87  //version with electric-field based charge asymmetry (combined compute and grad function)
88  __hostanddev__ void compute_or_grad_calc(int i, bool grad,
89  const double* nArr, vector3<const double*> DnArr, vector3<const double*> DphiArr, double* shape,
90  const double* A_shape, double* A_n, vector3<double*> A_Dn, vector3<double*> A_Dphi, double* A_pCavity,
91  const double nc, const double invSigmaSqrt2, const double pCavity)
92  { double n = nArr[i];
93  if(n<1e-8) { if(!grad) shape[i]=1.; return; }
94  //Regularized unit vector along Dn:
95  vector3<> Dn = loadVector(DnArr,i);
96  double normFac = 1./sqrt(Dn.length_squared() + 1e-4*nc*nc);
97  vector3<> e = Dn * normFac;
98  //Electric field along above unit vector, with saturation for stability:
99  vector3<> E = -loadVector(DphiArr,i);
100  double eDotE = dot(e,E);
101  const double x_eDotE = -fabs(pCavity);
102  double x = x_eDotE * eDotE;
103  double asymm=0., asymm_x=0.;
104  //modify cavity only for anion-like regions
105  if(x > 4.) { asymm = 1.; asymm_x = 0.; } //avoid Inf/Inf error
106  else if(x > 0.)
107  { double exp2x2 = exp(2.*x*x), den = 1./(1 + exp2x2);
108  asymm = (exp2x2 - 1.) * den; //tanh(x^2)
109  asymm_x = 8.*x * exp2x2 * den*den; //2x sech(x^2)
110  }
111  const double dlognMax = copysign(3., pCavity);
112  double comb = log(n/nc) - dlognMax*asymm;
113  if(!grad)
114  shape[i] = 0.5*erfc(invSigmaSqrt2*comb);
115  else
116  { double A_comb = (-invSigmaSqrt2/sqrt(M_PI)) * A_shape[i] * exp(-comb*comb*invSigmaSqrt2*invSigmaSqrt2);
117  A_n[i] += A_comb/n;
118  double A_x = A_comb*(-dlognMax)*asymm_x;
119  accumVector((A_x*x_eDotE*normFac) * (E - e*eDotE), A_Dn,i);
120  accumVector((A_x*x_eDotE*(-1.)) * e, A_Dphi,i);
121  A_pCavity[i] += A_x*(-copysign(1.,pCavity))*eDotE;
122  }
123  }
124 }
125 
126 namespace ShapeFunctionSGA13
127 {
128  __hostanddev__ void expandDensity_calc(int i, double alpha, const double* nBar, const double* DnBarSq, double* nEx, double* nEx_nBar, double* nEx_DnBarSq)
129  { double n = nBar[i], D2 = DnBarSq[i];
130  if(n < 1e-9) //Avoid numerical error in low density / gradient regions:
131  { nEx[i]=1e-9;
132  if(nEx_nBar) nEx_nBar[i]=0.;
133  if(nEx_DnBarSq) nEx_DnBarSq[i]=0.;
134  return;
135  }
136  double nInv = 1./n;
137  nEx[i] = alpha*n + D2*nInv;
138  if(nEx_nBar) { nEx_nBar[i] = alpha - D2*nInv*nInv; }
139  if(nEx_DnBarSq) { nEx_DnBarSq[i] = nInv; }
140  }
141 }
142 
143 //Cavity shape function and gradient for the SCCS models
144 namespace ShapeFunctionSCCS
145 {
146  __hostanddev__ void compute_calc(int i, const double* nCavity, double* shape,
147  const double rhoMin, const double rhoMax, const double epsBulk)
148  { double rho = nCavity[i];
149  if(rho >= rhoMax) { shape[i] = 0.; return; }
150  if(rho <= rhoMin) { shape[i] = 1.; return; }
151  const double logDen = log(rhoMax/rhoMin);
152  double f = log(rhoMax/rho)/logDen;
153  double t = f - sin(2*M_PI*f)/(2*M_PI);
154  shape[i] = (pow(epsBulk,t) - 1.)/(epsBulk - 1.);
155  }
156  __hostanddev__ void propagateGradient_calc(int i, const double* nCavity, const double* grad_shape, double* grad_nCavity,
157  const double rhoMin, const double rhoMax, const double epsBulk)
158  { double rho = nCavity[i];
159  if(rho >= rhoMax) return;
160  if(rho <= rhoMin) return;
161  const double logDen = log(rhoMax/rhoMin);
162  double f = log(rhoMax/rho)/logDen;
163  double f_rho = -1./(rho*logDen); //df/drho
164  double t = f - sin(2*M_PI*f)/(2*M_PI);
165  double t_f = 1. - cos(2*M_PI*f); //dt/df
166  double s_t = log(epsBulk) * pow(epsBulk,t)/(epsBulk - 1.); //dshape/dt
167  grad_nCavity[i] += grad_shape[i] * s_t * t_f * f_rho; //chain rule
168  }
169 }
170 
171 //------------- Helper classes for NonlinearPCM -------------
172 namespace NonlinearPCMeval
173 {
175  struct Screening
176  {
177  bool linear;
178  double NT, ZbyT, NZ;
179  double x0plus, x0minus, x0;
180 
181  Screening(bool linear, double T, double Nion, double Zion, double VhsPlus, double VhsMinus, double epsBulk); //epsBulk is used only for printing screening length
182 
183  #ifndef __in_a_cu_file__
184  inline double neutralityConstraint(const ScalarField& muPlus, const ScalarField& muMinus, const ScalarField& shape, double Qexp,
186  ScalarField* mu0_muPlus=0, ScalarField* mu0_muMinus=0, ScalarField* mu0_shape=0, double* mu0_Qexp=0)
187  {
188  if(linear)
189  { double Qsum = NZ * 2.*integral(shape);
190  double Qdiff = NZ * integral(shape*(muPlus+muMinus));
191  //Compute the constraint function and its derivatives w.r.t above moments:
192  double mu0 = -(Qexp + Qdiff)/Qsum;
193  double mu0_Qdiff = -1./Qsum;
194  double mu0_Qsum = (Qexp + Qdiff)/(Qsum*Qsum);
195  //Collect reuslt and optional gradients:
196  if(mu0_muPlus) *mu0_muPlus = (mu0_Qdiff * NZ) * shape;
197  if(mu0_muMinus) *mu0_muMinus = (mu0_Qdiff * NZ) * shape;
198  if(mu0_shape) *mu0_shape = NZ * (mu0_Qdiff*(muPlus+muMinus) + mu0_Qsum*2.);
199  if(mu0_Qexp) *mu0_Qexp = -1./Qsum;
200  return mu0;
201  }
202  else
203  { ScalarField etaPlus = exp(muPlus);
204  ScalarField etaMinus = exp(-muMinus);
205  double Qplus = +NZ * integral(shape * etaPlus);
206  double Qminus = -NZ * integral(shape * etaMinus);
207  //Compute the constraint function and its derivatives w.r.t above moments:
208  double mu0, mu0_Qplus, mu0_Qminus;
209  double disc = sqrt(Qexp*Qexp - 4.*Qplus*Qminus); //discriminant for quadratic
210  //Pick the numerically stable path (to avoid roundoff problems when |Qplus*Qminus| << Qexp^2):
211  if(Qexp<0)
212  { mu0 = log((disc-Qexp)/(2.*Qplus));
213  mu0_Qplus = -2.*Qminus/(disc*(disc-Qexp)) - 1./Qplus;
214  mu0_Qminus = -2.*Qplus/(disc*(disc-Qexp));
215  }
216  else
217  { mu0 = log(-2.*Qminus/(disc+Qexp));
218  mu0_Qplus = 2.*Qminus/(disc*(disc+Qexp));
219  mu0_Qminus = 2.*Qplus/(disc*(disc+Qexp)) + 1./Qminus;
220  }
221  //Collect result and optional gradients:
222  if(mu0_muPlus) *mu0_muPlus = (mu0_Qplus * NZ) * shape * etaPlus;
223  if(mu0_muMinus) *mu0_muMinus = (mu0_Qminus * NZ) * shape * etaMinus;
224  if(mu0_shape) *mu0_shape = NZ * (mu0_Qplus * etaPlus - mu0_Qminus * etaMinus);
225  if(mu0_Qexp) *mu0_Qexp = -1./disc;
226  return mu0;
227  }
228  }
229  #endif
230 
232  __hostanddev__ double fHS(double x, double& f_x) const
233  { if(x>=1.) { f_x = NAN; return NAN; }
234  double den = 1./(1-x), den0 = 1./(1-x0);
235  double comb = (x-x0)*den*den0, comb_x = den*den;
236  double prefac = (2./x0);
237  double f = prefac * comb*comb;
238  f_x = prefac * 2.*comb*comb_x;
239  return f;
240  }
241 
244  __hostanddev__ void compute(double muPlus, double muMinus, double& F, double& F_muPlus, double& F_muMinus, double& Rho, double& Rho_muPlus, double& Rho_muMinus) const
245  { if(linear)
246  { F = NT * 0.5*(muPlus*muPlus + muMinus*muMinus);
247  F_muPlus = NT * muPlus;
248  F_muMinus = NT * muMinus;
249  Rho = NZ * (muPlus + muMinus);
250  Rho_muPlus = NZ;
251  Rho_muMinus = NZ;
252  }
253  else
254  { double etaPlus = exp(muPlus), etaMinus=exp(-muMinus);
255  double x = x0plus*etaPlus + x0minus*etaMinus; //packing fraction
256  double f_x, f = fHS(x, f_x); //hard sphere free energy per particle
257  F = NT * (2. + etaPlus*(muPlus-1.) + etaMinus*(-muMinus-1.) + f);
258  F_muPlus = NT * etaPlus *(muPlus + f_x * x0plus);
259  F_muMinus = NT * etaMinus*(muMinus - f_x * x0minus);
260  Rho = NZ * (etaPlus - etaMinus);
261  Rho_muPlus = NZ * etaPlus;
262  Rho_muMinus = NZ * etaMinus;
263  }
264  }
265 
267  __hostanddev__ void freeEnergy_calc(size_t i, double mu0, const double* muPlus, const double* muMinus, const double* s, double* rho, double* A, double* A_muPlus, double* A_muMinus, double* A_s) const
268  { double F, F_muPlus, F_muMinus, Rho, Rho_muPlus, Rho_muMinus;
269  compute(muPlus[i]+mu0, muMinus[i]+mu0, F, F_muPlus, F_muMinus, Rho, Rho_muPlus, Rho_muMinus);
270  A[i] = s[i] * F;
271  A_muPlus[i] += s[i] * F_muPlus;
272  A_muMinus[i] += s[i] * F_muMinus;
273  if(A_s) A_s[i] += F;
274  rho[i] = s[i] * Rho;
275  }
276  void freeEnergy(size_t N, double mu0, const double* muPlus, const double* muMinus, const double* s, double* rho, double* A, double* A_muPlus, double* A_muMinus, double* A_s) const;
277  #ifdef GPU_ENABLED
278  void freeEnergy_gpu(size_t N, double mu0, const double* muPlus, const double* muMinus, const double* s, double* rho, double* A, double* A_muPlus, double* A_muMinus, double* A_s) const;
279  #endif
280 
282  __hostanddev__ void convertDerivative_calc(size_t i, double mu0, const double* muPlus, const double* muMinus, const double* s, const double* A_rho, double* A_muPlus, double* A_muMinus, double* A_s) const
283  { double F, F_muPlus, F_muMinus, Rho, Rho_muPlus, Rho_muMinus;
284  compute(muPlus[i]+mu0, muMinus[i]+mu0, F, F_muPlus, F_muMinus, Rho, Rho_muPlus, Rho_muMinus);
285  A_muPlus[i] += s[i] * Rho_muPlus * A_rho[i];
286  A_muMinus[i] += s[i] * Rho_muMinus * A_rho[i];
287  if(A_s) A_s[i] += Rho * A_rho[i];
288  }
289  void convertDerivative(size_t N, double mu0, const double* muPlus, const double* muMinus, const double* s, const double* A_rho, double* A_muPlus, double* A_muMinus, double* A_s) const;
290  #ifdef GPU_ENABLED
291  void convertDerivative_gpu(size_t N, double mu0, const double* muPlus, const double* muMinus, const double* s, const double* A_rho, double* A_muPlus, double* A_muMinus, double* A_s) const;
292  #endif
293 
295  __hostanddev__ double rootFunc(double x, double V) const
296  { double f_x; fHS(x, f_x); //hard sphere potential
297  return x - (x0plus*exp(-V-f_x*x0plus) + x0minus*exp(+V-f_x*x0minus));
298  }
299 
301  __hostanddev__ double x_from_V(double V) const
302  { double xLo = x0; while(rootFunc(xLo, V) > 0.) xLo *= 0.5;
303  double xHi = xLo; while(rootFunc(xHi, V) < 0.) xHi = 0.5*(xHi + 1.);
304  double x = 0.5*(xHi+xLo);
305  double dx = x*1e-13;
306  while(xHi-xLo > dx)
307  { if(rootFunc(x, V) < 0.)
308  xLo = x;
309  else
310  xHi = x;
311  x = 0.5*(xHi+xLo);
312  }
313  return x;
314  }
315 
317  __hostanddev__ void phiToState_calc(size_t i, const double* phi, const double* s, const RadialFunctionG& xLookup, bool setState, double* muPlus, double* muMinus, double* kappaSq) const
318  { double V = ZbyT * phi[i];
319  if(!setState)
320  { //Avoid V=0 in calculating kappaSq below
321  if(fabs(V) < -1e-7)
322  V = copysign(1e-7, V);
323  }
324  double twoCbrtV= 2.*pow(fabs(V), 1./3);
325  double Vmapped = copysign(twoCbrtV / (1. + sqrt(1. + twoCbrtV*twoCbrtV)), V);
326  double x = 1. - xLookup(1.+Vmapped);
327  double f_x; fHS(x, f_x); //hard sphere potential
328  double logEtaPlus = -V - f_x*x0plus;
329  double logEtaMinus = +V - f_x*x0minus;
330  if(setState)
331  { muPlus[i] = logEtaPlus;
332  muMinus[i] = -logEtaMinus;
333  }
334  else
335  kappaSq[i] = (4*M_PI)*s[i]*(NZ*ZbyT)*(exp(logEtaMinus) - exp(logEtaPlus))/V;
336  }
337  void phiToState(size_t N, const double* phi, const double* s, const RadialFunctionG& xLookup, bool setState, double* muPlus, double* muMinus, double* kappaSq) const;
338  #ifdef GPU_ENABLED
339  void phiToState_gpu(size_t N, const double* phi, const double* s, const RadialFunctionG& xLookup, bool setState, double* muPlus, double* muMinus, double* kappaSq) const;
340  #endif
341 
342  };
343 
345  struct Dielectric
346  {
347  bool linear;
348  double Np, pByT, NT;
349  double alpha, X;
350 
351  Dielectric(bool linear, double T, double Nmol, double pMol, double epsBulk, double epsInf);
352 
354  __hostanddev__ void calcFunctions(double eps, double& frac, double& frac_epsSqHlf, double& logsinch) const
355  { double epsSq = eps*eps;
356  if(linear)
357  { frac = 1.0/3;
358  frac_epsSqHlf = 0.;
359  logsinch = epsSq*(1.0/6);
360  }
361  else
362  { if(eps < 1e-1) //Use series expansions
363  { frac = 1.0/3 + epsSq*(-1.0/45 + epsSq*(2.0/945 + epsSq*(-1.0/4725)));
364  frac_epsSqHlf = -2.0/45 + epsSq*(8.0/945 + epsSq*(-6.0/4725));
365  logsinch = epsSq*(1.0/6 + epsSq*(-1.0/180 + epsSq*(1.0/2835)));
366  }
367  else
368  { frac = (eps/tanh(eps)-1)/epsSq;
369  frac_epsSqHlf = (2 - eps/tanh(eps) - pow(eps/sinh(eps),2))/(epsSq*epsSq);
370  logsinch = eps<20. ? log(sinh(eps)/eps) : eps - log(2.*eps);
371  }
372  }
373  }
374 
376  __hostanddev__ void compute(double epsSqHlf, double& F, double& F_epsSqHlf, double& ChiEff, double& ChiEff_epsSqHlf) const
377  { double epsSq = 2.*epsSqHlf, eps = sqrt(epsSq);
378  //----- Nonlinear functions of eps
379  double frac, frac_epsSqHlf, logsinch;
380  calcFunctions(eps, frac, frac_epsSqHlf, logsinch);
381  //----- Free energy and derivative
382  double screen = 1 - alpha*frac; //correlation screening factor = (pE/T) / eps where E is the real electric field
383  F = NT * (epsSq*(frac - 0.5*alpha*frac*frac + 0.5*X*screen*screen) - logsinch);
384  F_epsSqHlf = NT * (frac + X*screen + epsSq*frac_epsSqHlf*(1.-X*alpha)) * screen; //Simplified using logsinch_epsSqHlf = frac
385  //----- Effective susceptibility and derivative
386  ChiEff = Np * (frac + X*screen);
387  ChiEff_epsSqHlf = Np * frac_epsSqHlf * (1.-X*alpha);
388  }
389 
391  __hostanddev__ void freeEnergy_calc(size_t i, vector3<const double*> eps, const double* s, vector3<double*> p, double* A, vector3<double*> A_eps, double* A_s) const
392  { vector3<> epsVec = loadVector(eps, i);
393  double epsSqHlf = 0.5*epsVec.length_squared();
394  double F, F_epsSqHlf, ChiEff, ChiEff_epsSqHlf;
395  compute(epsSqHlf, F, F_epsSqHlf, ChiEff, ChiEff_epsSqHlf);
396  A[i] = F * s[i];
397  accumVector((F_epsSqHlf * s[i]) * epsVec, A_eps, i);
398  if(A_s) A_s[i] += F;
399  storeVector((ChiEff * s[i]) * epsVec, p, i);
400  }
401  void freeEnergy(size_t N, vector3<const double*> eps, const double* s, vector3<double*> p, double* A, vector3<double*> A_eps, double* A_s) const;
402  #ifdef GPU_ENABLED
403  void freeEnergy_gpu(size_t N, vector3<const double*> eps, const double* s, vector3<double*> p, double* A, vector3<double*> A_eps, double* A_s) const;
404  #endif
405 
407  __hostanddev__ void convertDerivative_calc(size_t i, vector3<const double*> eps, const double* s, vector3<const double*> A_p, vector3<double*> A_eps, double* A_s) const
408  { vector3<> epsVec = loadVector(eps, i);
409  double epsSqHlf = 0.5*epsVec.length_squared();
410  double F, F_epsSqHlf, ChiEff, ChiEff_epsSqHlf;
411  compute(epsSqHlf, F, F_epsSqHlf, ChiEff, ChiEff_epsSqHlf);
412  //Propagate derivatives:
413  vector3<> A_pVec = loadVector(A_p, i);
414  accumVector(s[i]*( ChiEff*A_pVec + ChiEff_epsSqHlf*epsVec*dot(A_pVec, epsVec)), A_eps, i);
415  if(A_s) A_s[i] += ChiEff * dot(A_pVec, epsVec);
416  }
417  void convertDerivative(size_t N, vector3<const double*> eps, const double* s, vector3<const double*> A_p, vector3<double*> A_eps, double* A_s) const;
418  #ifdef GPU_ENABLED
419  void convertDerivative_gpu(size_t N, vector3<const double*> eps, const double* s, vector3<const double*> A_p, vector3<double*> A_eps, double* A_s) const;
420  #endif
421 
423  __hostanddev__ double x_from_eps(double eps) const
424  { double frac, frac_epsSqHlf, logsinch;
425  calcFunctions(eps, frac, frac_epsSqHlf, logsinch);
426  return eps*(1. - alpha*frac);
427  }
428 
430  __hostanddev__ double eps_from_x(double x) const
431  { if(!x) return 0.;
432  double epsLo = x; while(x_from_eps(epsLo) > x) epsLo *= 0.95;
433  double epsHi = epsLo; while(x_from_eps(epsHi) < x) epsHi *= 1.05;
434  double eps = 0.5*(epsHi+epsLo);
435  double deps = eps*1e-13;
436  while(epsHi-epsLo > deps)
437  { if(x_from_eps(eps) < x)
438  epsLo = eps;
439  else
440  epsHi = eps;
441  eps = 0.5*(epsHi+epsLo);
442  }
443  return eps;
444  }
445 
447  __hostanddev__ void phiToState_calc(size_t i, vector3<const double*> Dphi, const double* s, const RadialFunctionG& gLookup, bool setState, vector3<double*> eps, double* epsilon) const
448  { vector3<> xVec = -pByT * loadVector(Dphi, i);
449  double x = xVec.length();
450  double g = gLookup(x/(1.+x));
451  if(setState)
452  storeVector(g * xVec, eps,i);
453  else
454  epsilon[i] = 1. + (4*M_PI)*s[i]*(Np*pByT)*((g-1.)/alpha + X);
455  }
456  void phiToState(size_t N, vector3<const double*> Dphi, const double* s, const RadialFunctionG& gLookup, bool setState, vector3<double*> eps, double* epsilon) const;
457  #ifdef GPU_ENABLED
458  void phiToState_gpu(size_t N, vector3<const double*> Dphi, const double* s, const RadialFunctionG& gLookup, bool setState, vector3<double*> eps, double* epsilon) const;
459  #endif
460  };
461 }
462 
464 #define FLUID_DUMP(object, suffix) \
465  filename = filenamePattern; \
466  filename.replace(filename.find("%s"), 2, suffix); \
467  logPrintf("Dumping '%s'... ", filename.c_str()); logFlush(); \
468  if(mpiUtil->isHead()) saveRawBinary(object, filename.c_str()); \
469  logPrintf("done.\n"); logFlush();
470 
471 #endif // JDFTX_ELECTRONIC_PCM_INTERNAL_H
NonlinearPCMeval::Screening::rootFunc
__hostanddev__ double rootFunc(double x, double V) const
Root function used for finding packing fraction x at a given dimensionless potential V = Z phi / T...
Definition: PCM_internal.h:295
pow
ScalarField pow(const ScalarField &, double alpha)
Elementwise power (preserve input)
NonlinearPCMeval::Dielectric::NT
double NT
N*p, p/T and N*T where N is molecular density and p is molecular dipole.
Definition: PCM_internal.h:348
NonlinearPCMeval::Dielectric::compute
__hostanddev__ void compute(double epsSqHlf, double &F, double &F_epsSqHlf, double &ChiEff, double &ChiEff_epsSqHlf) const
Compute the nonlinear functions in the free energy and effective susceptibility (p/eps) prior to scal...
Definition: PCM_internal.h:376
ShapeFunctionSCCS
Definition: PCM_internal.h:59
sqrt
ScalarField sqrt(const ScalarField &)
Elementwise square root (preserve input)
NonlinearPCMeval::Dielectric
Helper class for dielectric portion of NonlinearPCM.
Definition: PCM_internal.h:345
ShapeFunction
Definition: PCM_internal.h:33
NonlinearPCMeval::Dielectric::freeEnergy_calc
__hostanddev__ void freeEnergy_calc(size_t i, vector3< const double * > eps, const double *s, vector3< double * > p, double *A, vector3< double * > A_eps, double *A_s) const
Given shape function s and gradient of phi eps, compute polarization p, free energy density A and acc...
Definition: PCM_internal.h:391
dot
complex dot(const Tptr &X, const Tptr &Y)
Definition: Operators.h:196
NonlinearPCMeval::Screening::phiToState_calc
__hostanddev__ void phiToState_calc(size_t i, const double *phi, const double *s, const RadialFunctionG &xLookup, bool setState, double *muPlus, double *muMinus, double *kappaSq) const
Given shape function s and phi, calculate state mu&#39;s if setState=true or effective kappaSq if setStat...
Definition: PCM_internal.h:317
integral
double integral(const ScalarField &)
Integral in the unit cell (dV times sum())
NonlinearPCMeval::Dielectric::calcFunctions
__hostanddev__ void calcFunctions(double eps, double &frac, double &frac_epsSqHlf, double &logsinch) const
Calculate the various nonlinear functions of epsilon used in calculating the free energy and its deri...
Definition: PCM_internal.h:354
ScalarFieldTilde
std::shared_ptr< ScalarFieldTildeData > ScalarFieldTilde
A smart reference-counting pointer to ScalarFieldTildeData.
Definition: ScalarField.h:44
NonlinearPCMeval::Dielectric::phiToState_calc
__hostanddev__ void phiToState_calc(size_t i, vector3< const double * > Dphi, const double *s, const RadialFunctionG &gLookup, bool setState, vector3< double * > eps, double *epsilon) const
Given shape function s and gradient of phi Dphi, calculate state vector eps if setState=true or effec...
Definition: PCM_internal.h:447
NonlinearPCMeval::Screening::x_from_V
__hostanddev__ double x_from_V(double V) const
Calculate self-consistent packing fraction x at given dimensionless potential V = Z phi / T using a b...
Definition: PCM_internal.h:301
NonlinearPCMeval::Screening::linear
bool linear
whether ionic screening is linearized
Definition: PCM_internal.h:177
NonlinearPCMeval::Dielectric::eps_from_x
__hostanddev__ double eps_from_x(double x) const
Invert x_from_eps() using a bisection method. Note that x must be positive and finite.
Definition: PCM_internal.h:430
NonlinearPCMeval::Screening
Helper class for ionic screening portion of NonlinearPCM.
Definition: PCM_internal.h:175
storeVector
__hostanddev__ void storeVector(const vector3< scalar > &v, vector3< scalar * > &vArr, int i)
Store vector to a vector field.
Definition: vector3.h:171
accumVector
__hostanddev__ void accumVector(const vector3< scalar > &v, vector3< scalar * > &vArr, int i)
Accumulate vector onto a vector field.
Definition: vector3.h:175
log
ScalarField log(const ScalarField &)
Elementwise logarithm (preserve input)
ShapeFunctionSGA13
Definition: PCM_internal.h:53
NonlinearPCMeval::Dielectric::X
double X
dipole correlation factor and chi*T/p^2 where chi is the molecular susceptibility ...
Definition: PCM_internal.h:349
NonlinearPCMeval::Dielectric::x_from_eps
__hostanddev__ double x_from_eps(double eps) const
Calculate x = pMol E / T given eps.
Definition: PCM_internal.h:423
NonlinearPCMeval::Dielectric::linear
bool linear
whether dielectric is linearized
Definition: PCM_internal.h:347
loadVector
__hostanddev__ vector3< scalar > loadVector(const vector3< const scalar * > &vArr, int i)
Load vector from a constant vector field.
Definition: vector3.h:163
NonlinearPCMeval::Screening::NZ
double NZ
where T=temperature, N=bulk ionic concentration, Z=charge (all assumed +/- symmetric) ...
Definition: PCM_internal.h:178
NonlinearPCMeval::Screening::fHS
__hostanddev__ double fHS(double x, double &f_x) const
Hard sphere free energy per particle and derivative, where x is total packing fraction.
Definition: PCM_internal.h:232
RadialFunctionG
G-space radial function stored on a uniform grid (of |G|)
Definition: RadialFunction.h:28
exp
ScalarField exp(const ScalarField &)
Elementwise exponential (preserve input)
NonlinearPCMeval::Screening::compute
__hostanddev__ void compute(double muPlus, double muMinus, double &F, double &F_muPlus, double &F_muMinus, double &Rho, double &Rho_muPlus, double &Rho_muMinus) const
Definition: PCM_internal.h:244
NonlinearPCMeval
Definition: NonlinearPCM.h:29
ScalarField
std::shared_ptr< ScalarFieldData > ScalarField
A smart reference-counting pointer to ScalarFieldData.
Definition: ScalarField.h:40
NonlinearPCMeval::Screening::freeEnergy_calc
__hostanddev__ void freeEnergy_calc(size_t i, double mu0, const double *muPlus, const double *muMinus, const double *s, double *rho, double *A, double *A_muPlus, double *A_muMinus, double *A_s) const
Given shape function s and potential mu, compute induced charge rho, free energy density A and accumu...
Definition: PCM_internal.h:267
EnergyComponents.h
Represent components of the (free) energy.
NonlinearPCMeval::Screening::convertDerivative_calc
__hostanddev__ void convertDerivative_calc(size_t i, double mu0, const double *muPlus, const double *muMinus, const double *s, const double *A_rho, double *A_muPlus, double *A_muMinus, double *A_s) const
Propagate derivative A_rho and accumulate to A_mu and A_s.
Definition: PCM_internal.h:282
NonlinearPCMeval::Screening::x0
double x0
anion, cation and total packing fractions
Definition: PCM_internal.h:179
vector3
Generic 3-vector.
Definition: vector3.h:33
ShapeFunctionCANDLE
Definition: PCM_internal.h:42
NonlinearPCMeval::Dielectric::convertDerivative_calc
__hostanddev__ void convertDerivative_calc(size_t i, vector3< const double * > eps, const double *s, vector3< const double * > A_p, vector3< double * > A_eps, double *A_s) const
Propagate derivative A_p and accumulate to A_eps and A_s.
Definition: PCM_internal.h:407
Operators.h
Operators on ScalarField&#39;s and ScalarFieldTilde&#39;s.