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) 2016 Jeongnim Kim and QMCPACK developers.
6 //
7 // File developed by: Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
8 // Jeongnim Kim, jeongnim.kim@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 "Message/Communicate.h"
16 #include "KContainer.h"
17 #include "Utilities/qmc_common.h"
18 #include <map>
19
20 namespace qmcplusplus
21 {
UpdateKLists(ParticleLayout_t & lattice,RealType kc,bool useSphere)22 void KContainer::UpdateKLists(ParticleLayout_t& lattice, RealType kc, bool useSphere)
23 {
24 kcutoff = kc;
25 kcut2 = kc * kc;
26 if (kcutoff <= 0.0)
27 {
28 APP_ABORT(" Illegal cutoff for KContainer");
29 }
30 FindApproxMMax(lattice);
31 BuildKLists(lattice, useSphere);
32
33 app_log() << " KContainer initialised with cutoff " << kcutoff << std::endl;
34 app_log() << " # of K-shell = " << kshell.size() << std::endl;
35 app_log() << " # of K points = " << kpts.size() << std::endl;
36 app_log() << std::endl;
37 }
38
FindApproxMMax(ParticleLayout_t & lattice)39 void KContainer::FindApproxMMax(ParticleLayout_t& lattice)
40 {
41 //Estimate the size of the parallelpiped that encompasses a sphere of kcutoff.
42 //mmax is stored as integer translations of the reciprocal cell vectors.
43 //Does not require an orthorhombic cell.
44 /* Old method.
45 //2pi is not included in lattice.b
46 Matrix<RealType> mmat;
47 mmat.resize(3,3);
48 for(int j=0;j<3;j++)
49 for(int i=0;i<3;i++){
50 mmat[i][j] = 0.0;
51 for(int k=0;k<3;k++)
52 mmat[i][j] = mmat[i][j] + 4.0*M_PI*M_PI*lattice.b(k)[i]*lattice.b(j)[k];
53 }
54
55 TinyVector<RealType,3> x,temp;
56 RealType tempr;
57 for(int idim=0;idim<3;idim++){
58 int i = ((idim)%3);
59 int j = ((idim+1)%3);
60 int k = ((idim+2)%3);
61
62 x[i] = 1.0;
63 x[j] = (mmat[j][k]*mmat[k][i] - mmat[k][k]*mmat[i][j]);
64 x[j]/= (mmat[j][j]*mmat[k][k] - mmat[j][k]*mmat[j][k]);
65 x[k] = -(mmat[k][i] + mmat[j][k]*x[j])/mmat[k][k];
66
67 for(i=0;i<3;i++){
68 temp[i] = 0.0;
69 for(j=0;j<3;j++)
70 temp[i] += mmat[i][j]*x[j];
71 }
72
73 tempr = dot(x,temp);
74 mmax[idim] = static_cast<int>(sqrt(4.0*kcut2/tempr)) + 1;
75 }
76 */
77 // see rmm, Electronic Structure, p. 85 for details
78 for (int i = 0; i < DIM; i++)
79 mmax[i] = static_cast<int>(std::floor(std::sqrt(dot(lattice.a(i), lattice.a(i))) * kcutoff / (2 * M_PI))) + 1;
80
81
82 //overwrite the non-periodic directon to be zero
83 if (qmc_common.use_ewald)
84 {
85 app_log() << " Using Ewald sum for the slab " << std::endl;
86 #if OHMMS_DIM == 3
87 if (lattice.SuperCellEnum == SUPERCELL_SLAB)
88 mmax[2] = 0;
89 // if(lattice.SuperCellEnum == SUPERCELL_WIRE) mmax[1]=mmax[2]=0;
90 //#elif OHMMS_DIM==2
91 // if(lattice.SuperCellEnum == SUPERCELL_WIRE)
92 // mmax[1]=0;
93 #endif
94 }
95
96 mmax[DIM] = mmax[0];
97 for (int i = 1; i < DIM; ++i)
98 mmax[DIM] = std::max(mmax[i], mmax[DIM]);
99 }
100
BuildKLists(ParticleLayout_t & lattice,bool useSphere)101 void KContainer::BuildKLists(ParticleLayout_t& lattice, bool useSphere)
102 {
103 TinyVector<int, DIM + 1> TempActualMax;
104 TinyVector<int, DIM> kvec;
105 TinyVector<RealType, DIM> kvec_cart;
106 RealType modk2;
107 std::vector<TinyVector<int, DIM>> kpts_tmp;
108 VContainer_t kpts_cart_tmp;
109 SContainer_t ksq_tmp;
110 // reserve the space for memory efficiency
111 #if OHMMS_DIM == 3
112 if (useSphere)
113 {
114 //Loop over guesses for valid k-points.
115 for (int i = -mmax[0]; i <= mmax[0]; i++)
116 {
117 kvec[0] = i;
118 for (int j = -mmax[1]; j <= mmax[1]; j++)
119 {
120 kvec[1] = j;
121 for (int k = -mmax[2]; k <= mmax[2]; k++)
122 {
123 kvec[2] = k;
124 //Do not include k=0 in evaluations.
125 if (i == 0 && j == 0 && k == 0)
126 continue;
127 //Convert kvec to Cartesian
128 kvec_cart = lattice.k_cart(kvec);
129 //Find modk
130 modk2 = dot(kvec_cart, kvec_cart);
131 if (modk2 > kcut2)
132 continue; //Inside cutoff?
133 //This k-point should be added to the list
134 kpts_tmp.push_back(kvec);
135 kpts_cart_tmp.push_back(kvec_cart);
136 ksq_tmp.push_back(modk2);
137 //Update record of the allowed maximum translation.
138 for (int idim = 0; idim < 3; idim++)
139 if (std::abs(kvec[idim]) > TempActualMax[idim])
140 TempActualMax[idim] = std::abs(kvec[idim]);
141 }
142 }
143 }
144 }
145 else
146 {
147 // Loop over all k-points in the parallelpiped and add them to kcontainer
148 // note layout is for interfacing with fft, so for each dimension, the
149 // positive indexes come first then the negative indexes backwards
150 // e.g. 0, 1, .... mmax, -mmax+1, -mmax+2, ... -1
151 const int idimsize = mmax[0] * 2;
152 const int jdimsize = mmax[1] * 2;
153 const int kdimsize = mmax[2] * 2;
154 for (int i = 0; i < idimsize; i++)
155 {
156 kvec[0] = i;
157 if (kvec[0] > mmax[0])
158 kvec[0] -= idimsize;
159 for (int j = 0; j < jdimsize; j++)
160 {
161 kvec[1] = j;
162 if (kvec[1] > mmax[1])
163 kvec[1] -= jdimsize;
164 for (int k = 0; k < kdimsize; k++)
165 {
166 kvec[2] = k;
167 if (kvec[2] > mmax[2])
168 kvec[2] -= kdimsize;
169 // get cartesian location and modk2
170 kvec_cart = lattice.k_cart(kvec);
171 modk2 = dot(kvec_cart, kvec_cart);
172 // add k-point to lists
173 kpts_tmp.push_back(kvec);
174 kpts_cart_tmp.push_back(kvec_cart);
175 ksq_tmp.push_back(modk2);
176 }
177 }
178 }
179 // set allowed maximum translation
180 TempActualMax[0] = mmax[0];
181 TempActualMax[1] = mmax[1];
182 TempActualMax[2] = mmax[2];
183 }
184 #elif OHMMS_DIM == 2
185 if (useSphere)
186 {
187 //Loop over guesses for valid k-points.
188 for (int i = -mmax[0]; i <= mmax[0]; i++)
189 {
190 kvec[0] = i;
191 for (int j = -mmax[1]; j <= mmax[1]; j++)
192 {
193 kvec[1] = j;
194 //Do not include k=0 in evaluations.
195 if (i == 0 && j == 0)
196 continue;
197 //Convert kvec to Cartesian
198 kvec_cart = lattice.k_cart(kvec);
199 //Find modk
200 modk2 = dot(kvec_cart, kvec_cart);
201 if (modk2 > kcut2)
202 continue; //Inside cutoff?
203 //This k-point should be added to the list
204 kpts_tmp.push_back(kvec);
205 kpts_cart_tmp.push_back(kvec_cart);
206 ksq_tmp.push_back(modk2);
207 //Update record of the allowed maximum translation.
208 for (int idim = 0; idim < 3; idim++)
209 if (std::abs(kvec[idim]) > TempActualMax[idim])
210 TempActualMax[idim] = std::abs(kvec[idim]);
211 }
212 }
213 }
214 else
215 {
216 // Loop over all k-points in the parallelpiped and add them to kcontainer
217 // note layout is for interfacing with fft, so for each dimension, the
218 // positive indexes come first then the negative indexes backwards
219 // e.g. 0, 1, .... mmax, -mmax+1, -mmax+2, ... -1
220 const int idimsize = mmax[0] * 2;
221 const int jdimsize = mmax[1] * 2;
222 for (int i = 0; i < idimsize; i++)
223 {
224 kvec[0] = i;
225 if (kvec[0] > mmax[0])
226 kvec[0] -= idimsize;
227 for (int j = 0; j < jdimsize; j++)
228 {
229 kvec[1] = j;
230 if (kvec[1] > mmax[1])
231 kvec[1] -= jdimsize;
232 // get cartesian location and modk2
233 kvec_cart = lattice.k_cart(kvec);
234 modk2 = dot(kvec_cart, kvec_cart);
235 // add k-point to lists
236 kpts_tmp.push_back(kvec);
237 kpts_cart_tmp.push_back(kvec_cart);
238 ksq_tmp.push_back(modk2);
239 }
240 }
241 // set allowed maximum translation
242 TempActualMax[0] = mmax[0];
243 TempActualMax[1] = mmax[1];
244 }
245 //#elif OHMMS_DIM == 1
246 //add one-dimension
247 #else
248 #error "OHMMS_DIM != 2 || OHMMS_DIM != 3"
249 #endif
250 //Update a record of the number of k vectors
251 numk = kpts_tmp.size();
252 std::map<long long, std::vector<int>*> kpts_sorted;
253 //create the map: use simple integer with resolution of 0.00000001 in ksq
254 for (int ik = 0; ik < numk; ik++)
255 {
256 #ifdef MIXED_PRECISION
257 long long k_ind = static_cast<long long>(ksq_tmp[ik] * 1000);
258 #else
259 //This is a workaround for ewald bug (Issue #2105) for FULL PRECISION ONLY. Basically, 1e-7 is the resolution of |k|^2 for doubles,
260 //so we jack up the tolerance to match that.
261 long long k_ind = static_cast<long long>(ksq_tmp[ik] * 10000000);
262 #endif
263 std::map<long long, std::vector<int>*>::iterator it(kpts_sorted.find(k_ind));
264 if (it == kpts_sorted.end())
265 {
266 std::vector<int>* newSet = new std::vector<int>;
267 kpts_sorted[k_ind] = newSet;
268 newSet->push_back(ik);
269 }
270 else
271 {
272 (*it).second->push_back(ik);
273 }
274 }
275 std::map<long long, std::vector<int>*>::iterator it(kpts_sorted.begin());
276 kpts.resize(numk);
277 kpts_cart.resize(numk);
278 ksq.resize(numk);
279 kshell.resize(kpts_sorted.size() + 1, 0);
280 int ok = 0, ish = 0;
281 while (it != kpts_sorted.end())
282 {
283 std::vector<int>::iterator vit((*it).second->begin());
284 while (vit != (*it).second->end())
285 {
286 int ik = (*vit);
287 kpts[ok] = kpts_tmp[ik];
288 kpts_cart[ok] = kpts_cart_tmp[ik];
289 ksq[ok] = ksq_tmp[ik];
290 ++vit;
291 ++ok;
292 }
293 kshell[ish + 1] = kshell[ish] + (*it).second->size();
294 ++it;
295 ++ish;
296 }
297 it = kpts_sorted.begin();
298 std::map<long long, std::vector<int>*>::iterator e_it(kpts_sorted.end());
299 while (it != e_it)
300 {
301 delete it->second;
302 it++;
303 }
304 //Finished searching k-points. Copy list of maximum translations.
305 mmax[DIM] = 0;
306 for (int idim = 0; idim < DIM; idim++)
307 {
308 mmax[idim] = TempActualMax[idim];
309 mmax[DIM] = std::max(mmax[idim], mmax[DIM]);
310 //if(mmax[idim] > mmax[DIM]) mmax[DIM] = mmax[idim];
311 }
312 //Now fill the array that returns the index of -k when given the index of k.
313 minusk.resize(numk);
314
315 //Assigns a unique hash value to each kpoint.
316 auto getHashOfVec = [](const auto& inpv, int hashparam) -> long long {
317 long long hash = 0; // this will cause integral promotion below
318 for (int i = 0; i < inpv.Size; ++i)
319 hash += inpv[i] + hash * hashparam;
320 return hash;
321 };
322
323 // Create a map from the hash value for each k vector to the index
324 std::map<long long, int> hashToIndex;
325 for (int ki = 0; ki < numk; ki++)
326 {
327 hashToIndex[getHashOfVec(kpts[ki], numk)] = ki;
328 }
329 // Use the map to find the index of -k from the index of k
330 for (int ki = 0; ki < numk; ki++)
331 {
332 minusk[ki] = hashToIndex[getHashOfVec(-1 * kpts[ki], numk)];
333 }
334 }
335
336 } // namespace qmcplusplus
337