20 #ifndef JDFTX_FLUID_MIXEDFMT_INTERNAL_H 21 #define JDFTX_FLUID_MIXEDFMT_INTERNAL_H 24 #include <core/matrix3.h> 31 double Gsq = Gvec.length_squared();
32 mTilde.xy()[i] = minus_nTilde*Gvec.x()*Gvec.y();
33 mTilde.yz()[i] = minus_nTilde*Gvec.y()*Gvec.z();
34 mTilde.zx()[i] = minus_nTilde*Gvec.z()*Gvec.x();
35 mTilde.
xxr()[i] = minus_nTilde*(Gvec.x()*Gvec.x() - (1.0/3)*Gsq);
36 mTilde.
yyr()[i] = minus_nTilde*(Gvec.y()*Gvec.y() - (1.0/3)*Gsq);
41 { complex temp = complex(0,0);
44 double Gsq = Gvec.length_squared();
45 temp += grad_mTilde.xy()[i]*Gvec.x()*Gvec.y();
46 temp += grad_mTilde.yz()[i]*Gvec.y()*Gvec.z();
47 temp += grad_mTilde.zx()[i]*Gvec.z()*Gvec.x();
48 temp += grad_mTilde.
xxr()[i]*(Gvec.x()*Gvec.x() - (1.0/3)*Gsq);
49 temp += grad_mTilde.
yyr()[i]*(Gvec.y()*Gvec.y() - (1.0/3)*Gsq);
51 grad_nTilde[i] = -temp;
56 {
return 2*(m.xy()*v.x()*v.y() + m.yz()*v.y()*v.z() + m.zx()*v.z()*v.x())
57 + m.
xxr()*v.x()*v.x() + m.
yyr()*v.y()*v.y() - (m.
xxr()+m.
yyr())*v.z()*v.z();
60 __hostanddev__
void mul_vTmv_grad(
const double grad_mul,
const tensor3<>& m,
const vector3<>& v,
62 { grad_m.xy() += (2*grad_mul)*v.x()*v.y();
63 grad_m.yz() += (2*grad_mul)*v.y()*v.z();
64 grad_m.zx() += (2*grad_mul)*v.z()*v.x();
65 grad_m.
xxr() += grad_mul*(v.x()*v.x()-v.z()*v.z());
66 grad_m.
yyr() += grad_mul*(v.y()*v.y()-v.z()*v.z());
67 grad_v.x() += (2*grad_mul)*(m.xy()*v.y() + m.zx()*v.z() + m.
xxr()*v.x());
68 grad_v.y() += (2*grad_mul)*(m.xy()*v.x() + m.yz()*v.z() + m.
yyr()*v.y());
69 grad_v.z() += (2*grad_mul)*(m.yz()*v.y() + m.zx()*v.x()- (m.
xxr()+m.
yyr())*v.z());
74 __hostanddev__
double trace_cubed(
const tensor3<>& m)
75 {
return -3.0 * (-2.0*m.xy()*m.yz()*m.zx() +
pow(m.yz(),2)*m.
xxr() -
pow(m.xy(),2)*(m.
xxr()+m.
yyr()) + m.
yyr()*(
pow(m.zx(),2)+m.
xxr()*(m.
xxr()+m.
yyr())));
78 __hostanddev__
void trace_cubed_grad(
const double grad_trace,
const tensor3<>& m,
tensor3<>& grad_m)
79 { grad_m.xy() += (6*grad_trace)* (m.yz()*m.zx()+m.xy()*(m.
xxr()+m.
yyr()));
80 grad_m.yz() += (6*grad_trace)* (m.xy()*m.zx()-m.yz()*m.
xxr());
81 grad_m.zx() += (6*grad_trace)* (m.xy()*m.yz()-m.zx()*m.
yyr());
82 grad_m.
xxr() += (-3*grad_trace)* (-
pow(m.xy(),2)+
pow(m.yz(),2)+m.
yyr()*(2*m.
xxr()+m.
yyr()));
83 grad_m.
yyr() += (-3*grad_trace)* (-
pow(m.xy(),2)+
pow(m.zx(),2)+m.
xxr()*(m.
xxr()+2*m.
yyr()));
87 __hostanddev__
double WB_f2(
double x,
double& f2_x)
89 { f2_x = x*((2./9) + x*(3./18 + x*(4./30)));
90 return 1 + x*x*((1./9) + x*(1./18 + x*(1./30)));
93 { f2_x = (-1./3)*(x*(2+x) + 2*
log(1-x)) / (x*x);
94 return 1 + (1./3)*(2-x + 2*(1-x)*
log(1-x)/x);
97 __hostanddev__
double WB_f3(
double x,
double& f3_x)
99 { f3_x = -4./9 + x*(2./18 + x*(3./45 + x*(4./90)));
100 return 1 + x*(-4./9 + x*(1./18 + x*(1./45 + x*(1./90))));
103 { f3_x = 2*(1-x) * (x*(2+x) + 2*
log(1-x)) / (3*x*x*x);
104 return 1 - (x*(2+x*(-3+x*2)) + 2*(1-x)*(1-x)*
log(1-x)) / (3*x*x);
108 __hostanddev__
double phiFMT_calc(
int i,
109 const double *n0arr,
const double *n1arr,
const double *n2arr,
const double *n3arr,
111 double *grad_n0arr,
double *grad_n1arr,
double *grad_n2arr,
double *grad_n3arr,
114 double n0 = n0arr[i];
115 double n1 = n1arr[i];
116 double n2 = n2arr[i];
117 double n3 = n3arr[i];
118 if(n0<0. || n1<0. || n2<0. || n3<0.)
return 0.;
119 if(n3>=1.)
return NAN;
123 double tensorPart = mul_vTmv(n2m, n2v) - 0.5*trace_cubed(n2m);
125 double n1v_n2v =
dot(n1v, n2v);
126 double n2vsq = n2v.length_squared();
127 double pole = 1./(1-n3);
129 double f2prime, f2 = WB_f2(n3, f2prime), comb2 = (n1*n2-n1v_n2v);
130 double f3prime, f3 = WB_f3(n3, f3prime), comb3 = (n2*(n2*n2-3*n2vsq) + 9.*tensorPart);
131 double phi = n0*
log(pole) + pole*(f2*comb2 + pole*((1./((24*M_PI)))*f3*comb3));
132 double phi_n0 =
log(pole);
133 double phi_n1 = pole*(f2*n2);
134 double phi_n2 = pole*(f2*n1 + pole*((1./(8*M_PI))*f3*(n2*n2-n2vsq)));
135 double phi_n3 = pole*(n0 + pole*(f2*comb2 + pole*((1./(12*M_PI))*f3*comb3)))
136 + pole*(f2prime*comb2 + pole*((1./(24*M_PI))*f3prime*comb3));
138 vector3<> phi_n2v = (-pole)*(f2*n1v + (pole*f3*n2/(4*M_PI))*n2v );
139 double phi_tensorPart = pole*pole*f3*(9./(24*M_PI));
142 mul_vTmv_grad(phi_tensorPart, n2m, n2v, phi_n2m, phi_n2v);
143 trace_cubed_grad(-0.5*phi_tensorPart, n2m, phi_n2m);
147 grad_n3arr[i] += phi_n3;
148 grad_n2arr[i] += phi_n2;
149 grad_n1arr[i] += phi_n1;
150 grad_n0arr[i] += phi_n0;
154 __hostanddev__
double phiBond_calc(
int i,
double Rhm,
double scale,
156 double *grad_n0arr,
double *grad_n2arr,
double *grad_n3arr,
vector3<double*> grad_n2vArr)
158 double n0 = n0arr[i];
159 double n2 = n2arr[i];
160 double n3 = n3arr[i];
161 if(n0<0. || n2<0. || n3<0.)
return 0.;
162 double pole = 1.0/(1-n3);
165 double n2vsq = n2v.length_squared();
166 double zeta = (n2<=0 || n2vsq>n2*n2) ? 0.0 : (1-n2vsq/(n2*n2));
167 double zeta_n2 = zeta ? 2*n2vsq/(n2*n2*n2) : 0.0;
170 double gContact = pole*(1 + zeta * pole*(Rhm*n2 + pole*(2.0/9)*
pow(Rhm*n2,2) ));
171 double gContact_n2 = zeta *
pow(pole,2)*(Rhm + pole*(4.0/9)*
pow(Rhm,2)*n2 );
172 double gContact_n3 =
pow(pole,2)*(1 + zeta * pole*(2*Rhm*n2 + pole*(6.0/9)*
pow(Rhm*n2,2) ));
173 double gContact_zeta =
pow(pole,2) * (Rhm*n2 + pole*(2.0/9)*
pow(Rhm*n2,2) );
175 double phi = -scale*n0*
log(gContact);
176 double phi_n0 = -scale*
log(gContact);
177 double phi_n2 = -scale*n0*gContact_n2/gContact;
178 double phi_n3 = -scale*n0*gContact_n3/gContact;
179 double phi_zeta = -scale*n0*gContact_zeta/gContact;
181 grad_n0arr[i] += phi_n0;
182 grad_n2arr[i] += phi_n2 + phi_zeta * zeta_n2;
183 grad_n3arr[i] += phi_n3;
188 #endif // JDFTX_FLUID_MIXEDFMT_INTERNAL_H ScalarField pow(const ScalarField &, double alpha)
Elementwise power (preserve input)
__hostanddev__ tensor3< scalar > loadTensor(const tensor3< const scalar * > &tArr, int i)
Load tensor from a constant tensor field.
Definition: tensor3.h:55
complex dot(const Tptr &X, const Tptr &Y)
Definition: Operators.h:196
__hostanddev__ scalar & xxr()
xxr = x^2 - r^2/3
Definition: tensor3.h:42
__hostanddev__ void accumVector(const vector3< scalar > &v, vector3< scalar * > &vArr, int i)
Accumulate vector onto a vector field.
Definition: vector3.h:175
__hostanddev__ scalar & yyr()
yyr = y^2-r^2/3
Definition: tensor3.h:43
ScalarField log(const ScalarField &)
Elementwise logarithm (preserve input)
Symmetric traceless rank-2 tensor in 3D.
Definition: tensor3.h:32
__hostanddev__ void accumTensor(const tensor3< scalar > &t, tensor3< scalar * > &tArr, int i)
Accumulate tensor onto a tensor field.
Definition: tensor3.h:67
__hostanddev__ vector3< scalar > loadVector(const vector3< const scalar * > &vArr, int i)
Load vector from a constant vector field.
Definition: vector3.h:163
Complex number (need to define our own because we need operators for gpu code as well) ...
Definition: scalar.h:49