1 //////////////////////////////////////////////////////////////////////////////////////
2 // This file is distributed under the University of Illinois/NCSA Open Source License.
3 // See LICENSE file in top directory for details.
4 //
5 // Copyright (c) 2016 Jeongnim Kim and QMCPACK developers.
6 //
7 // File developed by: Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign
8 // Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
9 // Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
10 // Raymond Clay III, j.k.rofling@gmail.com, Lawrence Livermore National Laboratory
11 // Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
12 //
13 // File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
14 //////////////////////////////////////////////////////////////////////////////////////
15
16
17 /** @file LRHandlerTemp.h
18 * @brief Define a LRHandler with two template parameters
19 */
20 #include "LRCoulombSingleton.h"
21 #if OHMMS_DIM == 3
22 #include "LongRange/EwaldHandler.h"
23 #include "LongRange/EwaldHandler3D.h"
24 #elif OHMMS_DIM == 2
25 #include "LongRange/TwoDEwaldHandler.h"
26 #endif
27 #include <numeric>
28 namespace qmcplusplus
29 {
30 //initialization of the static data
31 std::unique_ptr<LRCoulombSingleton::LRHandlerType> LRCoulombSingleton::CoulombHandler;
32 std::unique_ptr<LRCoulombSingleton::LRHandlerType> LRCoulombSingleton::CoulombDerivHandler;
33 LRCoulombSingleton::lr_type LRCoulombSingleton::this_lr_type = ESLER;
34 /** CoulombFunctor
35 *
36 * An example for a Func for LRHandlerTemp. Four member functions have to be provided
37 * - reset(T volume) : reset the normalization factor
38 * - operator()(T r, T rinv): return a value of the original function, e.g., 1.0/r
39 * - Fk(T k, T rc)
40 * - Xk(T k, T rc)
41 */
42 #if OHMMS_DIM == 3
43 template<class T = double>
44 struct CoulombFunctor
45 {
46 T NormFactor;
CoulombFunctorqmcplusplus::CoulombFunctor47 inline CoulombFunctor() {}
resetqmcplusplus::CoulombFunctor48 void reset(ParticleSet& ref) { NormFactor = 4.0 * M_PI / ref.LRBox.Volume; }
resetqmcplusplus::CoulombFunctor49 void reset(ParticleSet& ref, T rs) { NormFactor = 4.0 * M_PI / ref.LRBox.Volume; }
operator ()qmcplusplus::CoulombFunctor50 inline T operator()(T r, T rinv) { return rinv; }
dfqmcplusplus::CoulombFunctor51 inline T df(T r) { return -1.0 / (r * r); }
df2qmcplusplus::CoulombFunctor52 inline T df2(T r) { return 2.0 / (r * r * r); }
Vkqmcplusplus::CoulombFunctor53 inline T Vk(T k) { return NormFactor / (k * k); }
54
Xk_dkqmcplusplus::CoulombFunctor55 inline T Xk_dk(T k) { return 0.0; }
Fkqmcplusplus::CoulombFunctor56 inline T Fk(T k, T rc) { return NormFactor / (k * k) * std::cos(k * rc); }
Xkqmcplusplus::CoulombFunctor57 inline T Xk(T k, T rc) { return -NormFactor / (k * k) * std::cos(k * rc); }
58
dVk_dkqmcplusplus::CoulombFunctor59 inline T dVk_dk(T k) { return -2 * NormFactor / k / k / k; }
dFk_dkqmcplusplus::CoulombFunctor60 inline T dFk_dk(T k, T rc) { return -NormFactor / k / k * (2.0 / k * std::cos(k * rc) + rc * std::sin(k * rc)); }
61
dXk_dkqmcplusplus::CoulombFunctor62 inline T dXk_dk(T k, T rc) { return NormFactor / k / k * (2.0 / k * std::cos(k * rc) + rc * std::sin(k * rc)); }
63
integrate_r2qmcplusplus::CoulombFunctor64 inline T integrate_r2(T r) const { return 0.5 * r * r; }
65 };
66 #elif OHMMS_DIM == 2
67 template<class T = double>
68 struct CoulombFunctor
69 {
70 T NormFactor;
CoulombFunctorqmcplusplus::CoulombFunctor71 inline CoulombFunctor() {}
resetqmcplusplus::CoulombFunctor72 void reset(ParticleSet& ref) { NormFactor = 2.0 * M_PI / ref.LRBox.Volume; }
resetqmcplusplus::CoulombFunctor73 void reset(ParticleSet& ref, T rs) { NormFactor = 2.0 * M_PI / ref.LRBox.Volume; }
operator ()qmcplusplus::CoulombFunctor74 inline T operator()(T r, T rinv) { return rinv; }
dfqmcplusplus::CoulombFunctor75 inline T df(T r) { return -1.0 / (r * r); }
df2qmcplusplus::CoulombFunctor76 inline T df2(T r) { return 2 / (r * r * r); }
Fkqmcplusplus::CoulombFunctor77 inline T Fk(T k, T rc) { return NormFactor / k * std::cos(k * rc); }
Xkqmcplusplus::CoulombFunctor78 inline T Xk(T k, T rc) { return -NormFactor / k * std::cos(k * rc); }
79
80
integrate_r2qmcplusplus::CoulombFunctor81 inline T integrate_r2(T r) const { return 0.5 * r * r; }
82 };
83 #endif
84
85 template<class T = double>
86 struct PseudoCoulombFunctor
87 {
88 //typedef OneDimCubicSpline<T> RadialFunctorType;
89 typedef OneDimLinearSpline<T> RadialFunctorType;
90 RadialFunctorType& radFunc;
91 T NormFactor;
PseudoCoulombFunctorqmcplusplus::PseudoCoulombFunctor92 inline PseudoCoulombFunctor(RadialFunctorType& rfunc) : radFunc(rfunc) {}
resetqmcplusplus::PseudoCoulombFunctor93 void reset(ParticleSet& ref) { NormFactor = 4.0 * M_PI / ref.LRBox.Volume; }
operator ()qmcplusplus::PseudoCoulombFunctor94 inline T operator()(T r, T rinv) { return radFunc.splint(r); }
dfqmcplusplus::PseudoCoulombFunctor95 inline T df(T r)
96 {
97 T du, d2u;
98 radFunc.splint(r, du, d2u);
99 return du;
100 }
Fkqmcplusplus::PseudoCoulombFunctor101 inline T Fk(T k, T rc) { return NormFactor / (k * k) * std::cos(k * rc); }
Xkqmcplusplus::PseudoCoulombFunctor102 inline T Xk(T k, T rc) { return -NormFactor / (k * k) * std::cos(k * rc); }
integrate_r2qmcplusplus::PseudoCoulombFunctor103 inline T integrate_r2(T r) const
104 {
105 //fix this to return the integration
106 return 0.5 * r * r;
107 }
108 };
109
110
getHandler(ParticleSet & ref)111 std::unique_ptr<LRCoulombSingleton::LRHandlerType> LRCoulombSingleton::getHandler(ParticleSet& ref)
112 {
113 if (CoulombHandler == 0)
114 {
115 #if OHMMS_DIM == 3
116 if (ref.SK->SuperCellEnum == SUPERCELL_SLAB)
117 {
118 app_log() << "\n Creating CoulombHandler using quasi-2D Ewald method for the slab. " << std::endl;
119 CoulombHandler = std::make_unique<EwaldHandler>(ref);
120 }
121 else //if(ref.LRBox.SuperCellEnum == SUPERCELL_BULK)
122 {
123 if (this_lr_type == ESLER)
124 {
125 app_log() << "\n Creating CoulombHandler with the Esler Optimized Breakup. " << std::endl;
126 CoulombHandler = std::make_unique<LRHandlerTemp<CoulombFunctor<mRealType>, LPQHIBasis>>(ref);
127 }
128 else if (this_lr_type == EWALD)
129 {
130 app_log() << "\n Creating CoulombHandler with the 3D Ewald Breakup. " << std::endl;
131 CoulombHandler = std::make_unique<EwaldHandler3D>(ref);
132 }
133 else if (this_lr_type == NATOLI)
134 {
135 app_log() << "\n Creating CoulombHandler with the Natoli Optimized Breakup. " << std::endl;
136 CoulombHandler = std::make_unique<LRHandlerSRCoulomb<CoulombFunctor<mRealType>, LPQHISRCoulombBasis>>(ref);
137 }
138 else
139 {
140 APP_ABORT("\n Long range breakup method not recognized.\n");
141 }
142 }
143 // else if(ref.LRBox.SuperCellEnum == SUPERCELL_SLAB)
144 // {
145 // app_log() << "\n Creating CoulombHandler using quasi-2D Ewald method for the slab. " << std::endl;
146 // CoulombHandler= new EwaldHandler(ref);
147 // }
148 #elif OHMMS_DIM == 2
149 app_log() << "\n Creating CoulombHandler using 2D Ewald method. " << std::endl;
150 CoulombHandler = std::make_unique<TwoDEwaldHandler>(ref);
151 #endif
152 CoulombHandler->initBreakup(ref);
153 return std::unique_ptr<LRHandlerType>(CoulombHandler->makeClone(ref));
154 }
155 else
156 {
157 app_log() << " Clone CoulombHandler. " << std::endl;
158 return std::unique_ptr<LRHandlerType>(CoulombHandler->makeClone(ref));
159 }
160 }
161
getDerivHandler(ParticleSet & ref)162 std::unique_ptr<LRCoulombSingleton::LRHandlerType> LRCoulombSingleton::getDerivHandler(ParticleSet& ref)
163 {
164 #if OHMMS_DIM != 3
165 APP_ABORT("energy derivative implemented for 3D only");
166 #endif
167 //APP_ABORT("SR Coulomb Basis Handler has cloning issues. Stress also has some kinks");
168 if (CoulombDerivHandler == 0)
169 {
170 if (this_lr_type == EWALD)
171 {
172 app_log() << "\n Creating CoulombDerivHandler with the 3D Ewald Breakup. " << std::endl;
173 CoulombDerivHandler = std::make_unique<EwaldHandler3D>(ref);
174 }
175 else if (this_lr_type == NATOLI)
176 {
177 app_log() << "\n Creating CoulombDerivHandler with the Natoli Optimized Breakup. " << std::endl;
178 CoulombDerivHandler = std::make_unique<LRHandlerSRCoulomb<CoulombFunctor<mRealType>, LPQHISRCoulombBasis>>(ref);
179 }
180 else if (this_lr_type == ESLER)
181 {
182 APP_ABORT("\n Derivatives are not supported with Esler Optimized Breakup.\n");
183 }
184 else
185 {
186 APP_ABORT("\n Long range breakup method for derivatives not recognized.\n");
187 }
188 CoulombDerivHandler->initBreakup(ref);
189 return std::unique_ptr<LRHandlerType>(CoulombDerivHandler->makeClone(ref));
190 }
191 else
192 {
193 app_log() << " Clone CoulombDerivHandler. " << std::endl;
194 return std::unique_ptr<LRHandlerType>(CoulombDerivHandler->makeClone(ref));
195 }
196 }
197
198 template<typename T>
createSpline4RbyVs_temp(LRHandlerBase * aLR,T rcut,LinearGrid<T> * agrid)199 OneDimCubicSpline<T>* createSpline4RbyVs_temp(LRHandlerBase* aLR, T rcut, LinearGrid<T>* agrid)
200 {
201 typedef OneDimCubicSpline<T> func_type;
202 if (agrid == nullptr)
203 {
204 agrid = new LinearGrid<T>;
205 agrid->set(0.0, rcut, 1001);
206 }
207 int ng = agrid->size();
208 std::vector<T> v(ng);
209 T r = (*agrid)[0];
210 //check if the first point is not zero
211 v[0] = (r > std::numeric_limits<T>::epsilon()) ? r * aLR->evaluate(r, 1.0 / r) : 0.0;
212 for (int ig = 1; ig < ng - 1; ig++)
213 {
214 r = (*agrid)[ig];
215 v[ig] = r * aLR->evaluate(r, 1.0 / r);
216 }
217 v[0] = 2.0 * v[1] - v[2];
218 v[ng - 1] = 0.0;
219 func_type* V0 = new func_type(agrid, v);
220 T deriv = (v[1] - v[0]) / ((*agrid)[1] - (*agrid)[0]);
221 V0->spline(0, deriv, ng - 1, 0.0);
222 return V0;
223 }
224
225 template<typename T>
createSpline4RbyVsDeriv_temp(LRHandlerBase * aLR,T rcut,LinearGrid<T> * agrid)226 OneDimCubicSpline<T>* createSpline4RbyVsDeriv_temp(LRHandlerBase* aLR, T rcut, LinearGrid<T>* agrid)
227 {
228 typedef OneDimCubicSpline<T> func_type;
229 if (agrid == nullptr)
230 {
231 agrid = new LinearGrid<T>;
232 agrid->set(0.0, rcut, 1001);
233 }
234 int ng = agrid->size();
235 std::vector<T> v(ng);
236 T r = (*agrid)[0];
237 //check if the first point is not zero
238 v[0] = (r > std::numeric_limits<T>::epsilon()) ? r * aLR->evaluate(r, 1.0 / r) : 0.0;
239 T v_val(0.0);
240 T v_deriv(0.0);
241
242 for (int ig = 1; ig < ng - 1; ig++)
243 {
244 r = (*agrid)[ig];
245 v_val = aLR->evaluate(r, 1.0 / r);
246 v_deriv = aLR->srDf(r, 1 / r);
247
248 v[ig] = v_val + v_deriv * r;
249 }
250 v[0] = 2.0 * v[1] - v[2];
251 v[ng - 1] = 0.0;
252 func_type* dV0 = new func_type(agrid, v);
253 T deriv = (v[1] - v[0]) / ((*agrid)[1] - (*agrid)[0]);
254 dV0->spline(0, deriv, ng - 1, 0.0);
255 return dV0;
256 }
257
258
createSpline4RbyVs(LRHandlerType * aLR,mRealType rcut,GridType * agrid)259 LRCoulombSingleton::RadFunctorType* LRCoulombSingleton::createSpline4RbyVs(LRHandlerType* aLR,
260 mRealType rcut,
261 GridType* agrid)
262 {
263 return createSpline4RbyVs_temp(aLR, static_cast<pRealType>(rcut), agrid);
264 }
265
createSpline4RbyVsDeriv(LRHandlerType * aLR,mRealType rcut,GridType * agrid)266 LRCoulombSingleton::RadFunctorType* LRCoulombSingleton::createSpline4RbyVsDeriv(LRHandlerType* aLR,
267 mRealType rcut,
268 GridType* agrid)
269 {
270 return createSpline4RbyVsDeriv_temp(aLR, static_cast<pRealType>(rcut), agrid);
271 }
272
273
274 } // namespace qmcplusplus
275