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