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: Raymond Clay III, j.k.rofling@gmail.com, Lawrence Livermore National Laboratory
8 //
9 // File created by: Raymond Clay III, j.k.rofling@gmail.com, Lawrence Livermore National Laboratory
10 //////////////////////////////////////////////////////////////////////////////////////
11
12
13 #include "EwaldHandler3D.h"
14
15 namespace qmcplusplus
16 {
initBreakup(ParticleSet & ref)17 void EwaldHandler3D::initBreakup(ParticleSet& ref)
18 {
19 SuperCellEnum = ref.Lattice.SuperCellEnum;
20 LR_rc = ref.Lattice.LR_rc;
21 LR_kc = ref.Lattice.LR_kc;
22
23 // Sigma=3.5;
24 //We provide two means of choosing sigma here...
25 //
26 //This condition on Sigma is based on the real space cutoff of the potential at r_c for the potential.
27 //while(erfc(Sigma)/LR_rc>1e-10)
28 //
29 //This condition on Sigma is based on the magnitude of the force at r_c for the potential.
30 // while( (erfc(Sigma)+std::exp(-Sigma*Sigma)*2*Sigma/std::sqrt(M_PI))/LR_rc/LR_rc>1e-14)
31 // {
32 // Sigma+=0.1;
33 // }
34 //
35 // app_log() << " EwaldHandler3D Sigma/LR_rc = " << Sigma ;
36 // Sigma/=ref.Lattice.LR_rc;
37
38 //This heuristic for choosing Sigma is from the 1992 Natoli Ceperley Optimized Breakup Paper.
39 Sigma = std::sqrt(LR_kc / (2.0 * LR_rc));
40 app_log() << " Sigma=" << Sigma << std::endl;
41 Volume = ref.Lattice.Volume;
42 PreFactors = 0.0;
43 fillFk(ref.SK->KLists);
44 fillYkgstrain(ref.SK->KLists);
45 filldFk_dk(ref.SK->KLists);
46 }
47
EwaldHandler3D(const EwaldHandler3D & aLR,ParticleSet & ref)48 EwaldHandler3D::EwaldHandler3D(const EwaldHandler3D& aLR, ParticleSet& ref)
49 : LRHandlerBase(aLR), Sigma(aLR.Sigma), Volume(aLR.Volume), Area(aLR.Area), PreFactors(aLR.PreFactors)
50 {
51 SuperCellEnum = aLR.SuperCellEnum;
52 }
53
fillFk(KContainer & KList)54 void EwaldHandler3D::fillFk(KContainer& KList)
55 {
56 Fk.resize(KList.kpts_cart.size());
57 Fkg.resize(KList.kpts_cart.size());
58 const std::vector<int>& kshell(KList.kshell);
59
60 MaxKshell = kshell.size() - 1;
61
62 Fk_symm.resize(MaxKshell);
63 kMag.resize(MaxKshell);
64 mRealType kgauss = 1.0 / (4 * Sigma * Sigma);
65 mRealType knorm = 4 * M_PI / Volume;
66 for (int ks = 0, ki = 0; ks < Fk_symm.size(); ks++)
67 {
68 mRealType t2e = KList.ksq[ki] * kgauss;
69 mRealType uk = knorm * std::exp(-t2e) / KList.ksq[ki];
70 Fk_symm[ks] = uk;
71 while (ki < KList.kshell[ks + 1] && ki < Fk.size())
72 Fk[ki++] = uk;
73 }
74
75 for (int ki = 0; ki < Fk.size(); ki++)
76 Fkg[ki] = Fk[ki];
77
78 PreFactors[3] = 0.0;
79 app_log().flush();
80 }
81
evaluate_vlr_k(mRealType k)82 EwaldHandler3D::mRealType EwaldHandler3D::evaluate_vlr_k(mRealType k)
83 {
84 mRealType kgauss = 1.0 / (4 * Sigma * Sigma);
85 mRealType knorm = 4 * M_PI / Volume;
86 mRealType k2 = k * k;
87 return knorm * std::exp(-k2 * kgauss) / k2;
88 }
89
90 } // namespace qmcplusplus
91