1 #ifndef AFQMC_READWFN_HPP
2 #define AFQMC_READWFN_HPP
3
4 #include <cstdlib>
5 #include <iostream>
6 #include <fstream>
7 #include <vector>
8 #include <string>
9 #include <ctype.h>
10
11 #include "Utilities/SimpleParser.h"
12 #include "readWfn.h"
13
14 #include "AFQMC/config.h"
15 #include "AFQMC/Matrix/csr_matrix.hpp"
16 #include "AFQMC/Matrix/csr_hdf5_readers.hpp"
17 #include "AFQMC/Matrix/csr_matrix_construct.hpp"
18
19 namespace qmcplusplus
20 {
21 namespace
22 {
read_header(std::ifstream & in,std::string & type,int & wfn_type,bool & fullMOMat,bool & Cstyle,int & ndet)23 void read_header(std::ifstream& in, std::string& type, int& wfn_type, bool& fullMOMat, bool& Cstyle, int& ndet)
24 {
25 std::vector<std::string> words;
26 getwords(words, in);
27 do
28 {
29 if (words.size() == 0)
30 app_error() << "Format error in ASCII integral file. End of file in header. \n";
31 for (std::vector<std::string>::iterator it = words.begin(); it != words.end(); it++)
32 {
33 if (*it == "&FCI")
34 {
35 // do nothing
36 }
37 else if (*it == "Type" || *it == "TYPE" || *it == "type")
38 {
39 if (it + 1 == words.end())
40 {
41 app_error() << "Format error in ASCII integral file. NAEB \n";
42 APP_ABORT("Format error in ASCII integral file. NAEB \n");
43 }
44 type = *(it + 1);
45 it++;
46 }
47 else if (*it == "NCI" || *it == "nci")
48 {
49 if (it + 1 == words.end())
50 {
51 app_error() << "Format error in ASCII integral file. NETOT \n";
52 APP_ABORT("Format error in ASCII integral file. NETOT \n");
53 }
54 ndet = atoi((++it)->c_str());
55 }
56 else if (*it == "UHF" || *it == "GHF")
57 {
58 if (it + 1 == words.end())
59 {
60 app_error() << "Format error in ASCII integral file. UHF/GHF \n";
61 APP_ABORT("Format error in ASCII integral file. UHF/GHF \n");
62 }
63 wfn_type = atoi((++it)->c_str());
64 switch (wfn_type)
65 {
66 case 0: {
67 app_log() << "Reading a RHF-type trial wave-function. \n";
68 break;
69 }
70 case 1: {
71 app_log() << "Reading a UHF-type trial wave-function. \n";
72 break;
73 }
74 case 2: {
75 app_log() << "Reading a GHF-type trial wave-function. \n";
76 break;
77 }
78 default: {
79 app_error() << "Unknown wave-function type in AFQMC/Utilities/readWfn.hpp: " << wfn_type << std::endl;
80 APP_ABORT("Unknown wave-function type in AFQMC/Utilities/readWfn.hpp \n");
81 }
82 }
83 }
84 else if (*it == "FullMO" || *it == "FULLMO")
85 {
86 fullMOMat = true;
87 }
88 else if (*it == "CMajor")
89 {
90 Cstyle = false;
91 }
92 }
93 getwords(words, in);
94 if (words.size() == 0)
95 app_error() << "Format error in ASCII integral file. End of file in header. \n";
96 } while ((words[0].find(std::string("/")) == std::string::npos &&
97 words[0].find(std::string("&END")) == std::string::npos));
98 }
99
skip_determinant(std::ifstream & in,bool Cstyle,bool fullMOMat,int NMO,int NAEA)100 void skip_determinant(std::ifstream& in, bool Cstyle, bool fullMOMat, int NMO, int NAEA)
101 {
102 int nread;
103 afqmc::ComplexType dummy;
104 if (Cstyle)
105 {
106 nread = fullMOMat ? NMO : NAEA;
107 for (int i = 0; i < NMO; i++)
108 for (int j = 0; j < nread; j++)
109 {
110 in >> dummy;
111 if (in.fail())
112 {
113 app_error() << "Problems reading ASCII file in readWfn.hpp. \n";
114 app_error() << i << " " << j << std::endl;
115 in.close();
116 APP_ABORT("Problems reading ASCII file in readWfn.hpp. \n");
117 }
118 }
119 }
120 else
121 {
122 nread = fullMOMat ? NMO : NAEA;
123 for (int j = 0; j < nread; j++)
124 for (int i = 0; i < NMO; i++)
125 {
126 in >> dummy;
127 if (in.fail())
128 {
129 app_error() << "Problems reading ASCII file in readWfn.hpp. \n";
130 app_error() << i << " " << j << std::endl;
131 in.close();
132 APP_ABORT("Problems reading wfn in readWfn.hpp.\n");
133 }
134 }
135 }
136 }
137
138 template<class Mat>
read_mat(std::ifstream & in,Mat && OrbMat,bool Cstyle,bool fullMOMat,int NMO,int NAEA)139 void read_mat(std::ifstream& in, Mat&& OrbMat, bool Cstyle, bool fullMOMat, int NMO, int NAEA)
140 {
141 int nread;
142 afqmc::ComplexType dummy;
143 if (Cstyle)
144 {
145 nread = fullMOMat ? NMO : NAEA;
146 for (int i = 0; i < NMO; i++)
147 for (int j = 0; j < nread; j++)
148 {
149 in >> dummy;
150 if (j < NAEA)
151 OrbMat[i][j] = dummy;
152 if (in.fail())
153 {
154 app_error() << "Problems reading ASCII file in readWfn.hpp. \n";
155 app_error() << i << " " << j << std::endl;
156 in.close();
157 APP_ABORT("Problems reading ASCII file in readWfn.hpp. \n");
158 }
159 }
160 }
161 else
162 {
163 nread = fullMOMat ? NMO : NAEA;
164 for (int j = 0; j < nread; j++)
165 for (int i = 0; i < NMO; i++)
166 {
167 in >> dummy;
168 if (j < NAEA)
169 OrbMat[i][j] = dummy;
170 if (in.fail())
171 {
172 app_error() << "Problems reading ASCII file in readWfn.hpp. \n";
173 app_error() << i << " " << j << std::endl;
174 in.close();
175 APP_ABORT("Problems reading wfn in readWfn.hpp.\n");
176 }
177 }
178 }
179 }
180
181 } // namespace
182
183 namespace afqmc
184 {
getWalkerType(std::string filename)185 WALKER_TYPES getWalkerType(std::string filename)
186 {
187 std::ifstream in;
188 in.open(filename.c_str());
189 if (in.fail())
190 {
191 app_error() << "Problems opening file: " << filename << std::endl;
192 APP_ABORT("Problems opening ASCII integral file. \n");
193 }
194
195 bool fullMOMat = false;
196 bool Cstyle = true;
197 int wfn_type = 0;
198 int ndet_in_file = -1;
199 std::string type;
200
201 read_header(in, type, wfn_type, fullMOMat, Cstyle, ndet_in_file);
202 in.close();
203
204 if (wfn_type == 0)
205 return CLOSED;
206 else if (wfn_type == 1)
207 return COLLINEAR;
208 else if (wfn_type == 2)
209 return NONCOLLINEAR;
210 else
211 return UNDEFINED_WALKER_TYPE;
212 }
213
getWalkerTypeHDF5(std::string filename,std::string type)214 WALKER_TYPES getWalkerTypeHDF5(std::string filename, std::string type)
215 {
216 hdf_archive dump;
217 if (!dump.open(filename, H5F_ACC_RDONLY))
218 {
219 std::cerr << " Error opening wavefunction file in read_info_from_wfn. \n";
220 APP_ABORT("");
221 }
222 if (!dump.push("Wavefunction", false))
223 {
224 std::cerr << " Error in getWalkerTypeHDF5: Group Wavefunction found. \n";
225 APP_ABORT("");
226 }
227 if (!dump.push(type, false))
228 {
229 std::cerr << " Error in getWalkerTypeHDF5: Group " << type << " not found. \n";
230 APP_ABORT("");
231 }
232
233 std::vector<int> Idata(5);
234 if (!dump.readEntry(Idata, "dims"))
235 {
236 std::cerr << " Error in getWalkerTypeHDF5: Problems reading dims. \n";
237 APP_ABORT("");
238 }
239
240 dump.pop();
241 int wfn_type = Idata[3];
242
243 if (wfn_type == 1)
244 return CLOSED;
245 else if (wfn_type == 2)
246 return COLLINEAR;
247 else if (wfn_type == 3)
248 return NONCOLLINEAR;
249 else
250 return UNDEFINED_WALKER_TYPE;
251 }
252
getWfnType(std::ifstream & in)253 std::string getWfnType(std::ifstream& in)
254 {
255 in.clear();
256 in.seekg(0, std::ios::beg);
257
258 bool fullMOMat = false;
259 bool Cstyle = true;
260 int wfn_type = 0;
261 int ndet_in_file = -1;
262 std::string type;
263 read_header(in, type, wfn_type, fullMOMat, Cstyle, ndet_in_file);
264
265 return type;
266 }
267
268 /*
269 * Reads ndets from the ascii file.
270 */
read_general_wavefunction(std::ifstream & in,int & ndets,WALKER_TYPES walker_type,boost::mpi3::shared_communicator & comm,int NMO,int NAEA,int NAEB,std::vector<PsiT_Matrix> & PsiT,std::vector<ComplexType> & ci)271 void read_general_wavefunction(std::ifstream& in,
272 int& ndets,
273 WALKER_TYPES walker_type,
274 boost::mpi3::shared_communicator& comm,
275 int NMO,
276 int NAEA,
277 int NAEB,
278 std::vector<PsiT_Matrix>& PsiT,
279 std::vector<ComplexType>& ci)
280 {
281 in.clear();
282 in.seekg(0, std::ios::beg);
283
284 assert(walker_type != UNDEFINED_WALKER_TYPE);
285 std::string type;
286 bool fullMOMat = false;
287 bool Cstyle = true;
288 int wfn_type = 0;
289 int ndet_in_file = -1;
290 int NEL = NAEA;
291 if (walker_type != CLOSED)
292 NEL += NAEB;
293
294 /*
295 * type:
296 * - matrix: Slater matrix for all terms in the expansion
297 */
298 read_header(in, type, wfn_type, fullMOMat, Cstyle, ndet_in_file);
299
300 std::transform(type.begin(), type.end(), type.begin(), (int (*)(int))tolower);
301
302 /*
303 * Expected order of inputs and tags:
304 * Coefficients:
305 * Determinant:
306 */
307
308 if (ndets <= 0)
309 ndets = ndet_in_file;
310
311 if (ndet_in_file < ndets)
312 APP_ABORT("Error: Requesting too many determinants from wfn file.\n");
313 if (type != "matrix")
314 APP_ABORT("Error: Expects type=matrix in read_general_wavefunction.\n");
315 if (wfn_type == 1 && walker_type == CLOSED)
316 APP_ABORT("Error in read_wavefunction: walker_type < wfn_type. \n");
317 if (wfn_type == 2 && (walker_type == CLOSED || walker_type == COLLINEAR))
318 APP_ABORT("Error in read_wavefunction: walker_type < wfn_type. \n");
319
320 if (walker_type == COLLINEAR && (NAEA == NAEB && wfn_type == 0))
321 app_log() << " MESSAGE: Using walker_type=colinear with a closed-shell wfn with NAEA==NAEB. \n"
322 << " Consider doing a closed shell calculation ( walker_type=closed in WalkerSet)\n";
323
324 ci.reserve(ndets);
325
326 std::string tag;
327 in >> tag;
328 if (tag != "Coefficients:")
329 APP_ABORT(" Error: Expecting Coefficients: tag in wavefunction file. \n");
330 ComplexType dummy;
331 for (int i = 0; i < ndet_in_file; i++)
332 {
333 in >> dummy;
334 if (i < ndets)
335 ci.emplace_back(dummy);
336 }
337
338 if (type == "matrix")
339 {
340 PsiT.reserve((walker_type != COLLINEAR) ? ndets : 2 * ndets);
341
342 if (wfn_type == 0)
343 {
344 boost::multi::array<ComplexType, 2> OrbMat({NMO, NAEA});
345 for (int i = 0, q = 0; i < ndets; i++)
346 {
347 if (comm.rank() == 0)
348 {
349 in >> tag >> q;
350 if (tag != "Determinant:" || q != i + 1)
351 APP_ABORT(" Error: Expecting Determinant: # tag in wavefunction file. \n");
352 read_mat(in, OrbMat, Cstyle, fullMOMat, NMO, NAEA);
353 }
354 if (walker_type == CLOSED)
355 {
356 PsiT.emplace_back(csr::shm::construct_csr_matrix_single_input<PsiT_Matrix>(OrbMat, 1e-8, 'H', comm));
357 }
358 else if (walker_type == COLLINEAR)
359 {
360 PsiT.emplace_back(csr::shm::construct_csr_matrix_single_input<PsiT_Matrix>(OrbMat, 1e-8, 'H', comm));
361 PsiT.emplace_back(
362 csr::shm::construct_csr_matrix_single_input<PsiT_Matrix>(OrbMat(OrbMat.extension(0), {0, NAEB}), 1e-8,
363 'H', comm));
364 }
365 else if (walker_type == NONCOLLINEAR)
366 {
367 APP_ABORT(" Error in readWfn: wfn_type==closed with walker_type=noncollinear.\n");
368 /*
369 boost::multi::array<ComplexType,2> Mat({2*NMO,NAEA+NAEB});
370 if(comm.rank()==0) {
371 std::fill_n(Mat.origin(),2*NMO*(NAEA+NAEB),ComplexType(0.0));
372 Mat({0,NMO},{0,NAEA}) = OrbMat;
373 Mat({NMO,2*NMO},{NAEA,NAEA+NAEB}) = OrbMat({0,NMO},{0,NAEB});
374 }
375 PsiT.emplace_back(csr::shm::construct_csr_matrix_single_input<PsiT_Matrix>(
376 OrbMat,1e-8,'H',comm));
377 */
378 }
379 }
380 }
381 else if (wfn_type == 1)
382 {
383 if (walker_type != COLLINEAR)
384 APP_ABORT(" Error in readWfn: wfn_type==collinear with walker_type!=collinear.\n");
385 boost::multi::array<ComplexType, 2> OrbMat({NMO, NAEA});
386 for (int i = 0, q = 0; i < ndets; i++)
387 {
388 if (comm.rank() == 0)
389 {
390 in >> tag >> q;
391 if (tag != "Determinant:" || q != i + 1)
392 APP_ABORT(" Error: Expecting Determinant: # tag in wavefunction file. \n");
393 read_mat(in, OrbMat, Cstyle, fullMOMat, NMO, NAEA);
394 }
395
396 PsiT.emplace_back(csr::shm::construct_csr_matrix_single_input<PsiT_Matrix>(OrbMat, 1e-8, 'H', comm));
397 if (comm.rank() == 0)
398 read_mat(in, OrbMat(OrbMat.extension(0), {0, NAEB}), Cstyle, fullMOMat, NMO, NAEB);
399 PsiT.emplace_back(
400 csr::shm::construct_csr_matrix_single_input<PsiT_Matrix>(OrbMat(OrbMat.extension(0), {0, NAEB}), 1e-8, 'H',
401 comm));
402 }
403 }
404 else if (wfn_type == 2)
405 {
406 if (walker_type != NONCOLLINEAR)
407 APP_ABORT(" Error in readWfn: wfn_type==collinear with walker_type!=collinear.\n");
408 if (NAEB != 0)
409 APP_ABORT(" Error in readWfn: walker_type==collinear with NAEB!=0.\n");
410 boost::multi::array<ComplexType, 2> OrbMat({2 * NMO, NAEA});
411 for (int i = 0, q = 0; i < ndets; i++)
412 {
413 if (comm.rank() == 0)
414 {
415 in >> tag >> q;
416 if (tag != "Determinant:" || q != i + 1)
417 APP_ABORT(" Error: Expecting Determinant: # tag in wavefunction file. \n");
418 read_mat(in, OrbMat, Cstyle, fullMOMat, 2 * NMO, NAEA);
419 }
420 PsiT.emplace_back(csr::shm::construct_csr_matrix_single_input<PsiT_Matrix>(OrbMat, 1e-8, 'H', comm));
421 }
422 } //type
423 }
424 else
425 {
426 APP_ABORT("Error: Unknown wavefunction type in file. Expected type matrix .\n");
427 }
428 }
429
430 // only root reads
431 /*
432 * The number of terms returned by this routine in PsiT depends on wfn_type, NOT on walker_type.
433 */
read_ph_wavefunction(std::ifstream & in,int & ndets,WALKER_TYPES walker_type,boost::mpi3::shared_communicator & comm,int NMO,int NAEA,int NAEB,std::vector<PsiT_Matrix> & PsiT)434 ph_excitations<int, ComplexType> read_ph_wavefunction(std::ifstream& in,
435 int& ndets,
436 WALKER_TYPES walker_type,
437 boost::mpi3::shared_communicator& comm,
438 int NMO,
439 int NAEA,
440 int NAEB,
441 std::vector<PsiT_Matrix>& PsiT)
442 {
443 if (comm.root())
444 {
445 in.clear();
446 in.seekg(0, std::ios::beg);
447 }
448
449 assert(walker_type != UNDEFINED_WALKER_TYPE);
450 bool fullMOMat = false;
451 bool Cstyle = true;
452 int wfn_type = 0;
453 int ndet_in_file = -1;
454 int NEL = NAEA;
455 bool mixed = false;
456 std::string type;
457 if (walker_type != CLOSED)
458 NEL += NAEB;
459
460 /*
461 * Expected order of inputs and tags:
462 * Reference:
463 * Configurations:
464 */
465
466 /*
467 * type:
468 * - occ: All determinants are specified with occupation numbers
469 * - mixed: mixed representation. Reference determinant in matrix form,
470 * determinant list (including reference) in occupation numbers.
471 *
472 * wfn_type:
473 * - 0: excitations out of a RHF reference
474 * NOTE: Does not mean perfect pairing, means excitations from a single reference
475 * - 1: excitations out of a UHF reference
476 */
477 if (comm.root())
478 {
479 read_header(in, type, wfn_type, fullMOMat, Cstyle, ndet_in_file);
480
481 std::transform(type.begin(), type.end(), type.begin(), (int (*)(int))tolower);
482
483
484 if (walker_type != COLLINEAR)
485 APP_ABORT(" Error: walker_type!=COLLINEAR not yet implemented in read_ph_wavefunction.\n");
486
487 if (ndet_in_file < ndets)
488 APP_ABORT("Error: Requesting too many determinants from wfn file.\n");
489 if (type != "occ" && type != "mixed")
490 APP_ABORT("Error: Expect type=occ or type==mixed in read_ph_wavefunction.\n");
491 if (wfn_type == 1 && walker_type == CLOSED)
492 APP_ABORT("Error in read_wavefunction: walker_type < wfn_type. \n");
493 if (wfn_type == 2 && (walker_type == CLOSED || walker_type == COLLINEAR))
494 APP_ABORT("Error in read_wavefunction: walker_type < wfn_type. \n");
495
496 if (type == "mixed")
497 mixed = true;
498
499 comm.broadcast_n(&mixed, 1, 0);
500 comm.broadcast_n(&ndet_in_file, 1, 0);
501 comm.broadcast_n(&wfn_type, 1, 0);
502 comm.broadcast_n(&fullMOMat, 1, 0);
503 }
504 else
505 {
506 comm.broadcast_n(&mixed, 1, 0);
507 comm.broadcast_n(&ndet_in_file, 1, 0);
508 comm.broadcast_n(&wfn_type, 1, 0);
509 comm.broadcast_n(&fullMOMat, 1, 0);
510 }
511
512 if (ndets <= 0)
513 ndets = ndet_in_file;
514
515 if (mixed)
516 { // read reference
517 int nmo_ = (walker_type == NONCOLLINEAR ? 2 * NMO : NMO);
518 if (not comm.root())
519 nmo_ = 0; // only root reads matrices
520 if (not fullMOMat)
521 APP_ABORT("Error: Wavefunction type mixed requires fullMOMat=true.\n");
522 PsiT.reserve((wfn_type != 1) ? 1 : 2);
523
524 boost::multi::array<ComplexType, 2> OrbMat({nmo_, nmo_});
525 if (comm.root())
526 {
527 std::string tag;
528 in >> tag;
529 if (tag != "Reference:")
530 APP_ABORT(" Error: Expecting Reference tag in wavefunction file. \n");
531 read_mat(in, OrbMat, Cstyle, fullMOMat, nmo_, nmo_);
532 }
533 PsiT.emplace_back(csr::shm::construct_csr_matrix_single_input<PsiT_Matrix>(OrbMat, 1e-8, 'H', comm));
534 if (wfn_type == 1)
535 {
536 if (comm.root())
537 read_mat(in, OrbMat, Cstyle, fullMOMat, NMO, NMO);
538 PsiT.emplace_back(csr::shm::construct_csr_matrix_single_input<PsiT_Matrix>(OrbMat, 1e-8, 'H', comm));
539 }
540 }
541
542 if (comm.root())
543 {
544 std::string tag;
545 in >> tag;
546 if (tag != "Configurations:")
547 {
548 app_error() << " tag: " << tag << std::endl;
549 APP_ABORT(" Error: Expecting Configurations: tag in wavefunction file. \n");
550 }
551 }
552
553 ComplexType ci;
554 // count number of k-particle excitations
555 // counts[0] has special meaning, it must be equal to NAEA+NAEB.
556 std::vector<size_t> counts_alpha(NAEA + 1);
557 std::vector<size_t> counts_beta(NAEB + 1);
558 // ugly but need dynamic memory allocation
559 std::vector<std::vector<int>> unique_alpha(NAEA + 1);
560 std::vector<std::vector<int>> unique_beta(NAEB + 1);
561 // reference configuration, taken as the first one right now
562 std::vector<int> refa;
563 std::vector<int> refb;
564 // space to read configurations
565 std::vector<int> confg;
566 // space for excitation string identifying the current configuration
567 std::vector<int> exct;
568 // record file position to come back
569 std::vector<int> Iwork; // work arrays for permutation calculation
570 std::streampos start;
571 if (comm.root())
572 {
573 confg.reserve(NAEA);
574 Iwork.resize(2 * NAEA);
575 exct.reserve(2 * NAEA);
576 start = in.tellg();
577 for (int i = 0; i < ndets; i++)
578 {
579 in >> ci;
580 if (in.fail())
581 APP_ABORT(" Error: Reading wfn file.\n");
582 // alpha
583 confg.clear();
584 for (int k = 0, q = 0; k < NAEA; k++)
585 {
586 in >> q;
587 if (in.fail())
588 APP_ABORT(" Error: Reading wfn file.\n");
589 if (q < 1 || q > NMO)
590 APP_ABORT("Error: Bad occupation number in wavefunction file. \n");
591 confg.emplace_back(q - 1);
592 }
593 if (i == 0)
594 {
595 refa = confg;
596 }
597 else
598 {
599 int np = get_excitation_number(true, refa, confg, exct, ci, Iwork);
600 push_excitation(exct, unique_alpha[np]);
601 }
602 // beta
603 confg.clear();
604 for (int k = 0, q = 0; k < NAEB; k++)
605 {
606 in >> q;
607 if (in.fail())
608 APP_ABORT(" Error: Reading wfn file.\n");
609 if (q <= NMO || q > 2 * NMO)
610 APP_ABORT("Error: Bad occupation number in wavefunction file. \n");
611 confg.emplace_back(q - 1);
612 }
613 if (i == 0)
614 {
615 refb = confg;
616 }
617 else
618 {
619 int np = get_excitation_number(true, refb, confg, exct, ci, Iwork);
620 push_excitation(exct, unique_beta[np]);
621 }
622 }
623 // now that we have all unique configurations, count
624 for (int i = 1; i <= NAEA; i++)
625 counts_alpha[i] = unique_alpha[i].size();
626 for (int i = 1; i <= NAEB; i++)
627 counts_beta[i] = unique_beta[i].size();
628 }
629 comm.broadcast_n(counts_alpha.begin(), counts_alpha.size());
630 comm.broadcast_n(counts_beta.begin(), counts_beta.size());
631 // using int for now, but should move to short later when everything works well
632 // ph_struct stores the reference configuration on the index [0]
633 ph_excitations<int, ComplexType> ph_struct(ndets, NAEA, NAEB, counts_alpha, counts_beta, shared_allocator<int>(comm));
634
635 if (comm.root())
636 {
637 in.clear();
638 in.seekg(start);
639 std::map<int, int> refa2loc;
640 for (int i = 0; i < NAEA; i++)
641 refa2loc[refa[i]] = i;
642 std::map<int, int> refb2loc;
643 for (int i = 0; i < NAEB; i++)
644 refb2loc[refb[i]] = i;
645 // add reference
646 ph_struct.add_reference(refa, refb);
647 // add unique configurations
648 // alpha
649 for (int n = 1; n < unique_alpha.size(); n++)
650 for (std::vector<int>::iterator it = unique_alpha[n].begin(); it < unique_alpha[n].end(); it += (2 * n))
651 ph_struct.add_alpha(n, it);
652 // beta
653 for (int n = 1; n < unique_beta.size(); n++)
654 for (std::vector<int>::iterator it = unique_beta[n].begin(); it < unique_beta[n].end(); it += (2 * n))
655 ph_struct.add_beta(n, it);
656 // read configurations
657 int alpha_index;
658 int beta_index;
659 int np;
660 for (int i = 0; i < ndets; i++)
661 {
662 in >> ci;
663 if (in.fail())
664 APP_ABORT(" Error: Reading wfn file.\n");
665 confg.clear();
666 for (int k = 0, q = 0; k < NAEA; k++)
667 {
668 in >> q;
669 if (in.fail())
670 APP_ABORT(" Error: Reading wfn file.\n");
671 if (q < 1 || q > NMO)
672 APP_ABORT("Error: Bad occupation number in wavefunction file. \n");
673 confg.emplace_back(q - 1);
674 }
675 np = get_excitation_number(true, refa, confg, exct, ci, Iwork);
676 alpha_index =
677 ((np == 0) ? (0)
678 : (find_excitation(exct, unique_alpha[np]) + ph_struct.number_of_unique_smaller_than(np)[0]));
679 confg.clear();
680 for (int k = 0, q = 0; k < NAEB; k++)
681 {
682 in >> q;
683 if (in.fail())
684 APP_ABORT(" Error: Reading wfn file.\n");
685 if (q <= NMO || q > 2 * NMO)
686 APP_ABORT("Error: Bad occupation number in wavefunction file. \n");
687 confg.emplace_back(q - 1);
688 }
689 np = get_excitation_number(true, refb, confg, exct, ci, Iwork);
690 beta_index =
691 ((np == 0) ? (0) : (find_excitation(exct, unique_beta[np]) + ph_struct.number_of_unique_smaller_than(np)[1]));
692 ph_struct.add_configuration(alpha_index, beta_index, ci);
693 }
694 }
695 comm.barrier();
696 return ph_struct;
697 }
698
read_ph_wavefunction_hdf(hdf_archive & dump,std::vector<ComplexType> & ci_coeff,std::vector<int> & occs,int & ndets,WALKER_TYPES walker_type,boost::mpi3::shared_communicator & comm,int NMO,int NAEA,int NAEB,std::vector<PsiT_Matrix> & PsiT,std::string & type)699 void read_ph_wavefunction_hdf(hdf_archive& dump,
700 std::vector<ComplexType>& ci_coeff,
701 std::vector<int>& occs,
702 int& ndets,
703 WALKER_TYPES walker_type,
704 boost::mpi3::shared_communicator& comm,
705 int NMO,
706 int NAEA,
707 int NAEB,
708 std::vector<PsiT_Matrix>& PsiT,
709 std::string& type)
710 {
711 using Alloc = shared_allocator<ComplexType>;
712 assert(walker_type != UNDEFINED_WALKER_TYPE);
713 int wfn_type = 0;
714 int NEL = NAEA;
715 bool mixed = false;
716 if (walker_type != CLOSED)
717 NEL += NAEB;
718
719 /*
720 * Expected order of inputs and tags:
721 * Reference:
722 * Configurations:
723 */
724
725 /*
726 * type:
727 * - occ: All determinants are specified with occupation numbers
728 *
729 * wfn_type:
730 * - 0: excitations out of a RHF reference
731 * NOTE: Does not mean perfect pairing, means excitations from a single reference
732 * - 1: excitations out of a UHF reference
733 */
734 WALKER_TYPES wtype;
735 getCommonInput(dump, NMO, NAEA, NAEB, ndets, ci_coeff, wtype, comm.root());
736 if (walker_type != wtype && NAEA != NAEB)
737 APP_ABORT(" Error: When Different walker_type between hdf5 and xml inputs when NAEA!=NAEB (not allowed). \n");
738 if (walker_type != COLLINEAR)
739 APP_ABORT(" Error: walker_type!=COLLINEAR not yet implemented in read_ph_wavefunction.\n");
740
741 int type_;
742 if (!dump.readEntry(type_, "type"))
743 {
744 app_error() << " Error in WavefunctionFactory::fromHDF5(): Problems reading type. \n";
745 APP_ABORT("");
746 }
747 if (type_ == 0)
748 type = "occ";
749 else
750 type = "mixed";
751 if (type == "mixed")
752 mixed = true;
753
754 if (mixed)
755 { // read reference
756 int nmo_ = (walker_type == NONCOLLINEAR ? 2 * NMO : NMO);
757 if (not comm.root())
758 nmo_ = 0; // only root reads matrices
759 PsiT.reserve((wfn_type != 1) ? 1 : 2);
760
761 if (!dump.push(std::string("PsiT_") + std::to_string(0), false))
762 {
763 app_error() << " Error in WavefunctionFactory: Group PsiT not found. \n";
764 APP_ABORT("");
765 }
766 PsiT.emplace_back(csr_hdf5::HDF2CSR<PsiT_Matrix, Alloc>(dump, comm));
767 dump.pop();
768 if (wfn_type == 1)
769 {
770 if (wtype == CLOSED)
771 {
772 if (!dump.push(std::string("PsiT_") + std::to_string(0), false))
773 {
774 app_error() << " Error in WavefunctionFactory: Group PsiT not found. \n";
775 APP_ABORT("");
776 }
777 }
778 else if (wtype == COLLINEAR)
779 {
780 if (!dump.push(std::string("PsiT_") + std::to_string(1), false))
781 {
782 app_error() << " Error in WavefunctionFactory: Group PsiT not found. \n";
783 APP_ABORT("");
784 }
785 }
786 PsiT.emplace_back(csr_hdf5::HDF2CSR<PsiT_Matrix, Alloc>(dump, comm));
787 dump.pop();
788 }
789 }
790 if (!dump.readEntry(occs, "occs"))
791 APP_ABORT("Error reading occs array.\n");
792 comm.barrier();
793 }
794
build_ph_struct(std::vector<ComplexType> ci_coeff,boost::multi::array_ref<int,2> & occs,int ndets,boost::mpi3::shared_communicator & comm,int NMO,int NAEA,int NAEB)795 ph_excitations<int, ComplexType> build_ph_struct(std::vector<ComplexType> ci_coeff,
796 boost::multi::array_ref<int, 2>& occs,
797 int ndets,
798 boost::mpi3::shared_communicator& comm,
799 int NMO,
800 int NAEA,
801 int NAEB)
802 {
803 using Alloc = shared_allocator<ComplexType>;
804
805 ComplexType ci;
806 // count number of k-particle excitations
807 // counts[0] has special meaning, it must be equal to NAEA+NAEB.
808 std::vector<size_t> counts_alpha(NAEA + 1);
809 std::vector<size_t> counts_beta(NAEB + 1);
810 // ugly but need dynamic memory allocation
811 std::vector<std::vector<int>> unique_alpha(NAEA + 1);
812 std::vector<std::vector<int>> unique_beta(NAEB + 1);
813 // reference configuration, taken as the first one right now
814 std::vector<int> refa;
815 std::vector<int> refb;
816 // space to read configurations
817 std::vector<int> confg;
818 // space for excitation string identifying the current configuration
819 std::vector<int> exct;
820 // record file position to come back
821 std::vector<int> Iwork; // work arrays for permutation calculation
822 std::streampos start;
823 if (comm.root())
824 {
825 confg.reserve(NAEA);
826 Iwork.resize(2 * NAEA);
827 exct.reserve(2 * NAEA);
828 for (int i = 0; i < ndets; i++)
829 {
830 ci = ci_coeff[i];
831 // alpha
832 confg.clear();
833 for (int k = 0, q = 0; k < NAEA; k++)
834 {
835 q = occs[i][k];
836 if (q < 0 || q >= NMO)
837 APP_ABORT("Error: Bad occupation number " << q << " in determinant " << i << " in wavefunction file. \n");
838 confg.emplace_back(q);
839 }
840 if (i == 0)
841 {
842 refa = confg;
843 }
844 else
845 {
846 int np = get_excitation_number(true, refa, confg, exct, ci, Iwork);
847 push_excitation(exct, unique_alpha[np]);
848 }
849 // beta
850 confg.clear();
851 for (int k = 0, q = 0; k < NAEB; k++)
852 {
853 q = occs[i][NAEA + k];
854 if (q < NMO || q >= 2 * NMO)
855 APP_ABORT("Error: Bad occupation number " << q << " in determinant " << i << " in wavefunction file. \n");
856 confg.emplace_back(q);
857 }
858 if (i == 0)
859 {
860 refb = confg;
861 }
862 else
863 {
864 int np = get_excitation_number(true, refb, confg, exct, ci, Iwork);
865 push_excitation(exct, unique_beta[np]);
866 }
867 }
868 // now that we have all unique configurations, count
869 for (int i = 1; i <= NAEA; i++)
870 counts_alpha[i] = unique_alpha[i].size();
871 for (int i = 1; i <= NAEB; i++)
872 counts_beta[i] = unique_beta[i].size();
873 }
874 comm.broadcast_n(counts_alpha.begin(), counts_alpha.size());
875 comm.broadcast_n(counts_beta.begin(), counts_beta.size());
876 // using int for now, but should move to short later when everything works well
877 // ph_struct stores the reference configuration on the index [0]
878 ph_excitations<int, ComplexType> ph_struct(ndets, NAEA, NAEB, counts_alpha, counts_beta, shared_allocator<int>(comm));
879
880 if (comm.root())
881 {
882 std::map<int, int> refa2loc;
883 for (int i = 0; i < NAEA; i++)
884 refa2loc[refa[i]] = i;
885 std::map<int, int> refb2loc;
886 for (int i = 0; i < NAEB; i++)
887 refb2loc[refb[i]] = i;
888 // add reference
889 ph_struct.add_reference(refa, refb);
890 // add unique configurations
891 // alpha
892 for (int n = 1; n < unique_alpha.size(); n++)
893 for (std::vector<int>::iterator it = unique_alpha[n].begin(); it < unique_alpha[n].end(); it += (2 * n))
894 ph_struct.add_alpha(n, it);
895 // beta
896 for (int n = 1; n < unique_beta.size(); n++)
897 for (std::vector<int>::iterator it = unique_beta[n].begin(); it < unique_beta[n].end(); it += (2 * n))
898 ph_struct.add_beta(n, it);
899 // read configurations
900 int alpha_index;
901 int beta_index;
902 int np;
903 for (int i = 0; i < ndets; i++)
904 {
905 ci = ci_coeff[i];
906 confg.clear();
907 for (int k = 0, q = 0; k < NAEA; k++)
908 {
909 q = occs[i][k];
910 if (q < 0 || q >= NMO)
911 APP_ABORT("Error: Bad occupation number " << q << " in determinant " << i << " in wavefunction file. \n");
912 confg.emplace_back(q);
913 }
914 np = get_excitation_number(true, refa, confg, exct, ci, Iwork);
915 alpha_index =
916 ((np == 0) ? (0)
917 : (find_excitation(exct, unique_alpha[np]) + ph_struct.number_of_unique_smaller_than(np)[0]));
918 confg.clear();
919 for (int k = 0, q = 0; k < NAEB; k++)
920 {
921 q = occs[i][NAEA + k];
922 if (q < NMO || q >= 2 * NMO)
923 APP_ABORT("Error: Bad occupation number " << q << " in determinant " << i << " in wavefunction file. \n");
924 confg.emplace_back(q);
925 }
926 np = get_excitation_number(true, refb, confg, exct, ci, Iwork);
927 beta_index =
928 ((np == 0) ? (0) : (find_excitation(exct, unique_beta[np]) + ph_struct.number_of_unique_smaller_than(np)[1]));
929 ph_struct.add_configuration(alpha_index, beta_index, ci);
930 }
931 }
932 comm.barrier();
933 return ph_struct;
934 }
935
936 /*
937 * Read trial wavefunction information from file.
938 */
getCommonInput(hdf_archive & dump,int NMO,int NAEA,int NAEB,int & ndets_to_read,std::vector<ComplexType> & ci,WALKER_TYPES & walker_type,bool root)939 void getCommonInput(hdf_archive& dump,
940 int NMO,
941 int NAEA,
942 int NAEB,
943 int& ndets_to_read,
944 std::vector<ComplexType>& ci,
945 WALKER_TYPES& walker_type,
946 bool root)
947 {
948 // check for consistency in parameters
949 std::vector<int> dims(5);
950 if (!dump.readEntry(dims, "dims"))
951 {
952 app_error() << " Error in getCommonInput(): Problems reading dims. \n";
953 APP_ABORT("");
954 }
955 if (NMO != dims[0])
956 {
957 app_error() << " Error in getCommonInput(): Inconsistent NMO . \n";
958 APP_ABORT("");
959 }
960 if (NAEA != dims[1])
961 {
962 app_error() << " Error in getCommonInput(): Inconsistent NAEA. \n";
963 APP_ABORT("");
964 }
965 if (NAEB != dims[2])
966 {
967 app_error() << " Error in getCommonInput(): Inconsistent NAEB. \n";
968 APP_ABORT("");
969 }
970 walker_type = afqmc::initWALKER_TYPES(dims[3]);
971 // just read walker_type, to allow flexibility
972 //if(walker_type != dims[3]) {
973 // app_error()<<" Error in getCommonInput(): Inconsistent walker_type. \n";
974 // APP_ABORT("");
975 //}
976 if (ndets_to_read < 1)
977 ndets_to_read = dims[4];
978 app_log() << " - Number of determinants in trial wavefunction: " << ndets_to_read << "\n";
979 if (ndets_to_read > dims[4])
980 {
981 app_error() << " Error in getCommonInput(): Inconsistent ndets_to_read. \n";
982 APP_ABORT("");
983 }
984 ci.resize(ndets_to_read);
985 if (!dump.readEntry(ci, "ci_coeffs"))
986 {
987 app_error() << " Error in getCommonInput(): Problems reading ci_coeffs. \n";
988 APP_ABORT("");
989 }
990 app_log() << " - Coefficient of first determinant: " << ci[0] << "\n";
991 }
992
993
994 // modify for multideterminant case based on type
readWfn(std::string fileName,boost::multi::array<ComplexType,3> & OrbMat,int NMO,int NAEA,int NAEB,int det)995 int readWfn(std::string fileName, boost::multi::array<ComplexType, 3>& OrbMat, int NMO, int NAEA, int NAEB, int det)
996 {
997 std::ifstream in;
998 in.open(fileName.c_str());
999 if (in.fail())
1000 {
1001 app_error() << "Problems opening ASCII integral file: " << fileName << std::endl;
1002 APP_ABORT("Problems opening ASCII integral file. \n");
1003 }
1004
1005 bool fullMOMat = false;
1006 bool Cstyle = true;
1007 int wfn_type;
1008 int ndet = 1;
1009 std::string type;
1010
1011 read_header(in, type, wfn_type, fullMOMat, Cstyle, ndet);
1012
1013 if (ndet != 1)
1014 APP_ABORT("Error: readWfn is for single determinant wave functions. \n");
1015 if (type != "matrix")
1016 APP_ABORT("Error: Only type=matrix accepted in readWfn. \n");
1017
1018 ComplexType dummy;
1019 std::string tag;
1020 in >> tag >> dummy;
1021 if (tag != "Coefficients:")
1022 APP_ABORT(" Error: Expecting Coefficients: tag in wavefunction file. \n");
1023
1024 int q;
1025 in >> tag >> q;
1026 if (tag != "Determinant:" || q != 1)
1027 APP_ABORT(" Error: Expecting Determinant: 1 tag in wavefunction file. \n");
1028
1029 if (wfn_type == 0)
1030 {
1031 OrbMat.reextent({1, NMO, NAEA});
1032 read_mat(in, OrbMat[0], Cstyle, fullMOMat, NMO, NAEA);
1033 }
1034 else if (wfn_type == 1)
1035 {
1036 OrbMat.reextent({2, NMO, NAEA});
1037 read_mat(in, OrbMat[0], Cstyle, fullMOMat, NMO, NAEA);
1038 read_mat(in, OrbMat[1], Cstyle, fullMOMat, NMO, NAEB);
1039 }
1040 else if (wfn_type == 2)
1041 {
1042 OrbMat.reextent({1, 2 * NMO, NAEA + NAEB});
1043 read_mat(in, OrbMat[0], Cstyle, fullMOMat, 2 * NMO, NAEA + NAEB);
1044 } //type
1045
1046 in.close();
1047
1048 return wfn_type;
1049 }
1050
1051 } // namespace afqmc
1052
1053 } // namespace qmcplusplus
1054 #endif
1055