20 #ifndef JDFTX_ELECTRONIC_PCM_INTERNAL_H 21 #define JDFTX_ELECTRONIC_PCM_INTERNAL_H 24 #include <electronic/RadialFunction.h> 27 #ifndef __in_a_cu_file__ 46 ScalarField& shape,
double nc,
double sigma,
double pCavity);
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;
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)));
88 __hostanddev__
void compute_or_grad_calc(
int i,
bool grad,
91 const double nc,
const double invSigmaSqrt2,
const double pCavity)
93 if(n<1e-8) {
if(!grad) shape[i]=1.;
return; }
96 double normFac = 1./
sqrt(Dn.length_squared() + 1e-4*nc*nc);
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.;
105 if(x > 4.) { asymm = 1.; asymm_x = 0.; }
107 {
double exp2x2 =
exp(2.*x*x), den = 1./(1 + exp2x2);
108 asymm = (exp2x2 - 1.) * den;
109 asymm_x = 8.*x * exp2x2 * den*den;
111 const double dlognMax = copysign(3., pCavity);
112 double comb =
log(n/nc) - dlognMax*asymm;
114 shape[i] = 0.5*erfc(invSigmaSqrt2*comb);
116 {
double A_comb = (-invSigmaSqrt2/
sqrt(M_PI)) * A_shape[i] *
exp(-comb*comb*invSigmaSqrt2*invSigmaSqrt2);
118 double A_x = A_comb*(-dlognMax)*asymm_x;
119 accumVector((A_x*x_eDotE*normFac) * (E - e*eDotE), A_Dn,i);
121 A_pCavity[i] += A_x*(-copysign(1.,pCavity))*eDotE;
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];
132 if(nEx_nBar) nEx_nBar[i]=0.;
133 if(nEx_DnBarSq) nEx_DnBarSq[i]=0.;
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; }
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.);
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);
164 double t = f - sin(2*M_PI*f)/(2*M_PI);
165 double t_f = 1. - cos(2*M_PI*f);
166 double s_t =
log(epsBulk) *
pow(epsBulk,t)/(epsBulk - 1.);
167 grad_nCavity[i] += grad_shape[i] * s_t * t_f * f_rho;
179 double x0plus, x0minus,
x0;
181 Screening(
bool linear,
double T,
double Nion,
double Zion,
double VhsPlus,
double VhsMinus,
double epsBulk);
183 #ifndef __in_a_cu_file__ 189 {
double Qsum = NZ * 2.*
integral(shape);
190 double Qdiff = NZ *
integral(shape*(muPlus+muMinus));
192 double mu0 = -(Qexp + Qdiff)/Qsum;
193 double mu0_Qdiff = -1./Qsum;
194 double mu0_Qsum = (Qexp + Qdiff)/(Qsum*Qsum);
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;
205 double Qplus = +NZ *
integral(shape * etaPlus);
206 double Qminus = -NZ *
integral(shape * etaMinus);
208 double mu0, mu0_Qplus, mu0_Qminus;
209 double disc =
sqrt(Qexp*Qexp - 4.*Qplus*Qminus);
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));
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;
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;
232 __hostanddev__
double fHS(
double xIn,
double& f_xIn)
const 233 {
double x = xIn, x_xIn = 1.;
235 {
double xInInv = 1./xIn;
237 x_xIn = 0.25*xInInv*xInInv;
239 double den = 1./(1-x), den0 = 1./(1-x0);
240 double comb = (x-x0)*den*den0, comb_x = den*den;
241 double prefac = (2./x0);
242 double f = prefac * comb*comb;
243 f_xIn = prefac * 2.*comb*comb_x * x_xIn;
249 __hostanddev__
void compute(
double muPlus,
double muMinus,
double& F,
double& F_muPlus,
double& F_muMinus,
double& Rho,
double& Rho_muPlus,
double& Rho_muMinus)
const 251 { F = NT * 0.5*(muPlus*muPlus + muMinus*muMinus);
252 F_muPlus = NT * muPlus;
253 F_muMinus = NT * muMinus;
254 Rho = NZ * (muPlus + muMinus);
259 {
double etaPlus =
exp(muPlus), etaMinus=
exp(-muMinus);
260 double x = x0plus*etaPlus + x0minus*etaMinus;
261 double f_x, f = fHS(x, f_x);
262 F = NT * (2. + etaPlus*(muPlus-1.) + etaMinus*(-muMinus-1.) + f);
263 F_muPlus = NT * etaPlus *(muPlus + f_x * x0plus);
264 F_muMinus = NT * etaMinus*(muMinus - f_x * x0minus);
265 Rho = NZ * (etaPlus - etaMinus);
266 Rho_muPlus = NZ * etaPlus;
267 Rho_muMinus = NZ * etaMinus;
272 __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 273 {
double F, F_muPlus, F_muMinus, Rho, Rho_muPlus, Rho_muMinus;
274 compute(muPlus[i]+mu0, muMinus[i]+mu0, F, F_muPlus, F_muMinus, Rho, Rho_muPlus, Rho_muMinus);
276 A_muPlus[i] += s[i] * F_muPlus;
277 A_muMinus[i] += s[i] * F_muMinus;
281 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;
283 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;
287 __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 288 {
double F, F_muPlus, F_muMinus, Rho, Rho_muPlus, Rho_muMinus;
289 compute(muPlus[i]+mu0, muMinus[i]+mu0, F, F_muPlus, F_muMinus, Rho, Rho_muPlus, Rho_muMinus);
290 A_muPlus[i] += s[i] * Rho_muPlus * A_rho[i];
291 A_muMinus[i] += s[i] * Rho_muMinus * A_rho[i];
292 if(A_s) A_s[i] += Rho * A_rho[i];
294 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;
296 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;
300 __hostanddev__
double rootFunc(
double x,
double V)
const 301 {
double f_x; fHS(x, f_x);
302 return x - (x0plus*
exp(-V-f_x*x0plus) + x0minus*
exp(+V-f_x*x0minus));
307 {
double xLo = x0;
while(rootFunc(xLo, V) > 0.) xLo *= 0.5;
308 double xHi = xLo;
while(rootFunc(xHi, V) < 0.) xHi *= 2.;
309 double x = 0.5*(xHi+xLo);
312 {
if(rootFunc(x, V) < 0.)
322 __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 323 {
double V = ZbyT * phi[i];
327 V = copysign(1e-7, V);
329 double twoCbrtV= 2.*
pow(fabs(V), 1./3);
330 double Vmapped = copysign(twoCbrtV / (1. +
sqrt(1. + twoCbrtV*twoCbrtV)), V);
331 double xMapped = xLookup(1.+Vmapped);
332 double x = 1./xMapped - 1.;
333 double f_x; fHS(x, f_x);
334 double logEtaPlus = -V - f_x*x0plus;
335 double logEtaMinus = +V - f_x*x0minus;
337 { muPlus[i] = logEtaPlus;
338 muMinus[i] = -logEtaMinus;
341 kappaSq[i] = (4*M_PI)*s[i]*(NZ*ZbyT)*(
exp(logEtaMinus) -
exp(logEtaPlus))/V;
343 void phiToState(
size_t N,
const double* phi,
const double* s,
const RadialFunctionG& xLookup,
bool setState,
double* muPlus,
double* muMinus,
double* kappaSq)
const;
345 void phiToState_gpu(
size_t N,
const double* phi,
const double* s,
const RadialFunctionG& xLookup,
bool setState,
double* muPlus,
double* muMinus,
double* kappaSq)
const;
357 Dielectric(
bool linear,
double T,
double Nmol,
double pMol,
double epsBulk,
double epsInf);
360 __hostanddev__
void calcFunctions(
double eps,
double& frac,
double& frac_epsSqHlf,
double& logsinch)
const 361 {
double epsSq = eps*eps;
365 logsinch = epsSq*(1.0/6);
369 { frac = 1.0/3 + epsSq*(-1.0/45 + epsSq*(2.0/945 + epsSq*(-1.0/4725)));
370 frac_epsSqHlf = -2.0/45 + epsSq*(8.0/945 + epsSq*(-6.0/4725));
371 logsinch = epsSq*(1.0/6 + epsSq*(-1.0/180 + epsSq*(1.0/2835)));
374 { frac = (eps/tanh(eps)-1)/epsSq;
375 frac_epsSqHlf = (2 - eps/tanh(eps) -
pow(eps/sinh(eps),2))/(epsSq*epsSq);
376 logsinch = eps<20. ?
log(sinh(eps)/eps) : eps -
log(2.*eps);
382 __hostanddev__
void compute(
double epsSqHlf,
double& F,
double& F_epsSqHlf,
double& ChiEff,
double& ChiEff_epsSqHlf)
const 383 {
double epsSq = 2.*epsSqHlf, eps =
sqrt(epsSq);
385 double frac, frac_epsSqHlf, logsinch;
386 calcFunctions(eps, frac, frac_epsSqHlf, logsinch);
388 double screen = 1 - alpha*frac;
389 F = NT * (epsSq*(frac - 0.5*alpha*frac*frac + 0.5*X*screen*screen) - logsinch);
390 F_epsSqHlf = NT * (frac + X*screen + epsSq*frac_epsSqHlf*(1.-X*alpha)) * screen;
392 ChiEff = Np * (frac + X*screen);
393 ChiEff_epsSqHlf = Np * frac_epsSqHlf * (1.-X*alpha);
399 double epsSqHlf = 0.5*epsVec.length_squared();
400 double F, F_epsSqHlf, ChiEff, ChiEff_epsSqHlf;
401 compute(epsSqHlf, F, F_epsSqHlf, ChiEff, ChiEff_epsSqHlf);
403 accumVector((F_epsSqHlf * s[i]) * epsVec, A_eps, i);
415 double epsSqHlf = 0.5*epsVec.length_squared();
416 double F, F_epsSqHlf, ChiEff, ChiEff_epsSqHlf;
417 compute(epsSqHlf, F, F_epsSqHlf, ChiEff, ChiEff_epsSqHlf);
420 accumVector(s[i]*( ChiEff*A_pVec + ChiEff_epsSqHlf*epsVec*
dot(A_pVec, epsVec)), A_eps, i);
421 if(A_s) A_s[i] += ChiEff *
dot(A_pVec, epsVec);
430 {
double frac, frac_epsSqHlf, logsinch;
431 calcFunctions(eps, frac, frac_epsSqHlf, logsinch);
432 return eps*(1. - alpha*frac);
438 double epsLo = x;
while(x_from_eps(epsLo) > x) epsLo *= 0.95;
439 double epsHi = epsLo;
while(x_from_eps(epsHi) < x) epsHi *= 1.05;
440 double eps = 0.5*(epsHi+epsLo);
441 double deps = eps*1e-13;
442 while(epsHi-epsLo > deps)
443 {
if(x_from_eps(eps) < x)
447 eps = 0.5*(epsHi+epsLo);
455 double x = xVec.length();
456 double g = gLookup(x/(1.+x));
460 epsilon[i] = 1. + (4*M_PI)*s[i]*(Np*pByT)*((g-1.)/alpha + X);
470 #define FLUID_DUMP(object, suffix) \ 471 filename = filenamePattern; \ 472 filename.replace(filename.find("%s"), 2, suffix); \ 473 logPrintf("Dumping '%s'... ", filename.c_str()); logFlush(); \ 474 if(mpiUtil->isHead()) saveRawBinary(object, filename.c_str()); \ 475 logPrintf("done.\n"); logFlush(); 477 #endif // JDFTX_ELECTRONIC_PCM_INTERNAL_H __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:300
ScalarField pow(const ScalarField &, double alpha)
Elementwise power (preserve input)
double NT
N*p, p/T and N*T where N is molecular density and p is molecular dipole.
Definition: PCM_internal.h:354
__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:382
Definition: PCM_internal.h:59
ScalarField sqrt(const ScalarField &)
Elementwise square root (preserve input)
Helper class for dielectric portion of NonlinearPCM.
Definition: PCM_internal.h:351
Definition: PCM_internal.h:33
__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:397
complex dot(const Tptr &X, const Tptr &Y)
Definition: Operators.h:196
__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's if setState=true or effective kappaSq if setStat...
Definition: PCM_internal.h:322
double integral(const ScalarField &)
Integral in the unit cell (dV times sum())
__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:360
std::shared_ptr< ScalarFieldTildeData > ScalarFieldTilde
A smart reference-counting pointer to ScalarFieldTildeData.
Definition: ScalarField.h:45
__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:453
__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:306
bool linear
whether ionic screening is linearized
Definition: PCM_internal.h:177
__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:436
__hostanddev__ double fHS(double xIn, double &f_xIn) const
Hard sphere free energy per particle and derivative, where x is total packing fraction.
Definition: PCM_internal.h:232
Helper class for ionic screening portion of NonlinearPCM.
Definition: PCM_internal.h:175
__hostanddev__ void storeVector(const vector3< scalar > &v, vector3< scalar * > &vArr, int i)
Store vector to a vector field.
Definition: vector3.h:171
__hostanddev__ void accumVector(const vector3< scalar > &v, vector3< scalar * > &vArr, int i)
Accumulate vector onto a vector field.
Definition: vector3.h:175
ScalarField log(const ScalarField &)
Elementwise logarithm (preserve input)
Definition: PCM_internal.h:53
double X
dipole correlation factor and chi*T/p^2 where chi is the molecular susceptibility ...
Definition: PCM_internal.h:355
__hostanddev__ double x_from_eps(double eps) const
Calculate x = pMol E / T given eps.
Definition: PCM_internal.h:429
bool linear
whether dielectric is linearized
Definition: PCM_internal.h:353
__hostanddev__ vector3< scalar > loadVector(const vector3< const scalar * > &vArr, int i)
Load vector from a constant vector field.
Definition: vector3.h:163
double NZ
where T=temperature, N=bulk ionic concentration, Z=charge (all assumed +/- symmetric) ...
Definition: PCM_internal.h:178
G-space radial function stored on a uniform grid (of |G|)
Definition: RadialFunction.h:28
ScalarField exp(const ScalarField &)
Elementwise exponential (preserve input)
__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:249
Definition: NonlinearPCM.h:29
std::shared_ptr< ScalarFieldData > ScalarField
A smart reference-counting pointer to ScalarFieldData.
Definition: ScalarField.h:41
__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:272
Represent components of the (free) energy.
__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:287
double x0
anion, cation and total packing fractions
Definition: PCM_internal.h:179
Generic 3-vector.
Definition: vector3.h:33
Definition: PCM_internal.h:42
__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:413
Operators on ScalarField's and ScalarFieldTilde's.