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: Jordan E. Vincent, University of Illinois at Urbana-Champaign
8 // Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
9 // Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
10 // Miguel Morales, moralessilva2@llnl.gov, Lawrence Livermore National Laboratory
11 // Anouar Benali, benali@anl.gov, Argonne National Laboratory
12 // Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
13 //
14 // File created by: Jordan E. Vincent, University of Illinois at Urbana-Champaign
15 //////////////////////////////////////////////////////////////////////////////////////
16
17
18 #include "QMCGaussianParserBase.h"
19 #include "ParticleIO/XMLParticleIO.h"
20 #include "Numerics/HDFSTLAttrib.h"
21 #include "OhmmsData/HDFStringAttrib.h"
22 #include <iterator>
23 #include <algorithm>
24 #include <numeric>
25 #include "hdf/hdf_archive.h"
26 #include <set>
27 #include <map>
28 #include <sstream>
29 #include <bitset>
30 #include <iomanip>
31
32
33 //std::vector<std::string> QMCGaussianParserBase::IonName;
34 const int OhmmsAsciiParser::bufferSize;
35 std::map<int, std::string> QMCGaussianParserBase::IonName;
36 std::vector<std::string> QMCGaussianParserBase::gShellType;
37 std::vector<int> QMCGaussianParserBase::gShellID;
38
QMCGaussianParserBase()39 QMCGaussianParserBase::QMCGaussianParserBase()
40 : multideterminant(false),
41 multidetH5(false),
42 BohrUnit(true),
43 SpinRestricted(false),
44 Periodicity(false),
45 UseHDF5(false),
46 PBC(false),
47 production(false),
48 zeroCI(false),
49 orderByExcitation(false),
50 addJastrow(true),
51 addJastrow3Body(false),
52 ECP(false),
53 debug(false),
54 Structure(false),
55 DoCusp(false),
56 FixValence(false),
57 singledetH5(false),
58 optDetCoeffs(false),
59 usingCSF(false),
60 NumberOfAtoms(0),
61 NumberOfEls(0),
62 target_state(0),
63 SpinMultiplicity(0),
64 NumberOfAlpha(0),
65 NumberOfBeta(0),
66 SizeOfBasisSet(0),
67 numMO(0),
68 readNO(0),
69 readGuess(0),
70 numMO2print(-1),
71 ci_size(0),
72 ci_nca(0),
73 ci_ncb(0),
74 ci_nea(0),
75 ci_neb(0),
76 ci_nstates(0),
77 NbKpts(0),
78 nbexcitedstates(0),
79 ci_threshold(1e-20),
80 Title("sample"),
81 basisType("Gaussian"),
82 basisName("generic"),
83 Normalized("no"),
84 CurrentCenter(""),
85 outputFile(""),
86 angular_type("spherical"),
87 h5file(""),
88 multih5file(""),
89 WFS_name("wfj"),
90 CodeName(""),
91 gShell(0),
92 gNumber(0),
93 gBound(0),
94 Occ_alpha(0),
95 Occ_beta(0),
96 gridPtr(0),
97 X(0),
98 Y(0),
99 Z(0)
100 {}
101
QMCGaussianParserBase(int argc,char ** argv)102 QMCGaussianParserBase::QMCGaussianParserBase(int argc, char** argv)
103 : multideterminant(false),
104 multidetH5(false),
105 BohrUnit(true),
106 SpinRestricted(false),
107 Periodicity(false),
108 UseHDF5(false),
109 PBC(false),
110 production(false),
111 zeroCI(false),
112 orderByExcitation(false),
113 addJastrow(true),
114 addJastrow3Body(false),
115 ECP(false),
116 debug(false),
117 Structure(false),
118 DoCusp(false),
119 FixValence(false),
120 singledetH5(false),
121 optDetCoeffs(false),
122 usingCSF(false),
123 NumberOfAtoms(0),
124 NumberOfEls(0),
125 target_state(0),
126 SpinMultiplicity(0),
127 NumberOfAlpha(0),
128 NumberOfBeta(0),
129 SizeOfBasisSet(0),
130 numMO(0),
131 readNO(0),
132 readGuess(0),
133 numMO2print(-1),
134 ci_size(0),
135 ci_nca(0),
136 ci_ncb(0),
137 ci_nea(0),
138 ci_neb(0),
139 ci_nstates(0),
140 NbKpts(0),
141 nbexcitedstates(0),
142 ci_threshold(1e-20),
143 Title("sample"),
144 basisType("Gaussian"),
145 basisName("generic"),
146 Normalized("no"),
147 CurrentCenter(""),
148 outputFile(""),
149 angular_type("spherical"),
150 h5file(""),
151 multih5file(""),
152 WFS_name("wfj"),
153 CodeName(""),
154 gShell(0),
155 gNumber(0),
156 gBound(0),
157 Occ_alpha(0),
158 Occ_beta(0),
159 gridPtr(0),
160 X(0),
161 Y(0),
162 Z(0)
163 {
164 IonChargeIndex = IonSystem.getSpeciesSet().addAttribute("charge");
165 ValenceChargeIndex = IonSystem.getSpeciesSet().addAttribute("valence");
166 AtomicNumberIndex = IonSystem.getSpeciesSet().addAttribute("atomicnumber");
167 std::cout << "Index of ion charge " << IonChargeIndex << std::endl;
168 std::cout << "Index of valence charge " << ValenceChargeIndex << std::endl;
169 Image.resize(3);
170 Image[0] = 8;
171 Image[1] = 8;
172 Image[2] = 8;
173 createGridNode(argc, argv);
174 }
175
init()176 void QMCGaussianParserBase::init()
177 {
178 IonName[1] = "H";
179 IonName[2] = "He";
180 IonName[3] = "Li";
181 IonName[4] = "Be";
182 IonName[5] = "B";
183 IonName[6] = "C";
184 IonName[7] = "N";
185 IonName[8] = "O";
186 IonName[9] = "F";
187 IonName[10] = "Ne";
188 IonName[11] = "Na";
189 IonName[12] = "Mg";
190 IonName[13] = "Al";
191 IonName[14] = "Si";
192 IonName[15] = "P";
193 IonName[16] = "S";
194 IonName[17] = "Cl";
195 IonName[18] = "Ar";
196 IonName[19] = "K";
197 IonName[20] = "Ca";
198 IonName[21] = "Sc";
199 IonName[22] = "Ti";
200 IonName[23] = "V";
201 IonName[24] = "Cr";
202 IonName[25] = "Mn";
203 IonName[26] = "Fe";
204 IonName[27] = "Co";
205 IonName[28] = "Ni";
206 IonName[29] = "Cu";
207 IonName[30] = "Zn";
208 IonName[31] = "Ga";
209 IonName[32] = "Ge";
210 IonName[33] = "As";
211 IonName[34] = "Se";
212 IonName[35] = "Br";
213 IonName[36] = "Kr";
214 IonName[37] = "Rb";
215 IonName[38] = "Sr";
216 IonName[39] = "Y";
217 IonName[40] = "Zr";
218 IonName[41] = "Nb";
219 IonName[42] = "Mo";
220 IonName[43] = "Tc";
221 IonName[44] = "Ru";
222 IonName[45] = "Rh";
223 IonName[46] = "Pd";
224 IonName[47] = "Ag";
225 IonName[48] = "Cd";
226 IonName[49] = "In";
227 IonName[50] = "Sn";
228 IonName[51] = "Sb";
229 IonName[52] = "Te";
230 IonName[53] = "I";
231 IonName[54] = "Xe";
232 IonName[55] = "Cs";
233 IonName[56] = "Ba";
234 IonName[57] = "La";
235 IonName[58] = "Ce";
236 IonName[59] = "Pr";
237 IonName[60] = "Nd";
238 IonName[61] = "Pm";
239 IonName[62] = "Sm";
240 IonName[63] = "Eu";
241 IonName[64] = "Gd";
242 IonName[65] = "Tb";
243 IonName[66] = "Dy";
244 IonName[67] = "Ho";
245 IonName[68] = "Er";
246 IonName[69] = "Tm";
247 IonName[70] = "Yb";
248 IonName[71] = "Lu";
249 IonName[72] = "Hf";
250 IonName[73] = "Ta";
251 IonName[74] = "W";
252 IonName[75] = "Re";
253 IonName[76] = "Os";
254 IonName[77] = "Ir";
255 IonName[78] = "Pt";
256 IonName[79] = "Au";
257 IonName[80] = "Hg";
258 IonName[81] = "Tl";
259 IonName[82] = "Pb";
260 IonName[83] = "Bi";
261 IonName[84] = "Po";
262 IonName[85] = "At";
263 IonName[86] = "Rn";
264 IonName[87] = "Fr";
265 IonName[88] = "Ra";
266 IonName[89] = "Ac";
267 IonName[90] = "Th";
268 IonName[91] = "Pa";
269 IonName[92] = "U";
270 IonName[93] = "Np";
271 gShellType.resize(10);
272 gShellType[1] = "s";
273 gShellType[2] = "sp";
274 gShellType[3] = "p";
275 gShellType[4] = "d";
276 gShellType[5] = "f";
277 gShellType[6] = "g";
278 gShellType[7] = "h";
279 gShellType[8] = "h1";
280 gShellType[9] = "h2";
281 gShellID.resize(10);
282 gShellID[1] = 0;
283 gShellID[2] = 0;
284 gShellID[3] = 1; //gShellID[4]=2; gShellID[5]=3; gShellID[6]=4; gShellID[7]=5;
285 for (int i = 4, l = 2; i < gShellID.size(); ++i, ++l)
286 gShellID[i] = l;
287 }
288
setOccupationNumbers()289 void QMCGaussianParserBase::setOccupationNumbers()
290 {
291 int ds = SpinMultiplicity - 1;
292 NumberOfBeta = (NumberOfEls - ds) / 2;
293 NumberOfAlpha = NumberOfEls - NumberOfBeta;
294 if (!SpinRestricted)
295 //UHF
296 {
297 std::multimap<value_type, int> e;
298 for (int i = 0; i < numMO; i++)
299 e.insert(std::pair<value_type, int>(EigVal_alpha[i], 0));
300 for (int i = 0; i < numMO; i++)
301 e.insert(std::pair<value_type, int>(EigVal_beta[i], 1));
302 int n = 0;
303 std::multimap<value_type, int>::iterator it(e.begin());
304 LOGMSG("Unrestricted HF. Sorted eigen values")
305 while (n < NumberOfEls && it != e.end())
306 {
307 LOGMSG(n << " " << (*it).first << " " << (*it).second)
308 ++it;
309 ++n;
310 }
311 }
312 //}
313 LOGMSG("Number of alpha electrons " << NumberOfAlpha)
314 LOGMSG("Number of beta electrons " << NumberOfBeta)
315 Occ_alpha.resize(numMO, 0);
316 Occ_beta.resize(numMO, 0);
317 for (int i = 0; i < NumberOfAlpha; i++)
318 Occ_alpha[i] = 1;
319 for (int i = 0; i < NumberOfBeta; i++)
320 Occ_beta[i] = 1;
321 }
322
createElectronSet(const std::string & ion_tag)323 xmlNodePtr QMCGaussianParserBase::createElectronSet(const std::string& ion_tag)
324 {
325 ParticleSet els;
326 els.setName("e");
327 std::vector<int> nel(2);
328 nel[0] = NumberOfAlpha;
329 nel[1] = NumberOfBeta;
330 els.create(nel);
331 int iu = els.getSpeciesSet().addSpecies("u");
332 int id = els.getSpeciesSet().addSpecies("d");
333 int ic = els.getSpeciesSet().addAttribute("charge");
334 els.getSpeciesSet()(ic, iu) = -1;
335 els.getSpeciesSet()(ic, id) = -1;
336
337 xmlNodePtr cur = xmlNewNode(NULL, (const xmlChar*)"particleset");
338 xmlNewProp(cur, (const xmlChar*)"name", (const xmlChar*)els.getName().c_str());
339 xmlNewProp(cur, (const xmlChar*)"random", (const xmlChar*)"yes");
340 xmlNewProp(cur, (const xmlChar*)"randomsrc", (const xmlChar*)ion_tag.c_str());
341
342 //Electron up
343 {
344 std::ostringstream ng;
345 ng << els.last(0) - els.first(0);
346 xmlNodePtr g = xmlNewNode(NULL, (const xmlChar*)"group");
347 xmlNewProp(g, (const xmlChar*)"name", (const xmlChar*)"u");
348 xmlNewProp(g, (const xmlChar*)"size", (const xmlChar*)ng.str().c_str());
349 xmlNodePtr p = xmlNewTextChild(g, NULL, (const xmlChar*)"parameter", (const xmlChar*)"-1");
350 xmlNewProp(p, (const xmlChar*)"name", (const xmlChar*)"charge");
351 xmlAddChild(cur, g);
352 }
353 //Electron dn
354 {
355 std::ostringstream ng;
356 ng << els.last(1) - els.first(1);
357 xmlNodePtr g = xmlNewNode(NULL, (const xmlChar*)"group");
358 xmlNewProp(g, (const xmlChar*)"name", (const xmlChar*)"d");
359 xmlNewProp(g, (const xmlChar*)"size", (const xmlChar*)ng.str().c_str());
360 xmlNodePtr p = xmlNewTextChild(g, NULL, (const xmlChar*)"parameter", (const xmlChar*)"-1");
361 xmlNewProp(p, (const xmlChar*)"name", (const xmlChar*)"charge");
362 xmlAddChild(cur, g);
363 }
364 return cur;
365 }
366
367
createCell()368 xmlNodePtr QMCGaussianParserBase::createCell()
369 {
370 xmlNodePtr cur = xmlNewNode(NULL, (const xmlChar*)"simulationcell");
371 std::ostringstream vec;
372 vec.setf(std::ios::scientific, std::ios::floatfield);
373 vec.setf(std::ios::right, std::ios::adjustfield);
374 vec.precision(14);
375 vec << "\n";
376 vec << std::setw(22) << X[0] << std::setw(22) << X[1] << std::setw(22) << X[2] << std::setw(22) << "\n";
377 vec << std::setw(22) << Y[0] << std::setw(22) << Y[1] << std::setw(22) << Y[2] << std::setw(22) << "\n";
378 vec << std::setw(22) << Z[0] << std::setw(22) << Z[1] << std::setw(22) << Z[2] << std::setw(22) << "\n";
379
380 xmlNodePtr LatVec = xmlNewTextChild(cur, NULL, (const xmlChar*)"parameter", (const xmlChar*)vec.str().c_str());
381 xmlNewProp(LatVec, (const xmlChar*)"name", (const xmlChar*)"lattice");
382 xmlAddChild(cur, LatVec);
383
384 xmlNodePtr bconds = xmlNewTextChild(cur, NULL, (const xmlChar*)"parameter", (const xmlChar*)"p p p");
385 xmlNewProp(bconds, (const xmlChar*)"name", (const xmlChar*)"bconds");
386 xmlAddChild(cur, bconds);
387
388 xmlNodePtr LRDim = xmlNewTextChild(cur, NULL, (const xmlChar*)"parameter", (const xmlChar*)"15");
389 xmlNewProp(LRDim, (const xmlChar*)"name", (const xmlChar*)"LR_dim_cutoff");
390 xmlAddChild(cur, LRDim);
391
392 return cur;
393 }
createIonSet()394 xmlNodePtr QMCGaussianParserBase::createIonSet()
395 {
396 const double ang_to_bohr = 1.889725989;
397 if (!BohrUnit)
398 {
399 IonSystem.R *= ang_to_bohr;
400 }
401 double CoreTable[] = {
402 0, /* index zero*/
403 1, 2, /*H He */
404 2, 2, 2, 2, 2, 2, 2, 10, /*Li-Ne*/
405 10, 10, 10, 10, 10, 10, 10, 18, /*Na-Ar*/
406 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, /*N-Zn*/
407 28, 28, 28, 28, 28, 36, /*Ga-Kr*/
408 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, /*Rb-Cd*/
409 46, 46, 46, 46, 46, 54 /*In-Xe*/
410 };
411 SpeciesSet& ionSpecies(IonSystem.getSpeciesSet());
412 for (int i = 0; i < ionSpecies.getTotalNum(); i++)
413 {
414 int z = static_cast<int>(ionSpecies(AtomicNumberIndex, i));
415 double valence = ionSpecies(IonChargeIndex, i);
416 if (valence > CoreTable[z] && FixValence != true)
417 valence -= CoreTable[z];
418 ionSpecies(ValenceChargeIndex, i) = valence;
419 ionSpecies(AtomicNumberIndex, i) = z;
420 }
421 XMLSaveParticle o(IonSystem);
422 if (UseHDF5)
423 {
424 hdf_archive hout;
425 hout.open(h5file.c_str(), H5F_ACC_RDWR);
426 hout.push("atoms", true);
427 hout.write(NumberOfAtoms, "number_of_atoms");
428 auto nbspecies = ionSpecies.getTotalNum();
429 hout.write(nbspecies, "number_of_species");
430 Matrix<double> Position(NumberOfAtoms, 3);
431 std::vector<int> speciesID;
432 speciesID.resize(NumberOfAtoms);
433 for (auto i = 0; i < NumberOfAtoms; i++)
434 {
435 Position[i][0] = IonSystem.R[i][0];
436 Position[i][1] = IonSystem.R[i][1];
437 Position[i][2] = IonSystem.R[i][2];
438 speciesID[i] = IonSystem.GroupID[i];
439 }
440
441 hout.write(speciesID, "species_ids");
442 hout.write(Position, "positions");
443
444 for (auto i = 0; i < nbspecies; i++)
445 {
446 std::ostringstream SpecieID;
447 SpecieID << "species_" << i;
448 hout.push(SpecieID.str().c_str(), true);
449 int z = static_cast<int>(ionSpecies(AtomicNumberIndex, i));
450 double valence = ionSpecies(IonChargeIndex, i);
451 if (valence > CoreTable[z] && FixValence != true)
452 valence -= CoreTable[z];
453 hout.write(z, "charge");
454 hout.write(z, "atomic_number");
455 hout.write(valence, "core");
456 hout.write(IonName[static_cast<int>(ionSpecies(AtomicNumberIndex, i))], "name");
457 hout.pop();
458 }
459
460 hout.close();
461 }
462 return o.createNode(Periodicity);
463 }
464
465
createBasisSet()466 xmlNodePtr QMCGaussianParserBase::createBasisSet()
467 {
468 xmlNodePtr bset = xmlNewNode(NULL, (const xmlChar*)"basisset");
469 xmlNewProp(bset, (const xmlChar*)"name", (const xmlChar*)"LCAOBSet");
470 xmlNodePtr cur = NULL;
471 std::map<int, int> species;
472 int gtot = 0;
473 if (!debug)
474 for (int iat = 0; iat < NumberOfAtoms; iat++)
475 {
476 int itype = IonSystem.GroupID[iat];
477 int ng = 0;
478 std::map<int, int>::iterator it = species.find(itype);
479 if (it == species.end())
480 {
481 for (int ig = gBound[iat]; ig < gBound[iat + 1]; ig++)
482 {
483 ng += gNumber[ig];
484 }
485 species[itype] = ng;
486 if (cur)
487 {
488 cur = xmlAddSibling(cur, createCenter(iat, gtot));
489 }
490 else
491 {
492 cur = xmlAddChild(bset, createCenter(iat, gtot));
493 }
494 }
495 else
496 {
497 ng = (*it).second;
498 }
499 gtot += ng;
500 }
501 return bset;
502 }
503
504
createBasisSetWithHDF5()505 xmlNodePtr QMCGaussianParserBase::createBasisSetWithHDF5()
506 {
507 int counter = 0;
508
509 xmlNodePtr bset = xmlNewNode(NULL, (const xmlChar*)"basisset");
510 hdf_archive hout;
511 hout.open(h5file.c_str(), H5F_ACC_RDWR);
512 hout.push("basisset", true);
513 std::string BasisSetName("LCAOBSet");
514 hout.write(BasisSetName, "name");
515
516 std::map<int, int> species;
517 int gtot = 0;
518 for (int iat = 0; iat < NumberOfAtoms; iat++)
519 {
520 int itype = IonSystem.GroupID[iat];
521 int ng = 0;
522 std::map<int, int>::iterator it = species.find(itype);
523 if (it == species.end())
524 {
525 for (int ig = gBound[iat]; ig < gBound[iat + 1]; ig++)
526 {
527 ng += gNumber[ig];
528 }
529 species[itype] = ng;
530 createCenterH5(iat, gtot, counter);
531 counter++;
532 }
533 else
534 {
535 ng = (*it).second;
536 }
537 gtot += ng;
538 }
539 hout.write(counter, "NbElements");
540 hout.close();
541 return bset;
542 }
543
544
createDeterminantSetWithHDF5()545 xmlNodePtr QMCGaussianParserBase::createDeterminantSetWithHDF5()
546 {
547 setOccupationNumbers();
548 xmlNodePtr slaterdet = xmlNewNode(NULL, (const xmlChar*)"slaterdeterminant");
549 std::ostringstream up_size, down_size, b_size;
550 up_size << NumberOfAlpha;
551 down_size << NumberOfBeta;
552 b_size << numMO;
553 //create a determinant Up
554 xmlNodePtr udet = xmlNewNode(NULL, (const xmlChar*)"determinant");
555 xmlNewProp(udet, (const xmlChar*)"id", (const xmlChar*)"updet");
556 xmlNewProp(udet, (const xmlChar*)"size", (const xmlChar*)up_size.str().c_str());
557 if (DoCusp == true)
558 xmlNewProp(udet, (const xmlChar*)"cuspInfo", (const xmlChar*)"../updet.cuspInfo.xml");
559
560 //add occupation
561 xmlNodePtr occ_data = xmlNewNode(NULL, (const xmlChar*)"occupation");
562 xmlNewProp(occ_data, (const xmlChar*)"mode", (const xmlChar*)"ground");
563 xmlAddChild(udet, occ_data);
564 //add coefficients
565 xmlNodePtr coeff_data = xmlNewNode(NULL, (const xmlChar*)"coefficient");
566 xmlNewProp(coeff_data, (const xmlChar*)"size", (const xmlChar*)b_size.str().c_str());
567 xmlNewProp(coeff_data, (const xmlChar*)"spindataset", (const xmlChar*)"0");
568 xmlAddChild(udet, coeff_data);
569 //add udet to slaterdet
570 xmlNodePtr cur = xmlAddChild(slaterdet, udet);
571
572 hdf_archive hout;
573 hout.open(h5file.c_str(), H5F_ACC_RDWR);
574 hout.push("Super_Twist", true);
575
576 Matrix<double> Ctemp(numMO, SizeOfBasisSet);
577
578 int n = 0;
579 for (int i = 0; i < numMO; i++)
580 for (int j = 0; j < SizeOfBasisSet; j++)
581 {
582 Ctemp[i][j] = EigVec[n];
583 n++;
584 }
585
586 hout.write(Ctemp, "eigenset_0");
587
588 xmlNodePtr ddet;
589 if (SpinRestricted)
590 {
591 ddet = xmlCopyNode(udet, 1);
592 xmlSetProp(ddet, (const xmlChar*)"id", (const xmlChar*)"downdet");
593 xmlSetProp(ddet, (const xmlChar*)"size", (const xmlChar*)down_size.str().c_str());
594 if (DoCusp == true)
595 xmlSetProp(ddet, (const xmlChar*)"cuspInfo", (const xmlChar*)"../downdet.cuspInfo.xml");
596 }
597 else
598 {
599 ddet = xmlCopyNode(udet, 2);
600 xmlSetProp(ddet, (const xmlChar*)"id", (const xmlChar*)"downdet");
601 xmlSetProp(ddet, (const xmlChar*)"size", (const xmlChar*)down_size.str().c_str());
602 if (DoCusp == true)
603 xmlSetProp(ddet, (const xmlChar*)"cuspInfo", (const xmlChar*)"../downdet.cuspInfo.xml");
604 xmlNodePtr o = xmlAddChild(ddet, xmlCopyNode(occ_data, 1));
605 xmlNodePtr c = xmlCopyNode(coeff_data, 1);
606 xmlSetProp(c, (const xmlChar*)"spindataset", (const xmlChar*)"1");
607 o = xmlAddSibling(o, c);
608
609 n = numMO * SizeOfBasisSet;
610 for (int i = 0; i < numMO; i++)
611 for (int j = 0; j < SizeOfBasisSet; j++)
612 {
613 Ctemp[i][j] = EigVec[n];
614 n++;
615 }
616
617 hout.write(Ctemp, "eigenset_1");
618 }
619 cur = xmlAddSibling(cur, ddet);
620 hout.close();
621 return slaterdet;
622 }
623
624
PrepareDeterminantSetFromHDF5()625 xmlNodePtr QMCGaussianParserBase::PrepareDeterminantSetFromHDF5()
626 {
627 setOccupationNumbers();
628 xmlNodePtr slaterdet = xmlNewNode(NULL, (const xmlChar*)"slaterdeterminant");
629 std::ostringstream up_size, down_size, b_size;
630 up_size << NumberOfAlpha;
631 down_size << NumberOfBeta;
632 b_size << numMO;
633 //create a determinant Up
634 xmlNodePtr udet = xmlNewNode(NULL, (const xmlChar*)"determinant");
635 xmlNewProp(udet, (const xmlChar*)"id", (const xmlChar*)"updet");
636 xmlNewProp(udet, (const xmlChar*)"size", (const xmlChar*)up_size.str().c_str());
637 if (DoCusp == true)
638 xmlNewProp(udet, (const xmlChar*)"cuspInfo", (const xmlChar*)"../updet.cuspInfo.xml");
639
640 //add occupation
641 xmlNodePtr occ_data = xmlNewNode(NULL, (const xmlChar*)"occupation");
642 xmlNewProp(occ_data, (const xmlChar*)"mode", (const xmlChar*)"ground");
643 xmlAddChild(udet, occ_data);
644 //add coefficients
645 xmlNodePtr coeff_data = xmlNewNode(NULL, (const xmlChar*)"coefficient");
646 xmlNewProp(coeff_data, (const xmlChar*)"size", (const xmlChar*)b_size.str().c_str());
647 xmlNewProp(coeff_data, (const xmlChar*)"spindataset", (const xmlChar*)"0");
648 xmlAddChild(udet, coeff_data);
649 //add udet to slaterdet
650 xmlNodePtr cur = xmlAddChild(slaterdet, udet);
651
652 xmlNodePtr ddet;
653 if (SpinRestricted)
654 {
655 ddet = xmlCopyNode(udet, 1);
656 xmlSetProp(ddet, (const xmlChar*)"id", (const xmlChar*)"downdet");
657 xmlSetProp(ddet, (const xmlChar*)"size", (const xmlChar*)down_size.str().c_str());
658 if (DoCusp == true)
659 xmlSetProp(ddet, (const xmlChar*)"cuspInfo", (const xmlChar*)"../downdet.cuspInfo.xml");
660 }
661 else
662 {
663 ddet = xmlCopyNode(udet, 2);
664 xmlSetProp(ddet, (const xmlChar*)"id", (const xmlChar*)"downdet");
665 xmlSetProp(ddet, (const xmlChar*)"size", (const xmlChar*)down_size.str().c_str());
666 if (DoCusp == true)
667 xmlSetProp(ddet, (const xmlChar*)"cuspInfo", (const xmlChar*)"../downdet.cuspInfo.xml");
668 xmlNodePtr o = xmlAddChild(ddet, xmlCopyNode(occ_data, 1));
669 xmlNodePtr c = xmlCopyNode(coeff_data, 1);
670 xmlSetProp(c, (const xmlChar*)"spindataset", (const xmlChar*)"1");
671 o = xmlAddSibling(o, c);
672 }
673 cur = xmlAddSibling(cur, ddet);
674 return slaterdet;
675 }
676
createDeterminantSet()677 xmlNodePtr QMCGaussianParserBase::createDeterminantSet()
678 {
679 setOccupationNumbers();
680 xmlNodePtr slaterdet = xmlNewNode(NULL, (const xmlChar*)"slaterdeterminant");
681 //check spin-dependent properties
682 std::ostringstream up_size, down_size, b_size, occ;
683 up_size << NumberOfAlpha;
684 down_size << NumberOfBeta;
685 b_size << numMO;
686 //create a determinant Up
687 xmlNodePtr adet = xmlNewNode(NULL, (const xmlChar*)"determinant");
688 xmlNewProp(adet, (const xmlChar*)"id", (const xmlChar*)"updet");
689 xmlNewProp(adet, (const xmlChar*)"size", (const xmlChar*)up_size.str().c_str());
690 if (DoCusp == true)
691 xmlNewProp(adet, (const xmlChar*)"cuspInfo", (const xmlChar*)"../updet.cuspInfo.xml");
692
693 xmlNodePtr occ_data = xmlNewNode(NULL, (const xmlChar*)"occupation");
694 xmlNewProp(occ_data, (const xmlChar*)"mode", (const xmlChar*)"ground");
695 xmlAddChild(adet, occ_data);
696 int btot = numMO * SizeOfBasisSet;
697 int n = btot / 4, b = 0;
698 int dn = btot - n * 4;
699 std::ostringstream eig;
700 eig.setf(std::ios::scientific, std::ios::floatfield);
701 eig.setf(std::ios::right, std::ios::adjustfield);
702 eig.precision(14);
703 eig << "\n";
704 for (int k = 0; k < n; k++)
705 {
706 eig << std::setw(22) << EigVec[b] << std::setw(22) << EigVec[b + 1] << std::setw(22) << EigVec[b + 2]
707 << std::setw(22) << EigVec[b + 3] << "\n";
708 b += 4;
709 }
710 for (int k = 0; k < dn; k++)
711 {
712 eig << std::setw(22) << EigVec[b++];
713 }
714 if (dn)
715 eig << std::endl;
716 xmlNodePtr det_data = xmlNewTextChild(adet, NULL, (const xmlChar*)"coefficient", (const xmlChar*)eig.str().c_str());
717 xmlNewProp(det_data, (const xmlChar*)"size", (const xmlChar*)b_size.str().c_str());
718 xmlNewProp(det_data, (const xmlChar*)"id", (const xmlChar*)"updetC");
719 xmlNodePtr cur = xmlAddChild(slaterdet, adet);
720 adet = xmlNewNode(NULL, (const xmlChar*)"determinant");
721 xmlNewProp(adet, (const xmlChar*)"id", (const xmlChar*)"downdet");
722 xmlNewProp(adet, (const xmlChar*)"size", (const xmlChar*)down_size.str().c_str());
723 if (DoCusp == true)
724 xmlNewProp(adet, (const xmlChar*)"cuspInfo", (const xmlChar*)"../downdet.cuspInfo.xml");
725
726 {
727 occ_data = xmlNewNode(NULL, (const xmlChar*)"occupation");
728 xmlNewProp(occ_data, (const xmlChar*)"mode", (const xmlChar*)"ground");
729 xmlAddChild(adet, occ_data);
730 std::ostringstream eigD;
731 eigD.setf(std::ios::scientific, std::ios::floatfield);
732 eigD.setf(std::ios::right, std::ios::adjustfield);
733 eigD.precision(14);
734 eigD << "\n";
735 b = numMO * SizeOfBasisSet;
736 for (int k = 0; k < n; k++)
737 {
738 eigD << std::setw(22) << EigVec[b] << std::setw(22) << EigVec[b + 1] << std::setw(22) << EigVec[b + 2]
739 << std::setw(22) << EigVec[b + 3] << "\n";
740 b += 4;
741 }
742 for (int k = 0; k < dn; k++)
743 {
744 eigD << std::setw(22) << EigVec[b++];
745 }
746 if (dn)
747 eigD << std::endl;
748 if (SpinRestricted)
749 det_data = xmlNewTextChild(adet, NULL, (const xmlChar*)"coefficient", (const xmlChar*)eig.str().c_str());
750 else
751 det_data = xmlNewTextChild(adet, NULL, (const xmlChar*)"coefficient", (const xmlChar*)eigD.str().c_str());
752 xmlNewProp(det_data, (const xmlChar*)"size", (const xmlChar*)b_size.str().c_str());
753 xmlNewProp(det_data, (const xmlChar*)"id", (const xmlChar*)"downdetC");
754 }
755 cur = xmlAddSibling(cur, adet);
756 return slaterdet;
757 }
758
createSPOSets(xmlNodePtr spoUP,xmlNodePtr spoDN)759 void QMCGaussianParserBase::createSPOSets(xmlNodePtr spoUP, xmlNodePtr spoDN)
760 {
761 setOccupationNumbers();
762 std::ostringstream up_size, down_size, b_size, occ, nstates_alpha, nstates_beta;
763 up_size << NumberOfAlpha;
764 down_size << NumberOfBeta;
765 b_size << numMO;
766 nstates_alpha << ci_nstates + ci_nca;
767 nstates_beta << ci_nstates + ci_ncb;
768
769 xmlNewProp(spoUP, (const xmlChar*)"name", (const xmlChar*)"spo-up");
770 xmlNewProp(spoDN, (const xmlChar*)"name", (const xmlChar*)"spo-dn");
771 xmlNewProp(spoUP, (const xmlChar*)"size", (const xmlChar*)nstates_alpha.str().c_str());
772 xmlNewProp(spoDN, (const xmlChar*)"size", (const xmlChar*)nstates_beta.str().c_str());
773 if (DoCusp == true)
774 {
775 xmlNewProp(spoUP, (const xmlChar*)"cuspInfo", (const xmlChar*)"../spo-up.cuspInfo.xml");
776 xmlNewProp(spoDN, (const xmlChar*)"cuspInfo", (const xmlChar*)"../spo-dn.cuspInfo.xml");
777 }
778 xmlNodePtr occ_data = xmlNewNode(NULL, (const xmlChar*)"occupation");
779 xmlNewProp(occ_data, (const xmlChar*)"mode", (const xmlChar*)"ground");
780 xmlAddChild(spoUP, occ_data);
781 int btot = numMO * SizeOfBasisSet;
782 int n = btot / 4, b = 0;
783 int dn = btot - n * 4;
784 std::ostringstream eig;
785 eig.setf(std::ios::scientific, std::ios::floatfield);
786 eig.setf(std::ios::right, std::ios::adjustfield);
787 eig.precision(14);
788 eig << "\n";
789 for (int k = 0; k < n; k++)
790 {
791 eig << std::setw(22) << EigVec[b] << std::setw(22) << EigVec[b + 1] << std::setw(22) << EigVec[b + 2]
792 << std::setw(22) << EigVec[b + 3] << "\n";
793 b += 4;
794 }
795 for (int k = 0; k < dn; k++)
796 {
797 eig << std::setw(22) << EigVec[b++];
798 }
799 if (dn)
800 eig << std::endl;
801 xmlNodePtr det_data = xmlNewTextChild(spoUP, NULL, (const xmlChar*)"coefficient", (const xmlChar*)eig.str().c_str());
802 xmlNewProp(det_data, (const xmlChar*)"size", (const xmlChar*)b_size.str().c_str());
803 xmlNewProp(det_data, (const xmlChar*)"id", (const xmlChar*)"updetC");
804 {
805 occ_data = xmlNewNode(NULL, (const xmlChar*)"occupation");
806 xmlNewProp(occ_data, (const xmlChar*)"mode", (const xmlChar*)"ground");
807 xmlAddChild(spoDN, occ_data);
808 std::ostringstream eigD;
809 eigD.setf(std::ios::scientific, std::ios::floatfield);
810 eigD.setf(std::ios::right, std::ios::adjustfield);
811 eigD.precision(14);
812 eigD << "\n";
813 b = numMO * SizeOfBasisSet;
814 for (int k = 0; k < n; k++)
815 {
816 eigD << std::setw(22) << EigVec[b] << std::setw(22) << EigVec[b + 1] << std::setw(22) << EigVec[b + 2]
817 << std::setw(22) << EigVec[b + 3] << "\n";
818 b += 4;
819 }
820 for (int k = 0; k < dn; k++)
821 {
822 eigD << std::setw(22) << EigVec[b++];
823 }
824 if (dn)
825 eigD << std::endl;
826 if (SpinRestricted)
827 det_data = xmlNewTextChild(spoDN, NULL, (const xmlChar*)"coefficient", (const xmlChar*)eig.str().c_str());
828 else
829 det_data = xmlNewTextChild(spoDN, NULL, (const xmlChar*)"coefficient", (const xmlChar*)eigD.str().c_str());
830 xmlNewProp(det_data, (const xmlChar*)"size", (const xmlChar*)b_size.str().c_str());
831 xmlNewProp(det_data, (const xmlChar*)"id", (const xmlChar*)"downdetC");
832 }
833 }
834
createSPOSetsH5(xmlNodePtr spoUP,xmlNodePtr spoDN)835 void QMCGaussianParserBase::createSPOSetsH5(xmlNodePtr spoUP, xmlNodePtr spoDN)
836 {
837 setOccupationNumbers();
838 Matrix<double> Ctemp(numMO, SizeOfBasisSet);
839 int n = 0;
840 hdf_archive hout;
841 hout.open(h5file.c_str(), H5F_ACC_RDWR);
842 hout.push("Super_Twist", true);
843
844 std::ostringstream up_size, down_size, b_size, occ, nstates_alpha, nstates_beta;
845 up_size << NumberOfAlpha;
846 down_size << NumberOfBeta;
847 b_size << numMO;
848 nstates_alpha << ci_nstates + ci_nca;
849 ;
850 nstates_beta << ci_nstates + ci_ncb;
851
852 //create a spoUp
853 xmlNewProp(spoUP, (const xmlChar*)"name", (const xmlChar*)"spo-up");
854 xmlNewProp(spoUP, (const xmlChar*)"size", (const xmlChar*)nstates_alpha.str().c_str());
855
856 //create a spoDN
857 xmlNewProp(spoDN, (const xmlChar*)"name", (const xmlChar*)"spo-dn");
858 xmlNewProp(spoDN, (const xmlChar*)"size", (const xmlChar*)nstates_beta.str().c_str());
859
860
861 if (DoCusp == true)
862 {
863 xmlNewProp(spoUP, (const xmlChar*)"cuspInfo", (const xmlChar*)"../spo-up.cuspInfo.xml");
864 xmlNewProp(spoDN, (const xmlChar*)"cuspInfo", (const xmlChar*)"../spo-dn.cuspInfo.xml");
865 }
866
867
868 //add occupation UP
869 xmlNodePtr occ_data = xmlNewNode(NULL, (const xmlChar*)"occupation");
870 xmlNewProp(occ_data, (const xmlChar*)"mode", (const xmlChar*)"ground");
871 xmlAddChild(spoUP, occ_data);
872
873
874 //add coefficients
875 xmlNodePtr coeff_data = xmlNewNode(NULL, (const xmlChar*)"coefficient");
876 xmlNewProp(coeff_data, (const xmlChar*)"size", (const xmlChar*)b_size.str().c_str());
877 xmlNewProp(coeff_data, (const xmlChar*)"spindataset", (const xmlChar*)"0");
878 xmlAddChild(spoUP, coeff_data);
879
880
881 for (int i = 0; i < numMO; i++)
882 for (int j = 0; j < SizeOfBasisSet; j++)
883 {
884 Ctemp[i][j] = EigVec[n];
885 n++;
886 }
887
888
889 hout.write(Ctemp, "eigenset_0");
890
891
892 //add occupation DN
893 occ_data = xmlNewNode(NULL, (const xmlChar*)"occupation");
894 xmlNewProp(occ_data, (const xmlChar*)"mode", (const xmlChar*)"ground");
895 xmlAddChild(spoDN, occ_data);
896
897 coeff_data = xmlNewNode(NULL, (const xmlChar*)"coefficient");
898 xmlNewProp(coeff_data, (const xmlChar*)"size", (const xmlChar*)b_size.str().c_str());
899 if (SpinRestricted)
900 xmlNewProp(coeff_data, (const xmlChar*)"spindataset", (const xmlChar*)"0");
901 else
902 {
903 xmlNewProp(coeff_data, (const xmlChar*)"spindataset", (const xmlChar*)"1");
904 n = numMO * SizeOfBasisSet;
905 for (int i = 0; i < numMO; i++)
906 for (int j = 0; j < SizeOfBasisSet; j++)
907 {
908 Ctemp[i][j] = EigVec[n];
909 n++;
910 }
911 hout.write(Ctemp, "eigenset_1");
912 }
913 xmlAddChild(spoDN, coeff_data);
914
915 hout.close();
916 }
917
createMultiDeterminantSetCIHDF5()918 xmlNodePtr QMCGaussianParserBase::createMultiDeterminantSetCIHDF5()
919 {
920 xmlNodePtr multislaterdet = xmlNewNode(NULL, (const xmlChar*)"multideterminant");
921 if (optDetCoeffs)
922 xmlNewProp(multislaterdet, (const xmlChar*)"optimize", (const xmlChar*)"yes");
923 else
924 xmlNewProp(multislaterdet, (const xmlChar*)"optimize", (const xmlChar*)"no");
925 xmlNewProp(multislaterdet, (const xmlChar*)"spo_up", (const xmlChar*)"spo-up");
926 xmlNewProp(multislaterdet, (const xmlChar*)"spo_dn", (const xmlChar*)"spo-dn");
927 xmlNodePtr detlist = xmlNewNode(NULL, (const xmlChar*)"detlist");
928 std::ostringstream nstates, cisize, cinca, cincb, cinea, cineb, ci_thr;
929 cisize << ci_size;
930 nstates << ci_nstates;
931 cinca << ci_nca;
932 cincb << ci_ncb;
933 cinea << ci_nea;
934 cineb << ci_neb;
935 ci_thr << ci_threshold;
936 xmlNewProp(detlist, (const xmlChar*)"size", (const xmlChar*)cisize.str().c_str());
937 xmlNewProp(detlist, (const xmlChar*)"type", (const xmlChar*)"DETS");
938 xmlNewProp(detlist, (const xmlChar*)"nca", (const xmlChar*)cinca.str().c_str());
939 xmlNewProp(detlist, (const xmlChar*)"ncb", (const xmlChar*)cincb.str().c_str());
940 xmlNewProp(detlist, (const xmlChar*)"nea", (const xmlChar*)cinea.str().c_str());
941 xmlNewProp(detlist, (const xmlChar*)"neb", (const xmlChar*)cineb.str().c_str());
942 xmlNewProp(detlist, (const xmlChar*)"nstates", (const xmlChar*)nstates.str().c_str());
943 xmlNewProp(detlist, (const xmlChar*)"cutoff", (const xmlChar*)ci_thr.str().c_str());
944 xmlNewProp(detlist, (const xmlChar*)"href", (const xmlChar*)h5file.c_str());
945 if (CIcoeff.size() == 0)
946 {
947 std::cerr << " CI configuration list is empty. \n";
948 exit(101);
949 }
950 if (CIcoeff.size() != CIalpha.size() || CIcoeff.size() != CIbeta.size())
951 {
952 std::cerr << " Problem with CI configuration lists. \n";
953 exit(102);
954 }
955 /// 64 bit fixed width integer
956 const unsigned bit_kind = 64;
957 static_assert(bit_kind == sizeof(int64_t) * 8, "Must be 64 bit fixed width integer");
958 static_assert(bit_kind == sizeof(unsigned long long) * 8, "Must be 64 bit fixed width integer");
959 /// the number of 64 bit integers which represent the binary string for occupation
960 int N_int;
961 N_int = int(ci_nstates / bit_kind);
962 if (ci_nstates % bit_kind > 0)
963 N_int += 1;
964
965 hdf_archive hout;
966 hout.open(h5file.c_str(), H5F_ACC_RDWR);
967 hout.push("MultiDet", true);
968 hout.write(ci_size, "NbDet");
969 hout.write(ci_nstates, "nstate");
970 hout.write(CIcoeff, "Coeff");
971 hout.write(N_int, "Nbits");
972 hout.write(nbexcitedstates, "nexcitedstate");
973
974 Matrix<int64_t> tempAlpha(ci_size, N_int);
975 Matrix<int64_t> tempBeta(ci_size, N_int);
976 for (int i = 0; i < CIcoeff.size(); i++)
977 {
978 std::string loc_alpha = CIalpha[i].substr(0, ci_nstates);
979 std::string loc_beta = CIbeta[i].substr(0, ci_nstates);
980 std::string BiteSizeStringAlpha;
981 std::string BiteSizeStringBeta;
982 BiteSizeStringAlpha.resize(N_int * bit_kind);
983 BiteSizeStringBeta.resize(N_int * bit_kind);
984
985 for (std::size_t n = 0; n < (N_int * bit_kind); n++)
986 {
987 BiteSizeStringAlpha[n] = '0';
988 BiteSizeStringBeta[n] = '0';
989 }
990
991 for (std::size_t n = 0; n < ci_nstates; n++)
992 {
993 BiteSizeStringAlpha[n] = loc_alpha[n];
994 BiteSizeStringBeta[n] = loc_beta[n];
995 }
996
997 std::size_t offset = 0;
998 for (std::size_t l = 0; l < N_int; l++)
999 {
1000 offset = bit_kind * l;
1001 int64_t Val;
1002 std::string Var_alpha, Var_beta;
1003 Var_alpha.resize(bit_kind);
1004 Var_beta.resize(bit_kind);
1005
1006 for (auto j = 0; j < bit_kind; j++)
1007 {
1008 Var_alpha[j] = BiteSizeStringAlpha[j + offset];
1009 Var_beta[j] = BiteSizeStringBeta[j + offset];
1010 }
1011 std::reverse(Var_alpha.begin(), Var_alpha.end());
1012 std::reverse(Var_beta.begin(), Var_beta.end());
1013 std::bitset<bit_kind> bit_alpha(Var_alpha);
1014 std::bitset<bit_kind> bit_beta(Var_beta);
1015 tempAlpha[i][l] = bit_alpha.to_ullong();
1016 tempBeta[i][l] = bit_beta.to_ullong();
1017 }
1018 }
1019 hout.write(tempAlpha, "CI_Alpha");
1020 hout.write(tempBeta, "CI_Beta");
1021
1022 hout.pop();
1023 xmlAddChild(multislaterdet, detlist);
1024 hout.close();
1025 return multislaterdet;
1026 }
1027
createMultiDeterminantSet()1028 xmlNodePtr QMCGaussianParserBase::createMultiDeterminantSet()
1029 {
1030 xmlNodePtr multislaterdet = xmlNewNode(NULL, (const xmlChar*)"multideterminant");
1031 if (optDetCoeffs)
1032 xmlNewProp(multislaterdet, (const xmlChar*)"optimize", (const xmlChar*)"yes");
1033 else
1034 xmlNewProp(multislaterdet, (const xmlChar*)"optimize", (const xmlChar*)"no");
1035 xmlNewProp(multislaterdet, (const xmlChar*)"spo_up", (const xmlChar*)"spo-up");
1036 xmlNewProp(multislaterdet, (const xmlChar*)"spo_dn", (const xmlChar*)"spo-dn");
1037 if (usingCSF)
1038 {
1039 xmlNodePtr detlist = xmlNewNode(NULL, (const xmlChar*)"detlist");
1040 std::ostringstream nstates, cisize, cinca, cincb, cinea, cineb, ci_thr;
1041 cisize << ci_size;
1042 nstates << ci_nstates;
1043 cinca << ci_nca;
1044 cincb << ci_ncb;
1045 cinea << ci_nea;
1046 cineb << ci_neb;
1047 ci_thr << ci_threshold;
1048 xmlNewProp(detlist, (const xmlChar*)"size", (const xmlChar*)cisize.str().c_str());
1049 xmlNewProp(detlist, (const xmlChar*)"type", (const xmlChar*)"CSF");
1050 xmlNewProp(detlist, (const xmlChar*)"nca", (const xmlChar*)cinca.str().c_str());
1051 xmlNewProp(detlist, (const xmlChar*)"ncb", (const xmlChar*)cincb.str().c_str());
1052 xmlNewProp(detlist, (const xmlChar*)"nea", (const xmlChar*)cinea.str().c_str());
1053 xmlNewProp(detlist, (const xmlChar*)"neb", (const xmlChar*)cineb.str().c_str());
1054 xmlNewProp(detlist, (const xmlChar*)"nstates", (const xmlChar*)nstates.str().c_str());
1055 xmlNewProp(detlist, (const xmlChar*)"cutoff", (const xmlChar*)ci_thr.str().c_str());
1056 CIexcitLVL.clear();
1057 for (int i = 0; i < CSFocc.size(); i++)
1058 {
1059 CIexcitLVL.push_back(numberOfExcitationsCSF(CSFocc[i]));
1060 //cout<<CSFocc[i] <<" " <<CIexcitLVL.back() << std::endl;
1061 }
1062 // order dets according to ci coeff
1063 std::vector<std::pair<double, int>> order;
1064 if (orderByExcitation)
1065 {
1066 std::cout << "Ordering csfs by excitation level. \n";
1067 int maxE = *max_element(CIexcitLVL.begin(), CIexcitLVL.end());
1068 std::vector<int> pos(maxE);
1069 int ip1, ip2, cnt = 0, cnt2;
1070 // add by excitations, and do partial sorts of the list
1071 // messy but I dont want std::pair< std::pair<> > types right now
1072 for (int i = maxE; i >= 0; i--)
1073 {
1074 ip1 = ip2 = cnt;
1075 cnt2 = 0;
1076 for (int k = 0; k < CIexcitLVL.size(); k++)
1077 {
1078 if (CIexcitLVL[k] == i)
1079 {
1080 std::pair<double, int> cic(std::abs(coeff2csf[k].second), k);
1081 order.push_back(cic);
1082 cnt2++;
1083 cnt++;
1084 }
1085 }
1086 if (cnt2 > 0)
1087 sort(order.begin() + ip1, order.end());
1088 }
1089 }
1090 else
1091 {
1092 for (int i = 0; i < coeff2csf.size(); i++)
1093 {
1094 std::pair<double, int> cic(std::abs(coeff2csf[i].second), i);
1095 order.push_back(cic);
1096 }
1097 sort(order.begin(), order.end());
1098 }
1099 std::vector<std::pair<double, int>>::reverse_iterator it(order.rbegin());
1100 std::vector<std::pair<double, int>>::reverse_iterator last(order.rend());
1101 int iv = 0;
1102 while (it != last)
1103 {
1104 int nq = (*it).second;
1105 xmlNodePtr csf = xmlNewNode(NULL, (const xmlChar*)"csf");
1106 std::ostringstream qc_coeff;
1107 qc_coeff << coeff2csf[nq].second;
1108 std::ostringstream coeff;
1109 std::ostringstream exct;
1110 exct << CIexcitLVL[nq];
1111 if (zeroCI && iv == 0)
1112 {
1113 coeff << 1.0;
1114 }
1115 else if (zeroCI && iv > 0)
1116 {
1117 coeff << 0.0;
1118 }
1119 else
1120 {
1121 coeff << coeff2csf[nq].second;
1122 }
1123 std::ostringstream tag;
1124 tag << "CSFcoeff_" << iv;
1125 xmlNewProp(csf, (const xmlChar*)"id", (const xmlChar*)tag.str().c_str());
1126 xmlNewProp(csf, (const xmlChar*)"exctLvl", (const xmlChar*)exct.str().c_str());
1127 xmlNewProp(csf, (const xmlChar*)"coeff", (const xmlChar*)coeff.str().c_str());
1128 xmlNewProp(csf, (const xmlChar*)"qchem_coeff", (const xmlChar*)qc_coeff.str().c_str());
1129 xmlNewProp(csf, (const xmlChar*)"occ", (const xmlChar*)CSFocc[nq].substr(0, ci_nstates).c_str());
1130 for (int i = 0; i < CSFexpansion[nq].size(); i++)
1131 {
1132 xmlNodePtr ci = xmlNewNode(NULL, (const xmlChar*)"det");
1133 std::ostringstream coeff0;
1134 coeff0 << CSFexpansion[nq][i];
1135 std::ostringstream tag0;
1136 tag0 << "csf_" << iv << "-" << i;
1137 xmlNewProp(ci, (const xmlChar*)"id", (const xmlChar*)tag0.str().c_str());
1138 xmlNewProp(ci, (const xmlChar*)"coeff", (const xmlChar*)coeff0.str().c_str());
1139 xmlNewProp(ci, (const xmlChar*)"alpha", (const xmlChar*)CSFalpha[nq][i].substr(0, ci_nstates).c_str());
1140 xmlNewProp(ci, (const xmlChar*)"beta", (const xmlChar*)CSFbeta[nq][i].substr(0, ci_nstates).c_str());
1141 xmlAddChild(csf, ci);
1142 }
1143 xmlAddChild(detlist, csf);
1144 it++;
1145 iv++;
1146 }
1147 xmlAddChild(multislaterdet, detlist);
1148 }
1149 else
1150 // usingCSF
1151 {
1152 xmlNodePtr detlist = xmlNewNode(NULL, (const xmlChar*)"detlist");
1153 ci_size = 0;
1154 for (int i = 0; i < CIcoeff.size(); i++)
1155 if (std::abs(CIcoeff[i]) > ci_threshold)
1156 ci_size++;
1157 std::ostringstream nstates, cisize, cinca, cincb, cinea, cineb, ci_thr;
1158 cisize << ci_size;
1159 nstates << ci_nstates;
1160 cinca << ci_nca;
1161 cincb << ci_ncb;
1162 cinea << ci_nea;
1163 cineb << ci_neb;
1164 ci_thr << ci_threshold;
1165 xmlNewProp(detlist, (const xmlChar*)"size", (const xmlChar*)cisize.str().c_str());
1166 xmlNewProp(detlist, (const xmlChar*)"type", (const xmlChar*)"DETS");
1167 xmlNewProp(detlist, (const xmlChar*)"nca", (const xmlChar*)cinca.str().c_str());
1168 xmlNewProp(detlist, (const xmlChar*)"ncb", (const xmlChar*)cincb.str().c_str());
1169 xmlNewProp(detlist, (const xmlChar*)"nea", (const xmlChar*)cinea.str().c_str());
1170 xmlNewProp(detlist, (const xmlChar*)"neb", (const xmlChar*)cineb.str().c_str());
1171 xmlNewProp(detlist, (const xmlChar*)"nstates", (const xmlChar*)nstates.str().c_str());
1172 xmlNewProp(detlist, (const xmlChar*)"cutoff", (const xmlChar*)ci_thr.str().c_str());
1173 if (CIcoeff.size() == 0)
1174 {
1175 std::cerr << " CI configuration list is empty. \n";
1176 exit(101);
1177 }
1178 if (CIcoeff.size() != CIalpha.size() || CIcoeff.size() != CIbeta.size())
1179 {
1180 std::cerr << " Problem with CI configuration lists. \n";
1181 exit(102);
1182 }
1183 int iv = 0;
1184 for (int i = 0; i < CIcoeff.size(); i++)
1185 {
1186 if (std::abs(CIcoeff[i]) > ci_threshold)
1187 {
1188 xmlNodePtr ci = xmlNewNode(NULL, (const xmlChar*)"ci");
1189 std::ostringstream coeff;
1190 std::ostringstream qc_coeff;
1191 qc_coeff << CIcoeff[i];
1192 if (zeroCI && i == 0)
1193 {
1194 coeff << 1.0;
1195 }
1196 else if (zeroCI && i > 0)
1197 {
1198 coeff << 0.0;
1199 }
1200 else
1201 {
1202 coeff << CIcoeff[i];
1203 }
1204 std::ostringstream tag;
1205 tag << "CIcoeff_" << iv++;
1206 xmlNewProp(ci, (const xmlChar*)"id", (const xmlChar*)tag.str().c_str());
1207 xmlNewProp(ci, (const xmlChar*)"coeff", (const xmlChar*)coeff.str().c_str());
1208 xmlNewProp(ci, (const xmlChar*)"qc_coeff", (const xmlChar*)qc_coeff.str().c_str());
1209 xmlNewProp(ci, (const xmlChar*)"alpha", (const xmlChar*)CIalpha[i].substr(0, ci_nstates).c_str());
1210 xmlNewProp(ci, (const xmlChar*)"beta", (const xmlChar*)CIbeta[i].substr(0, ci_nstates).c_str());
1211 xmlAddChild(detlist, ci);
1212 }
1213 }
1214 xmlAddChild(multislaterdet, detlist);
1215 } //usingCSF
1216 return multislaterdet;
1217 }
1218
createCenterH5(int iat,int off_,int numelem)1219 void QMCGaussianParserBase::createCenterH5(int iat, int off_, int numelem)
1220 {
1221 CurrentCenter = GroupName[iat];
1222 int numbasisgroups(0);
1223 std::stringstream tempcenter;
1224 std::string CenterID;
1225 std::string expandYlm("Gamess");
1226 tempcenter << CurrentCenter << "";
1227 CenterID = tempcenter.str();
1228 std::stringstream tempElem;
1229 std::string ElemID0 = "atomicBasisSet", ElemID;
1230 tempElem << ElemID0 << numelem;
1231 ElemID = tempElem.str();
1232
1233 hdf_archive hout;
1234 hout.open(h5file.c_str(), H5F_ACC_RDWR);
1235 hout.push("basisset");
1236 hout.push(ElemID.c_str(), true);
1237 hout.write(basisName, "name");
1238 hout.write(angular_type, "angular");
1239 hout.write(CenterID, "elementType");
1240 hout.write(expandYlm, "expandYlm");
1241 hout.write(Normalized, "normalized");
1242
1243 double ValgridFirst(1.e-6), ValgridLast(1.e2);
1244 int Valnpts(1001);
1245 std::string gridType("log");
1246 double gridFirst(ValgridFirst);
1247 double gridLast(ValgridLast);
1248 int gridSize(Valnpts);
1249
1250 hout.write(gridType, "grid_type");
1251 hout.write(gridFirst, "grid_ri");
1252 hout.write(gridLast, "grid_rf");
1253 hout.write(gridSize, "grid_npts");
1254
1255 for (int ig = gBound[iat], n = 0; ig < gBound[iat + 1]; ig++, n++)
1256 {
1257 createShellH5(n, ig, off_, numelem);
1258 off_ += gNumber[ig];
1259 numbasisgroups = n + 1;
1260 }
1261
1262 hout.write(numbasisgroups, "NbBasisGroups");
1263 hout.close();
1264 }
1265
1266
createCenter(int iat,int off_)1267 xmlNodePtr QMCGaussianParserBase::createCenter(int iat, int off_)
1268 {
1269 CurrentCenter = GroupName[iat];
1270
1271 xmlNodePtr abasis = xmlNewNode(NULL, (const xmlChar*)"atomicBasisSet");
1272 xmlNewProp(abasis, (const xmlChar*)"name", (const xmlChar*)basisName.c_str());
1273 //xmlNewProp(abasis,(const xmlChar*)"angular",(const xmlChar*)"spherical");
1274 xmlNewProp(abasis, (const xmlChar*)"angular", (const xmlChar*)angular_type.c_str());
1275 xmlNewProp(abasis, (const xmlChar*)"type", (const xmlChar*)basisType.c_str());
1276 xmlNewProp(abasis, (const xmlChar*)"elementType", (const xmlChar*)CurrentCenter.c_str());
1277 xmlNewProp(abasis, (const xmlChar*)"normalized", (const xmlChar*)Normalized.c_str());
1278 xmlAddChild(abasis, xmlCopyNode(gridPtr, 1));
1279 for (int ig = gBound[iat], n = 0; ig < gBound[iat + 1]; ig++, n++)
1280 {
1281 createShell(n, ig, off_, abasis);
1282 off_ += gNumber[ig];
1283 }
1284 return abasis;
1285 }
1286
1287
createShellH5(int n,int ig,int off_,int numelem)1288 void QMCGaussianParserBase::createShellH5(int n, int ig, int off_, int numelem)
1289 {
1290 int gid(gShell[ig]);
1291 int ng(gNumber[ig]);
1292
1293
1294 char l_name[4], n_name[4], a_name[32];
1295 sprintf(a_name, "%s%d%d", CurrentCenter.c_str(), n, gShellID[gid]);
1296 sprintf(l_name, "%d", gShellID[gid]);
1297 sprintf(n_name, "%d", n);
1298
1299 std::string aa_name(a_name);
1300 std::string an_name(n_name);
1301 std::string al_name(l_name);
1302 std::string at_name("Gaussian");
1303 std::string basisGroupID = "basisGroup" + an_name;
1304
1305 std::stringstream tempElem;
1306 std::string ElemID0 = "atomicBasisSet", ElemID;
1307 tempElem << ElemID0 << numelem;
1308 ElemID = tempElem.str();
1309
1310 hdf_archive hout;
1311 hout.open(h5file.c_str(), H5F_ACC_RDWR);
1312 hout.push("basisset");
1313 hout.push(ElemID.c_str());
1314 hout.push(basisGroupID.c_str(), true);
1315 hout.write(aa_name, "rid");
1316 hout.write(n, "n");
1317 hout.write(gShellID[gid], "l");
1318 hout.write(at_name, "type");
1319 hout.write(ng, "NbRadFunc");
1320
1321 hout.push("radfunctions", true);
1322
1323 if (gid == 2)
1324 {
1325 std::string basisGroupID2 = "basisGroup2" + aa_name;
1326 std::string valac("1");
1327 hout.push(basisGroupID2.c_str(), true);
1328 hout.write(aa_name, "rid");
1329 hout.write(an_name, "n");
1330 hout.write(valac, "l");
1331 hout.write(at_name, "type");
1332 hout.write(ng, "NbRadFunc");
1333 hout.push("radfunctions2", true);
1334 hout.pop();
1335 hout.pop();
1336 }
1337
1338 for (int ig = 0, i = off_; ig < ng; ig++, i++)
1339 {
1340 std::stringstream tempdata;
1341 std::string dataradID0 = "DataRad", dataradID;
1342 tempdata << dataradID0 << ig;
1343 dataradID = tempdata.str();
1344 hout.push(dataradID.c_str(), true);
1345 hout.write(gExp[i], "exponent");
1346 hout.write(gC0[i], "contraction");
1347 hout.pop();
1348 if (gid == 2)
1349 {
1350 std::string basisGroupID2 = "basisGroup2" + aa_name;
1351 hout.push(basisGroupID2.c_str());
1352 hout.push("radfunctions2");
1353 std::stringstream tempdata2;
1354 std::string datarad2ID0 = "DataRad2", datarad2ID;
1355 tempdata2 << datarad2ID0 << ig;
1356 datarad2ID = tempdata2.str();
1357 hout.push(datarad2ID.c_str(), true);
1358 hout.write(gExp[i], "exponent");
1359 hout.write(gC1[i], "contraction");
1360 hout.pop();
1361 hout.pop();
1362 }
1363 }
1364 hout.close();
1365 }
1366
1367
createShell(int n,int ig,int off_,xmlNodePtr abasis)1368 void QMCGaussianParserBase::createShell(int n, int ig, int off_, xmlNodePtr abasis)
1369 {
1370 int gid(gShell[ig]);
1371 int ng(gNumber[ig]);
1372 xmlNodePtr ag = xmlNewNode(NULL, (const xmlChar*)"basisGroup");
1373 xmlNodePtr ag1 = 0;
1374 char l_name[4], n_name[4], a_name[32];
1375 sprintf(a_name, "%s%d%d", CurrentCenter.c_str(), n, gShellID[gid]);
1376 sprintf(l_name, "%d", gShellID[gid]);
1377 sprintf(n_name, "%d", n);
1378 xmlNewProp(ag, (const xmlChar*)"rid", (const xmlChar*)a_name);
1379 xmlNewProp(ag, (const xmlChar*)"n", (const xmlChar*)n_name);
1380 xmlNewProp(ag, (const xmlChar*)"l", (const xmlChar*)l_name);
1381 xmlNewProp(ag, (const xmlChar*)"type", (const xmlChar*)"Gaussian");
1382 if (gid == 2)
1383 {
1384 sprintf(a_name, "%s%d1", CurrentCenter.c_str(), n);
1385 ag1 = xmlNewNode(NULL, (const xmlChar*)"basisGroup");
1386 xmlNewProp(ag1, (const xmlChar*)"rid", (const xmlChar*)a_name);
1387 xmlNewProp(ag1, (const xmlChar*)"n", (const xmlChar*)n_name);
1388 xmlNewProp(ag1, (const xmlChar*)"l", (const xmlChar*)"1");
1389 xmlNewProp(ag1, (const xmlChar*)"type", (const xmlChar*)"Gaussian");
1390 }
1391 for (int ig = 0, i = off_; ig < ng; ig++, i++)
1392 {
1393 std::ostringstream a, b, c;
1394 a.setf(std::ios::scientific, std::ios::floatfield);
1395 b.setf(std::ios::scientific, std::ios::floatfield);
1396 a.precision(12);
1397 b.precision(12);
1398 a << gExp[i];
1399 b << gC0[i];
1400 xmlNodePtr anode = xmlNewNode(NULL, (const xmlChar*)"radfunc");
1401 xmlNewProp(anode, (const xmlChar*)"exponent", (const xmlChar*)a.str().c_str());
1402 xmlNewProp(anode, (const xmlChar*)"contraction", (const xmlChar*)b.str().c_str());
1403 xmlAddChild(ag, anode);
1404 if (gid == 2)
1405 {
1406 c.setf(std::ios::scientific, std::ios::floatfield);
1407 c.precision(12);
1408 c << gC1[i];
1409 anode = xmlNewNode(NULL, (const xmlChar*)"radfunc");
1410 xmlNewProp(anode, (const xmlChar*)"exponent", (const xmlChar*)a.str().c_str());
1411 xmlNewProp(anode, (const xmlChar*)"contraction", (const xmlChar*)c.str().c_str());
1412 xmlAddChild(ag1, anode);
1413 }
1414 }
1415 xmlAddChild(abasis, ag);
1416 if (gid == 2)
1417 xmlAddChild(abasis, ag1);
1418 }
1419
1420
createJ3()1421 xmlNodePtr QMCGaussianParserBase::createJ3()
1422 {
1423 xmlNodePtr j3 = xmlNewNode(NULL, (const xmlChar*)"jastrow");
1424 xmlNewProp(j3, (const xmlChar*)"name", (const xmlChar*)"J3");
1425 xmlNewProp(j3, (const xmlChar*)"type", (const xmlChar*)"eeI");
1426 xmlNewProp(j3, (const xmlChar*)"function", (const xmlChar*)"polynomial");
1427 xmlNewProp(j3, (const xmlChar*)"source", (const xmlChar*)"ion0");
1428 xmlNewProp(j3, (const xmlChar*)"print", (const xmlChar*)"yes");
1429 SpeciesSet& ionSpecies(IonSystem.getSpeciesSet());
1430 for (int i = 0; i < ionSpecies.getTotalNum(); i++)
1431 {
1432 xmlNodePtr uuc = xmlNewNode(NULL, (const xmlChar*)"correlation");
1433 xmlNewProp(uuc, (const xmlChar*)"ispecies", (const xmlChar*)ionSpecies.speciesName[i].c_str());
1434 xmlNewProp(uuc, (const xmlChar*)"especies", (const xmlChar*)"u");
1435 xmlNewProp(uuc, (const xmlChar*)"isize", (const xmlChar*)"3");
1436 xmlNewProp(uuc, (const xmlChar*)"esize", (const xmlChar*)"3");
1437 if (!PBC)
1438 xmlNewProp(uuc, (const xmlChar*)"rcut", (const xmlChar*)"5");
1439
1440 xmlNodePtr a = xmlNewTextChild(uuc, NULL, (const xmlChar*)"coefficients", (const xmlChar*)"\n ");
1441 std::ostringstream o1;
1442 o1 << "uu" << ionSpecies.speciesName[i];
1443 xmlNewProp(a, (const xmlChar*)"id", (const xmlChar*)o1.str().c_str());
1444 xmlNewProp(a, (const xmlChar*)"type", (const xmlChar*)"Array");
1445 xmlNewProp(a, (const xmlChar*)"optimize", (const xmlChar*)"yes");
1446 xmlAddChild(j3, uuc);
1447
1448 xmlNodePtr udc = xmlNewNode(NULL, (const xmlChar*)"correlation");
1449 xmlNewProp(udc, (const xmlChar*)"ispecies", (const xmlChar*)ionSpecies.speciesName[i].c_str());
1450 xmlNewProp(udc, (const xmlChar*)"especies1", (const xmlChar*)"u");
1451 xmlNewProp(udc, (const xmlChar*)"especies2", (const xmlChar*)"d");
1452 xmlNewProp(udc, (const xmlChar*)"isize", (const xmlChar*)"3");
1453 xmlNewProp(udc, (const xmlChar*)"esize", (const xmlChar*)"3");
1454 if (!PBC)
1455 xmlNewProp(udc, (const xmlChar*)"rcut", (const xmlChar*)"5");
1456
1457 xmlNodePtr b = xmlNewTextChild(udc, NULL, (const xmlChar*)"coefficients", (const xmlChar*)"\n ");
1458 std::ostringstream o2;
1459 o2 << "ud" << ionSpecies.speciesName[i];
1460 xmlNewProp(b, (const xmlChar*)"id", (const xmlChar*)o2.str().c_str());
1461 xmlNewProp(b, (const xmlChar*)"type", (const xmlChar*)"Array");
1462 xmlNewProp(b, (const xmlChar*)"optimize", (const xmlChar*)"yes");
1463 xmlAddChild(j3, udc);
1464 }
1465 return j3;
1466 }
1467
createJ2()1468 xmlNodePtr QMCGaussianParserBase::createJ2()
1469 {
1470 xmlNodePtr j2 = xmlNewNode(NULL, (const xmlChar*)"jastrow");
1471 xmlNewProp(j2, (const xmlChar*)"name", (const xmlChar*)"J2");
1472 xmlNewProp(j2, (const xmlChar*)"type", (const xmlChar*)"Two-Body");
1473 xmlNewProp(j2, (const xmlChar*)"function", (const xmlChar*)"Bspline");
1474 xmlNewProp(j2, (const xmlChar*)"print", (const xmlChar*)"yes");
1475 if (NumberOfAlpha > 1 || NumberOfBeta > 1)
1476 {
1477 xmlNodePtr uu = xmlNewNode(NULL, (const xmlChar*)"correlation");
1478 if (!PBC)
1479 xmlNewProp(uu, (const xmlChar*)"rcut", (const xmlChar*)"10");
1480 xmlNewProp(uu, (const xmlChar*)"size", (const xmlChar*)"10");
1481 xmlNewProp(uu, (const xmlChar*)"speciesA", (const xmlChar*)"u");
1482 xmlNewProp(uu, (const xmlChar*)"speciesB", (const xmlChar*)"u");
1483 xmlNodePtr a = xmlNewTextChild(uu, NULL, (const xmlChar*)"coefficients",
1484 (const xmlChar*)"0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0");
1485 xmlNewProp(a, (const xmlChar*)"id", (const xmlChar*)"uu");
1486 xmlNewProp(a, (const xmlChar*)"type", (const xmlChar*)"Array");
1487 xmlAddChild(j2, uu);
1488 }
1489 if (NumberOfAlpha > 0 && NumberOfBeta > 0)
1490 {
1491 xmlNodePtr uu = xmlNewNode(NULL, (const xmlChar*)"correlation");
1492 if (!PBC)
1493 xmlNewProp(uu, (const xmlChar*)"rcut", (const xmlChar*)"10");
1494 xmlNewProp(uu, (const xmlChar*)"size", (const xmlChar*)"10");
1495 xmlNewProp(uu, (const xmlChar*)"speciesA", (const xmlChar*)"u");
1496 xmlNewProp(uu, (const xmlChar*)"speciesB", (const xmlChar*)"d");
1497 //xmlNodePtr a = xmlNewNode(NULL,(const xmlChar*)"coefficients");
1498 xmlNodePtr a = xmlNewTextChild(uu, NULL, (const xmlChar*)"coefficients",
1499 (const xmlChar*)"0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0");
1500 xmlNewProp(a, (const xmlChar*)"id", (const xmlChar*)"ud");
1501 xmlNewProp(a, (const xmlChar*)"type", (const xmlChar*)"Array");
1502 // xmlAddChild(uu,a);
1503 xmlAddChild(j2, uu);
1504 }
1505 return j2;
1506 }
1507
createJ1()1508 xmlNodePtr QMCGaussianParserBase::createJ1()
1509 {
1510 xmlNodePtr j1 = xmlNewNode(NULL, (const xmlChar*)"jastrow");
1511 xmlNewProp(j1, (const xmlChar*)"name", (const xmlChar*)"J1");
1512 xmlNewProp(j1, (const xmlChar*)"type", (const xmlChar*)"One-Body");
1513 xmlNewProp(j1, (const xmlChar*)"function", (const xmlChar*)"Bspline");
1514 xmlNewProp(j1, (const xmlChar*)"source", (const xmlChar*)"ion0");
1515 xmlNewProp(j1, (const xmlChar*)"print", (const xmlChar*)"yes");
1516 SpeciesSet& ionSpecies(IonSystem.getSpeciesSet());
1517 for (int i = 0; i < ionSpecies.getTotalNum(); i++)
1518 {
1519 xmlNodePtr c = xmlNewNode(NULL, (const xmlChar*)"correlation");
1520 if (!PBC)
1521 xmlNewProp(c, (const xmlChar*)"rcut", (const xmlChar*)"10");
1522 xmlNewProp(c, (const xmlChar*)"size", (const xmlChar*)"10");
1523 xmlNewProp(c, (const xmlChar*)"cusp", (const xmlChar*)"0");
1524 xmlNewProp(c, (const xmlChar*)"elementType", (const xmlChar*)ionSpecies.speciesName[i].c_str());
1525 xmlNodePtr a = xmlNewTextChild(c, NULL, (const xmlChar*)"coefficients",
1526 (const xmlChar*)"0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0");
1527 std::ostringstream o;
1528 o << 'e' << ionSpecies.speciesName[i];
1529 xmlNewProp(a, (const xmlChar*)"id", (const xmlChar*)o.str().c_str());
1530 xmlNewProp(a, (const xmlChar*)"type", (const xmlChar*)"Array");
1531 xmlAddChild(j1, c);
1532 }
1533 return j1;
1534 }
1535
createGridNode(int argc,char ** argv)1536 void QMCGaussianParserBase::createGridNode(int argc, char** argv)
1537 {
1538 gridPtr = xmlNewNode(NULL, (const xmlChar*)"grid");
1539 std::string gridType("log");
1540 std::string gridFirst("1.e-6");
1541 std::string gridLast("1.e2");
1542 std::string gridSize("1001");
1543 int iargc = 0;
1544 while (iargc < argc)
1545 {
1546 std::string a(argv[iargc]);
1547 if (a == "-gridtype")
1548 {
1549 gridType = argv[++iargc];
1550 }
1551 else if (a == "-frst")
1552 {
1553 gridFirst = argv[++iargc];
1554 }
1555 else if (a == "-last")
1556 {
1557 gridLast = argv[++iargc];
1558 }
1559 else if (a == "-size")
1560 {
1561 gridSize = argv[++iargc];
1562 }
1563 else if (a == "-numMO")
1564 {
1565 numMO2print = atoi(argv[++iargc]);
1566 }
1567 ++iargc;
1568 }
1569 xmlNewProp(gridPtr, (const xmlChar*)"type", (const xmlChar*)gridType.c_str());
1570 xmlNewProp(gridPtr, (const xmlChar*)"ri", (const xmlChar*)gridFirst.c_str());
1571 xmlNewProp(gridPtr, (const xmlChar*)"rf", (const xmlChar*)gridLast.c_str());
1572 xmlNewProp(gridPtr, (const xmlChar*)"npts", (const xmlChar*)gridSize.c_str());
1573 }
1574
dump(const std::string & psi_tag,const std::string & ion_tag)1575 void QMCGaussianParserBase::dump(const std::string& psi_tag, const std::string& ion_tag)
1576 {
1577 std::cout << " QMCGaussianParserBase::dump " << std::endl;
1578 if (!Structure)
1579 {
1580 //if (UseHDF5 || multidetH5)
1581 if (UseHDF5)
1582 {
1583 bool IsComplex = false;
1584 hdf_archive hout;
1585 hout.create(h5file.c_str(), H5F_ACC_TRUNC);
1586 hout.push("PBC", true);
1587 hout.write(PBC, "PBC");
1588 hout.pop();
1589 //Adding generic code name to the H5 file.
1590 std::string CodeName("generic");
1591 hout.push("application", true);
1592 hout.write(CodeName, "code");
1593 hout.pop();
1594 hout.push("parameters", true);
1595 hout.write(ECP, "ECP");
1596 //Assumes MO-Coeff always real as this path is only for molecules and for generating stand alone H5file.
1597 hout.write(IsComplex, "IsComplex");
1598 hout.write(multideterminant, "Multidet");
1599 hout.write(NumberOfAlpha, "NbAlpha");
1600 hout.write(NumberOfBeta, "NbBeta");
1601 hout.write(NumberOfEls, "NbTotElec");
1602 hout.write(SpinRestricted, "SpinRestricted");
1603 hout.write(BohrUnit, "Unit");
1604 hout.write(numMO, "numAO");
1605 hout.write(numMO, "numMO");
1606 hout.write(SpinMultiplicity, "spin");
1607 hout.pop();
1608 hout.close();
1609 }
1610
1611 xmlDocPtr doc_p = xmlNewDoc((const xmlChar*)"1.0");
1612 xmlNodePtr qm_root_p = xmlNewNode(NULL, BAD_CAST "qmcsystem");
1613 if (PBC)
1614 {
1615 app_log()
1616 << "ABORT::THIS IS NOT SUPPOSED TO HAPPEN. PBC are ON but you are not in an HDF5 path. Contact developers"
1617 << std::endl;
1618 exit(0);
1619 // xmlAddChild(qm_root_p, createCell());
1620 }
1621 xmlAddChild(qm_root_p, createIonSet());
1622 xmlAddChild(qm_root_p, createElectronSet(ion_tag));
1623 xmlDocSetRootElement(doc_p, qm_root_p);
1624 std::string fname = Title + ".structure.xml";
1625 xmlSaveFormatFile(fname.c_str(), doc_p, 1);
1626 xmlFreeDoc(doc_p);
1627 Structure = true;
1628 }
1629 xmlDocPtr doc = xmlNewDoc((const xmlChar*)"1.0");
1630 xmlNodePtr qm_root = xmlNewNode(NULL, BAD_CAST "qmcsystem");
1631 {
1632 //wavefunction
1633 xmlNodePtr wfPtr = xmlNewNode(NULL, (const xmlChar*)"wavefunction");
1634 xmlNewProp(wfPtr, (const xmlChar*)"name", (const xmlChar*)psi_tag.c_str());
1635 xmlNewProp(wfPtr, (const xmlChar*)"target", (const xmlChar*)"e");
1636 {
1637 xmlNodePtr detPtr = xmlNewNode(NULL, (const xmlChar*)"determinantset");
1638 xmlNewProp(detPtr, (const xmlChar*)"type", (const xmlChar*)"MolecularOrbital");
1639 xmlNewProp(detPtr, (const xmlChar*)"name", (const xmlChar*)"LCAOBSet");
1640 xmlNewProp(detPtr, (const xmlChar*)"source", (const xmlChar*)ion_tag.c_str());
1641 xmlNewProp(detPtr, (const xmlChar*)"transform", (const xmlChar*)"yes");
1642
1643 if (DoCusp == true)
1644 xmlNewProp(detPtr, (const xmlChar*)"cuspCorrection", (const xmlChar*)"yes");
1645 if (UseHDF5 || singledetH5)
1646 xmlNewProp(detPtr, (const xmlChar*)"href", (const xmlChar*)h5file.c_str());
1647 //BASISSET
1648 {
1649 if (UseHDF5)
1650 xmlNodePtr bsetPtr = createBasisSetWithHDF5();
1651 else
1652 {
1653 if (!singledetH5)
1654 {
1655 xmlNodePtr bsetPtr = createBasisSet();
1656 xmlAddChild(detPtr, bsetPtr);
1657 }
1658 }
1659 if (multideterminant)
1660 {
1661 xmlNodePtr spoupPtr = xmlNewNode(NULL, (const xmlChar*)"sposet");
1662 xmlNodePtr spodnPtr = xmlNewNode(NULL, (const xmlChar*)"sposet");
1663 xmlNewProp(spoupPtr, (const xmlChar*)"basisset", (const xmlChar*)"LCAOBSet");
1664 xmlNewProp(spodnPtr, (const xmlChar*)"basisset", (const xmlChar*)"LCAOBSet");
1665 if (multidetH5)
1666 {
1667 PrepareSPOSetsFromH5(spoupPtr, spodnPtr);
1668 xmlAddChild(detPtr, spoupPtr);
1669 xmlAddChild(detPtr, spodnPtr);
1670 xmlNodePtr multislaterdetPtr = NULL;
1671 multislaterdetPtr = createMultiDeterminantSetFromH5();
1672 xmlAddChild(detPtr, multislaterdetPtr);
1673 }
1674 else
1675 {
1676 if (UseHDF5)
1677 {
1678 createSPOSetsH5(spoupPtr, spodnPtr);
1679 xmlAddChild(detPtr, spoupPtr);
1680 xmlAddChild(detPtr, spodnPtr);
1681 xmlNodePtr multislaterdetPtr = NULL;
1682 if (usingCSF)
1683 {
1684 app_log() << "Warning: CSF in HDF5 not implemented. Will attempt to revert multideterminant to xml. It "
1685 "is recommended to verify input for accuracy or avoid using -hdf5"
1686 << std::endl;
1687 multislaterdetPtr = createMultiDeterminantSet();
1688 }
1689 else
1690 multislaterdetPtr = createMultiDeterminantSetCIHDF5();
1691 xmlAddChild(detPtr, multislaterdetPtr);
1692 }
1693 else
1694 {
1695 createSPOSets(spoupPtr, spodnPtr);
1696 xmlAddChild(detPtr, spoupPtr);
1697 xmlAddChild(detPtr, spodnPtr);
1698 xmlNodePtr multislaterdetPtr = NULL;
1699 multislaterdetPtr = createMultiDeterminantSet();
1700 xmlAddChild(detPtr, multislaterdetPtr);
1701 }
1702 }
1703 }
1704 else
1705 {
1706 xmlNodePtr slaterdetPtr = NULL;
1707 if (UseHDF5)
1708 {
1709 slaterdetPtr = createDeterminantSetWithHDF5();
1710 }
1711 else
1712 {
1713 if (singledetH5)
1714 slaterdetPtr = PrepareDeterminantSetFromHDF5();
1715 else
1716 slaterdetPtr = createDeterminantSet();
1717 }
1718 xmlAddChild(detPtr, slaterdetPtr);
1719 }
1720 }
1721 xmlAddChild(wfPtr, detPtr);
1722 if (addJastrow)
1723 {
1724 std::cout << "Adding Two-Body and One-Body jastrows with rcut=\"10\" and size=\"10\"" << std::endl;
1725 if (NumberOfEls > 1)
1726 {
1727 xmlAddChild(wfPtr, createJ2());
1728 }
1729 xmlAddChild(wfPtr, createJ1());
1730 if (NumberOfEls > 1)
1731 {
1732 std::cout << "Adding Three-Body jastrows with rcut=\"5\"" << std::endl;
1733 xmlAddChild(wfPtr, createJ3());
1734 }
1735 }
1736 }
1737 xmlAddChild(qm_root, wfPtr);
1738 }
1739 xmlDocSetRootElement(doc, qm_root);
1740 std::string fname = Title + ".wf" + WFS_name + ".xml";
1741 xmlSaveFormatFile(fname.c_str(), doc, 1);
1742 xmlFreeDoc(doc);
1743 if (numMO * SizeOfBasisSet >= 4000 && !UseHDF5)
1744 if (!singledetH5)
1745 std::cout << "Consider using HDF5 via -hdf5 for higher performance and smaller wavefunction files" << std::endl;
1746 }
1747
dumpPBC(const std::string & psi_tag,const std::string & ion_tag)1748 void QMCGaussianParserBase::dumpPBC(const std::string& psi_tag, const std::string& ion_tag)
1749 {
1750 std::cout << " QMCGaussianParserBase::dumpPBC " << std::endl;
1751 if (!Structure)
1752 {
1753 xmlDocPtr doc_p = xmlNewDoc((const xmlChar*)"1.0");
1754 xmlNodePtr qm_root_p = xmlNewNode(NULL, BAD_CAST "qmcsystem");
1755 if (PBC)
1756 xmlAddChild(qm_root_p, createCell());
1757 xmlAddChild(qm_root_p, createIonSet());
1758 xmlAddChild(qm_root_p, createElectronSet(ion_tag));
1759 xmlDocSetRootElement(doc_p, qm_root_p);
1760 std::string fname = Title + ".structure.xml";
1761 xmlSaveFormatFile(fname.c_str(), doc_p, 1);
1762 xmlFreeDoc(doc_p);
1763 Structure = true;
1764 }
1765 xmlDocPtr doc = xmlNewDoc((const xmlChar*)"1.0");
1766 xmlNodePtr qm_root = xmlNewNode(NULL, BAD_CAST "qmcsystem");
1767 {
1768 //wavefunction
1769 xmlNodePtr wfPtr = xmlNewNode(NULL, (const xmlChar*)"wavefunction");
1770 xmlNewProp(wfPtr, (const xmlChar*)"name", (const xmlChar*)psi_tag.c_str());
1771 xmlNewProp(wfPtr, (const xmlChar*)"target", (const xmlChar*)"e");
1772 {
1773 xmlNodePtr detPtr = xmlNewNode(NULL, (const xmlChar*)"determinantset");
1774 xmlNewProp(detPtr, (const xmlChar*)"type", (const xmlChar*)"MolecularOrbital");
1775 xmlNewProp(detPtr, (const xmlChar*)"name", (const xmlChar*)"LCAOBSet");
1776 xmlNewProp(detPtr, (const xmlChar*)"source", (const xmlChar*)ion_tag.c_str());
1777 xmlNewProp(detPtr, (const xmlChar*)"transform", (const xmlChar*)"yes");
1778
1779 std::stringstream ss;
1780 ss << std::setprecision(10) << STwist_Coord[0] << " " << std::setprecision(10) << STwist_Coord[1] << " "
1781 << std::setprecision(10) << STwist_Coord[2];
1782 xmlNewProp(detPtr, (const xmlChar*)"twist", (const xmlChar*)(ss.str()).c_str());
1783
1784 if (DoCusp == true)
1785 xmlNewProp(detPtr, (const xmlChar*)"cuspCorrection", (const xmlChar*)"yes");
1786
1787 xmlNewProp(detPtr, (const xmlChar*)"href", (const xmlChar*)h5file.c_str());
1788
1789 std::stringstream sss;
1790 sss << Image[0] << " " << Image[1] << " " << Image[2];
1791 xmlNewProp(detPtr, (const xmlChar*)"PBCimages", (const xmlChar*)(sss.str()).c_str());
1792
1793 {
1794 if (multideterminant)
1795 {
1796 xmlNodePtr spoupPtr = xmlNewNode(NULL, (const xmlChar*)"sposet");
1797 xmlNodePtr spodnPtr = xmlNewNode(NULL, (const xmlChar*)"sposet");
1798 xmlNewProp(spoupPtr, (const xmlChar*)"basisset", (const xmlChar*)"LCAOBSet");
1799 xmlNewProp(spodnPtr, (const xmlChar*)"basisset", (const xmlChar*)"LCAOBSet");
1800
1801 PrepareSPOSetsFromH5(spoupPtr, spodnPtr);
1802 xmlAddChild(detPtr, spoupPtr);
1803 xmlAddChild(detPtr, spodnPtr);
1804 xmlNodePtr multislaterdetPtr = NULL;
1805
1806 multislaterdetPtr = createMultiDeterminantSetFromH5();
1807
1808
1809 xmlAddChild(detPtr, multislaterdetPtr);
1810 }
1811 else
1812 {
1813 xmlNodePtr slaterdetPtr = NULL;
1814 slaterdetPtr = PrepareDeterminantSetFromHDF5();
1815 xmlAddChild(detPtr, slaterdetPtr);
1816 }
1817 }
1818 xmlAddChild(wfPtr, detPtr);
1819 if (addJastrow)
1820 {
1821 std::cout << "Adding Two-Body and One-Body jastrows with rcut=\"10\" and size=\"10\"" << std::endl;
1822 if (NumberOfEls > 1)
1823 {
1824 xmlAddChild(wfPtr, createJ2());
1825 }
1826 xmlAddChild(wfPtr, createJ1());
1827 if (NumberOfEls > 1)
1828 {
1829 std::cout << "Adding Three-Body jastrows with rcut=\"5\"" << std::endl;
1830 xmlAddChild(wfPtr, createJ3());
1831 }
1832 }
1833 }
1834 xmlAddChild(qm_root, wfPtr);
1835 }
1836 xmlDocSetRootElement(doc, qm_root);
1837 std::string fname = Title + ".wf" + WFS_name + ".xml";
1838 xmlSaveFormatFile(fname.c_str(), doc, 1);
1839 xmlFreeDoc(doc);
1840 }
1841
1842
dumpStdInputProd(const std::string & psi_tag,const std::string & ion_tag)1843 void QMCGaussianParserBase::dumpStdInputProd(const std::string& psi_tag, const std::string& ion_tag)
1844 {
1845 std::cout << " Generating production input file designed for large calculations." << std::endl;
1846 std::cout << " Modify according to the accuracy you would like to achieve. " << std::endl;
1847
1848 std::string fname = Title + ".qmc.in-wf" + WFS_name + ".xml";
1849
1850 xmlDocPtr doc_input = xmlNewDoc((const xmlChar*)"1.0");
1851 xmlNodePtr qm_root_input = xmlNewNode(NULL, BAD_CAST "simulation");
1852
1853 ///Adding Project id
1854 {
1855 xmlNodePtr project = xmlNewNode(NULL, (const xmlChar*)"project");
1856 xmlNewProp(project, (const xmlChar*)"id", (const xmlChar*)Title.c_str());
1857 xmlNewProp(project, (const xmlChar*)"series", (const xmlChar*)"0");
1858 xmlAddChild(qm_root_input, project);
1859 }
1860
1861 ///Adding Link to Partcle Set and Wave function
1862 {
1863 std::string Ptclname = Title + ".structure.xml";
1864 xmlNodePtr ptcl = xmlNewNode(NULL, (const xmlChar*)"include");
1865 xmlNewProp(ptcl, (const xmlChar*)"href", (const xmlChar*)Ptclname.c_str());
1866 xmlAddChild(qm_root_input, ptcl);
1867
1868 std::string Wfsname = Title + ".wf" + WFS_name + ".xml";
1869 xmlNodePtr wfs = xmlNewNode(NULL, (const xmlChar*)"include");
1870 xmlNewProp(wfs, (const xmlChar*)"href", (const xmlChar*)Wfsname.c_str());
1871 xmlAddChild(qm_root_input, wfs);
1872 }
1873
1874 ///Adding Hamiltonian
1875 {
1876 xmlAddChild(qm_root_input, createHamiltonian(ion_tag, psi_tag));
1877 }
1878 ///Adding Optimization Block based on One Shift Only
1879 if (addJastrow)
1880 {
1881 ///Adding a VMC Block to help equilibrate
1882 xmlNodePtr initvmc = xmlNewNode(NULL, (const xmlChar*)"qmc");
1883 xmlNewProp(initvmc, (const xmlChar*)"method", (const xmlChar*)"vmc");
1884 xmlNewProp(initvmc, (const xmlChar*)"move", (const xmlChar*)"pbyp");
1885 xmlNewProp(initvmc, (const xmlChar*)"checkpoint", (const xmlChar*)"-1");
1886 //xmlNewProp(initvmc,(const xmlChar*)"gpu", (const xmlChar*)"no");
1887 {
1888 xmlNodePtr estimator = xmlNewNode(NULL, (const xmlChar*)"estimator");
1889 xmlNewProp(estimator, (const xmlChar*)"name", (const xmlChar*)"LocalEnergy");
1890 xmlNewProp(estimator, (const xmlChar*)"hdf5", (const xmlChar*)"no");
1891 xmlAddChild(initvmc, estimator);
1892
1893 xmlAddChild(initvmc, parameter(initvmc, "walkers", "1"));
1894 xmlAddChild(initvmc, parameter(initvmc, "samplesperthread", "1"));
1895 xmlAddChild(initvmc, parameter(initvmc, "stepsbetweensamples", "10"));
1896 xmlAddChild(initvmc, parameter(initvmc, "substeps", "5"));
1897 xmlAddChild(initvmc, parameter(initvmc, "warmupSteps", "20"));
1898 xmlAddChild(initvmc, parameter(initvmc, "blocks", "10"));
1899 xmlAddChild(initvmc, parameter(initvmc, "timestep", "0.5"));
1900 xmlAddChild(initvmc, parameter(initvmc, "usedrift", "no"));
1901 }
1902 xmlAddChild(qm_root_input, initvmc);
1903
1904
1905 ///Adding First loop of Cheap optimization blocks
1906 xmlNodePtr loopopt1 = xmlNewNode(NULL, (const xmlChar*)"loop");
1907 xmlNewProp(loopopt1, (const xmlChar*)"max", (const xmlChar*)"4");
1908 {
1909 xmlNodePtr initopt = xmlNewNode(NULL, (const xmlChar*)"qmc");
1910 xmlNewProp(initopt, (const xmlChar*)"method", (const xmlChar*)"linear");
1911 xmlNewProp(initopt, (const xmlChar*)"move", (const xmlChar*)"pbyp");
1912 xmlNewProp(initopt, (const xmlChar*)"checkpoint", (const xmlChar*)"-1");
1913 //xmlNewProp(initopt,(const xmlChar*)"gpu", (const xmlChar*)"no");
1914 {
1915 xmlNodePtr estimator = xmlNewNode(NULL, (const xmlChar*)"estimator");
1916 xmlNewProp(estimator, (const xmlChar*)"name", (const xmlChar*)"LocalEnergy");
1917 xmlNewProp(estimator, (const xmlChar*)"hdf5", (const xmlChar*)"no");
1918 xmlAddChild(initopt, estimator);
1919
1920 xmlAddChild(initopt, parameter(initopt, "blocks", "20"));
1921 xmlAddChild(initopt, parameter(initopt, "warmupSteps", "2"));
1922 xmlAddChild(initopt, parameter(initopt, "timestep", "0.5"));
1923 xmlAddChild(initopt, parameter(initopt, "walkers", "1"));
1924 xmlAddChild(initopt, parameter(initopt, "samples", "80000"));
1925 xmlAddChild(initopt, parameter(initopt, "substeps", "5"));
1926 xmlAddChild(initopt, parameter(initopt, "usedrift", "no"));
1927 xmlAddChild(initopt, parameter(initopt, "MinMethod", "OneShiftOnly"));
1928 xmlAddChild(initopt, parameter(initopt, "minwalkers", "0.1"));
1929 }
1930 xmlAddChild(loopopt1, initopt);
1931 }
1932 xmlAddChild(qm_root_input, loopopt1);
1933
1934 ///Adding loop for optimization blocks
1935 xmlNodePtr loopopt = xmlNewNode(NULL, (const xmlChar*)"loop");
1936 xmlNewProp(loopopt, (const xmlChar*)"max", (const xmlChar*)"10");
1937 {
1938 xmlNodePtr initopt = xmlNewNode(NULL, (const xmlChar*)"qmc");
1939 xmlNewProp(initopt, (const xmlChar*)"method", (const xmlChar*)"linear");
1940 xmlNewProp(initopt, (const xmlChar*)"move", (const xmlChar*)"pbyp");
1941 xmlNewProp(initopt, (const xmlChar*)"checkpoint", (const xmlChar*)"-1");
1942 //xmlNewProp(initopt,(const xmlChar*)"gpu", (const xmlChar*)"no");
1943 {
1944 xmlNodePtr estimator = xmlNewNode(NULL, (const xmlChar*)"estimator");
1945 xmlNewProp(estimator, (const xmlChar*)"name", (const xmlChar*)"LocalEnergy");
1946 xmlNewProp(estimator, (const xmlChar*)"hdf5", (const xmlChar*)"no");
1947 xmlAddChild(initopt, estimator);
1948
1949 xmlAddChild(initopt, parameter(initopt, "blocks", "40"));
1950 xmlAddChild(initopt, parameter(initopt, "warmupSteps", "5"));
1951 xmlAddChild(initopt, parameter(initopt, "timestep", "0.5"));
1952 xmlAddChild(initopt, parameter(initopt, "walkers", "1"));
1953 xmlAddChild(initopt, parameter(initopt, "samples", "160000"));
1954 xmlAddChild(initopt, parameter(initopt, "substeps", "5"));
1955 xmlAddChild(initopt, parameter(initopt, "usedrift", "no"));
1956 xmlAddChild(initopt, parameter(initopt, "MinMethod", "OneShiftOnly"));
1957 xmlAddChild(initopt, parameter(initopt, "minwalkers", "0.5"));
1958 }
1959 xmlAddChild(loopopt, initopt);
1960 }
1961 xmlAddChild(qm_root_input, loopopt);
1962 }
1963
1964
1965 ///Adding a VMC Block to the Input
1966 xmlNodePtr vmc = xmlNewNode(NULL, (const xmlChar*)"qmc");
1967 xmlNewProp(vmc, (const xmlChar*)"method", (const xmlChar*)"vmc");
1968 xmlNewProp(vmc, (const xmlChar*)"move", (const xmlChar*)"pbyp");
1969 xmlNewProp(vmc, (const xmlChar*)"checkpoint", (const xmlChar*)"-1");
1970 //xmlNewProp(vmc,(const xmlChar*)"gpu", (const xmlChar*)"no");
1971 {
1972 xmlNodePtr estimator = xmlNewNode(NULL, (const xmlChar*)"estimator");
1973 xmlNewProp(estimator, (const xmlChar*)"name", (const xmlChar*)"LocalEnergy");
1974 xmlNewProp(estimator, (const xmlChar*)"hdf5", (const xmlChar*)"no");
1975 xmlAddChild(vmc, estimator);
1976
1977 xmlAddChild(vmc, parameter(vmc, "walkers", "1"));
1978 xmlAddChild(vmc, parameter(vmc, "samplesperthread", "1"));
1979 xmlAddChild(vmc, parameter(vmc, "stepsbetweensamples", "10"));
1980 xmlAddChild(vmc, parameter(vmc, "substeps", "30"));
1981 xmlAddChild(vmc, parameter(vmc, "warmupSteps", "100"));
1982 xmlAddChild(vmc, parameter(vmc, "blocks", "200"));
1983 xmlAddChild(vmc, parameter(vmc, "timestep", "0.1"));
1984 xmlAddChild(vmc, parameter(vmc, "usedrift", "no"));
1985 }
1986 xmlAddChild(qm_root_input, vmc);
1987
1988 ///Adding a DMC Block to the Input
1989 xmlNodePtr dmc = xmlNewNode(NULL, (const xmlChar*)"qmc");
1990 xmlNewProp(dmc, (const xmlChar*)"method", (const xmlChar*)"dmc");
1991 xmlNewProp(dmc, (const xmlChar*)"move", (const xmlChar*)"pbyp");
1992 xmlNewProp(dmc, (const xmlChar*)"checkpoint", (const xmlChar*)"20");
1993 //xmlNewProp(dmc,(const xmlChar*)"gpu", (const xmlChar*)"no");
1994 {
1995 xmlNodePtr estimator = xmlNewNode(NULL, (const xmlChar*)"estimator");
1996 xmlNewProp(estimator, (const xmlChar*)"name", (const xmlChar*)"LocalEnergy");
1997 xmlNewProp(estimator, (const xmlChar*)"hdf5", (const xmlChar*)"no");
1998 xmlAddChild(dmc, estimator);
1999
2000 xmlAddChild(dmc, parameter(dmc, "targetwalkers", "16000"));
2001 xmlAddChild(dmc, parameter(dmc, "reconfiguration", "no"));
2002 xmlAddChild(dmc, parameter(dmc, "warmupSteps", "100"));
2003 xmlAddChild(dmc, parameter(dmc, "timestep", "0.0005"));
2004 xmlAddChild(dmc, parameter(dmc, "steps", "30"));
2005 xmlAddChild(dmc, parameter(dmc, "blocks", "1000"));
2006 xmlAddChild(dmc, parameter(dmc, "nonlocalmoves", "v3"));
2007 }
2008 xmlAddChild(qm_root_input, dmc);
2009
2010 xmlDocSetRootElement(doc_input, qm_root_input);
2011 xmlSaveFormatFile(fname.c_str(), doc_input, 1);
2012 xmlFreeDoc(doc_input);
2013 }
2014
dumpStdInput(const std::string & psi_tag,const std::string & ion_tag)2015 void QMCGaussianParserBase::dumpStdInput(const std::string& psi_tag, const std::string& ion_tag)
2016 {
2017 std::cout << " Generating Standard Input file containing VMC, standard optmization, and DMC blocks." << std::endl;
2018 std::cout << " Modify according to the accuracy you would like to achieve. " << std::endl;
2019
2020 std::string fname = Title + ".qmc.in-wf" + WFS_name + ".xml";
2021
2022 xmlDocPtr doc_input = xmlNewDoc((const xmlChar*)"1.0");
2023 xmlNodePtr qm_root_input = xmlNewNode(NULL, BAD_CAST "simulation");
2024
2025 std::ostringstream Comment;
2026 ///Adding Project id
2027 {
2028 Comment.str("");
2029 Comment.clear();
2030 Comment << "\n \nExample QMCPACK input file produced by convert4qmc\n \nIt is recommend to start with only the "
2031 "initial VMC block and adjust\nparameters based on the measured energies, variance, and statistics.\n\n";
2032 xmlNodePtr MyComment = xmlNewComment((const xmlChar*)Comment.str().c_str());
2033 xmlAddChild(qm_root_input, MyComment);
2034 Comment.str("");
2035 Comment.clear();
2036 Comment << "Name and Series number of the project.";
2037 MyComment = xmlNewComment((const xmlChar*)Comment.str().c_str());
2038 xmlAddChild(qm_root_input, MyComment);
2039
2040 xmlNodePtr project = xmlNewNode(NULL, (const xmlChar*)"project");
2041 xmlNewProp(project, (const xmlChar*)"id", (const xmlChar*)Title.c_str());
2042 xmlNewProp(project, (const xmlChar*)"series", (const xmlChar*)"0");
2043 xmlAddChild(qm_root_input, project);
2044 }
2045
2046 ///Adding Link to Partcle Set and Wave function
2047 {
2048 Comment.str("");
2049 Comment.clear();
2050 Comment << "Link to the location of the Atomic Coordinates and the location of the Wavefunction.";
2051 xmlNodePtr MyComment = xmlNewComment((const xmlChar*)Comment.str().c_str());
2052 xmlAddChild(qm_root_input, MyComment);
2053 std::string Ptclname = Title + ".structure.xml";
2054 xmlNodePtr ptcl = xmlNewNode(NULL, (const xmlChar*)"include");
2055 xmlNewProp(ptcl, (const xmlChar*)"href", (const xmlChar*)Ptclname.c_str());
2056 xmlAddChild(qm_root_input, ptcl);
2057
2058 std::string Wfsname = Title + ".wf" + WFS_name + ".xml";
2059 xmlNodePtr wfs = xmlNewNode(NULL, (const xmlChar*)"include");
2060 xmlNewProp(wfs, (const xmlChar*)"href", (const xmlChar*)Wfsname.c_str());
2061 xmlAddChild(qm_root_input, wfs);
2062 }
2063
2064 ///Adding Hamiltonian
2065 {
2066 Comment.str("");
2067 Comment.clear();
2068 if (ECP == true)
2069 Comment << "Hamiltonian of the system. Default ECP filenames are assumed.";
2070 else
2071 Comment << "Hamiltonian of the system.\n";
2072
2073 xmlNodePtr MyComment = xmlNewComment((const xmlChar*)Comment.str().c_str());
2074 xmlAddChild(qm_root_input, MyComment);
2075 xmlAddChild(qm_root_input, createHamiltonian(ion_tag, psi_tag));
2076 }
2077
2078 {
2079 Comment.str("");
2080 Comment.clear();
2081 Comment << "\n \nExample initial VMC to measure initial energy and variance \n\n";
2082 xmlNodePtr MyComment = xmlNewComment((const xmlChar*)Comment.str().c_str());
2083 xmlAddChild(qm_root_input, MyComment);
2084 }
2085 xmlNodePtr initvmc = xmlNewNode(NULL, (const xmlChar*)"qmc");
2086 xmlNewProp(initvmc, (const xmlChar*)"method", (const xmlChar*)"vmc");
2087 xmlNewProp(initvmc, (const xmlChar*)"move", (const xmlChar*)"pbyp");
2088 xmlNewProp(initvmc, (const xmlChar*)"checkpoint", (const xmlChar*)"-1");
2089 //xmlNewProp(initvmc,(const xmlChar*)"gpu", (const xmlChar*)"no");
2090 {
2091 xmlNodePtr estimator = xmlNewNode(NULL, (const xmlChar*)"estimator");
2092 xmlNewProp(estimator, (const xmlChar*)"name", (const xmlChar*)"LocalEnergy");
2093 xmlNewProp(estimator, (const xmlChar*)"hdf5", (const xmlChar*)"no");
2094 xmlAddChild(initvmc, estimator);
2095
2096 xmlAddChild(initvmc, parameter(initvmc, "warmupSteps", "100"));
2097 xmlAddChild(initvmc, parameter(initvmc, "blocks", "20"));
2098 xmlAddChild(initvmc, parameter(initvmc, "steps", "50"));
2099 xmlAddChild(initvmc, parameter(initvmc, "substeps", "8"));
2100 xmlAddChild(initvmc, parameter(initvmc, "timestep", "0.5"));
2101 xmlAddChild(initvmc, parameter(initvmc, "usedrift", "no"));
2102 }
2103 xmlAddChild(qm_root_input, initvmc);
2104
2105
2106 if (addJastrow)
2107 {
2108 ///Adding First loop of Cheap optimization blocks
2109 {
2110 Comment.str("");
2111 Comment.clear();
2112 Comment << "\n \nExample initial VMC optimization \n \nNumber of steps required will be computed from total "
2113 "requested sample \ncount and total number of walkers \n\n";
2114 xmlNodePtr MyComment = xmlNewComment((const xmlChar*)Comment.str().c_str());
2115 xmlAddChild(qm_root_input, MyComment);
2116 }
2117 xmlNodePtr loopopt1 = xmlNewNode(NULL, (const xmlChar*)"loop");
2118 xmlNewProp(loopopt1, (const xmlChar*)"max", (const xmlChar*)"4");
2119 {
2120 xmlNodePtr initopt = xmlNewNode(NULL, (const xmlChar*)"qmc");
2121 xmlNewProp(initopt, (const xmlChar*)"method", (const xmlChar*)"linear");
2122 xmlNewProp(initopt, (const xmlChar*)"move", (const xmlChar*)"pbyp");
2123 xmlNewProp(initopt, (const xmlChar*)"checkpoint", (const xmlChar*)"-1");
2124 {
2125 xmlNodePtr estimator = xmlNewNode(NULL, (const xmlChar*)"estimator");
2126 xmlNewProp(estimator, (const xmlChar*)"name", (const xmlChar*)"LocalEnergy");
2127 xmlNewProp(estimator, (const xmlChar*)"hdf5", (const xmlChar*)"no");
2128 xmlAddChild(initopt, estimator);
2129
2130 xmlAddChild(initopt, parameter(initopt, "warmupSteps", "100"));
2131 xmlAddChild(initopt, parameter(initopt, "blocks", "20"));
2132 xmlAddChild(initopt, parameter(initopt, "timestep", "0.5"));
2133 xmlAddChild(initopt, parameter(initopt, "walkers", "1"));
2134 xmlAddChild(initopt, parameter(initopt, "samples", "16000"));
2135 xmlAddChild(initopt, parameter(initopt, "substeps", "4"));
2136 xmlAddChild(initopt, parameter(initopt, "usedrift", "no"));
2137 xmlAddChild(initopt, parameter(initopt, "MinMethod", "OneShiftOnly"));
2138 xmlAddChild(initopt, parameter(initopt, "minwalkers", "0.0001"));
2139 }
2140 xmlAddChild(loopopt1, initopt);
2141 }
2142 xmlAddChild(qm_root_input, loopopt1);
2143
2144 ///Adding loop for optimization blocks
2145 {
2146 Comment.str("");
2147 Comment.clear();
2148 Comment << "\n \nExample follow-up VMC optimization using more samples for greater accuracy\n\n";
2149 xmlNodePtr MyComment = xmlNewComment((const xmlChar*)Comment.str().c_str());
2150 xmlAddChild(qm_root_input, MyComment);
2151 }
2152 xmlNodePtr loopopt = xmlNewNode(NULL, (const xmlChar*)"loop");
2153 xmlNewProp(loopopt, (const xmlChar*)"max", (const xmlChar*)"10");
2154 {
2155 xmlNodePtr initopt = xmlNewNode(NULL, (const xmlChar*)"qmc");
2156 xmlNewProp(initopt, (const xmlChar*)"method", (const xmlChar*)"linear");
2157 xmlNewProp(initopt, (const xmlChar*)"move", (const xmlChar*)"pbyp");
2158 xmlNewProp(initopt, (const xmlChar*)"checkpoint", (const xmlChar*)"-1");
2159 //xmlNewProp(initopt,(const xmlChar*)"gpu", (const xmlChar*)"no");
2160 {
2161 xmlNodePtr estimator = xmlNewNode(NULL, (const xmlChar*)"estimator");
2162 xmlNewProp(estimator, (const xmlChar*)"name", (const xmlChar*)"LocalEnergy");
2163 xmlNewProp(estimator, (const xmlChar*)"hdf5", (const xmlChar*)"no");
2164 xmlAddChild(initopt, estimator);
2165
2166 xmlAddChild(initopt, parameter(initopt, "warmupSteps", "100"));
2167 xmlAddChild(initopt, parameter(initopt, "blocks", "20"));
2168 xmlAddChild(initopt, parameter(initopt, "timestep", "0.5"));
2169 xmlAddChild(initopt, parameter(initopt, "walkers", "1"));
2170 xmlAddChild(initopt, parameter(initopt, "samples", "64000"));
2171 xmlAddChild(initopt, parameter(initopt, "substeps", "4"));
2172 xmlAddChild(initopt, parameter(initopt, "usedrift", "no"));
2173 xmlAddChild(initopt, parameter(initopt, "MinMethod", "OneShiftOnly"));
2174 xmlAddChild(initopt, parameter(initopt, "minwalkers", "0.3"));
2175 }
2176 xmlAddChild(loopopt, initopt);
2177 }
2178 xmlAddChild(qm_root_input, loopopt);
2179
2180
2181 ///Adding a VMC Block to the Input
2182 {
2183 Comment.str("");
2184 Comment.clear();
2185 Comment
2186 << "\n\nProduction VMC and DMC\n\nExamine the results of the optimization before running these blocks.\ne.g. "
2187 "Choose the best optimized jastrow from all obtained, put in \nwavefunction file, do not reoptimize.\n\n";
2188 xmlNodePtr MyComment = xmlNewComment((const xmlChar*)Comment.str().c_str());
2189 xmlAddChild(qm_root_input, MyComment);
2190 }
2191 xmlNodePtr vmc = xmlNewNode(NULL, (const xmlChar*)"qmc");
2192 xmlNewProp(vmc, (const xmlChar*)"method", (const xmlChar*)"vmc");
2193 xmlNewProp(vmc, (const xmlChar*)"move", (const xmlChar*)"pbyp");
2194 xmlNewProp(vmc, (const xmlChar*)"checkpoint", (const xmlChar*)"-1");
2195 //xmlNewProp(vmc,(const xmlChar*)"gpu", (const xmlChar*)"no");
2196 {
2197 std::ostringstream Comment;
2198 xmlNodePtr estimator = xmlNewNode(NULL, (const xmlChar*)"estimator");
2199 xmlNewProp(estimator, (const xmlChar*)"name", (const xmlChar*)"LocalEnergy");
2200 xmlNewProp(estimator, (const xmlChar*)"hdf5", (const xmlChar*)"no");
2201 xmlAddChild(vmc, estimator);
2202
2203 xmlAddChild(vmc, parameter(vmc, "warmupSteps", "100"));
2204 xmlAddChild(vmc, parameter(vmc, "blocks", "200"));
2205 xmlAddChild(vmc, parameter(vmc, "steps", "50"));
2206 xmlAddChild(vmc, parameter(vmc, "substeps", "8"));
2207 xmlAddChild(vmc, parameter(vmc, "timestep", "0.5"));
2208 xmlAddChild(vmc, parameter(vmc, "usedrift", "no"));
2209 Comment.str("");
2210 Comment.clear();
2211 Comment << "Sample count should match targetwalker count for DMC. Will be obtained from all nodes.";
2212 xmlNodePtr MyComment = xmlNewComment((const xmlChar*)Comment.str().c_str());
2213 xmlAddChild(vmc, MyComment);
2214 xmlAddChild(vmc, parameter(vmc, "samples", "16000"));
2215 }
2216 xmlAddChild(qm_root_input, vmc);
2217
2218 ///Adding a DMC Block to the Input
2219 xmlNodePtr dmc = xmlNewNode(NULL, (const xmlChar*)"qmc");
2220 xmlNewProp(dmc, (const xmlChar*)"method", (const xmlChar*)"dmc");
2221 xmlNewProp(dmc, (const xmlChar*)"move", (const xmlChar*)"pbyp");
2222 xmlNewProp(dmc, (const xmlChar*)"checkpoint", (const xmlChar*)"20");
2223 //xmlNewProp(dmc,(const xmlChar*)"gpu", (const xmlChar*)"no");
2224 {
2225 xmlNodePtr estimator = xmlNewNode(NULL, (const xmlChar*)"estimator");
2226 xmlNewProp(estimator, (const xmlChar*)"name", (const xmlChar*)"LocalEnergy");
2227 xmlNewProp(estimator, (const xmlChar*)"hdf5", (const xmlChar*)"no");
2228 xmlAddChild(dmc, estimator);
2229
2230 xmlAddChild(dmc, parameter(dmc, "targetwalkers", "16000"));
2231 xmlAddChild(dmc, parameter(dmc, "reconfiguration", "no"));
2232 xmlAddChild(dmc, parameter(dmc, "warmupSteps", "100"));
2233 xmlAddChild(dmc, parameter(dmc, "timestep", "0.0005"));
2234 xmlAddChild(dmc, parameter(dmc, "steps", "30"));
2235 xmlAddChild(dmc, parameter(dmc, "blocks", "1000"));
2236 xmlAddChild(dmc, parameter(dmc, "nonlocalmoves", "v3"));
2237 }
2238 xmlAddChild(qm_root_input, dmc);
2239 }
2240
2241 xmlDocSetRootElement(doc_input, qm_root_input);
2242 xmlSaveFormatFile(fname.c_str(), doc_input, 1);
2243 xmlFreeDoc(doc_input);
2244 }
2245
createHamiltonian(const std::string & ion_tag,const std::string & psi_tag)2246 xmlNodePtr QMCGaussianParserBase::createHamiltonian(const std::string& ion_tag, const std::string& psi_tag)
2247 {
2248 xmlNodePtr hamPtr = xmlNewNode(NULL, (const xmlChar*)"hamiltonian");
2249 xmlNewProp(hamPtr, (const xmlChar*)"name", (const xmlChar*)"h0");
2250 xmlNewProp(hamPtr, (const xmlChar*)"type", (const xmlChar*)"generic");
2251 xmlNewProp(hamPtr, (const xmlChar*)"target", (const xmlChar*)"e");
2252
2253 {
2254 xmlNodePtr pairpot1 = xmlNewNode(NULL, (const xmlChar*)"pairpot");
2255 xmlNewProp(pairpot1, (const xmlChar*)"name", (const xmlChar*)"ElecElec");
2256 xmlNewProp(pairpot1, (const xmlChar*)"type", (const xmlChar*)"coulomb");
2257 xmlNewProp(pairpot1, (const xmlChar*)"source", (const xmlChar*)"e");
2258 xmlNewProp(pairpot1, (const xmlChar*)"target", (const xmlChar*)"e");
2259 xmlNewProp(pairpot1, (const xmlChar*)"physical", (const xmlChar*)"true");
2260 xmlAddChild(hamPtr, pairpot1);
2261
2262 xmlNodePtr pairpot2 = xmlNewNode(NULL, (const xmlChar*)"pairpot");
2263 xmlNewProp(pairpot2, (const xmlChar*)"name", (const xmlChar*)"IonIon");
2264 xmlNewProp(pairpot2, (const xmlChar*)"type", (const xmlChar*)"coulomb");
2265 xmlNewProp(pairpot2, (const xmlChar*)"source", (const xmlChar*)ion_tag.c_str());
2266 xmlNewProp(pairpot2, (const xmlChar*)"target", (const xmlChar*)ion_tag.c_str());
2267 xmlAddChild(hamPtr, pairpot2);
2268 if (ECP == false)
2269 {
2270 xmlNodePtr pairpot3 = xmlNewNode(NULL, (const xmlChar*)"pairpot");
2271 xmlNewProp(pairpot3, (const xmlChar*)"name", (const xmlChar*)"IonElec");
2272 xmlNewProp(pairpot3, (const xmlChar*)"type", (const xmlChar*)"coulomb");
2273 xmlNewProp(pairpot3, (const xmlChar*)"source", (const xmlChar*)ion_tag.c_str());
2274 xmlNewProp(pairpot3, (const xmlChar*)"target", (const xmlChar*)"e");
2275 xmlAddChild(hamPtr, pairpot3);
2276 }
2277 else
2278 {
2279 std::cout << "Hamiltonian using ECP for Electron Ion=" << ECP << std::endl;
2280 xmlNodePtr pairpot3 = xmlNewNode(NULL, (const xmlChar*)"pairpot");
2281 xmlNewProp(pairpot3, (const xmlChar*)"name", (const xmlChar*)"PseudoPot");
2282 xmlNewProp(pairpot3, (const xmlChar*)"type", (const xmlChar*)"pseudo");
2283 xmlNewProp(pairpot3, (const xmlChar*)"source", (const xmlChar*)ion_tag.c_str());
2284 xmlNewProp(pairpot3, (const xmlChar*)"wavefunction", (const xmlChar*)psi_tag.c_str());
2285 xmlNewProp(pairpot3, (const xmlChar*)"format", (const xmlChar*)"xml");
2286 {
2287 std::vector<std::string> AtomNames(GroupName);
2288 sort(AtomNames.begin(), AtomNames.end());
2289 AtomNames.erase(unique(AtomNames.begin(), AtomNames.end()), AtomNames.end());
2290
2291 for (int iat = 0; iat < AtomNames.size(); iat++)
2292 {
2293 std::string PPname = AtomNames[iat] + ".qmcpp.xml";
2294 xmlNodePtr a = xmlNewNode(NULL, (const xmlChar*)"pseudo");
2295 xmlNewProp(a, (const xmlChar*)"elementType", (const xmlChar*)AtomNames[iat].c_str());
2296 xmlNewProp(a, (const xmlChar*)"href", (const xmlChar*)PPname.c_str());
2297 xmlAddChild(pairpot3, a);
2298 }
2299 }
2300 xmlAddChild(hamPtr, pairpot3);
2301 }
2302 }
2303 return hamPtr;
2304 }
2305
numberOfExcitationsCSF(std::string & occ)2306 int QMCGaussianParserBase::numberOfExcitationsCSF(std::string& occ)
2307 {
2308 int res = 0;
2309 for (int i = ci_neb; i < ci_nstates; i++)
2310 {
2311 if (i < ci_nea && occ[i] == '2')
2312 {
2313 res++; //excitation into singly occupied alpha states in the reference
2314 }
2315 else
2316 {
2317 if (occ[i] == '1')
2318 res++;
2319 else if (occ[i] == '2')
2320 res += 2;
2321 }
2322 }
2323 return res;
2324 }
2325
2326
parameter(xmlNodePtr Parent,std::string Mypara,std::string a)2327 xmlNodePtr QMCGaussianParserBase::parameter(xmlNodePtr Parent, std::string Mypara, std::string a)
2328 {
2329 xmlNodePtr e = xmlNewTextChild(Parent, NULL, (const xmlChar*)"parameter", (const xmlChar*)a.c_str());
2330 xmlNewProp(e, (const xmlChar*)"name", (const xmlChar*)Mypara.c_str());
2331 return e;
2332 }
2333
2334
PrepareSPOSetsFromH5(xmlNodePtr spoUP,xmlNodePtr spoDN)2335 void QMCGaussianParserBase::PrepareSPOSetsFromH5(xmlNodePtr spoUP, xmlNodePtr spoDN)
2336 {
2337 setOccupationNumbers();
2338 std::ostringstream up_size, down_size, b_size, occ, nstates_alpha, nstates_beta;
2339 up_size << NumberOfAlpha;
2340 down_size << NumberOfBeta;
2341 b_size << numMO;
2342 nstates_alpha << ci_nstates + ci_nca;
2343 ;
2344 nstates_beta << ci_nstates + ci_ncb;
2345
2346 //create a spoUp
2347 xmlNewProp(spoUP, (const xmlChar*)"name", (const xmlChar*)"spo-up");
2348 xmlNewProp(spoUP, (const xmlChar*)"size", (const xmlChar*)nstates_alpha.str().c_str());
2349
2350 //create a spoDN
2351 xmlNewProp(spoDN, (const xmlChar*)"name", (const xmlChar*)"spo-dn");
2352 xmlNewProp(spoDN, (const xmlChar*)"size", (const xmlChar*)nstates_beta.str().c_str());
2353
2354
2355 if (DoCusp == true)
2356 {
2357 xmlNewProp(spoUP, (const xmlChar*)"cuspInfo", (const xmlChar*)"../spo-up.cuspInfo.xml");
2358 xmlNewProp(spoDN, (const xmlChar*)"cuspInfo", (const xmlChar*)"../spo-dn.cuspInfo.xml");
2359 }
2360
2361
2362 //add occupation UP
2363 xmlNodePtr occ_data = xmlNewNode(NULL, (const xmlChar*)"occupation");
2364 xmlNewProp(occ_data, (const xmlChar*)"mode", (const xmlChar*)"ground");
2365 xmlAddChild(spoUP, occ_data);
2366
2367
2368 //add coefficients
2369 xmlNodePtr coeff_data = xmlNewNode(NULL, (const xmlChar*)"coefficient");
2370 xmlNewProp(coeff_data, (const xmlChar*)"size", (const xmlChar*)b_size.str().c_str());
2371 xmlNewProp(coeff_data, (const xmlChar*)"spindataset", (const xmlChar*)"0");
2372 xmlAddChild(spoUP, coeff_data);
2373
2374
2375 //add occupation DN
2376 occ_data = xmlNewNode(NULL, (const xmlChar*)"occupation");
2377 xmlNewProp(occ_data, (const xmlChar*)"mode", (const xmlChar*)"ground");
2378 xmlAddChild(spoDN, occ_data);
2379
2380 coeff_data = xmlNewNode(NULL, (const xmlChar*)"coefficient");
2381 xmlNewProp(coeff_data, (const xmlChar*)"size", (const xmlChar*)b_size.str().c_str());
2382 if (SpinRestricted)
2383 xmlNewProp(coeff_data, (const xmlChar*)"spindataset", (const xmlChar*)"0");
2384 else
2385 xmlNewProp(coeff_data, (const xmlChar*)"spindataset", (const xmlChar*)"1");
2386 xmlAddChild(spoDN, coeff_data);
2387 }
createMultiDeterminantSetFromH5()2388 xmlNodePtr QMCGaussianParserBase::createMultiDeterminantSetFromH5()
2389 {
2390 xmlNodePtr multislaterdet = xmlNewNode(NULL, (const xmlChar*)"multideterminant");
2391 if (optDetCoeffs)
2392 xmlNewProp(multislaterdet, (const xmlChar*)"optimize", (const xmlChar*)"yes");
2393 else
2394 xmlNewProp(multislaterdet, (const xmlChar*)"optimize", (const xmlChar*)"no");
2395 xmlNewProp(multislaterdet, (const xmlChar*)"spo_up", (const xmlChar*)"spo-up");
2396 xmlNewProp(multislaterdet, (const xmlChar*)"spo_dn", (const xmlChar*)"spo-dn");
2397 xmlNodePtr detlist = xmlNewNode(NULL, (const xmlChar*)"detlist");
2398 std::ostringstream nstates, cisize, cinca, cincb, cinea, cineb, ci_thr;
2399 cisize << ci_size;
2400 nstates << ci_nstates;
2401 cinca << ci_nca;
2402 cincb << ci_ncb;
2403 cinea << ci_nea;
2404 cineb << ci_neb;
2405 ci_thr << ci_threshold;
2406 xmlNewProp(detlist, (const xmlChar*)"size", (const xmlChar*)cisize.str().c_str());
2407 xmlNewProp(detlist, (const xmlChar*)"type", (const xmlChar*)"DETS");
2408 xmlNewProp(detlist, (const xmlChar*)"nca", (const xmlChar*)cinca.str().c_str());
2409 xmlNewProp(detlist, (const xmlChar*)"ncb", (const xmlChar*)cincb.str().c_str());
2410 xmlNewProp(detlist, (const xmlChar*)"nea", (const xmlChar*)cinea.str().c_str());
2411 xmlNewProp(detlist, (const xmlChar*)"neb", (const xmlChar*)cineb.str().c_str());
2412 xmlNewProp(detlist, (const xmlChar*)"nstates", (const xmlChar*)nstates.str().c_str());
2413 xmlNewProp(detlist, (const xmlChar*)"cutoff", (const xmlChar*)ci_thr.str().c_str());
2414 if (nbexcitedstates >= 1)
2415 {
2416 app_log() << "WARNING!! THE HDF5 Contains CI coefficients for " << nbexcitedstates
2417 << ". By default, the ground state coefficients will be loaded ( ext_level=0). If you want to evaluate "
2418 "an excited for which the coefficients are stored in the HDF5 file, modify the value of ext_level "
2419 "Using [1," << nbexcitedstates<< "]" <<std::endl;
2420 xmlNewProp(detlist, (const xmlChar*)"ext_level", (const xmlChar*)"0");
2421 }
2422 if (!debug)
2423 xmlNewProp(detlist, (const xmlChar*)"href", (const xmlChar*)multih5file.c_str());
2424 if (CIcoeff.size() == 0)
2425 {
2426 std::cerr << " CI configuration list is empty. \n";
2427 exit(101);
2428 }
2429 if (CIcoeff.size() != CIalpha.size() || CIcoeff.size() != CIbeta.size())
2430 {
2431 std::cerr << " Problem with CI configuration lists. \n";
2432 exit(102);
2433 }
2434 int iv = 0;
2435 xmlAddChild(multislaterdet, detlist);
2436 return multislaterdet;
2437 }
2438