1 //////////////////////////////////////////////////////////////////
2 // (c) Copyright 2008- by Jeongnim Kim
3 //////////////////////////////////////////////////////////////////
4 //////////////////////////////////////////////////////////////////
5 // National Center for Supercomputing Applications &
6 // Materials Computation Center
7 // University of Illinois, Urbana-Champaign
8 // Urbana, IL 61801
9 // e-mail: jnkim@ncsa.uiuc.edu
10 //
11 // Supported by
12 // National Center for Supercomputing Applications, UIUC
13 // Materials Computation Center, UIUC
14 //////////////////////////////////////////////////////////////////
15 // -*- C++ -*-
16 #include "SkAllEstimator.h"
17 #include "LongRange/StructFact.h"
18 #include "Utilities/IteratorUtility.h"
19 #include "OhmmsData/AttributeSet.h"
20 #include "Configuration.h"
21 #include <fstream>
22 #include <sstream>
23 #include <string>
24 namespace qmcplusplus
25 {
SkAllEstimator(ParticleSet & source,ParticleSet & target)26 SkAllEstimator::SkAllEstimator(ParticleSet& source, ParticleSet& target)
27 {
28 elns = ⌖
29 ions = &source;
30 NumeSpecies = elns->getSpeciesSet().getTotalNum();
31 NumIonSpecies = ions->getSpeciesSet().getTotalNum();
32 UpdateMode.set(COLLECTABLE, 1);
33
34 NumK = source.SK->KLists.numk;
35 OneOverN = 1.0 / static_cast<RealType>(source.getTotalNum());
36 Kshell = source.SK->KLists.kshell;
37 MaxKshell = Kshell.size() - 1;
38
39 #if defined(USE_REAL_STRUCT_FACTOR)
40 RhokTot_r.resize(NumK);
41 RhokTot_i.resize(NumK);
42
43 #else
44 RhokTot.resize(NumK);
45 #endif
46 //for values, we are including e-e structure factor, and e-Ion. So a total of NumIonSpecies+1 structure factors.
47 //+2 for the real and imaginary parts of rho_k^e
48 values.resize(3 * NumK);
49 Kmag.resize(MaxKshell);
50 OneOverDnk.resize(MaxKshell);
51 for (int ks = 0; ks < MaxKshell; ks++)
52 {
53 Kmag[ks] = std::sqrt(source.SK->KLists.ksq[Kshell[ks]]);
54 OneOverDnk[ks] = 1.0 / static_cast<RealType>(Kshell[ks + 1] - Kshell[ks]);
55 }
56 hdf5_out = false;
57 }
58
resetTargetParticleSet(ParticleSet & P)59 void SkAllEstimator::resetTargetParticleSet(ParticleSet& P) { elns = &P; }
60
61
evaluateIonIon()62 void SkAllEstimator::evaluateIonIon()
63 {
64 std::stringstream ss;
65 ss << " kx ky kz";
66
67 std::ofstream skfile;
68 std::stringstream filebuffer;
69 skfile.open("rhok_IonIon.dat");
70 for (int i = 0; i < NumIonSpecies; i++)
71 ss << " rho_" << i << "_r"
72 << " rho_" << i << "_i";
73
74 filebuffer << ss.str() << std::endl;
75
76 for (int k = 0; k < NumK; k++)
77 {
78 PosType kvec = ions->SK->KLists.kpts_cart[k];
79
80 filebuffer << kvec;
81 for (int i = 0; i < NumIonSpecies; i++)
82 {
83 double rho_i(0.0);
84 double rho_r(0.0);
85
86 #if defined(USE_REAL_STRUCT_FACTOR)
87 rho_r = ions->SK->rhok_r[i][k];
88 rho_i = ions->SK->rhok_i[i][k];
89 #else
90 rho_r = ions->SK->rhok[i][k].real();
91 rho_i = ions->SK->rhok[i][k].imag();
92 #endif
93 filebuffer << " " << rho_r << " " << rho_i;
94 }
95
96 filebuffer << std::endl;
97 }
98
99 skfile << filebuffer.str();
100
101 skfile.close();
102 }
103
104
evaluate(ParticleSet & P)105 SkAllEstimator::Return_t SkAllEstimator::evaluate(ParticleSet& P)
106 {
107 // Segfaults unless setHistories() is called prior
108 RealType w = tWalker->Weight;
109 #if defined(USE_REAL_STRUCT_FACTOR)
110 //sum over species
111 std::copy(P.SK->rhok_r[0], P.SK->rhok_r[0] + NumK, RhokTot_r.begin());
112 std::copy(P.SK->rhok_i[0], P.SK->rhok_i[0] + NumK, RhokTot_i.begin());
113 for (int i = 1; i < NumeSpecies; ++i)
114 accumulate_elements(P.SK->rhok_r[i], P.SK->rhok_r[i] + NumK, RhokTot_r.begin());
115 for (int i = 1; i < NumeSpecies; ++i)
116 accumulate_elements(P.SK->rhok_i[i], P.SK->rhok_i[i] + NumK, RhokTot_i.begin());
117
118 for (int k = 0; k < NumK; k++)
119 values[k] = w * (RhokTot_r[k] * RhokTot_r[k] + RhokTot_i[k] * RhokTot_i[k]);
120
121 // for (int ionSpec=0; ionSpec<NumIonSpecies; ionSpec++)
122 // {
123 // for(int k=0; k<NumK; k++)
124 // {
125 // RealType rhok_A_r(ions->SK->rhok_r[ionSpec][k]), rhok_A_i(ions->SK->rhok_i[ionSpec][k]);
126 // values[(ionSpec+1)*NumK+k]=RhokTot_r[k]*rhok_A_r+RhokTot_i[k]*rhok_A_i;
127 // }
128 // }
129
130 if (hdf5_out)
131 {
132 int kc = myIndex;
133 for (int k = 0; k < NumK; k++, kc++)
134 P.Collectables[kc] += values[k];
135 for (int k = 0; k < NumK; k++, kc++)
136 P.Collectables[kc] += w * RhokTot_r[k];
137 for (int k = 0; k < NumK; k++, kc++)
138 P.Collectables[kc] += w * RhokTot_i[k];
139 }
140 else
141 {
142 for (int k = 0; k < NumK; k++)
143 {
144 values[NumK + 2 * k] = w * RhokTot_r[k];
145 values[NumK + 2 * k + 1] = w * RhokTot_i[k];
146 }
147 }
148 #else // Is this path ever touched?
149 //sum over species
150 std::copy(P.SK->rhok[0], P.SK->rhok[0] + NumK, RhokTot.begin());
151 for (int i = 1; i < NumeSpecies; ++i)
152 accumulate_elements(P.SK->rhok[i], P.SK->rhok[i] + NumK, RhokTot.begin());
153 for (int k = 0; k < NumK; k++)
154 values[k] = w * (rhok[k].real() * rhok[k].real() + rhok[k].imag() * rhok[k].imag());
155
156 // for (int ionSpec=0; ionSpec<NumIonSpecies; ionSpec++)
157 // {
158 // for(int k=0; k<NumK; k++)
159 // {
160 // RealType rhok_A_r(ions->SK->rhok[ionSpec][k].real()), rhok_A_i(ions->SK->rhok[ionSpec][k].imag());
161 // values[(ionSpec+1)*NumK+k]=rhok[k].real()*rhok_A_r+rho_k[k].imag()*rhok_A_i;
162 // }
163 // }
164 if (hdf5_out)
165 {
166 int kc = myIndex;
167 for (int k = 0; k < NumK; k++, kc++)
168 P.Collectables[kc] += values[k];
169 for (int k = 0; k < NumK; k++, kc++)
170 P.Collectables[kc] += w * rhok[k].real();
171 for (int k = 0; k < NumK; k++, kc++)
172 P.Collectables[kc] += w * rhok[k].imag();
173 }
174 else
175 {
176 for (int k = 0, count = 0; k < NumK; k++)
177 {
178 value[NumK + count] = w * rhok[k].real();
179 count++;
180 value[NumK + count] = w * rhok[k].imag();
181 count++;
182 }
183 }
184 #endif
185 return 0.0;
186 }
187
addObservables(PropertySetType & plist,BufferType & collectables)188 void SkAllEstimator::addObservables(PropertySetType& plist, BufferType& collectables)
189 {
190 if (hdf5_out)
191 {
192 myIndex = collectables.size();
193 std::vector<RealType> tmp(3 * NumK); // space for e-e modulus, e real, e imag
194 collectables.add(tmp.begin(), tmp.end());
195 }
196 else
197 {
198 myIndex = plist.size();
199 //First the electron structure factor
200 for (int i = 0; i < NumK; i++)
201 {
202 std::stringstream sstr;
203 sstr << "rhok_e_e_" << i;
204 int id = plist.add(sstr.str());
205 }
206 //Now the e-Ion structure factors. IonIon are dumped to file, and not evaluated.
207 // for (int ionSpec=0; ionSpec<NumIonSpecies; ionSpec++)
208 // {
209 // for (int i=0; i<NumK; i++)
210 // {
211 // std::stringstream sstr;
212 // sstr << "rhok_e_" <<ionSpec<<"_"<<i;
213 // int id=plist.add(sstr.str());
214 // }
215 // }
216
217 for (int k = 0; k < NumK; k++)
218 {
219 std::stringstream sstr1, sstr2;
220 sstr1 << "rhok_e_r_" << k;
221 sstr2 << "rhok_e_i_" << k;
222 int id = plist.add(sstr1.str());
223 int id2 = plist.add(sstr2.str());
224 }
225 }
226 }
227
addObservables(PropertySetType & plist)228 void SkAllEstimator::addObservables(PropertySetType& plist)
229 {
230 myIndex = plist.size();
231 //First the electron structure factor
232 for (int i = 0; i < NumK; i++)
233 {
234 std::stringstream sstr;
235 sstr << "rhok_e_e_" << i;
236 int id = plist.add(sstr.str());
237 }
238 //Now the e-Ion structure factors. IonIon are dumped to file, and not evaluated.
239 // for (int ionSpec=0; ionSpec<NumIonSpecies; ionSpec++)
240 // {
241 // for (int i=0; i<NumK; i++)
242 // {
243 // std::stringstream sstr;
244 // sstr << "rhok_e_" <<ionSpec<<"_"<<i;
245 // int id=plist.add(sstr.str());
246 // }
247 // }
248 for (int k = 0; k < NumK; k++)
249 {
250 std::stringstream sstr1, sstr2;
251 sstr1 << "rhok_e_r_" << k;
252 sstr2 << "rhok_e_i_" << k;
253 int id = plist.add(sstr1.str());
254 int id2 = plist.add(sstr2.str());
255 }
256 }
257
setObservables(PropertySetType & plist)258 void SkAllEstimator::setObservables(PropertySetType& plist)
259 {
260 if (!hdf5_out)
261 std::copy(values.begin(), values.end(), plist.begin() + myIndex);
262 }
263
setParticlePropertyList(PropertySetType & plist,int offset)264 void SkAllEstimator::setParticlePropertyList(PropertySetType& plist, int offset)
265 {
266 if (!hdf5_out)
267 std::copy(values.begin(), values.end(), plist.begin() + myIndex + offset);
268 }
269
270
registerCollectables(std::vector<observable_helper * > & h5desc,hid_t gid) const271 void SkAllEstimator::registerCollectables(std::vector<observable_helper*>& h5desc, hid_t gid) const
272 {
273 if (hdf5_out)
274 {
275 // Create HDF group in stat.h5 with SkAllEstimator's name
276 hid_t sgid = H5Gcreate(gid, myName.c_str(), 0);
277
278 // Add k-point information
279 observable_helper* oh = new observable_helper("kpoints");
280 oh->open(sgid); // add to SkAll hdf group
281 oh->addProperty(const_cast<std::vector<PosType>&>(ions->SK->KLists.kpts_cart), "value");
282 h5desc.push_back(oh);
283
284 // Add electron-electron S(k)
285 std::vector<int> ng(1);
286 ng[0] = NumK;
287 // modulus
288 oh = new observable_helper("rhok_e_e");
289 oh->set_dimensions(ng, myIndex);
290 oh->open(sgid); // add to SkAll hdf group
291 h5desc.push_back(oh);
292 // real part
293 oh = new observable_helper("rhok_e_r");
294 oh->set_dimensions(ng, myIndex + NumK);
295 oh->open(sgid); // add to SkAll hdf group
296 h5desc.push_back(oh);
297 // imaginary part
298 oh = new observable_helper("rhok_e_i");
299 oh->set_dimensions(ng, myIndex + 2 * NumK);
300 oh->open(sgid); // add to SkAll hdf group
301 h5desc.push_back(oh);
302 }
303 }
304
put(xmlNodePtr cur)305 bool SkAllEstimator::put(xmlNodePtr cur)
306 {
307 OhmmsAttributeSet pAttrib;
308 std::string hdf5_flag = "no";
309 std::string write_ionion_flag = "no";
310
311 pAttrib.add(hdf5_flag, "hdf5");
312 pAttrib.add(write_ionion_flag, "writeionion");
313 pAttrib.put(cur);
314 if (hdf5_flag == "yes")
315 hdf5_out = true;
316 else
317 hdf5_out = false;
318 app_log() << "hdf5_out= " << hdf5_out << "\n";
319
320 if (write_ionion_flag == "Yes" || write_ionion_flag == "yes")
321 {
322 app_log() << "SkAll: evaluateIonIon()\n";
323 evaluateIonIon();
324 }
325 return true;
326 }
327
328
329 // Display some stuff
get(std::ostream & os) const330 bool SkAllEstimator::get(std::ostream& os) const
331 {
332 app_log() << std::fixed;
333 app_log() << "\n";
334 app_log() << " SkAllEstimator Settings:\n";
335 app_log() << " ========================\n";
336 app_log() << " NumeSpecies " << std::setw(10) << std::setprecision(4) << NumeSpecies << "\n";
337 app_log() << " NumIonSpecies " << std::setw(10) << std::setprecision(4) << NumIonSpecies << "\n";
338 app_log() << " NumK " << std::setw(10) << std::setprecision(4) << NumK << "\n";
339 app_log() << " MaxKshell " << std::setw(10) << std::setprecision(4) << MaxKshell << "\n";
340
341 return true;
342 }
343
makeClone(ParticleSet & qp,TrialWaveFunction & psi)344 OperatorBase* SkAllEstimator::makeClone(ParticleSet& qp, TrialWaveFunction& psi)
345 {
346 SkAllEstimator* myclone = new SkAllEstimator(*this);
347 myclone->hdf5_out = hdf5_out;
348 myclone->myIndex = myIndex;
349 return myclone;
350 }
351 } // namespace qmcplusplus
352