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) 2020 QMCPACK developers.
6 //
7 // File developed by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
8 //                    Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
9 //                    Raymond Clay III, j.k.rofling@gmail.com, Lawrence Livermore National Laboratory
10 //                    Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
11 //
12 // File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
13 //////////////////////////////////////////////////////////////////////////////////////
14 
15 
16 /** @file LRHandlerBase.h
17  * @brief Define LRHandlerBase and DummyLRHandler<typename Func>
18  */
19 #ifndef QMCPLUSPLUS_LRHANLDERBASE_AND_DUMMY_H
20 #define QMCPLUSPLUS_LRHANLDERBASE_AND_DUMMY_H
21 
22 #include "coulomb_types.h"
23 #include "LongRange/StructFact.h"
24 
25 namespace qmcplusplus
26 {
27 /** base class for LRHandlerTemp<FUNC,BASIS> and DummyLRHanlder<typename Func>
28  */
29 struct LRHandlerBase
30 {
31   DECLARE_COULOMB_TYPES
32 
33   /// Maxkimum Kshell for the given Kc
34   int MaxKshell;
35   /// Maximum k cutoff
36   mRealType LR_kc;
37   /// Maximum r cutoff
38   mRealType LR_rc;
39   ///Fourier component for all the k-point
40   Vector<mRealType> Fk;
41   ///Fourier component of the LR part, fit to optimize the gradients.
42   Vector<mRealType> Fkg;
43   ///Fourier component of the LR part of strain tensor, by optimized breakup.
44   std::vector<SymTensor<mRealType, OHMMS_DIM>> dFk_dstrain;
45   ///Vector of df_k/dk, fit as to optimize strains.
46   Vector<mRealType> Fkgstrain;
47   ///Fourier component for each k-shell
48   Vector<mRealType> Fk_symm;
49 
50   ///Fourier component for each k-shell
51   ///Coefficient
52   std::vector<mRealType> coefs;
53   ///Coefficient for gradient fit.
54   std::vector<mRealType> gcoefs;
55   ///Coefficient for strain fit.
56   std::vector<mRealType> gstraincoefs;
57 
58   virtual mRealType evaluate_vlr_k(mRealType k) = 0;
59 
60 
61   //constructor
LRHandlerBaseLRHandlerBase62   explicit LRHandlerBase(mRealType kc) : MaxKshell(0), LR_kc(kc), LR_rc(0), ClassName("LRHandlerBase") {}
63 
64   // virtual destructor
~LRHandlerBaseLRHandlerBase65   virtual ~LRHandlerBase() {}
66 
67   //return r cutoff
get_rcLRHandlerBase68   inline mRealType get_rc() const { return LR_rc; }
69   //return k cutoff
get_kcLRHandlerBase70   inline mRealType get_kc() const { return LR_kc; }
71 
72   /** evaluate \f$\sum_k F_{k} \rho^1_{-{\bf k} \rho^2_{\bf k}\f$
73    * @param kshell degeneracies of the vectors
74    * @param rk1 starting address of \f$\rho^1_{{\bf k}\f$
75    * @param rk2 starting address of \f$\rho^2_{{\bf k}\f$
76    *
77    * Valid for the strictly ordered k and \f$F_{k}\f$.
78    */
evaluateLRHandlerBase79   inline mRealType evaluate(const std::vector<int>& kshell,
80                             const pComplexType* restrict rk1,
81                             const pComplexType* restrict rk2)
82   {
83     mRealType vk = 0.0;
84     for (int ks = 0, ki = 0; ks < MaxKshell; ks++)
85     {
86       mRealType u = 0;
87       for (; ki < kshell[ks + 1]; ki++, rk1++, rk2++)
88         u += ((*rk1).real() * (*rk2).real() + (*rk1).imag() * (*rk2).imag());
89       vk += Fk_symm[ks] * u;
90     }
91     return vk;
92   }
93 
evaluate_w_skLRHandlerBase94   inline mRealType evaluate_w_sk(const std::vector<int>& kshell, const pRealType* restrict sk)
95   {
96     mRealType vk = 0.0;
97     for (int ks = 0, ki = 0; ks < MaxKshell; ks++)
98     {
99       mRealType u = 0;
100       for (; ki < kshell[ks + 1]; ki++)
101         u += (*sk++);
102       vk += Fk_symm[ks] * u;
103     }
104     return vk;
105   }
106 
evaluateLRHandlerBase107   inline mRealType evaluate(const std::vector<int>& kshell,
108                             const pRealType* restrict rk1_r,
109                             const pRealType* restrict rk1_i,
110                             const pRealType* restrict rk2_r,
111                             const pRealType* restrict rk2_i)
112   {
113     mRealType vk = 0.0;
114     for (int ks = 0, ki = 0; ks < MaxKshell; ks++)
115     {
116       mRealType u = 0;
117       for (; ki < kshell[ks + 1]; ki++)
118         u += ((*rk1_r++) * (*rk2_r++) + (*rk1_i++) * (*rk2_i++));
119       vk += Fk_symm[ks] * u;
120     }
121     return vk;
122   }
123 
124   /** Evaluate the long-range potential with the open BC for the D-1 direction */
evaluate_slabLRHandlerBase125   virtual mRealType evaluate_slab(pRealType z,
126                                   const std::vector<int>& kshell,
127                                   const pComplexType* restrict eikr_i,
128                                   const pComplexType* restrict eikr_j)
129   {
130     return 0.0;
131   }
132 
evaluateLRHandlerBase133   inline mRealType evaluate(const std::vector<int>& kshell, int iat, const pComplexType* restrict rk2, ParticleSet& P)
134   {
135     mRealType vk = 0.0;
136 #if !defined(USE_REAL_STRUCT_FACTOR)
137     for (int ks = 0, ki = 0; ks < MaxKshell; ks++)
138     {
139       mRealType u = 0;
140       for (; ki < kshell[ks + 1]; ki++, rk2++)
141       {
142         pComplexType eikr = P.SK->eikr(iat, ki);
143         u += eikr.real() * (*rk2).real() + eikr.imag() * (*rk2).imag();
144       }
145       vk += Fk_symm[ks] * u;
146     }
147 #endif
148     return vk;
149   }
150 
evaluateLRHandlerBase151   inline mRealType evaluate(const std::vector<int>& kshell,
152                             int iat,
153                             const pRealType* restrict rk2_r,
154                             const pRealType* restrict rk2_i,
155                             ParticleSet& P)
156   {
157     mRealType vk = 0.0;
158 #if defined(USE_REAL_STRUCT_FACTOR)
159     const pRealType* restrict eikr_r = P.SK->eikr_r[iat];
160     const pRealType* restrict eikr_i = P.SK->eikr_i[iat];
161     for (int ks = 0, ki = 0; ks < MaxKshell; ks++)
162     {
163       mRealType u = 0;
164       for (; ki < kshell[ks + 1]; ki++)
165         u += eikr_r[ki] * (*rk2_r++) + eikr_i[ki] * (*rk2_i++);
166       vk += Fk_symm[ks] * u;
167     }
168 #endif
169     return vk;
170   }
171 
172   /** evaluate \f$\sum_k F_{k} \rho^1_{-{\bf k} \rho^2_{\bf k}\f$
173    * and \f$\sum_k F_{k} \rho^1_{-{\bf k} \rho^2_{\bf k}\f$
174    * @param kshell degeneracies of the vectors
175    * @param rk1 starting address of \f$\rho^1_{{\bf k}\f$
176    * @param rk2 starting address of \f$\rho^2_{{\bf k}\f$
177    *
178    * Valid for the strictly ordered k and \f$F_{k}\f$.
179    */
evaluateGradLRHandlerBase180   inline void evaluateGrad(const ParticleSet& A,
181                            const ParticleSet& B,
182                            int specB,
183                            std::vector<pRealType>& Zat,
184                            std::vector<TinyVector<pRealType, OHMMS_DIM>>& grad1)
185   {
186 #if !defined(USE_REAL_STRUCT_FACTOR)
187     const Matrix<pComplexType>& e2ikrA = A.SK->eikr;
188     const pComplexType* rhokB          = B.SK->rhok[specB];
189     const std::vector<PosType>& kpts   = A.SK->KLists.kpts_cart;
190     for (int ki = 0; ki < Fk.size(); ki++)
191     {
192       PosType k = kpts[ki];
193       for (int iat = 0; iat < Zat.size(); iat++)
194       {
195         grad1[iat] -= Zat[iat] * k * Fkg[ki] *
196             (e2ikrA(iat, ki).real() * rhokB[ki].imag() - e2ikrA(iat, ki).imag() * rhokB[ki].real());
197       }
198     }
199 #else
200 
201     const Matrix<pRealType>& e2ikrA_r = A.SK->eikr_r;
202     const Matrix<pRealType>& e2ikrA_i = A.SK->eikr_i;
203     const pRealType* rhokB_r          = B.SK->rhok_r[specB];
204     const pRealType* rhokB_i          = B.SK->rhok_i[specB];
205     const std::vector<PosType>& kpts  = A.SK->KLists.kpts_cart;
206     for (int ki = 0; ki < Fk.size(); ki++)
207     {
208       PosType k = kpts[ki];
209       for (int iat = 0; iat < Zat.size(); iat++)
210       {
211         grad1[iat] -= Zat[iat] * k * Fkg[ki] * (e2ikrA_r(iat, ki) * rhokB_i[ki] - e2ikrA_i(iat, ki) * rhokB_r[ki]);
212       }
213     }
214 #endif
215   }
216 
217   ///FIX_PRECISION
evaluateStressLRHandlerBase218   inline SymTensor<pRealType, OHMMS_DIM> evaluateStress(const std::vector<int>& kshell,
219                                                         const pRealType* rhokA_r,
220                                                         const pRealType* rhokA_i,
221                                                         const pRealType* rhokB_r,
222                                                         const pRealType* rhokB_i)
223   {
224     SymTensor<pRealType, OHMMS_DIM> stress;
225     for (int ki = 0; ki < dFk_dstrain.size(); ki++)
226     {
227       stress += (rhokA_r[ki] * rhokB_r[ki] + rhokA_i[ki] * rhokB_i[ki]) * dFk_dstrain[ki];
228     }
229 
230     return stress;
231   }
232 
233   ///FIX_PRECISION
evaluateStressLRHandlerBase234   inline SymTensor<pRealType, OHMMS_DIM> evaluateStress(const std::vector<int>& kshell,
235                                                         const pComplexType* rhokA,
236                                                         const pComplexType* rhokB)
237   {
238     SymTensor<pRealType, OHMMS_DIM> stress;
239     for (int ki = 0; ki < dFk_dstrain.size(); ki++)
240     {
241       stress += (rhokA[ki].real() * rhokB[ki].real() + rhokA[ki].imag() * rhokB[ki].imag()) * dFk_dstrain[ki];
242     }
243     return stress;
244   }
245 
246   /** evaluate \f$ v_{s}(k=0) = \frac{4\pi}{V}\int_0^{r_c} r^2 v_s(r) dr \f$
247    */
evaluateSR_k0LRHandlerBase248   virtual mRealType evaluateSR_k0() { return 0.0; }
249   /** evaluate \f$ v_s(r=0) \f$ for the self-interaction term
250    */
evaluateLR_r0LRHandlerBase251   virtual mRealType evaluateLR_r0() { return 0.0; }
252 
253   ///These functions return the strain derivatives of all corresponding quantities
254   /// in total energy.  See documentation (forthcoming).
evaluateLR_r0_dstrainLRHandlerBase255   virtual SymTensor<mRealType, OHMMS_DIM> evaluateLR_r0_dstrain() { return 0; };
evaluateSR_k0_dstrainLRHandlerBase256   virtual SymTensor<mRealType, OHMMS_DIM> evaluateSR_k0_dstrain() { return 0; };
evaluateLR_dstrainLRHandlerBase257   virtual SymTensor<mRealType, OHMMS_DIM> evaluateLR_dstrain(TinyVector<pRealType, OHMMS_DIM> k, pRealType kmag)
258   {
259     return 0;
260   };
evaluateSR_dstrainLRHandlerBase261   virtual SymTensor<mRealType, OHMMS_DIM> evaluateSR_dstrain(TinyVector<pRealType, OHMMS_DIM> r, pRealType rmag)
262   {
263     return 0;
264   };
265 
266   virtual void initBreakup(ParticleSet& ref)              = 0;
267   virtual void Breakup(ParticleSet& ref, mRealType rs_in) = 0;
268   virtual void resetTargetParticleSet(ParticleSet& ref)   = 0;
269 
270   virtual mRealType evaluate(mRealType r, mRealType rinv) = 0;
271   virtual mRealType evaluateLR(mRealType r)               = 0;
272   virtual mRealType srDf(mRealType r, mRealType rinv)     = 0;
273 
lrDfLRHandlerBase274   virtual mRealType lrDf(mRealType r)
275   {
276     APP_ABORT("Error: lrDf(r) is not implemented in " + ClassName + "\n");
277     return 0.0;
278   };
279 
280   /** make clone */
281   virtual LRHandlerBase* makeClone(ParticleSet& ref) = 0;
282 
283 protected:
284   std::string ClassName;
285 };
286 
287 /** LRHandler without breakup.
288  *
289  * The template parameter Func should impelement operator()(kk) which
290  * returns the Fourier component of a long-range function. Here kk
291  * is \f$|{\bf k}|^2\f$.
292  */
293 template<class Func>
294 struct DummyLRHandler : public LRHandlerBase
295 {
296   Func myFunc;
297 
LRHandlerBaseDummyLRHandler298   DummyLRHandler(mRealType kc = -1.0) : LRHandlerBase(kc) {}
299 
DummyLRHandlerDummyLRHandler300   DummyLRHandler(const DummyLRHandler& aLR) : LRHandlerBase(aLR) {}
301 
initBreakupDummyLRHandler302   void initBreakup(ParticleSet& ref) override
303   {
304     mRealType norm = 4.0 * M_PI / ref.Lattice.Volume;
305     mRealType kcsq = LR_kc * LR_kc;
306     KContainer& KList(ref.SK->KLists);
307     int maxshell = KList.kshell.size() - 1;
308     const KContainer::SContainer_t& kk(KList.ksq);
309     int ksh = 0, ik = 0;
310     while (ksh < maxshell)
311     {
312       if (kk[ik] > kcsq)
313         break; //exit
314       ik = KList.kshell[++ksh];
315     }
316     MaxKshell = ksh;
317     Fk_symm.resize(MaxKshell);
318     Fk.resize(KList.kpts_cart.size());
319     for (ksh = 0, ik = 0; ksh < MaxKshell; ksh++, ik++)
320     {
321       mRealType v  = norm * myFunc(kk[KList.kshell[ksh]]); //rpa=u0/kk[ik];
322       Fk_symm[ksh] = v;
323       for (; ik < KList.kshell[ksh + 1]; ik++)
324         Fk[ik] = v;
325     }
326   }
327 
evaluate_vlr_kDummyLRHandler328   mRealType evaluate_vlr_k(mRealType k) override { return 0.0; }
evaluateDummyLRHandler329   mRealType evaluate(mRealType r, mRealType rinv) override { return 0.0; }
evaluateLRDummyLRHandler330   mRealType evaluateLR(mRealType r) override { return 0.0; }
srDfDummyLRHandler331   mRealType srDf(mRealType r, mRealType rinv) override { return 0.0; }
BreakupDummyLRHandler332   void Breakup(ParticleSet& ref, mRealType rs_in) override {}
resetTargetParticleSetDummyLRHandler333   void resetTargetParticleSet(ParticleSet& ref) override {}
makeCloneDummyLRHandler334   virtual LRHandlerBase* makeClone(ParticleSet& ref) override { return new DummyLRHandler<Func>(LR_kc); }
335 };
336 
337 } // namespace qmcplusplus
338 #endif
339