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 // Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory 10 // 11 // File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign 12 ////////////////////////////////////////////////////////////////////////////////////// 13 14 15 /** @file LRHandlerTemp.h 16 * @brief Define a LRHandler with two template parameters 17 */ 18 #ifndef QMCPLUSPLUS_LRRPAHANLDERTEMP_H 19 #define QMCPLUSPLUS_LRRPAHANLDERTEMP_H 20 21 #include "LongRange/LRHandlerBase.h" 22 #include "LongRange/LPQHIBasis.h" 23 #include "LongRange/LRBreakup.h" 24 25 namespace qmcplusplus 26 { 27 /* Templated LRHandler class 28 * 29 * LRHandlerTemp<Func,BreakupBasis> is a modification of LRHandler 30 * and a derived class from LRHanlderBase. 31 * The first template parameter Func is a generic functor, e.g., CoulombFunctor. 32 * The second template parameter is a BreakupBasis and the default is set to LPQHIBasis. 33 * LRHandlerBase is introduced to enable run-time options. See RPAContstraints.h 34 */ 35 template<class Func, class BreakupBasis = LPQHIBasis> 36 struct LRRPAHandlerTemp : public LRHandlerBase 37 { 38 DECLARE_COULOMB_TYPES 39 40 //Typedef for the lattice-type. 41 typedef ParticleSet::ParticleLayout_t ParticleLayout_t; 42 typedef BreakupBasis BreakupBasisType; 43 44 bool FirstTime; 45 mRealType rs; 46 mRealType kc; 47 BreakupBasisType Basis; //This needs a Lattice for the constructor... 48 Func myFunc; 49 50 //Constructor LRHandlerBaseLRRPAHandlerTemp51 LRRPAHandlerTemp(ParticleSet& ref, mRealType kc_in = -1.0) : LRHandlerBase(kc_in), FirstTime(true), Basis(ref.Lattice) 52 { 53 LRHandlerBase::ClassName = "LRRPAHandlerTemp"; 54 myFunc.reset(ref); 55 } 56 57 //LRHandlerTemp(ParticleSet& ref, mRealType rs, mRealType kc=-1.0): LRHandlerBase(kc), Basis(ref.Lattice) 58 //{ 59 // myFunc.reset(ref,rs); 60 //} 61 62 /** "copy" constructor 63 * @param aLR LRHandlerTemp 64 * @param ref Particleset 65 * 66 * Copy the content of aLR 67 * References to ParticleSet or ParticleLayoutout_t are not copied. 68 */ LRRPAHandlerTempLRRPAHandlerTemp69 LRRPAHandlerTemp(const LRRPAHandlerTemp& aLR, ParticleSet& ref) 70 : LRHandlerBase(aLR), FirstTime(true), Basis(aLR.Basis, ref.Lattice) 71 { 72 myFunc.reset(ref); 73 fillFk(ref.SK->KLists); 74 } 75 makeCloneLRRPAHandlerTemp76 LRHandlerBase* makeClone(ParticleSet& ref) override { return new LRRPAHandlerTemp<Func, BreakupBasis>(*this, ref); } 77 initBreakupLRRPAHandlerTemp78 void initBreakup(ParticleSet& ref) override 79 { 80 InitBreakup(ref.Lattice, 1); 81 fillFk(ref.SK->KLists); 82 LR_rc = Basis.get_rc(); 83 } 84 BreakupLRRPAHandlerTemp85 void Breakup(ParticleSet& ref, mRealType rs_ext) override 86 { 87 //ref.Lattice.Volume=ref.getTotalNum()*4.0*M_PI/3.0*rs*rs*rs; 88 if (rs_ext > 0) 89 rs = rs_ext; 90 else 91 rs = std::pow(3.0 / 4.0 / M_PI * ref.Lattice.Volume / static_cast<mRealType>(ref.getTotalNum()), 1.0 / 3.0); 92 myFunc.reset(ref, rs); 93 InitBreakup(ref.Lattice, 1); 94 fillFk(ref.SK->KLists); 95 LR_rc = Basis.get_rc(); 96 } 97 resetTargetParticleSetLRRPAHandlerTemp98 void resetTargetParticleSet(ParticleSet& ref) override { myFunc.reset(ref); } 99 resetTargetParticleSetLRRPAHandlerTemp100 void resetTargetParticleSet(ParticleSet& ref, mRealType rs) { myFunc.reset(ref, rs); } 101 evaluateLRRPAHandlerTemp102 inline mRealType evaluate(mRealType r, mRealType rinv) override 103 { 104 mRealType v = 0.0; 105 for (int n = 0; n < coefs.size(); n++) 106 v += coefs[n] * Basis.h(n, r); 107 return v; 108 } 109 110 /** evaluate the first derivative of the short range part at r 111 * 112 * @param r radius 113 * @param rinv 1/r 114 */ srDfLRRPAHandlerTemp115 inline mRealType srDf(mRealType r, mRealType rinv) override 116 { 117 mRealType df = 0.0; 118 //mRealType df = myFunc.df(r, rinv); 119 for (int n = 0; n < coefs.size(); n++) 120 df += coefs[n] * Basis.dh_dr(n, r); 121 return df; 122 } 123 124 /** evaluate the contribution from the long-range part for for spline 125 */ evaluateLRLRRPAHandlerTemp126 inline mRealType evaluateLR(mRealType r) override 127 { 128 mRealType vk = 0.0; 129 return vk; 130 // for(int n=0; n<coefs.size(); n++) v -= coefs[n]*Basis.h(n,r); 131 } 132 133 /** evaluate \f$\sum_k F_{k} \rho^1_{-{\bf k} \rho^2_{\bf k}\f$ 134 * @param kshell degeneracies of the vectors 135 * @param rk1 starting address of \f$\rho^1_{{\bf k}\f$ 136 * @param rk2 starting address of \f$\rho^2_{{\bf k}\f$ 137 * 138 * Valid for the strictly ordered k and \f$F_{k}\f$. 139 */ evaluateLRRPAHandlerTemp140 inline mRealType evaluate(const std::vector<int>& kshell, 141 const pComplexType* restrict rk1, 142 const pComplexType* restrict rk2) 143 { 144 mRealType vk = 0.0; 145 for (int ks = 0, ki = 0; ks < MaxKshell; ks++) 146 { 147 mRealType u = 0; 148 for (; ki < kshell[ks + 1]; ki++, rk1++, rk2++) 149 u += ((*rk1).real() * (*rk2).real() + (*rk1).imag() * (*rk2).imag()); 150 vk += Fk_symm[ks] * u; 151 } 152 //for(int ki=0; ki<Fk.size(); ki++) { 153 // //vk += (rk1[ki]*rk2[minusk[ki]]).real()*Fk[ki]; 154 // vk += (rk1[ki].real()*rk2[ki].real()+rk1[ki].imag()*rk2[ki].imag())*Fk[ki]; 155 //} //ki 156 return vk; 157 } 158 159 // use what is put in fillFk. Multiplies evalFk by -1 evaluate_vlr_kLRRPAHandlerTemp160 inline mRealType evaluate_vlr_k(mRealType k) override { return -1.0 * evalFk(k); } 161 162 private: evalFkLRRPAHandlerTemp163 inline mRealType evalFk(mRealType k) 164 { 165 //FatK = 4.0*M_PI/(Basis.get_CellVolume()*k*k)* std::cos(k*Basis.get_rc()); 166 mRealType FatK = myFunc.Fk(k, Basis.get_rc()); 167 for (int n = 0; n < Basis.NumBasisElem(); n++) 168 FatK += coefs[n] * Basis.c(n, k); 169 return FatK; 170 } 171 evalXkLRRPAHandlerTemp172 inline mRealType evalXk(mRealType k) 173 { 174 //mRealType FatK; 175 //FatK = -4.0*M_PI/(Basis.get_CellVolume()*k*k)* std::cos(k*Basis.get_rc()); 176 //return (FatK); 177 return myFunc.Xk(k, Basis.get_rc()); 178 } 179 180 /** Initialise the basis and coefficients for the long-range beakup. 181 * 182 * We loocally create a breakup handler and pass in the basis 183 * that has been initialised here. We then discard the handler, leaving 184 * basis and coefs in a usable state. 185 * This method can be re-called later if lattice changes shape. 186 */ InitBreakupLRRPAHandlerTemp187 void InitBreakup(ParticleLayout_t& ref, int NumFunctions) 188 { 189 //First we send the new Lattice to the Basis, in case it has been updated. 190 Basis.set_Lattice(ref); 191 //Compute RC from box-size - in constructor? 192 //No here...need update if box changes 193 int NumKnots(15); 194 Basis.set_NumKnots(NumKnots); 195 Basis.set_rc(ref.LR_rc); 196 //Initialise the breakup - pass in basis. 197 LRBreakup<BreakupBasis> breakuphandler(Basis); 198 // std::cout<<" finding kc: "<<ref.LR_kc<<" , "<<LR_kc<< std::endl; 199 //Find size of basis from cutoffs 200 kc = (LR_kc < 0) ? ref.LR_kc : LR_kc; 201 //mRealType kc(ref.LR_kc); //User cutoff parameter... 202 //kcut is the cutoff for switching to approximate k-point degeneracies for 203 //better performance in making the breakup. A good bet is 30*K-spacing so that 204 //there are 30 "boxes" in each direction that are treated with exact degeneracies. 205 //Assume orthorhombic cell just for deriving this cutoff - should be insensitive. 206 //K-Spacing = (kpt_vol)**1/3 = 2*pi/(cellvol**1/3) 207 mRealType kcut = (30.0) * 2 * M_PI * std::pow(Basis.get_CellVolume(), -1.0 / 3.0); 208 //Use 3000/LMax here...==6000/rc for non-ortho cells 209 mRealType kmax(6000.0 / ref.LR_rc); 210 // std::cout<<"K_STATS !!! "<<kcut<<" "<<kmax<<std::endl; 211 MaxKshell = static_cast<int>(breakuphandler.SetupKVecs(kc, kcut, kmax)); 212 if (FirstTime) 213 { 214 app_log() << " finding kc: " << ref.LR_kc << " , " << LR_kc << std::endl; 215 app_log() << " LRBreakp parameter Kc =" << kc << std::endl; 216 app_log() << " Continuum approximation in k = [" << kcut << "," << kmax << ")" << std::endl; 217 FirstTime = false; 218 } 219 //Set up x_k 220 //This is the FT of -V(r) from r_c to infinity. 221 //This is the only data that the breakup handler needs to do the breakup. 222 //We temporarily store it in Fk, which is replaced with the full FT (0->inf) 223 //of V_l(r) after the breakup has been done. 224 fillXk(breakuphandler.KList); 225 //Allocate the space for the coefficients. 226 coefs.resize(Basis.NumBasisElem()); //This must be after SetupKVecs. 227 breakuphandler.DoBreakup(Fk.data(), coefs.data()); //Fill array of coefficients. 228 } 229 fillXkLRRPAHandlerTemp230 void fillXk(std::vector<TinyVector<mRealType, 2>>& KList) 231 { 232 Fk.resize(KList.size()); 233 for (int ki = 0; ki < KList.size(); ki++) 234 { 235 mRealType k = KList[ki][0]; 236 Fk[ki] = evalXk(k); //Call derived fn. 237 } 238 } 239 fillFkLRRPAHandlerTemp240 void fillFk(KContainer& KList) 241 { 242 Fk.resize(KList.kpts_cart.size()); 243 const std::vector<int>& kshell(KList.kshell); 244 if (MaxKshell >= kshell.size()) 245 MaxKshell = kshell.size() - 1; 246 Fk_symm.resize(MaxKshell); 247 // std::cout<<"Filling FK :"<<std::endl; 248 for (int ks = 0, ki = 0; ks < Fk_symm.size(); ks++) 249 { 250 mRealType k = std::pow(KList.ksq[ki], 0.5); 251 mRealType uk = -1.0 * evalFk(k); 252 Fk_symm[ks] = uk; 253 // std::cout<<uk<<std::endl; 254 while (ki < KList.kshell[ks + 1] && ki < Fk.size()) 255 Fk[ki++] = uk; 256 } 257 //for(int ki=0; ki<KList.kpts_cart.size(); ki++){ 258 // mRealType k=dot(KList.kpts_cart[ki],KList.kpts_cart[ki]); 259 // k=std::sqrt(k); 260 // Fk[ki] = evalFk(k); //Call derived fn. 261 //} 262 } 263 }; 264 } // namespace qmcplusplus 265 #endif 266