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 #include "EwaldHandler.h"
16
17 namespace qmcplusplus
18 {
initBreakup(ParticleSet & ref)19 void EwaldHandler::initBreakup(ParticleSet& ref)
20 {
21 SuperCellEnum = ref.Lattice.SuperCellEnum;
22 LR_rc = ref.Lattice.LR_rc;
23 LR_kc = ref.Lattice.LR_kc;
24 Sigma = 3.5;
25 //determine the sigma
26 while (erfc(Sigma) / LR_rc > 1e-10)
27 {
28 Sigma += 0.1;
29 }
30 app_log() << " EwaldHandler Sigma/LR_rc = " << Sigma;
31 Sigma /= ref.Lattice.LR_rc;
32 app_log() << " Sigma=" << Sigma << std::endl;
33 Volume = ref.Lattice.Volume;
34 PreFactors = 0.0;
35 //See A.18
36 if (SuperCellEnum == SUPERCELL_SLAB)
37 {
38 Area = std::abs(ref.Lattice.R(0, 0) * ref.Lattice.R(1, 1) - ref.Lattice.R(0, 1) * ref.Lattice.R(1, 0));
39 PreFactors[0] = 2.0 * M_PI / Area; //\f$ \frac{2\pi}{A}\f$
40 PreFactors[1] = 2.0 * std::sqrt(M_PI) / (Sigma * Area); //\f$ \frac{2\pi}{A}\frac{1}{Sigma\pi}\f$
41 PreFactors[2] = 1.0 / (2 * Sigma); //Used for the k-dependent term
42 }
43 fillFk(ref.SK->KLists);
44 }
45
EwaldHandler(const EwaldHandler & aLR,ParticleSet & ref)46 EwaldHandler::EwaldHandler(const EwaldHandler& aLR, ParticleSet& ref)
47 : LRHandlerBase(aLR), Sigma(aLR.Sigma), Volume(aLR.Volume), Area(aLR.Area), PreFactors(aLR.PreFactors)
48 {
49 SuperCellEnum = aLR.SuperCellEnum;
50 if (SuperCellEnum == SUPERCELL_SLAB)
51 kMag = aLR.kMag;
52 }
53
fillFk(KContainer & KList)54 void EwaldHandler::fillFk(KContainer& KList)
55 {
56 Fk.resize(KList.kpts_cart.size());
57 const std::vector<int>& kshell(KList.kshell);
58 MaxKshell = kshell.size() - 1;
59 Fk_symm.resize(MaxKshell);
60 kMag.resize(MaxKshell);
61 if (SuperCellEnum == SUPERCELL_SLAB)
62 {
63 mRealType knorm = M_PI / Area;
64 mRealType ksum = 0.0;
65 mRealType oneovertwosigma = 1.0 / (2.0 * Sigma);
66 for (int ks = 0, ki = 0; ks < Fk_symm.size(); ks++)
67 {
68 kMag[ks] = std::sqrt(KList.ksq[ki]);
69 mRealType uk = knorm / kMag[ks]; //pi/(A*k)
70 Fk_symm[ks] = uk;
71 while (ki < KList.kshell[ks + 1] && ki < Fk.size())
72 Fk[ki++] = uk;
73 ksum += uk * erfc(kMag[ks] * oneovertwosigma) * (KList.kshell[ks + 1] - KList.kshell[ks]);
74 }
75 //correction for the self-energy returned by evaluateLR_r0
76 PreFactors[3] = 2.0 * (std::sqrt(M_PI) / (Area * Sigma) - ksum);
77 }
78 else
79 {
80 #if OHMMS_DIM == 2
81 mRealType kgauss = 1.0 / (4 * Sigma * Sigma);
82 mRealType knorm = 2 * M_PI / Volume;
83 for (int ks = 0, ki = 0; ks < Fk_symm.size(); ks++)
84 {
85 mRealType t2e = KList.ksq[ki] * kgauss;
86 mRealType uk = knorm * std::exp(-t2e) / KList.ksq[ki];
87 Fk_symm[ks] = uk;
88 while (ki < KList.kshell[ks + 1] && ki < Fk.size())
89 Fk[ki++] = uk;
90 }
91 PreFactors[3] = 0.0;
92 #elif OHMMS_DIM == 3
93 mRealType kgauss = 1.0 / (4 * Sigma * Sigma);
94 mRealType knorm = 4 * M_PI / Volume;
95 for (int ks = 0, ki = 0; ks < Fk_symm.size(); ks++)
96 {
97 mRealType t2e = KList.ksq[ki] * kgauss;
98 mRealType uk = knorm * std::exp(-t2e) / KList.ksq[ki];
99 Fk_symm[ks] = uk;
100 while (ki < KList.kshell[ks + 1] && ki < Fk.size())
101 Fk[ki++] = uk;
102 }
103 PreFactors[3] = 0.0;
104 #endif
105 }
106 app_log().flush();
107 }
108
evaluate_vlr_k(mRealType k)109 EwaldHandler::mRealType EwaldHandler::evaluate_vlr_k(mRealType k)
110 {
111 mRealType uk = 0.0;
112 if (SuperCellEnum == SUPERCELL_SLAB)
113 {
114 mRealType knorm = M_PI / Area;
115 uk = knorm / k; //pi/(A*k)
116 }
117 else
118 {
119 #if OHMMS_DIM == 2
120 mRealType kgauss = 1.0 / (4 * Sigma * Sigma);
121 mRealType knorm = 2 * M_PI / Volume;
122 mRealType k2 = k * k;
123 uk = knorm * std::exp(-k2 * kgauss) / k2;
124 #elif OHMMS_DIM == 3
125 mRealType kgauss = 1.0 / (4 * Sigma * Sigma);
126 mRealType knorm = 4 * M_PI / Volume;
127 mRealType k2 = k * k;
128 uk = knorm * std::exp(-k2 * kgauss) / k2;
129 }
130 #endif
131 return uk;
132 }
133
134 EwaldHandler::mRealType EwaldHandler::evaluate_slab(pRealType z, const std::vector<int>& kshell,
135 const pComplexType* restrict eikr_i,
136 const pComplexType* restrict eikr_j)
137 {
138 mRealType zp = z * Sigma;
139 mRealType vk = -SlabFunc0(z, zp);
140 //cout << "### SLAB " << z << " " << zp << std::endl;
141 for (int ks = 0, ki = 0; ks < MaxKshell; ks++)
142 {
143 mRealType u = 0; //\sum Real (e^ikr_i e^(-ikr_j))
144 for (; ki < kshell[ks + 1]; ki++, eikr_i++, eikr_j++)
145 u += ((*eikr_i).real() * (*eikr_j).real() + (*eikr_i).imag() * (*eikr_j).imag());
146 vk += u * Fk_symm[ks] * SlabFuncK(ks, z, zp);
147 }
148 return vk;
149 }
150 } // namespace qmcplusplus
151