1 /*-------------------------------------------------------------------
2 Copyright 2012 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_CORE_COULOMB_INTERNAL_H
21 #define JDFTX_CORE_COULOMB_INTERNAL_H
22
23 //! @addtogroup LongRange
24 //! @{
25
26 //! @file Coulomb_internal.h Shared inline functions / internal declarations for Coulomb framework
27
28 #include <core/matrix3.h>
29 #include <core/Spline.h>
30 #include <gsl/gsl_integration.h>
31
32 //Common citation for Coulomb truncation
33 #define wsTruncationPaper "R. Sundararaman and T.A. Arias, Phys. Rev. B 87, 165122 (2013)"
34
35 //Ion-margin error message
36 #define ionMarginMessage "Expand unit cell, or if absolutely sure, reduce coulomb-truncation-ion-margin.\n"
37
38 //! Periodic coulomb interaction (4 pi/G^2)
39 struct CoulombPeriodic_calc
operatorCoulombPeriodic_calc40 { __hostanddev__ double operator()(const vector3<int>& iG, const matrix3<>& GGT) const
41 { double Gsq = GGT.metric_length_squared(iG);
42 return Gsq ? (4*M_PI)/Gsq : 0.;
43 }
44 __hostanddev__ symmetricMatrix3<> latticeGradient(const vector3<int>& iG, const matrix3<>& GGT) const
45 { double Gsq = GGT.metric_length_squared(iG);
46 return (Gsq ? (8*M_PI)/(Gsq*Gsq) : 0.) * outer(vector3<>(iG));
47 }
48 };
49
50 //! Slab-truncated coulomb interaction
51 struct CoulombSlab_calc
52 { int iDir; double hlfL;
CoulombSlab_calcCoulombSlab_calc53 CoulombSlab_calc(int iDir, double hlfL) : iDir(iDir), hlfL(hlfL) {}
operatorCoulombSlab_calc54 __hostanddev__ double operator()(const vector3<int>& iG, const matrix3<>& GGT) const
55 { double Gsq = GGT.metric_length_squared(iG);
56 double Gplane = Gsq - GGT(iDir,iDir) * iG[iDir]*iG[iDir]; //G along the non-truncated directions
57 Gplane = Gplane>0. ? sqrt(Gplane) : 0.; //safe sqrt to prevent NaN from roundoff errors
58 return (4*M_PI) * (Gsq ? (1. - exp(-Gplane*hlfL) * cos(M_PI*iG[iDir]))/Gsq : -0.5*hlfL*hlfL);
59 }
latticeGradientCoulombSlab_calc60 __hostanddev__ symmetricMatrix3<> latticeGradient(const vector3<int>& iG, const matrix3<>& GGT) const
61 { double Gsq = GGT.metric_length_squared(iG);
62 if(Gsq)
63 { double GsqInv = 1./Gsq;
64 vector3<int> iGplane = iG; iGplane[iDir] = 0.; //iG along non-truncated directions
65 double Gplane = sqrt(GGT.metric_length_squared(iGplane));
66 double GplaneInv = Gplane ? 1./Gplane : 0.;
67 double expCosTerm = exp(-Gplane*hlfL) * cos(M_PI*iG[iDir]);
68 //Compute planar part:
69 double prefac1 = 2.*(1.-expCosTerm)*GsqInv;
70 symmetricMatrix3<> result = (prefac1 - expCosTerm*hlfL*GplaneInv) * outer(vector3<>(iGplane));
71 //Fix normal part:
72 ((double*)&result)[iDir] = prefac1*iG[iDir]*iG[iDir] + expCosTerm*hlfL*Gplane/GGT(iDir,iDir);
73 return (4*M_PI*GsqInv) * result;
74 }
75 else
76 { symmetricMatrix3<> result;
77 ((double*)&result)[iDir] = (-4*M_PI)*hlfL*hlfL/GGT(iDir,iDir);
78 return result;
79 }
80 }
81 };
82
83 //! Sphere-truncated coulomb interaction
84 struct CoulombSpherical_calc
85 { double Rc;
CoulombSpherical_calcCoulombSpherical_calc86 CoulombSpherical_calc(double Rc) : Rc(Rc) {}
operatorCoulombSpherical_calc87 __hostanddev__ double operator()(const vector3<int>& iG, const matrix3<>& GGT) const
88 { double Gsq = GGT.metric_length_squared(iG);
89 return Gsq ? (4*M_PI) * (1. - cos(Rc*sqrt(Gsq)))/Gsq : (2*M_PI)*Rc*Rc;
90 }
latticeGradientCoulombSpherical_calc91 __hostanddev__ symmetricMatrix3<> latticeGradient(const vector3<int>& iG, const matrix3<>& GGT) const
92 { double Gsq = GGT.metric_length_squared(iG);
93 double GRc = sqrt(Gsq)*Rc, cosGRc, sinGRc;
94 sincos(GRc, &sinGRc, &cosGRc);
95 return (Gsq ? (4*M_PI)*(2*(1.-cosGRc) - GRc*sinGRc)/(Gsq*Gsq) : 0.) * outer(vector3<>(iG));
96 }
97 };
98
99 //! Lattice derivative calculation for ionKernel used in embedded mode
100 struct CoulombIonKernel_calc
101 { double expFac;
CoulombIonKernel_calcCoulombIonKernel_calc102 CoulombIonKernel_calc(double ionWidth) : expFac(0.5*ionWidth*ionWidth) {}
latticeGradientCoulombIonKernel_calc103 __hostanddev__ symmetricMatrix3<> latticeGradient(const vector3<int>& iG, const matrix3<>& GGT) const
104 { double Gsq = GGT.metric_length_squared(iG);
105 return (Gsq ? (8*M_PI)*(1.-exp(-expFac*Gsq)*(1.+expFac*Gsq))/(Gsq*Gsq) : 0.) * outer(vector3<>(iG));
106 }
107 };
108
109
110 #ifdef GPU_ENABLED
111 void coulombAnalytic_gpu(vector3<int> S, const matrix3<>& GGT, const CoulombPeriodic_calc& calc, complex* data);
112 void coulombAnalytic_gpu(vector3<int> S, const matrix3<>& GGT, const CoulombSlab_calc& calc, complex* data);
113 void coulombAnalytic_gpu(vector3<int> S, const matrix3<>& GGT, const CoulombSpherical_calc& calc, complex* data);
114 void coulombAnalyticStress_gpu(vector3<int> S, const matrix3<>& GGT, const CoulombPeriodic_calc& calc, const complex* X, const complex* Y, symmetricMatrix3<>* grad_RRT);
115 void coulombAnalyticStress_gpu(vector3<int> S, const matrix3<>& GGT, const CoulombSlab_calc& calc, const complex* X, const complex* Y, symmetricMatrix3<>* grad_RRT);
116 void coulombAnalyticStress_gpu(vector3<int> S, const matrix3<>& GGT, const CoulombSpherical_calc& calc, const complex* X, const complex* Y, symmetricMatrix3<>* grad_RRT);
117 void coulombAnalyticStress_gpu(vector3<int> S, const matrix3<>& GGT, const CoulombIonKernel_calc& calc, const complex* X, const complex* Y, symmetricMatrix3<>* grad_RRT);
118 void coulombNumericalStress_gpu(vector3<int> S, const matrix3<>& GGT, const symmetricMatrix3<>* Vc_RRT, const complex* X, const complex* Y, symmetricMatrix3<>* grad_RRT);
119 #endif
120 void coulombAnalytic(vector3<int> S, const matrix3<>& GGT, const CoulombPeriodic_calc& calc, complex* data);
121 void coulombAnalytic(vector3<int> S, const matrix3<>& GGT, const CoulombSlab_calc& calc, complex* data);
122 void coulombAnalytic(vector3<int> S, const matrix3<>& GGT, const CoulombSpherical_calc& calc, complex* data);
123 void coulombAnalyticStress(vector3<int> S, const matrix3<>& GGT, const CoulombPeriodic_calc& calc, const complex* X, const complex* Y, symmetricMatrix3<>* grad_RRT);
124 void coulombAnalyticStress(vector3<int> S, const matrix3<>& GGT, const CoulombSlab_calc& calc, const complex* X, const complex* Y, symmetricMatrix3<>* grad_RRT);
125 void coulombAnalyticStress(vector3<int> S, const matrix3<>& GGT, const CoulombSpherical_calc& calc, const complex* X, const complex* Y, symmetricMatrix3<>* grad_RRT);
126 void coulombAnalyticStress(vector3<int> S, const matrix3<>& GGT, const CoulombIonKernel_calc& calc, const complex* X, const complex* Y, symmetricMatrix3<>* grad_RRT);
127 void coulombNumericalStress(vector3<int> S, const matrix3<>& GGT, const symmetricMatrix3<>* Vc_RRT, const complex* X, const complex* Y, symmetricMatrix3<>* grad_RRT);
128
129 //! Compute erf(x)/x (with x~0 handled properly)
erf_by_x(double x)130 __hostanddev__ double erf_by_x(double x)
131 { double xSq = x*x;
132 if(xSq<1e-6) return (1./sqrt(M_PI))*(2. + xSq*(-2./3 + 0.2*xSq));
133 else return erf(x)/x;
134 }
135
136 //! Compute (1/x) d(erf(x)/x)/dx (with x~0 handled properly)
erf_by_xPrime_by_x(double x)137 __hostanddev__ double erf_by_xPrime_by_x(double x)
138 { double xSq = x*x;
139 if(xSq<1e-6) return (2./sqrt(M_PI))*(-2./3 + 0.4*xSq);
140 else return ((2./sqrt(M_PI))*x*exp(-xSq) - erf(x))/(x*xSq);
141 }
142
143
144
145 //--------------- Special function for cylinder/wire modes ------------
146 // (implemented in CoulombWire.cpp)
147
148 //! Compute Cbar_k^sigma - the gaussian convolved cylindrical coulomb kernel - by numerical quadrature
149 struct Cbar
150 { Cbar();
151 ~Cbar();
152 double operator()(double k, double sigma, double rho, double rho0=1.); //!< Compute Cbar_k^sigma(rho)
153 private:
154 static const size_t maxIntervals = 1000; //!< Size of integration workspace
155 gsl_integration_workspace* iWS; //!< Integration workspace
156 static double integrandSmallRho(double t, void* params); //!< Integrand for rho < sigma
157 static double integrandLargeRho(double t, void* params); //!< Integrand for rho > sigma
158 };
159
160 //! Look-up table for Cbar_k^sigma(rho) for specific values of k and sigma
161 //! If prime=true, construct lookup table for -d(Cbar_k^sigma(rho))/dk instead
162 struct Cbar_k_sigma
163 { Cbar_k_sigma(double k, double sigma, double rhoMax, double rho0=1., bool prime=false);
164 //! Get value:
valueCbar_k_sigma165 inline double value(double rho) const
166 { double f = QuinticSpline::value(coeff.data(), drhoInv * rho);
167 return isLog ? exp(f) : f;
168 }
169 //! Get derivative:
derivCbar_k_sigma170 inline double deriv(double rho) const
171 { double fp = QuinticSpline::deriv(coeff.data(), drhoInv * rho) * drhoInv;
172 return isLog ? fp * value(rho) : fp;
173 }
174 private:
175 double drhoInv; bool isLog;
176 std::vector<double> coeff;
177 };
178
179
180 //---------------------- Exchange Kernels --------------------
181
182 //In each of the following functions, kSq is the square of the appropriate
183 //wave vector (includes reciprocal lattice vector and k-point difference),
184 //and will not be zero (the G=0 term is handled in the calling routine)
185
186 //! Radial fourier transform of erfc(omega r)/r (not valid at G=0)
erfcTilde(double Gsq,double omegaSq)187 __hostanddev__ double erfcTilde(double Gsq, double omegaSq)
188 { return (4*M_PI) * (omegaSq ? (1.-exp(-0.25*Gsq/omegaSq)) : 1.) / Gsq;
189 }
190
191
192 //! Periodic exchange
193 struct ExchangePeriodic_calc
operatorExchangePeriodic_calc194 { __hostanddev__ double operator()(double kSq) const
195 { return (4*M_PI) / kSq;
196 }
197 __hostanddev__ double latticeGradientPrefac(double kSq, double& grad_lnDetR) const
198 { return (8*M_PI) / (kSq*kSq);
199 }
200 };
201
202 //! Erfc-screened Periodic exchange
203 struct ExchangePeriodicScreened_calc
204 { double inv4omegaSq; //!< 1/(4 omega^2)
ExchangePeriodicScreened_calcExchangePeriodicScreened_calc205 ExchangePeriodicScreened_calc(double omega) : inv4omegaSq(0.25/(omega*omega)) {}
206
operatorExchangePeriodicScreened_calc207 __hostanddev__ double operator()(double kSq) const
208 { return (4*M_PI) * (1.-exp(-inv4omegaSq*kSq)) / kSq;
209 }
latticeGradientPrefacExchangePeriodicScreened_calc210 __hostanddev__ double latticeGradientPrefac(double kSq, double& grad_lnDetR) const
211 { return (8*M_PI) * (1.-exp(-inv4omegaSq*kSq)*(1.+inv4omegaSq*kSq)) / (kSq*kSq);
212 }
213 };
214
215 //! Spherical-truncated exchange
216 struct ExchangeSpherical_calc
217 { double Rc;
ExchangeSpherical_calcExchangeSpherical_calc218 ExchangeSpherical_calc(double Rc) : Rc(Rc) {}
219
operatorExchangeSpherical_calc220 __hostanddev__ double operator()(double kSq) const
221 { return (4*M_PI) * (1. - cos(Rc * sqrt(kSq))) / kSq;
222 }
latticeGradientPrefacExchangeSpherical_calc223 __hostanddev__ double latticeGradientPrefac(double kSq, double& grad_lnDetR) const
224 { double kR = sqrt(kSq)*Rc, coskR, sinkR;
225 sincos(kR, &sinkR, &coskR);
226 grad_lnDetR = (4*M_PI) * (kR*sinkR) / (3*kSq);
227 return (4*M_PI) * (2*(1.-coskR) - kR*sinkR)/(kSq*kSq);
228 }
229 };
230
231 //! Erfc-screened Spherical-truncated exchange
232 struct ExchangeSphericalScreened_calc
233 { double* coeff; //!< quintic spline coefficients
234 double dGinv; //!< inverse of coefficient spacing
235 size_t nSamples; //!< number of coefficients
ExchangeSphericalScreened_calcExchangeSphericalScreened_calc236 ExchangeSphericalScreened_calc() : coeff(0) {}
237
operatorExchangeSphericalScreened_calc238 __hostanddev__ double operator()(double kSq) const
239 { double t = dGinv * sqrt(kSq);
240 if(t >= nSamples) return 0.;
241 else return QuinticSpline::value(coeff, t);
242 }
243
244 //Additional members set for stress calculation:
245 double Rc; //!< truncation radius
246 double erfcOmegaRc_4piBy3; //!< (4 pi/3) erfc(omega*Rc), where omega is the screening parameter
247
latticeGradientPrefacExchangeSphericalScreened_calc248 __hostanddev__ double latticeGradientPrefac(double kSq, double& grad_lnDetR) const
249 { double k = sqrt(kSq);
250 grad_lnDetR = erfcOmegaRc_4piBy3 * sin(k*Rc) * Rc/k;
251 double t = dGinv * k;
252 if(t >= nSamples) return 0.;
253 else return QuinticSpline::deriv(coeff, t) * (-t/kSq);
254 }
255 };
256
257 //! Slab-truncated exchange
258 struct ExchangeSlab_calc
259 { int iDir; double hlfL;
260 double* coeff; double dGinv; size_t nSamples, nCoeff; //quintic-spline coefficients for screened mode
261 symmetricMatrix3<> zHatOuter;
ExchangeSlab_calcExchangeSlab_calc262 ExchangeSlab_calc() : coeff(0) {}
263
operatorExchangeSlab_calc264 __hostanddev__ double operator()(const vector3<int>& iG, const matrix3<>& GGT, const vector3<>& kDiff, double Vzero, double thresholdSq) const
265 { vector3<> g = iG + kDiff; //net G-vector in reciprocal lattice coordinates including k-point
266 double Gsq = GGT.metric_length_squared(g);
267 if(Gsq < thresholdSq)
268 return Vzero;
269 double Gplane = Gsq - GGT(iDir,iDir) * iG[iDir]*iG[iDir]; //G along the non-truncated directions (note kDiff[iDir]=0)
270 Gplane = Gplane>0. ? sqrt(Gplane) : 0.; //safe sqrt to prevent NaN from roundoff errors
271 double Vc = (4*M_PI) * (1. - exp(-Gplane*hlfL) * cos(M_PI*iG[iDir]))/Gsq; //Unscreened exchange (calculate analytically)
272 if(coeff)
273 { //Correct for Screened exchange using lookup table:
274 const double* coeffPlane = coeff + abs(iG[iDir]) * nCoeff; //coefficients for this plane
275 double t = dGinv * Gplane;
276 double prefac = iG[iDir] ? 1. : 1./Gplane;
277 if(t<nSamples) Vc += prefac * QuinticSpline::value(coeffPlane, t);
278 }
279 return Vc;
280 }
latticeGradientExchangeSlab_calc281 __hostanddev__ symmetricMatrix3<> latticeGradient(const vector3<int>& iG, const matrix3<>& G, const vector3<>& kDiff, double thresholdSq) const
282 { double Gsq = ((iG+kDiff)*G).length_squared();
283 if(Gsq > thresholdSq)
284 { double GsqInv = 1./Gsq;
285 vector3<> iGplane = iG+kDiff; iGplane[iDir] = 0.; //iG along non-truncated directions
286 vector3<> iGplaneCart = iGplane * G;
287 double GplaneSq = iGplaneCart.length_squared();
288 double GaxisSq = Gsq - GplaneSq;
289 double Gplane = sqrt(GplaneSq);
290 double GplaneInv = Gplane ? 1./Gplane : 0.;
291 double expCosTerm = exp(-Gplane*hlfL) * cos(M_PI*iG[iDir]);
292 //Planar part:
293 double Vc0 = 4*M_PI*GsqInv;
294 double prefac1 = 2.*(1.-expCosTerm)*GsqInv;
295 symmetricMatrix3<> result = (Vc0*(prefac1 - expCosTerm*hlfL*GplaneInv)) * outer(iGplaneCart);
296 //Normal part:
297 double result_zz = Vc0*(prefac1*GaxisSq + expCosTerm*hlfL*Gplane);
298 //Combine and return:
299 result += result_zz * zHatOuter;
300 return result;
301 }
302 else return symmetricMatrix3<>();
303 }
304 };
305
306
307 void exchangeAnalytic(vector3<int> S, const matrix3<>& GGT, const ExchangePeriodic_calc& calc, complex* data, const vector3<>& kDiff, double Vzero, double thresholdSq);
308 void exchangeAnalytic(vector3<int> S, const matrix3<>& GGT, const ExchangePeriodicScreened_calc& calc, complex* data, const vector3<>& kDiff, double Vzero, double thresholdSq);
309 void exchangeAnalytic(vector3<int> S, const matrix3<>& GGT, const ExchangeSpherical_calc& calc, complex* data, const vector3<>& kDiff, double Vzero, double thresholdSq);
310 void exchangeAnalytic(vector3<int> S, const matrix3<>& GGT, const ExchangeSphericalScreened_calc& calc, complex* data, const vector3<>& kDiff, double Vzero, double thresholdSq);
311 void exchangeAnalytic(vector3<int> S, const matrix3<>& GGT, const ExchangeSlab_calc& calc, complex* data, const vector3<>& kDiff, double Vzero, double thresholdSq);
312 #ifdef GPU_ENABLED
313 void exchangeAnalytic_gpu(vector3<int> S, const matrix3<>& GGT, const ExchangePeriodic_calc& calc, complex* data, const vector3<>& kDiff, double Vzero, double thresholdSq);
314 void exchangeAnalytic_gpu(vector3<int> S, const matrix3<>& GGT, const ExchangePeriodicScreened_calc& calc, complex* data, const vector3<>& kDiff, double Vzero, double thresholdSq);
315 void exchangeAnalytic_gpu(vector3<int> S, const matrix3<>& GGT, const ExchangeSpherical_calc& calc, complex* data, const vector3<>& kDiff, double Vzero, double thresholdSq);
316 void exchangeAnalytic_gpu(vector3<int> S, const matrix3<>& GGT, const ExchangeSphericalScreened_calc& calc, complex* data, const vector3<>& kDiff, double Vzero, double thresholdSq);
317 void exchangeAnalytic_gpu(vector3<int> S, const matrix3<>& GGT, const ExchangeSlab_calc& calc, complex* data, const vector3<>& kDiff, double Vzero, double thresholdSq);
318 #endif
319
320 void exchangeAnalyticStress(vector3<int> S, const matrix3<>& G, const ExchangePeriodic_calc& calc, const complex* X, symmetricMatrix3<>* grad_RRT, const vector3<>& kDiff, double thresholdSq);
321 void exchangeAnalyticStress(vector3<int> S, const matrix3<>& G, const ExchangePeriodicScreened_calc& calc, const complex* X, symmetricMatrix3<>* grad_RRT, const vector3<>& kDiff, double thresholdSq);
322 void exchangeAnalyticStress(vector3<int> S, const matrix3<>& G, const ExchangeSpherical_calc& calc, const complex* X, symmetricMatrix3<>* grad_RRT, const vector3<>& kDiff, double thresholdSq);
323 void exchangeAnalyticStress(vector3<int> S, const matrix3<>& G, const ExchangeSphericalScreened_calc& calc, const complex* X, symmetricMatrix3<>* grad_RRT, const vector3<>& kDiff, double thresholdSq);
324 void exchangeAnalyticStress(vector3<int> S, const matrix3<>& G, const ExchangeSlab_calc& calc, const complex* X, symmetricMatrix3<>* grad_RRT, const vector3<>& kDiff, double thresholdSq);
325 #ifdef GPU_ENABLED
326 void exchangeAnalyticStress_gpu(vector3<int> S, const matrix3<>& G, const ExchangePeriodic_calc& calc, const complex* X, symmetricMatrix3<>* grad_RRT, const vector3<>& kDiff, double thresholdSq);
327 void exchangeAnalyticStress_gpu(vector3<int> S, const matrix3<>& G, const ExchangePeriodicScreened_calc& calc, const complex* X, symmetricMatrix3<>* grad_RRT, const vector3<>& kDiff, double thresholdSq);
328 void exchangeAnalyticStress_gpu(vector3<int> S, const matrix3<>& G, const ExchangeSpherical_calc& calc, const complex* X, symmetricMatrix3<>* grad_RRT, const vector3<>& kDiff, double thresholdSq);
329 void exchangeAnalyticStress_gpu(vector3<int> S, const matrix3<>& G, const ExchangeSphericalScreened_calc& calc, const complex* X, symmetricMatrix3<>* grad_RRT, const vector3<>& kDiff, double thresholdSq);
330 void exchangeAnalyticStress_gpu(vector3<int> S, const matrix3<>& G, const ExchangeSlab_calc& calc, const complex* X, symmetricMatrix3<>* grad_RRT, const vector3<>& kDiff, double thresholdSq);
331 #endif
332
333
334 //! Multiply a complexScalarFieldTilde's data by a RealKernel (real-symmetry reduced)
multRealKernel_calc(size_t i,const vector3<int> & iG,const vector3<int> & S,const double * kernel,complex * data)335 __hostanddev__ void multRealKernel_calc(size_t i, const vector3<int>& iG,
336 const vector3<int>& S, const double* kernel, complex* data)
337 { //Compute index on the real kernel:
338 vector3<int> iGreal = iG;
339 if(iGreal[2]<0) iGreal = -iGreal; //inversion symmetry in G-space for real-kernels
340 if(iGreal[1]<0) iGreal[1] += S[1];
341 if(iGreal[0]<0) iGreal[0] += S[0];
342 size_t iReal = iGreal[2] + size_t(1+S[2]/2) * (iGreal[1] + S[1]*iGreal[0]);
343 //Multiply:
344 data[i] *= kernel[iReal];
345 }
346 void multRealKernel(vector3<int> S, const double* kernel, complex* data);
347 #ifdef GPU_ENABLED
348 void multRealKernel_gpu(vector3<int> S, const double* kernel, complex* data);
349 #endif
350
351 //! Multiply a complexScalarFieldTilde's data by a kernel sampled with offset and rotation by rot
multTransformedKernel_calc(size_t i,const vector3<int> & iG,const vector3<int> & S,const double * kernel,complex * data,const vector3<int> & offset)352 __hostanddev__ void multTransformedKernel_calc(size_t i, const vector3<int>& iG,
353 const vector3<int>& S, const double* kernel, complex* data, const vector3<int>& offset)
354 { vector3<int> iGkernel = (iG - offset); //Compute index on the real kernel
355 for(int k=0; k<3; k++) if(iGkernel[k]<0) iGkernel[k] += S[k]; //Reduce to [0,S-1) in each dimension
356 size_t iReal = iGkernel[2] + S[2]*size_t(iGkernel[1] + S[1]*iGkernel[0]); //net index into kernel
357 data[i] *= kernel[iReal];
358 }
359 void multTransformedKernel(vector3<int> S, const double* kernel, complex* data, const vector3<int>& offset);
360 #ifdef GPU_ENABLED
361 void multTransformedKernel_gpu(vector3<int> S, const double* kernel, complex* data, const vector3<int>& offset);
362 #endif
363
364
365 //! Compute stress corresponding to multRealKernel()
realKernelStress_calc(size_t i,const vector3<int> & iG,const vector3<int> & S,const symmetricMatrix3<> * kernel_RRT,const complex * X,symmetricMatrix3<> * grad_RRT)366 __hostanddev__ void realKernelStress_calc(size_t i, const vector3<int>& iG,
367 const vector3<int>& S, const symmetricMatrix3<>* kernel_RRT, const complex* X, symmetricMatrix3<>* grad_RRT)
368 { //Compute index on the real kernel:
369 vector3<int> iGreal = iG;
370 if(iGreal[2]<0) iGreal = -iGreal; //inversion symmetry in G-space for real-kernels
371 if(iGreal[1]<0) iGreal[1] += S[1];
372 if(iGreal[0]<0) iGreal[0] += S[0];
373 size_t iReal = iGreal[2] + size_t(1+S[2]/2) * (iGreal[1] + S[1]*iGreal[0]);
374 //Set stress contribution:
375 grad_RRT[i] = kernel_RRT[iReal] * X[i].norm();
376 }
377 void realKernelStress(vector3<int> S, const symmetricMatrix3<>* kernel_RRT, const complex* X, symmetricMatrix3<>* grad_RRT);
378 #ifdef GPU_ENABLED
379 void realKernelStress_gpu(vector3<int> S, const symmetricMatrix3<>* kernel_RRT, const complex* X, symmetricMatrix3<>* grad_RRT);
380 #endif
381
382
383 //! Compute stress corresponding to multTransformedKernel()
transformedKernelStress_calc(size_t i,const vector3<int> & iG,const vector3<int> & S,const symmetricMatrix3<> * kernel_RRT,const complex * X,symmetricMatrix3<> * grad_RRT,const vector3<int> & offset)384 __hostanddev__ void transformedKernelStress_calc(size_t i, const vector3<int>& iG,
385 const vector3<int>& S, const symmetricMatrix3<>* kernel_RRT, const complex* X, symmetricMatrix3<>* grad_RRT, const vector3<int>& offset)
386 { vector3<int> iGkernel = (iG - offset); //Compute index on the real kernel
387 for(int k=0; k<3; k++) if(iGkernel[k]<0) iGkernel[k] += S[k]; //Reduce to [0,S-1) in each dimension
388 size_t iReal = iGkernel[2] + S[2]*size_t(iGkernel[1] + S[1]*iGkernel[0]); //net index into kernel
389 //Set stress contribution:
390 grad_RRT[i] = kernel_RRT[iReal] * X[i].norm();
391 }
392 void transformedKernelStress(vector3<int> S, const symmetricMatrix3<>* kernel_RRT, const complex* X, symmetricMatrix3<>* grad_RRT, const vector3<int>& offset);
393 #ifdef GPU_ENABLED
394 void transformedKernelStress_gpu(vector3<int> S, const symmetricMatrix3<>* kernel_RRT, const complex* X, symmetricMatrix3<>* grad_RRT, const vector3<int>& offset);
395 #endif
396
397
398 //! @}
399 #endif // JDFTX_CORE_COULOMB_INTERNAL_H
400