From ac8075e89446ce19a05a09384d8506ea40f938eb Mon Sep 17 00:00:00 2001 From: camelto2 Date: Wed, 17 Feb 2021 14:15:13 -0700 Subject: [PATCH 01/16] multislater to allow arbitrary particle groups --- .../Fermion/SlaterDetBuilder.cpp | 49 +++++++++++++++---- 1 file changed, 39 insertions(+), 10 deletions(-) diff --git a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp index bb4f4f97d9..dd85c02e66 100644 --- a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp +++ b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp @@ -77,7 +77,7 @@ WaveFunctionComponent* SlaterDetBuilder::buildComponent(xmlNodePtr cur) ///save the current node xmlNodePtr curRoot = cur; xmlNodePtr BFnode = nullptr; - bool success = true; + bool success = true; std::string cname, tname; std::map spomap; bool multiDet = false; @@ -207,25 +207,30 @@ WaveFunctionComponent* SlaterDetBuilder::buildComponent(xmlNodePtr cur) std::string spo_alpha_name; std::string spo_beta_name; std::string fastAlg; + int nGroups = -1; OhmmsAttributeSet spoAttrib; spoAttrib.add(spo_alpha_name, "spo_up"); spoAttrib.add(spo_beta_name, "spo_dn"); spoAttrib.add(fastAlg, "Fast", {"", "yes", "no"}, TagStatus::DEPRECATED); spoAttrib.add(msd_algorithm, "algorithm", {"", "precomputed_table_method", "table_method", "all_determinants"}); + spoAttrib.add(nGroups, "ngroups"); spoAttrib.put(cur); if (msd_algorithm.empty()) { - if(fastAlg.empty()) + if (fastAlg.empty()) msd_algorithm = "precomputed_table_method"; - else if(fastAlg == "yes") + else if (fastAlg == "yes") msd_algorithm = "table_method"; else msd_algorithm = "all_determinants"; } - SPOSetPtr spo_alpha = sposet_builder_factory_.getSPOSet(spo_alpha_name); - SPOSetPtr spo_beta = sposet_builder_factory_.getSPOSet(spo_beta_name); + SPOSetPtr spo_alpha; + SPOSetPtr spo_beta; + std::vector SPOSetPtrs; + spo_alpha = sposet_builder_factory_.getSPOSet(spo_alpha_name); + spo_beta = sposet_builder_factory_.getSPOSet(spo_beta_name); std::unique_ptr spo_alpha_clone, spo_beta_clone; if (spo_alpha == nullptr) { @@ -245,6 +250,29 @@ WaveFunctionComponent* SlaterDetBuilder::buildComponent(xmlNodePtr cur) else spo_beta_clone.reset(spo_beta->makeClone()); + OhmmsAttributeSet grpAttrib; + std::vector spoNames(nGroups); + for (int grp = 0; grp < nGroups; grp++) + grpAttrib.add(spoNames[grp], "spo_" + std::to_string(grp)); + grpAttrib.put(cur); + app_log() << "GROUPS" << std::endl; + for (int grp = 0; grp < nGroups; grp++) + app_log() << spoNames[grp] << std::endl; + std::vector spos(nGroups, nullptr); + std::vector> spo_clones(nGroups); + for (int grp = 0; grp < nGroups; grp++) + { + spos[grp] = sposet_builder_factory_.getSPOSet(spoNames[grp]); + if (spos[grp] == nullptr) + { + app_error() << "In SlaterDetBuilder: SPOSet \"" << spoNames[grp] + << "\" is not found. Expected for MultiSlaterDeterminant.\n"; + abort(); + } + spo_clones[grp].reset(spos[grp]->makeClone()); + } + + if (msd_algorithm == "precomputed_table_method" || msd_algorithm == "table_method") { app_summary() << " Using Bryan's table method." << std::endl; @@ -253,11 +281,12 @@ WaveFunctionComponent* SlaterDetBuilder::buildComponent(xmlNodePtr cur) APP_ABORT("Backflow is not implemented with the table method."); } - std::vector> dets; - app_log() << " Creating base determinant (up) for MSD expansion. \n"; - dets.emplace_back(std::make_unique(std::move(spo_alpha_clone), 0)); - app_log() << " Creating base determinant (down) for MSD expansion. \n"; - dets.emplace_back(std::make_unique(std::move(spo_beta_clone), 1)); + std::vector> dets(nGroups); + for (int grp = 0; grp < nGroups; grp++) + { + app_log() << " Creating base determinant (" << grp << ") for MSD expansion. \n"; + dets.emplace_back(std::make_unique(std::move(spo_clones[grp]), grp)); + } if (msd_algorithm == "precomputed_table_method") { From 6992a4e969f564546b1fc8926fe41e397da8e8cc Mon Sep 17 00:00:00 2001 From: camelto2 Date: Wed, 17 Feb 2021 14:35:24 -0700 Subject: [PATCH 02/16] allow for old input for the time being --- .../Fermion/SlaterDetBuilder.cpp | 92 +++++++++++-------- 1 file changed, 54 insertions(+), 38 deletions(-) diff --git a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp index dd85c02e66..bb00d81667 100644 --- a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp +++ b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp @@ -226,52 +226,58 @@ WaveFunctionComponent* SlaterDetBuilder::buildComponent(xmlNodePtr cur) msd_algorithm = "all_determinants"; } + //old format SPOSetPtr spo_alpha; SPOSetPtr spo_beta; - std::vector SPOSetPtrs; - spo_alpha = sposet_builder_factory_.getSPOSet(spo_alpha_name); - spo_beta = sposet_builder_factory_.getSPOSet(spo_beta_name); std::unique_ptr spo_alpha_clone, spo_beta_clone; - if (spo_alpha == nullptr) - { - app_error() << "In SlaterDetBuilder: SPOSet \"" << spo_alpha_name - << "\" is not found. Expected for MultiSlaterDeterminant.\n"; - abort(); - } - else - spo_alpha_clone.reset(spo_alpha->makeClone()); - - if (spo_beta == nullptr) - { - app_error() << "In SlaterDetBuilder: SPOSet \"" << spo_beta_name - << "\" is not found. Expected for MultiSlaterDeterminant.\n"; - abort(); - } - else - spo_beta_clone.reset(spo_beta->makeClone()); - - OhmmsAttributeSet grpAttrib; + //new format std::vector spoNames(nGroups); - for (int grp = 0; grp < nGroups; grp++) - grpAttrib.add(spoNames[grp], "spo_" + std::to_string(grp)); - grpAttrib.put(cur); - app_log() << "GROUPS" << std::endl; - for (int grp = 0; grp < nGroups; grp++) - app_log() << spoNames[grp] << std::endl; std::vector spos(nGroups, nullptr); std::vector> spo_clones(nGroups); - for (int grp = 0; grp < nGroups; grp++) + if (nGroups == -1) { - spos[grp] = sposet_builder_factory_.getSPOSet(spoNames[grp]); - if (spos[grp] == nullptr) + spo_alpha = sposet_builder_factory_.getSPOSet(spo_alpha_name); + spo_beta = sposet_builder_factory_.getSPOSet(spo_beta_name); + if (spo_alpha == nullptr) { - app_error() << "In SlaterDetBuilder: SPOSet \"" << spoNames[grp] + app_error() << "In SlaterDetBuilder: SPOSet \"" << spo_alpha_name << "\" is not found. Expected for MultiSlaterDeterminant.\n"; abort(); } - spo_clones[grp].reset(spos[grp]->makeClone()); - } + else + spo_alpha_clone.reset(spo_alpha->makeClone()); + if (spo_beta == nullptr) + { + app_error() << "In SlaterDetBuilder: SPOSet \"" << spo_beta_name + << "\" is not found. Expected for MultiSlaterDeterminant.\n"; + abort(); + } + else + spo_beta_clone.reset(spo_beta->makeClone()); + } + else + { + app_warning() << "Using new group method" << std::endl; + OhmmsAttributeSet grpAttrib; + for (int grp = 0; grp < nGroups; grp++) + grpAttrib.add(spoNames[grp], "spo_" + std::to_string(grp)); + grpAttrib.put(cur); + app_log() << "GROUPS" << std::endl; + for (int grp = 0; grp < nGroups; grp++) + app_log() << spoNames[grp] << std::endl; + for (int grp = 0; grp < nGroups; grp++) + { + spos[grp] = sposet_builder_factory_.getSPOSet(spoNames[grp]); + if (spos[grp] == nullptr) + { + app_error() << "In SlaterDetBuilder: SPOSet \"" << spoNames[grp] + << "\" is not found. Expected for MultiSlaterDeterminant.\n"; + abort(); + } + spo_clones[grp].reset(spos[grp]->makeClone()); + } + } if (msd_algorithm == "precomputed_table_method" || msd_algorithm == "table_method") { @@ -281,11 +287,21 @@ WaveFunctionComponent* SlaterDetBuilder::buildComponent(xmlNodePtr cur) APP_ABORT("Backflow is not implemented with the table method."); } - std::vector> dets(nGroups); - for (int grp = 0; grp < nGroups; grp++) + std::vector> dets; + if (nGroups == -1) + { + app_log() << " Creating base determinant (up) for MSD expansion. \n"; + dets.emplace_back(std::make_unique(std::move(spo_alpha_clone), 0)); + app_log() << " Creating base determinant (down) for MSD expansion. \n"; + dets.emplace_back(std::make_unique(std::move(spo_beta_clone), 1)); + } + else { - app_log() << " Creating base determinant (" << grp << ") for MSD expansion. \n"; - dets.emplace_back(std::make_unique(std::move(spo_clones[grp]), grp)); + for (int grp = 0; grp < nGroups; grp++) + { + app_log() << " Creating base determinant (" << grp << ") for MSD expansion. \n"; + dets.emplace_back(std::make_unique(std::move(spo_clones[grp]), grp)); + } } if (msd_algorithm == "precomputed_table_method") From 68936f62b06adaa6ce89ee6271f13136d4f17725 Mon Sep 17 00:00:00 2001 From: camelto2 Date: Wed, 17 Feb 2021 15:27:30 -0700 Subject: [PATCH 03/16] generalize createMSDFast need to generalize readDetList in each case --- .../Fermion/SlaterDetBuilder.cpp | 810 ++++++++++++++++-- .../Fermion/SlaterDetBuilder.h | 19 + 2 files changed, 767 insertions(+), 62 deletions(-) diff --git a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp index bb00d81667..3b99880abd 100644 --- a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp +++ b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp @@ -586,7 +586,7 @@ bool SlaterDetBuilder::putDeterminant(xmlNodePtr cur, int spin_group) } bool SlaterDetBuilder::createMSDFast(std::vector>& Dets, - std::vector>& C2node, + std::vector>& C2nodes, std::vector& C, std::vector& CSFcoeff, std::vector& DetsPerCSF, @@ -598,12 +598,18 @@ bool SlaterDetBuilder::createMSDFast(std::vector uniqueConfg_up, uniqueConfg_dn; + + std::vector> uniqueConfgs; std::vector CItags; bool optimizeCI; - int nels_up = targetPtcl.groupsize(0); - int nels_dn = targetPtcl.groupsize(1); + + int nGroups = targetPtcl.groups(); + assert(nGroups == Dets.size()); + std::vector nptcls(nGroups); + for (int grp = 0; grp < nGroups; grp++) + nptcls[grp] = targetPtcl.groupsize(grp); + //Check id multideterminants are in HDF5 xmlNodePtr curTemp = cur, DetListNode = nullptr; @@ -622,51 +628,37 @@ bool SlaterDetBuilder::createMSDFast(std::vectorReferenceDeterminant = 0; // for now - Dets[0]->NumDets = uniqueConfg_up.size(); - std::vector& list_up = *Dets[0]->ciConfigList; - list_up.resize(uniqueConfg_up.size()); - for (int i = 0; i < list_up.size(); i++) + for (int grp = 0; grp < nGroups; grp++) { - list_up[i].occup.resize(nels_up); - int cnt = 0; - for (int k = 0; k < uniqueConfg_up[i].occup.size(); k++) - if (uniqueConfg_up[i].occup[k]) - list_up[i].occup[cnt++] = k; - if (cnt != nels_up) + // you should choose the det with highest weight for reference + Dets[grp]->ReferenceDeterminant = 0; // for now + Dets[grp]->NumDets = uniqueConfgs[grp].size(); + std::vector& list = *Dets[grp]->ciConfigList; + list.resize(uniqueConfgs[grp].size()); + for (int i = 0; i < list.size(); i++) { - APP_ABORT("Error in SlaterDetBuilder::createMSDFast, problems with ci configuration list. \n"); + list[i].occup.resize(nptcls[grp]); + int cnt = 0; + for (int k = 0; k < uniqueConfgs[grp][i].occup.size(); k++) + if (uniqueConfgs[grp][i].occup[k]) + list[i].occup[cnt++] = k; + if (cnt != nptcls[grp]) + { + APP_ABORT("Error in SlaterDetBuilder::createMSDFast for ptcl group " + << grp << ", problems with ci configuration list. \n"); + } } + Dets[grp]->set(targetPtcl.first(grp), nptcls[grp], Dets[grp]->Phi->getOrbitalSetSize()); } - Dets[0]->set(targetPtcl.first(0), nels_up, Dets[0]->Phi->getOrbitalSetSize()); - Dets[1]->ReferenceDeterminant = 0; // for now - Dets[1]->NumDets = uniqueConfg_dn.size(); - std::vector& list_dn = *Dets[1]->ciConfigList; - list_dn.resize(uniqueConfg_dn.size()); - for (int i = 0; i < list_dn.size(); i++) - { - list_dn[i].occup.resize(nels_dn); - int cnt = 0; - for (int k = 0; k < uniqueConfg_dn[i].occup.size(); k++) - if (uniqueConfg_dn[i].occup[k]) - list_dn[i].occup[cnt++] = k; - if (cnt != nels_dn) - { - APP_ABORT("Error in SlaterDetBuilder::createMSDFast, problems with ci configuration list. \n"); - } - } - Dets[1]->set(targetPtcl.first(1), nels_dn, Dets[1]->Phi->getOrbitalSetSize()); if (CSFcoeff.size() == 1) optimizeCI = false; if (optimizeCI) @@ -700,26 +692,37 @@ bool SlaterDetBuilder::createMSDFast(std::vectorOptimizable == true || Dets[1]->Optimizable == true) + bool any_optimizable = false; + for (int grp = 0; grp < nGroups; grp++) + { + if (Dets[grp]->Optimizable == true) + { + any_optimizable = true; + break; + } + } + if (any_optimizable) { - //safety checks for orbital optimization - if (Dets[0]->Optimizable != Dets[1]->Optimizable) - APP_ABORT("Optimizing the SPOSet of only one spin is not supported!\n"); + for (int grp = 0; grp < nGroups; grp++) + { + if (Dets[grp]->Optimizable != true) + APP_ABORT("Optimizing the SPOSet of only only species is not supported!\n"); + } if (usingCSF) APP_ABORT("Currently, Using CSF is not available with MSJ Orbital Optimization!\n"); - //checks that the hartree fock determinant is the first in the multislater expansion - for (int i = 0; i < nels_up; i++) + for (int grp = 0; grp < nGroups; grp++) { - if ((uniqueConfg_up[0].occup[i] != true) || (uniqueConfg_dn[0].occup[i] != true)) - APP_ABORT( - "The Hartree Fock Reference Determinant must be first in the Multi-Slater expansion for the input!\n"); + for (int i = 0; i < nptcls[grp]; i++) + { + if (uniqueConfgs[grp][0].occup[i] != true) + APP_ABORT( + "The Hartee Fock Reference Determinant must be the first in the Multi-Slater expansion for the input!\n"); + } } + app_warning() << "Unrestricted Orbital Optimization will be performed. Spin symmetry is not guaranteed to be " + "preserved!\n"; - app_log() << "WARNING: Unrestricted Orbital Optimization will be performed. Spin symmetry is not guaranteed to be " - "preserved!\n"; - - // mark the overall optimization flag Optimizable = true; } @@ -1298,18 +1301,700 @@ bool SlaterDetBuilder::readDetList(xmlNodePtr cur, return success; } +bool SlaterDetBuilder::readDetList(xmlNodePtr cur, + std::vector>& uniqueConfgs, + std::vector>& C2nodes, + std::vector& CItags, + std::vector& coeff, + bool& optimizeCI, + std::vector& nptcls, + std::vector& CSFcoeff, + std::vector& DetsPerCSF, + std::vector& CSFexpansion, + bool& usingCSF) +{ + APP_ABORT("readDetList for arbitrary ptcls groups not yet implemented.\n"); + /* + bool success = true; + uniqueConfg_up.clear(); + uniqueConfg_dn.clear(); + C2node_up.clear(); + C2node_dn.clear(); + CItags.clear(); + coeff.clear(); + CSFcoeff.clear(); + DetsPerCSF.clear(); + CSFexpansion.clear(); + std::vector confgList_up; + std::vector confgList_dn; + std::string optCI = "no"; + RealType cutoff = 0.0; + RealType zero_cutoff = 0.0; + OhmmsAttributeSet ciAttrib; + ciAttrib.add(optCI, "optimize"); + ciAttrib.add(optCI, "Optimize"); + ciAttrib.put(cur); + optimizeCI = (optCI == "yes"); + xmlNodePtr curRoot = cur, DetListNode = nullptr; + std::string cname, cname0; + cur = curRoot->children; + while (cur != NULL) //check the basis set + { + getNodeName(cname, cur); + if (cname == "detlist") + { + DetListNode = cur; + app_log() << "Found determinant list. \n"; + } + cur = cur->next; + } + size_t NCA = 0, NCB = 0; + size_t NEA = 0, NEB = 0; + size_t nstates = 0; + size_t ndets = 0; + size_t count = 0; + size_t cnt0 = 0; + std::string Dettype = "DETS"; + std::string CSFChoice = "qchem_coeff"; + OhmmsAttributeSet spoAttrib; + spoAttrib.add(NCA, "nca"); + spoAttrib.add(NCB, "ncb"); + spoAttrib.add(NEA, "nea"); + spoAttrib.add(NEB, "neb"); + spoAttrib.add(ndets, "size"); + spoAttrib.add(Dettype, "type"); + spoAttrib.add(nstates, "nstates"); + spoAttrib.add(cutoff, "cutoff"); + spoAttrib.add(zero_cutoff, "zero_cutoff"); + spoAttrib.add(zero_cutoff, "zerocutoff"); + spoAttrib.add(CSFChoice, "sortby"); + spoAttrib.put(DetListNode); + + if (ndets == 0) + { + APP_ABORT("size==0 in detlist is not allowed. Use slaterdeterminant in this case.\n"); + } -bool SlaterDetBuilder::readDetListH5(xmlNodePtr cur, - std::vector& uniqueConfg_up, - std::vector& uniqueConfg_dn, - std::vector& C2node_up, - std::vector& C2node_dn, - std::vector& CItags, - std::vector& coeff, - bool& optimizeCI, - int nels_up, - int nels_dn) + if (Dettype == "DETS" || Dettype == "Determinants") + usingCSF = false; + else if (Dettype == "CSF") + usingCSF = true; + else + { + APP_ABORT("Only allowed type in detlist is DETS or CSF.\n"); + } + + if (zero_cutoff > 0) + app_log() << " Initializing CI coeffs less than " << zero_cutoff << " to zero." << std::endl; + + if (NEA == 0) + NEA = nels_up - NCA; + else if (NEA != nels_up - NCA) + throw std::runtime_error("nea is not equal to nels_up - nca"); + if (NEB == 0) + NEB = nels_dn - NCB; + else if (NEB != nels_dn - NCB) + throw std::runtime_error("neb is not equal to nels_dn - ncb"); + + // mmorales: a little messy for now, clean up later + cur = DetListNode->children; + ci_configuration dummyC_alpha; + ci_configuration dummyC_beta; + dummyC_alpha.occup.resize(NCA + nstates, false); + for (size_t i = 0; i < NCA + NEA; i++) + dummyC_alpha.occup[i] = true; + dummyC_beta.occup.resize(NCB + nstates, false); + for (size_t i = 0; i < NCB + NEB; i++) + dummyC_beta.occup[i] = true; + RealType sumsq_qc = 0.0; + RealType sumsq = 0.0; + //app_log() <<"alpha reference: \n" < cutoff)) || ((CSFChoice == "coeff") && (std::abs(ci) < cutoff))) + { + cur = cur->next; + cnt0++; + continue; + } + cnt0++; + if (std::abs(qc_ci) < zero_cutoff) +#ifdef QMC_COMPLEX + ci_real = 0.0; +#else + ci = 0.0; +#endif + CSFcoeff.push_back(ci); + sumsq_qc += qc_ci * qc_ci; + DetsPerCSF.push_back(0); + CItags.push_back(tag); + count++; + xmlNodePtr csf = cur->children; + while (csf != NULL) + { + getNodeName(cname0, csf); + if (cname0 == "det") + { + std::string alpha, beta, tag0; + RealType coef = 0.0; + OhmmsAttributeSet detAttrib; + detAttrib.add(tag0, "id"); + detAttrib.add(coef, "coeff"); + detAttrib.add(beta, "beta"); + detAttrib.add(alpha, "alpha"); + detAttrib.put(csf); + size_t nq = 0; + if (alpha.size() < nstates) + { + std::cerr << "alpha: " << alpha << std::endl; + APP_ABORT("Found incorrect alpha determinant label. size < nca+nstates"); + } + for (size_t i = 0; i < nstates; i++) + { + if (alpha[i] != '0' && alpha[i] != '1') + { + std::cerr << alpha << std::endl; + APP_ABORT("Found incorrect determinant label."); + } + if (alpha[i] == '1') + nq++; + } + if (nq != NEA) + { + std::cerr << "alpha: " << alpha << std::endl; + APP_ABORT("Found incorrect alpha determinant label. noccup != nca+nea"); + } + nq = 0; + if (beta.size() < nstates) + { + std::cerr << "beta: " << beta << std::endl; + APP_ABORT("Found incorrect beta determinant label. size < ncb+nstates"); + } + for (size_t i = 0; i < nstates; i++) + { + if (beta[i] != '0' && beta[i] != '1') + { + std::cerr << beta << std::endl; + APP_ABORT("Found incorrect determinant label."); + } + if (beta[i] == '1') + nq++; + } + if (nq != NEB) + { + std::cerr << "beta: " << beta << std::endl; + APP_ABORT("Found incorrect beta determinant label. noccup != ncb+neb"); + } + //app_log() <<" " << std::endl; + DetsPerCSF.back()++; + CSFexpansion.push_back(coef); + coeff.push_back(coef * ci); + confgList_up.push_back(dummyC_alpha); + for (size_t i = 0; i < NCA; i++) + confgList_up.back().occup[i] = true; + for (size_t i = NCA; i < NCA + nstates; i++) + confgList_up.back().occup[i] = (alpha[i - NCA] == '1'); + confgList_dn.push_back(dummyC_beta); + for (size_t i = 0; i < NCB; i++) + confgList_dn.back().occup[i] = true; + for (size_t i = NCB; i < NCB + nstates; i++) + confgList_dn.back().occup[i] = (beta[i - NCB] == '1'); + } // if(name=="det") + csf = csf->next; + } // csf loop + if (DetsPerCSF.back() == 0) + { + APP_ABORT("Found empty CSF (no det blocks)."); + } + } // if (name == "csf") + cur = cur->next; + } + if (cnt0 != ndets) + { + std::cerr << "count, ndets: " << cnt0 << " " << ndets << std::endl; + APP_ABORT("Problems reading determinant ci_configurations. Found a number of determinants inconsistent with xml " + "file size parameter.\n"); + } + //if(!usingCSF) + // if(confgList_up.size() != ndets || confgList_dn.size() != ndets || coeff.size() != ndets) { + // APP_ABORT("Problems reading determinant ci_configurations."); + // } + C2node_up.resize(coeff.size()); + C2node_dn.resize(coeff.size()); + app_log() << "Found " << coeff.size() << " terms in the MSD expansion.\n"; + RealType sumsq = 0.0; + for (size_t i = 0; i < coeff.size(); i++) + sumsq += std::abs(coeff[i] * coeff[i]); + app_log() << "Norm of ci vector (sum of ci^2): " << sumsq << std::endl; + app_log() << "Norm of qchem ci vector (sum of qchem_ci^2): " << sumsq_qc << std::endl; + for (size_t i = 0; i < confgList_up.size(); i++) + { + bool found = false; + size_t k = -1; + for (size_t j = 0; j < uniqueConfg_up.size(); j++) + { + if (confgList_up[i] == uniqueConfg_up[j]) + { + found = true; + k = j; + break; + } + } + if (found) + { + C2node_up[i] = k; + } + else + { + uniqueConfg_up.push_back(confgList_up[i]); + C2node_up[i] = uniqueConfg_up.size() - 1; + } + } + for (size_t i = 0; i < confgList_dn.size(); i++) + { + bool found = false; + size_t k = -1; + for (size_t j = 0; j < uniqueConfg_dn.size(); j++) + { + if (confgList_dn[i] == uniqueConfg_dn[j]) + { + found = true; + k = j; + break; + } + } + if (found) + { + C2node_dn[i] = k; + } + else + { + uniqueConfg_dn.push_back(confgList_dn[i]); + C2node_dn[i] = uniqueConfg_dn.size() - 1; + } + } + } + else + { + app_log() << "Reading CI expansion." << std::endl; + + int cntup = 0; + int cntdn = 0; + std::unordered_map MyMapUp; + std::unordered_map MyMapDn; + while (cur != NULL) //check the basis set + { + getNodeName(cname, cur); + if (cname == "configuration" || cname == "ci") + { + RealType qc_ci = 0.0; + std::string alpha, beta, tag; + OhmmsAttributeSet confAttrib; +#ifdef QMC_COMPLEX + RealType ci_real = 0.0, ci_imag = 0.0; + confAttrib.add(ci_real, "coeff_real"); + confAttrib.add(ci_imag, "coeff_imag"); +#else + RealType ci = 0.0; + confAttrib.add(ci, "coeff"); +#endif + confAttrib.add(qc_ci, "qchem_coeff"); + confAttrib.add(alpha, "alpha"); + confAttrib.add(beta, "beta"); + confAttrib.add(tag, "id"); + confAttrib.put(cur); + +#ifdef QMC_COMPLEX + ValueType ci(ci_real, ci_imag); +#endif + + //Will always loop through the whole determinant set as no assumption on the order of the determinant is made + if (std::abs(ci) < cutoff) + { + cur = cur->next; + continue; + } + + for (size_t i = 0; i < nstates; i++) + { + if (alpha[i] != '0' && alpha[i] != '1') + { + std::cerr << alpha << std::endl; + APP_ABORT("Found incorrect determinant label."); + } + if (beta[i] != '0' && beta[i] != '1') + { + std::cerr << beta << std::endl; + APP_ABORT("Found incorrect determinant label."); + } + } + + if (alpha.size() < nstates) + { + std::cerr << "alpha: " << alpha << std::endl; + APP_ABORT("Found incorrect alpha determinant label. size < nca+nstates"); + } + if (beta.size() < nstates) + { + std::cerr << "beta: " << beta << std::endl; + APP_ABORT("Found incorrect beta determinant label. size < ncb+nstates"); + } + + coeff.push_back(ci); + CItags.push_back(tag); + + std::unordered_map::const_iterator gotup = MyMapUp.find(alpha); + + + if (gotup == MyMapUp.end()) + { + uniqueConfg_up.push_back(dummyC_alpha); + uniqueConfg_up.back().add_occupation(alpha); + C2node_up.push_back(cntup); + MyMapUp.insert(std::pair(alpha, cntup)); + cntup++; + } + else + { + C2node_up.push_back(gotup->second); + } + + std::unordered_map::const_iterator gotdn = MyMapDn.find(beta); + + if (gotdn == MyMapDn.end()) + { + uniqueConfg_dn.push_back(dummyC_beta); + uniqueConfg_dn.back().add_occupation(beta); + C2node_dn.push_back(cntdn); + MyMapDn.insert(std::pair(beta, cntdn)); + cntdn++; + } + else + { + C2node_dn.push_back(gotdn->second); + } + + //if (qc_ci == 0.0) + // qc_ci = ci; + + cnt0++; + sumsq_qc += qc_ci * qc_ci; + sumsq += std::abs(ci * ci); + } + cur = cur->next; + } + + app_log() << "Found " << coeff.size() << " terms in the MSD expansion.\n"; + app_log() << "Norm of ci vector (sum of ci^2): " << sumsq << std::endl; + app_log() << "Norm of qchem ci vector (sum of qchem_ci^2): " << sumsq_qc << std::endl; + + } //usingCSF + + app_log() << "Found " << uniqueConfg_up.size() << " unique up determinants.\n"; + app_log() << "Found " << uniqueConfg_dn.size() << " unique down determinants.\n"; + + return success; + */ +} + + +bool SlaterDetBuilder::readDetListH5(xmlNodePtr cur, + std::vector& uniqueConfg_up, + std::vector& uniqueConfg_dn, + std::vector& C2node_up, + std::vector& C2node_dn, + std::vector& CItags, + std::vector& coeff, + bool& optimizeCI, + int nels_up, + int nels_dn) +{ + bool success = true; + int extlevel(0); + uniqueConfg_up.clear(); + uniqueConfg_dn.clear(); + C2node_up.clear(); + C2node_dn.clear(); + CItags.clear(); + coeff.clear(); + std::string CICoeffH5path(""); + std::vector confgList_up; + std::vector confgList_dn; + std::vector CIcoeff; + std::vector ConfigTag; + std::string optCI = "no"; + RealType cutoff = 0.0; + OhmmsAttributeSet ciAttrib; + ciAttrib.add(optCI, "optimize"); + ciAttrib.put(cur); + optimizeCI = (optCI == "yes"); + xmlNodePtr curRoot = cur, DetListNode = nullptr; + std::string cname, cname0, multidetH5path; + cur = curRoot->children; + while (cur != NULL) //check the basis set + { + getNodeName(cname, cur); + if (cname == "detlist") + { + DetListNode = cur; + app_log() << "Found determinant list. \n"; + } + cur = cur->next; + } + + app_log() << " H5 code path implicitly assumes NCA = NCB = 0, " + << "NEA = " << nels_up << " and NEB = " << nels_dn << std::endl; + size_t nstates = 0; + size_t ndets = 0; + size_t H5_ndets, H5_nstates; + /// 64 bit fixed width integer + const unsigned bit_kind = 64; + static_assert(bit_kind == sizeof(int64_t) * 8, "Must be 64 bit fixed width integer"); + /// the number of 64 bit integers which represent the binary string for occupation + int N_int; + std::string Dettype = "DETS"; + ValueType sumsq = 0.0; + OhmmsAttributeSet spoAttrib; + spoAttrib.add(ndets, "size"); + spoAttrib.add(Dettype, "type"); + spoAttrib.add(nstates, "nstates"); + spoAttrib.add(extlevel, "ext_level"); + spoAttrib.add(cutoff, "cutoff"); + spoAttrib.add(multidetH5path, "href"); + spoAttrib.add(CICoeffH5path, "opt_coeffs"); + spoAttrib.put(DetListNode); + if (ndets == 0) + { + APP_ABORT("size==0 in detlist is not allowed. Use slaterdeterminant in this case.\n"); + } + if (Dettype != "DETS" && Dettype != "Determinants") + APP_ABORT("Reading from HDF5 is only enabled for CI DETS. Must be accessed through (type=\"DETS\") or " + "(type=\"Determinants\") .\n"); + app_log() << "Reading CI expansion from HDF5:" << multidetH5path << std::endl; + + hdf_archive hin; + if (!hin.open(multidetH5path.c_str(), H5F_ACC_RDONLY)) + { + std::cerr << "Could not open H5 file" << std::endl; + abort(); + } + + if (!hin.push("MultiDet")) + { + std::cerr << "Could not open Multidet Group in H5 file" << std::endl; + abort(); + } + + hin.read(H5_ndets, "NbDet"); + if (ndets != H5_ndets) + { + std::cerr << "Number of determinants in H5 file (" << H5_ndets << ") different from number of dets in XML (" + << ndets << ")" << std::endl; + abort(); + } + + hin.read(H5_nstates, "nstate"); + if (nstates == 0) + nstates = H5_nstates; + else if (nstates != H5_nstates) + { + std::cerr << "Number of states/orbitals in H5 file (" << H5_nstates + << ") different from number of states/orbitals in XML (" << nstates << ")" << std::endl; + abort(); + } + + hin.read(N_int, "Nbits"); + CIcoeff.resize(ndets); + ConfigTag.resize(ndets); + + readCoeffs(hin, CIcoeff, ndets, extlevel); + + ///IF OPTIMIZED COEFFICIENTS ARE PRESENT IN opt_coeffs Path + ///THEY ARE READ FROM DIFFERENT HDF5 the replace the previous coeff + ///It is important to still read all old coeffs and only replace the optimized ones + ///in order to keep coherence with the cutoff on the number of determinants + ///REMEMBER!! FIRST COEFF IS FIXED. THEREFORE WE DO NOT REPLACE IT!!! + if (CICoeffH5path != "") + { + int OptCiSize = 0; + std::vector CIcoeffopt; + hdf_archive coeffin; + if (!coeffin.open(CICoeffH5path.c_str(), H5F_ACC_RDONLY)) + { + std::cerr << "Could not open H5 file containing Optimized Coefficients" << std::endl; + abort(); + } + + if (!coeffin.push("MultiDet")) + { + std::cerr << "Could not open Multidet Group in H5 file" << std::endl; + abort(); + } + coeffin.read(OptCiSize, "NbDet"); + CIcoeffopt.resize(OptCiSize); + + readCoeffs(coeffin, CIcoeffopt, ndets, extlevel); + + coeffin.close(); + + for (int i = 0; i < OptCiSize; i++) + CIcoeff[i + 1] = CIcoeffopt[i]; + + app_log() << "The first " << OptCiSize + << " Optimized coefficients were substituted to the original set of coefficients." << std::endl; + } + + Matrix tempAlpha(ndets, N_int); + hin.read(tempAlpha, "CI_Alpha"); + + Matrix tempBeta(ndets, N_int); + hin.read(tempBeta, "CI_Beta"); + + std::string MyCIAlpha, MyCIBeta; + MyCIAlpha.resize(nstates); + MyCIBeta.resize(nstates); + + ci_configuration dummyC_alpha; + ci_configuration dummyC_beta; + + dummyC_alpha.occup.resize(nstates, false); + dummyC_beta.occup.resize(nstates, false); + + for (size_t i = 0; i < nstates; i++) + dummyC_alpha.occup[i] = true; + + for (size_t i = 0; i < nstates; i++) + dummyC_beta.occup[i] = true; + + hin.close(); + app_log() << " Done reading " << ndets << " CIs from H5!" << std::endl; + + int cntup = 0; + int cntdn = 0; + std::unordered_map MyMapUp; + std::unordered_map MyMapDn; + + app_log() << " Sorting unique CIs" << std::endl; + ///This loop will find all unique Determinants in and store them "unsorted" in a new container uniqueConfg_up + ///and uniqueConfg_dn. The sorting is not done here + for (int ni = 0; ni < ndets; ni++) + { + if (std::abs(CIcoeff[ni]) < cutoff) + continue; + + int j = 0; + for (int k = 0; k < N_int; k++) + { + int64_t a = tempAlpha[ni][k]; + std::bitset a2 = a; + + auto b = tempBeta[ni][k]; + auto b2 = std::bitset(b); + + for (int i = 0; i < bit_kind; i++) + { + if (j < nstates) + { + MyCIAlpha[j] = a2[i] ? '1' : '0'; + MyCIBeta[j] = b2[i] ? '1' : '0'; + j++; + } + } + } + coeff.push_back(CIcoeff[ni]); + std::ostringstream h5tag; + h5tag << "CIcoeff_" << ni; + CItags.push_back(h5tag.str()); + + std::unordered_map::const_iterator gotup = MyMapUp.find(MyCIAlpha); + + if (gotup == MyMapUp.end()) + { + uniqueConfg_up.push_back(dummyC_alpha); + uniqueConfg_up.back().add_occupation(MyCIAlpha); + C2node_up.push_back(cntup); + MyMapUp.insert(std::pair(MyCIAlpha, cntup)); + cntup++; + } + else + { + C2node_up.push_back(gotup->second); + } + + std::unordered_map::const_iterator gotdn = MyMapDn.find(MyCIBeta); + + if (gotdn == MyMapDn.end()) + { + uniqueConfg_dn.push_back(dummyC_beta); + uniqueConfg_dn.back().add_occupation(MyCIBeta); + C2node_dn.push_back(cntdn); + MyMapDn.insert(std::pair(MyCIBeta, cntdn)); + cntdn++; + } + else + { + C2node_dn.push_back(gotdn->second); + } + + sumsq += CIcoeff[ni] * CIcoeff[ni]; + } + + app_log() << " Done Sorting unique CIs" << std::endl; + app_log() << "Found " << coeff.size() << " terms in the MSD expansion.\n"; + app_log() << "Norm of ci vector (sum of ci^2): " << sumsq << std::endl; + + app_log() << "Found " << uniqueConfg_up.size() << " unique up determinants.\n"; + app_log() << "Found " << uniqueConfg_dn.size() << " unique down determinants.\n"; + + + return success; +} + +bool SlaterDetBuilder::readDetListH5(xmlNodePtr cur, + std::vector>& uniqueConfgs, + std::vector>& C2nodes, + std::vector& CItags, + std::vector& coeff, + bool& optimizeCI, + std::vector& nptcls) { + APP_ABORT("readDetListH5 not implemetned for arbitrary numbers of particle groups."); + /* bool success = true; int extlevel(0); uniqueConfg_up.clear(); @@ -1549,5 +2234,6 @@ bool SlaterDetBuilder::readDetListH5(xmlNodePtr cur, return success; + */ } } // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.h b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.h index 42b13c09bc..d5c87586c0 100644 --- a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.h +++ b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.h @@ -104,6 +104,17 @@ class SlaterDetBuilder : public WaveFunctionComponentBuilder std::vector& DetsPerCSF, std::vector& CSFexpansion, bool& usingCSF); + bool readDetList(xmlNodePtr cur, + std::vector>& uniqueConfgs, + std::vector>& C2nodes, + std::vector& CItags, + std::vector& coeff, + bool& optimizeCI, + std::vector& nptcls, + std::vector& CSFcoeff, + std::vector& DetsPerCSF, + std::vector& CSFexpansion, + bool& usingCSF); bool readDetListH5(xmlNodePtr cur, @@ -117,6 +128,14 @@ class SlaterDetBuilder : public WaveFunctionComponentBuilder int nels_up, int nels_dn); + bool readDetListH5(xmlNodePtr cur, + std::vector>& uniqueConfgs, + std::vector>& C2nodes, + std::vector& CItags, + std::vector& coeff, + bool& optimizeCI, + std::vector& nptcls); + // clang-format off template::value) && From 469aa988fd5bee2d8a7614392e77b987dacc0f3e Mon Sep 17 00:00:00 2001 From: camelto2 Date: Wed, 17 Feb 2021 17:05:27 -0700 Subject: [PATCH 04/16] readDetList generalized --- .../Fermion/SlaterDetBuilder.cpp | 308 +++++++----------- 1 file changed, 118 insertions(+), 190 deletions(-) diff --git a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp index 3b99880abd..56fc7d16d1 100644 --- a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp +++ b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp @@ -599,8 +599,6 @@ bool SlaterDetBuilder::createMSDFast(std::vector> uniqueConfgs; - std::vector CItags; bool optimizeCI; @@ -610,6 +608,9 @@ bool SlaterDetBuilder::createMSDFast(std::vector> uniqueConfgs(nGroups); + std::vector CItags; + //Check id multideterminants are in HDF5 xmlNodePtr curTemp = cur, DetListNode = nullptr; @@ -1313,20 +1314,20 @@ bool SlaterDetBuilder::readDetList(xmlNodePtr cur, std::vector& CSFexpansion, bool& usingCSF) { - APP_ABORT("readDetList for arbitrary ptcls groups not yet implemented.\n"); - /* bool success = true; - uniqueConfg_up.clear(); - uniqueConfg_dn.clear(); - C2node_up.clear(); - C2node_dn.clear(); + + int nGroups = uniqueConfgs.size(); + for (int grp = 0; grp < nGroups; grp++) + { + uniqueConfgs[grp].clear(); + C2nodes[grp].clear(); + } CItags.clear(); coeff.clear(); CSFcoeff.clear(); DetsPerCSF.clear(); CSFexpansion.clear(); - std::vector confgList_up; - std::vector confgList_dn; + std::vector> confgLists; std::string optCI = "no"; RealType cutoff = 0.0; RealType zero_cutoff = 0.0; @@ -1348,22 +1349,24 @@ bool SlaterDetBuilder::readDetList(xmlNodePtr cur, } cur = cur->next; } - size_t NCA = 0, NCB = 0; - size_t NEA = 0, NEB = 0; - size_t nstates = 0; + + std::vector NCs(nGroups); + std::vector NEs(nGroups); + std::vector nstates(nGroups); size_t ndets = 0; size_t count = 0; size_t cnt0 = 0; std::string Dettype = "DETS"; std::string CSFChoice = "qchem_coeff"; OhmmsAttributeSet spoAttrib; - spoAttrib.add(NCA, "nca"); - spoAttrib.add(NCB, "ncb"); - spoAttrib.add(NEA, "nea"); - spoAttrib.add(NEB, "neb"); + for (int grp = 0; grp < nGroups; grp++) + { + spoAttrib.add(NCs[grp], "nc" + std::to_string(grp)); + spoAttrib.add(NEs[grp], "ne" + std::to_string(grp)); + spoAttrib.add(nstates[grp], "nstates" + std::to_string(grp)); + } spoAttrib.add(ndets, "size"); spoAttrib.add(Dettype, "type"); - spoAttrib.add(nstates, "nstates"); spoAttrib.add(cutoff, "cutoff"); spoAttrib.add(zero_cutoff, "zero_cutoff"); spoAttrib.add(zero_cutoff, "zerocutoff"); @@ -1387,29 +1390,24 @@ bool SlaterDetBuilder::readDetList(xmlNodePtr cur, if (zero_cutoff > 0) app_log() << " Initializing CI coeffs less than " << zero_cutoff << " to zero." << std::endl; - if (NEA == 0) - NEA = nels_up - NCA; - else if (NEA != nels_up - NCA) - throw std::runtime_error("nea is not equal to nels_up - nca"); - if (NEB == 0) - NEB = nels_dn - NCB; - else if (NEB != nels_dn - NCB) - throw std::runtime_error("neb is not equal to nels_dn - ncb"); + for (int grp = 0; grp < nGroups; grp++) + { + if (NEs[grp] == 0) + NEs[grp] = nptcls[grp] - NCs[grp]; + else if (NEs[grp] != nptcls[grp] - NCs[grp]) + throw std::runtime_error("ne is not equal to n - nc for group " + std::to_string(grp)); + } - // mmorales: a little messy for now, clean up later cur = DetListNode->children; - ci_configuration dummyC_alpha; - ci_configuration dummyC_beta; - dummyC_alpha.occup.resize(NCA + nstates, false); - for (size_t i = 0; i < NCA + NEA; i++) - dummyC_alpha.occup[i] = true; - dummyC_beta.occup.resize(NCB + nstates, false); - for (size_t i = 0; i < NCB + NEB; i++) - dummyC_beta.occup[i] = true; + std::vector dummyCs(nGroups); + for (int grp = 0; grp < nGroups; grp++) + { + dummyCs[grp].occup.resize(NCs[grp] + nstates[grp], false); + for (size_t i = 0; i < NCs[grp] + NEs[grp]; i++) + dummyCs[grp].occup[i] = true; + } RealType sumsq_qc = 0.0; RealType sumsq = 0.0; - //app_log() <<"alpha reference: \n" < occs(nGroups); + std::string tag0; RealType coef = 0.0; OhmmsAttributeSet detAttrib; detAttrib.add(tag0, "id"); detAttrib.add(coef, "coeff"); - detAttrib.add(beta, "beta"); - detAttrib.add(alpha, "alpha"); + for (int grp = 0; grp < nGroups; grp++) + detAttrib.add(occs[grp], "occ" + std::to_string(grp)); detAttrib.put(csf); - size_t nq = 0; - if (alpha.size() < nstates) + for (int grp = 0; grp < nGroups; grp++) { - std::cerr << "alpha: " << alpha << std::endl; - APP_ABORT("Found incorrect alpha determinant label. size < nca+nstates"); - } - for (size_t i = 0; i < nstates; i++) - { - if (alpha[i] != '0' && alpha[i] != '1') + size_t nq = 0; + if (occs[grp].size() < nstates[grp]) { - std::cerr << alpha << std::endl; - APP_ABORT("Found incorrect determinant label."); + std::cerr << "occ" << grp << ": " << occs[grp] << std::endl; + APP_ABORT("Found incorrect group" + std::to_string(grp) + " determinant label. size < nc+nstates"); } - if (alpha[i] == '1') - nq++; - } - if (nq != NEA) - { - std::cerr << "alpha: " << alpha << std::endl; - APP_ABORT("Found incorrect alpha determinant label. noccup != nca+nea"); - } - nq = 0; - if (beta.size() < nstates) - { - std::cerr << "beta: " << beta << std::endl; - APP_ABORT("Found incorrect beta determinant label. size < ncb+nstates"); - } - for (size_t i = 0; i < nstates; i++) - { - if (beta[i] != '0' && beta[i] != '1') + for (size_t i = 0; i < nstates[grp]; i++) { - std::cerr << beta << std::endl; - APP_ABORT("Found incorrect determinant label."); + if (occs[grp][i] != '0' && occs[grp][i] != '1') + { + std::cerr << occs[grp] << std::endl; + APP_ABORT("Found incorrect determinant label."); + } + if (occs[grp][i] == '1') + nq++; + } + if (nq != NEs[grp]) + { + std::cerr << "occ" << grp << ": " << occs[grp] << std::endl; + APP_ABORT("Found incorrect group" + std::to_string(grp) + " determinant label. nocc != nc+ne"); } - if (beta[i] == '1') - nq++; - } - if (nq != NEB) - { - std::cerr << "beta: " << beta << std::endl; - APP_ABORT("Found incorrect beta determinant label. noccup != ncb+neb"); } - //app_log() <<" " << std::endl; DetsPerCSF.back()++; CSFexpansion.push_back(coef); coeff.push_back(coef * ci); - confgList_up.push_back(dummyC_alpha); - for (size_t i = 0; i < NCA; i++) - confgList_up.back().occup[i] = true; - for (size_t i = NCA; i < NCA + nstates; i++) - confgList_up.back().occup[i] = (alpha[i - NCA] == '1'); - confgList_dn.push_back(dummyC_beta); - for (size_t i = 0; i < NCB; i++) - confgList_dn.back().occup[i] = true; - for (size_t i = NCB; i < NCB + nstates; i++) - confgList_dn.back().occup[i] = (beta[i - NCB] == '1'); + for (int grp = 0; grp < nGroups; grp++) + { + confgLists[grp].push_back(dummyCs[grp]); + for (size_t i = 0; i < NCs[grp]; i++) + confgLists[grp].back().occup[i] = true; + for (size_t i = NCs[grp]; i < NCs[grp] + nstates[grp]; i++) + confgLists[grp].back().occup[i] = (occs[grp][i - NCs[grp]] == '1'); + } } // if(name=="det") csf = csf->next; } // csf loop @@ -1550,62 +1528,39 @@ bool SlaterDetBuilder::readDetList(xmlNodePtr cur, APP_ABORT("Problems reading determinant ci_configurations. Found a number of determinants inconsistent with xml " "file size parameter.\n"); } - //if(!usingCSF) - // if(confgList_up.size() != ndets || confgList_dn.size() != ndets || coeff.size() != ndets) { - // APP_ABORT("Problems reading determinant ci_configurations."); - // } - C2node_up.resize(coeff.size()); - C2node_dn.resize(coeff.size()); + + for (int grp = 0; grp < nGroups; grp++) + C2nodes[grp].resize(coeff.size()); app_log() << "Found " << coeff.size() << " terms in the MSD expansion.\n"; RealType sumsq = 0.0; for (size_t i = 0; i < coeff.size(); i++) sumsq += std::abs(coeff[i] * coeff[i]); app_log() << "Norm of ci vector (sum of ci^2): " << sumsq << std::endl; app_log() << "Norm of qchem ci vector (sum of qchem_ci^2): " << sumsq_qc << std::endl; - for (size_t i = 0; i < confgList_up.size(); i++) + for (int grp = 0; grp < nGroups; grp++) { - bool found = false; - size_t k = -1; - for (size_t j = 0; j < uniqueConfg_up.size(); j++) + for (size_t i = 0; i < confgLists[grp].size(); i++) { - if (confgList_up[i] == uniqueConfg_up[j]) + bool found = false; + size_t k = -1; + for (size_t j = 0; j < uniqueConfgs[grp].size(); j++) { - found = true; - k = j; - break; + if (confgLists[grp][i] == uniqueConfgs[grp][j]) + { + found = true; + k = j; + break; + } } - } - if (found) - { - C2node_up[i] = k; - } - else - { - uniqueConfg_up.push_back(confgList_up[i]); - C2node_up[i] = uniqueConfg_up.size() - 1; - } - } - for (size_t i = 0; i < confgList_dn.size(); i++) - { - bool found = false; - size_t k = -1; - for (size_t j = 0; j < uniqueConfg_dn.size(); j++) - { - if (confgList_dn[i] == uniqueConfg_dn[j]) + if (found) { - found = true; - k = j; - break; + C2nodes[grp][i] = k; + } + else + { + uniqueConfgs[grp].push_back(confgLists[grp][i]); + C2nodes[grp][i] = uniqueConfgs[grp].size() - 1; } - } - if (found) - { - C2node_dn[i] = k; - } - else - { - uniqueConfg_dn.push_back(confgList_dn[i]); - C2node_dn[i] = uniqueConfg_dn.size() - 1; } } } @@ -1613,17 +1568,16 @@ bool SlaterDetBuilder::readDetList(xmlNodePtr cur, { app_log() << "Reading CI expansion." << std::endl; - int cntup = 0; - int cntdn = 0; - std::unordered_map MyMapUp; - std::unordered_map MyMapDn; + std::vector cnts(nGroups, 0); + std::vector> MyMaps(nGroups); while (cur != NULL) //check the basis set { getNodeName(cname, cur); if (cname == "configuration" || cname == "ci") { RealType qc_ci = 0.0; - std::string alpha, beta, tag; + std::vector occs(nGroups); + std::string tag; OhmmsAttributeSet confAttrib; #ifdef QMC_COMPLEX RealType ci_real = 0.0, ci_imag = 0.0; @@ -1634,8 +1588,8 @@ bool SlaterDetBuilder::readDetList(xmlNodePtr cur, confAttrib.add(ci, "coeff"); #endif confAttrib.add(qc_ci, "qchem_coeff"); - confAttrib.add(alpha, "alpha"); - confAttrib.add(beta, "beta"); + for (int grp = 0; grp < nGroups; grp++) + confAttrib.add(occs[grp], "occ" + std::to_string(grp)); confAttrib.add(tag, "id"); confAttrib.put(cur); @@ -1650,68 +1604,43 @@ bool SlaterDetBuilder::readDetList(xmlNodePtr cur, continue; } - for (size_t i = 0; i < nstates; i++) + for (int grp = 0; grp < nGroups; grp++) { - if (alpha[i] != '0' && alpha[i] != '1') + for (size_t i = 0; i < nstates[grp]; i++) { - std::cerr << alpha << std::endl; - APP_ABORT("Found incorrect determinant label."); + if (occs[grp][i] != '0' && occs[grp][i] != '1') + { + std::cerr << occs[grp] << std::endl; + APP_ABORT("Found incorrect determinant label for group " + std::to_string(grp)); + } } - if (beta[i] != '0' && beta[i] != '1') + if (occs[grp].size() < nstates[grp]) { - std::cerr << beta << std::endl; - APP_ABORT("Found incorrect determinant label."); + std::cerr << "occ" << grp << ": " << occs[grp] << std::endl; + APP_ABORT("Found incorrect group" + std::to_string(grp) + " determinant label. size < nc+nstates"); } } - if (alpha.size() < nstates) - { - std::cerr << "alpha: " << alpha << std::endl; - APP_ABORT("Found incorrect alpha determinant label. size < nca+nstates"); - } - if (beta.size() < nstates) - { - std::cerr << "beta: " << beta << std::endl; - APP_ABORT("Found incorrect beta determinant label. size < ncb+nstates"); - } - coeff.push_back(ci); CItags.push_back(tag); - std::unordered_map::const_iterator gotup = MyMapUp.find(alpha); - - - if (gotup == MyMapUp.end()) - { - uniqueConfg_up.push_back(dummyC_alpha); - uniqueConfg_up.back().add_occupation(alpha); - C2node_up.push_back(cntup); - MyMapUp.insert(std::pair(alpha, cntup)); - cntup++; - } - else - { - C2node_up.push_back(gotup->second); - } - - std::unordered_map::const_iterator gotdn = MyMapDn.find(beta); - - if (gotdn == MyMapDn.end()) - { - uniqueConfg_dn.push_back(dummyC_beta); - uniqueConfg_dn.back().add_occupation(beta); - C2node_dn.push_back(cntdn); - MyMapDn.insert(std::pair(beta, cntdn)); - cntdn++; - } - else + for (int grp = 0; grp < nGroups; grp++) { - C2node_dn.push_back(gotdn->second); + std::unordered_map::const_iterator got = MyMaps[grp].find(occs[grp]); + if (got == MyMaps[grp].end()) + { + uniqueConfgs[grp].push_back(dummyCs[grp]); + uniqueConfgs[grp].back().add_occupation(occs[grp]); + C2nodes[grp].push_back(cnts[grp]); + MyMaps[grp].insert(std::pair(occs[grp], cnts[grp])); + cnts[grp]++; + } + else + { + C2nodes[grp].push_back(got->second); + } } - //if (qc_ci == 0.0) - // qc_ci = ci; - cnt0++; sumsq_qc += qc_ci * qc_ci; sumsq += std::abs(ci * ci); @@ -1725,11 +1654,10 @@ bool SlaterDetBuilder::readDetList(xmlNodePtr cur, } //usingCSF - app_log() << "Found " << uniqueConfg_up.size() << " unique up determinants.\n"; - app_log() << "Found " << uniqueConfg_dn.size() << " unique down determinants.\n"; + for (int grp = 0; grp < nGroups; grp++) + app_log() << "Found " << uniqueConfgs[grp].size() << " unique group" << grp << " determinants.\n"; return success; - */ } From 70874a0dc7fd47cdc7eb4b68fb9e6f3250f509b6 Mon Sep 17 00:00:00 2001 From: camelto2 Date: Wed, 17 Feb 2021 17:37:18 -0700 Subject: [PATCH 05/16] readDetListH5 --- .../Fermion/SlaterDetBuilder.cpp | 141 ++++++++---------- 1 file changed, 61 insertions(+), 80 deletions(-) diff --git a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp index 56fc7d16d1..a0d7ad1033 100644 --- a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp +++ b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp @@ -1352,7 +1352,7 @@ bool SlaterDetBuilder::readDetList(xmlNodePtr cur, std::vector NCs(nGroups); std::vector NEs(nGroups); - std::vector nstates(nGroups); + size_t nstates = 0; size_t ndets = 0; size_t count = 0; size_t cnt0 = 0; @@ -1363,9 +1363,9 @@ bool SlaterDetBuilder::readDetList(xmlNodePtr cur, { spoAttrib.add(NCs[grp], "nc" + std::to_string(grp)); spoAttrib.add(NEs[grp], "ne" + std::to_string(grp)); - spoAttrib.add(nstates[grp], "nstates" + std::to_string(grp)); } spoAttrib.add(ndets, "size"); + spoAttrib.add(nstates, "nstates"); spoAttrib.add(Dettype, "type"); spoAttrib.add(cutoff, "cutoff"); spoAttrib.add(zero_cutoff, "zero_cutoff"); @@ -1402,7 +1402,7 @@ bool SlaterDetBuilder::readDetList(xmlNodePtr cur, std::vector dummyCs(nGroups); for (int grp = 0; grp < nGroups; grp++) { - dummyCs[grp].occup.resize(NCs[grp] + nstates[grp], false); + dummyCs[grp].occup.resize(NCs[grp] + nstates, false); for (size_t i = 0; i < NCs[grp] + NEs[grp]; i++) dummyCs[grp].occup[i] = true; } @@ -1480,12 +1480,12 @@ bool SlaterDetBuilder::readDetList(xmlNodePtr cur, for (int grp = 0; grp < nGroups; grp++) { size_t nq = 0; - if (occs[grp].size() < nstates[grp]) + if (occs[grp].size() < nstates) { std::cerr << "occ" << grp << ": " << occs[grp] << std::endl; APP_ABORT("Found incorrect group" + std::to_string(grp) + " determinant label. size < nc+nstates"); } - for (size_t i = 0; i < nstates[grp]; i++) + for (size_t i = 0; i < nstates; i++) { if (occs[grp][i] != '0' && occs[grp][i] != '1') { @@ -1509,7 +1509,7 @@ bool SlaterDetBuilder::readDetList(xmlNodePtr cur, confgLists[grp].push_back(dummyCs[grp]); for (size_t i = 0; i < NCs[grp]; i++) confgLists[grp].back().occup[i] = true; - for (size_t i = NCs[grp]; i < NCs[grp] + nstates[grp]; i++) + for (size_t i = NCs[grp]; i < NCs[grp] + nstates; i++) confgLists[grp].back().occup[i] = (occs[grp][i - NCs[grp]] == '1'); } } // if(name=="det") @@ -1606,7 +1606,7 @@ bool SlaterDetBuilder::readDetList(xmlNodePtr cur, for (int grp = 0; grp < nGroups; grp++) { - for (size_t i = 0; i < nstates[grp]; i++) + for (size_t i = 0; i < nstates; i++) { if (occs[grp][i] != '0' && occs[grp][i] != '1') { @@ -1614,7 +1614,7 @@ bool SlaterDetBuilder::readDetList(xmlNodePtr cur, APP_ABORT("Found incorrect determinant label for group " + std::to_string(grp)); } } - if (occs[grp].size() < nstates[grp]) + if (occs[grp].size() < nstates) { std::cerr << "occ" << grp << ": " << occs[grp] << std::endl; APP_ABORT("Found incorrect group" + std::to_string(grp) + " determinant label. size < nc+nstates"); @@ -1921,19 +1921,18 @@ bool SlaterDetBuilder::readDetListH5(xmlNodePtr cur, bool& optimizeCI, std::vector& nptcls) { - APP_ABORT("readDetListH5 not implemetned for arbitrary numbers of particle groups."); - /* bool success = true; int extlevel(0); - uniqueConfg_up.clear(); - uniqueConfg_dn.clear(); - C2node_up.clear(); - C2node_dn.clear(); + int nGroups = uniqueConfgs.size(); + for (int grp = 0; grp < nGroups; grp++) + { + uniqueConfgs[grp].clear(); + C2nodes[grp].clear(); + } CItags.clear(); coeff.clear(); std::string CICoeffH5path(""); - std::vector confgList_up; - std::vector confgList_dn; + std::vector> confgLists(nGroups); std::vector CIcoeff; std::vector ConfigTag; std::string optCI = "no"; @@ -1956,8 +1955,10 @@ bool SlaterDetBuilder::readDetListH5(xmlNodePtr cur, cur = cur->next; } - app_log() << " H5 code path implicitly assumes NCA = NCB = 0, " - << "NEA = " << nels_up << " and NEB = " << nels_dn << std::endl; + app_log() << " H5 code path implicitly assumes NC0 = NC1 = ... = 0" << std::endl; + for (int grp = 0; grp < nGroups; grp++) + app_log() << "NE" << grp << " = " << nptcls[grp] << ", "; + app_log() << std::endl; size_t nstates = 0; size_t ndets = 0; size_t H5_ndets, H5_nstates; @@ -2058,35 +2059,29 @@ bool SlaterDetBuilder::readDetListH5(xmlNodePtr cur, << " Optimized coefficients were substituted to the original set of coefficients." << std::endl; } - Matrix tempAlpha(ndets, N_int); - hin.read(tempAlpha, "CI_Alpha"); - - Matrix tempBeta(ndets, N_int); - hin.read(tempBeta, "CI_Beta"); - - std::string MyCIAlpha, MyCIBeta; - MyCIAlpha.resize(nstates); - MyCIBeta.resize(nstates); - - ci_configuration dummyC_alpha; - ci_configuration dummyC_beta; - - dummyC_alpha.occup.resize(nstates, false); - dummyC_beta.occup.resize(nstates, false); - - for (size_t i = 0; i < nstates; i++) - dummyC_alpha.occup[i] = true; + std::vector> temps; + for (int grp = 0; grp < nGroups; grp++) + { + Matrix tmp(ndets, N_int); + temps.push_back(tmp); + hin.read(temps[grp], "CI_" + std::to_string(grp)); + } - for (size_t i = 0; i < nstates; i++) - dummyC_beta.occup[i] = true; + std::vector MyCIs(nGroups); + std::vector dummyCs(nGroups); + for (int grp = 0; grp < nGroups; grp++) + { + MyCIs[grp].resize(nstates); + dummyCs[grp].occup.resize(nstates, false); + for (size_t i = 0; i < nstates; i++) + dummyCs[grp].occup[i] = true; + } hin.close(); app_log() << " Done reading " << ndets << " CIs from H5!" << std::endl; - int cntup = 0; - int cntdn = 0; - std::unordered_map MyMapUp; - std::unordered_map MyMapDn; + std::vector cnts(nGroups, 0); + std::vector> MyMaps(nGroups); app_log() << " Sorting unique CIs" << std::endl; ///This loop will find all unique Determinants in and store them "unsorted" in a new container uniqueConfg_up @@ -2099,18 +2094,19 @@ bool SlaterDetBuilder::readDetListH5(xmlNodePtr cur, int j = 0; for (int k = 0; k < N_int; k++) { - int64_t a = tempAlpha[ni][k]; - std::bitset a2 = a; - - auto b = tempBeta[ni][k]; - auto b2 = std::bitset(b); + std::vector> a2s(nGroups); + for (int grp = 0; grp < nGroups; grp++) + { + int64_t a = temps[grp][ni][k]; + a2s[grp] = a; + } for (int i = 0; i < bit_kind; i++) { if (j < nstates) { - MyCIAlpha[j] = a2[i] ? '1' : '0'; - MyCIBeta[j] = b2[i] ? '1' : '0'; + for (int grp = 0; grp < nGroups; grp++) + MyCIs[grp][j] = a2s[grp][i] ? '1' : '0'; j++; } } @@ -2120,36 +2116,23 @@ bool SlaterDetBuilder::readDetListH5(xmlNodePtr cur, h5tag << "CIcoeff_" << ni; CItags.push_back(h5tag.str()); - std::unordered_map::const_iterator gotup = MyMapUp.find(MyCIAlpha); - - if (gotup == MyMapUp.end()) - { - uniqueConfg_up.push_back(dummyC_alpha); - uniqueConfg_up.back().add_occupation(MyCIAlpha); - C2node_up.push_back(cntup); - MyMapUp.insert(std::pair(MyCIAlpha, cntup)); - cntup++; - } - else + for (int grp = 0; grp < nGroups; grp++) { - C2node_up.push_back(gotup->second); - } + std::unordered_map::const_iterator got = MyMaps[grp].find(MyCIs[grp]); - std::unordered_map::const_iterator gotdn = MyMapDn.find(MyCIBeta); - - if (gotdn == MyMapDn.end()) - { - uniqueConfg_dn.push_back(dummyC_beta); - uniqueConfg_dn.back().add_occupation(MyCIBeta); - C2node_dn.push_back(cntdn); - MyMapDn.insert(std::pair(MyCIBeta, cntdn)); - cntdn++; - } - else - { - C2node_dn.push_back(gotdn->second); + if (got == MyMaps[grp].end()) + { + uniqueConfgs[grp].push_back(dummyCs[grp]); + uniqueConfgs[grp].back().add_occupation(MyCIs[grp]); + C2nodes[grp].push_back(cnts[grp]); + MyMaps[grp].insert(std::pair(MyCIs[grp], cnts[grp])); + cnts[grp]++; + } + else + { + C2nodes[grp].push_back(got->second); + } } - sumsq += CIcoeff[ni] * CIcoeff[ni]; } @@ -2157,11 +2140,9 @@ bool SlaterDetBuilder::readDetListH5(xmlNodePtr cur, app_log() << "Found " << coeff.size() << " terms in the MSD expansion.\n"; app_log() << "Norm of ci vector (sum of ci^2): " << sumsq << std::endl; - app_log() << "Found " << uniqueConfg_up.size() << " unique up determinants.\n"; - app_log() << "Found " << uniqueConfg_dn.size() << " unique down determinants.\n"; - + for (int grp = 0; grp < nGroups; grp++) + app_log() << "Found " << uniqueConfgs[grp].size() << " unique group" << grp << " determinants.\n"; return success; - */ } } // namespace qmcplusplus From f78c10e07204b50ef0f5a6ce3db1e0e10dc57a58 Mon Sep 17 00:00:00 2001 From: camelto2 Date: Thu, 18 Feb 2021 10:41:40 -0700 Subject: [PATCH 06/16] Properly guard against alpha,beta input style --- .../Fermion/SlaterDetBuilder.cpp | 47 ++++++++++++++----- .../Fermion/SlaterDetBuilder.h | 3 +- 2 files changed, 37 insertions(+), 13 deletions(-) diff --git a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp index a0d7ad1033..3c203a1c0e 100644 --- a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp +++ b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp @@ -231,9 +231,9 @@ WaveFunctionComponent* SlaterDetBuilder::buildComponent(xmlNodePtr cur) SPOSetPtr spo_beta; std::unique_ptr spo_alpha_clone, spo_beta_clone; //new format - std::vector spoNames(nGroups); - std::vector spos(nGroups, nullptr); - std::vector> spo_clones(nGroups); + std::vector spoNames; + std::vector spos; + std::vector> spo_clones; if (nGroups == -1) { spo_alpha = sposet_builder_factory_.getSPOSet(spo_alpha_name); @@ -259,6 +259,9 @@ WaveFunctionComponent* SlaterDetBuilder::buildComponent(xmlNodePtr cur) else { app_warning() << "Using new group method" << std::endl; + spoNames.resize(nGroups); + spos.resize(nGroups); + spo_clones.resize(nGroups); OhmmsAttributeSet grpAttrib; for (int grp = 0; grp < nGroups; grp++) grpAttrib.add(spoNames[grp], "spo_" + std::to_string(grp)); @@ -315,12 +318,21 @@ WaveFunctionComponent* SlaterDetBuilder::buildComponent(xmlNodePtr cur) multislaterdetfast_0 = new MultiSlaterDeterminantFast(targetPtcl, std::move(dets), false); } + bool ab_format = (nGroups == -1) ? true : false; + multislaterdetfast_0->initialize(); - success = createMSDFast(multislaterdetfast_0->Dets, *multislaterdetfast_0->C2node, *multislaterdetfast_0->C, - *multislaterdetfast_0->CSFcoeff, *multislaterdetfast_0->DetsPerCSF, - *multislaterdetfast_0->CSFexpansion, multislaterdetfast_0->usingCSF, - *multislaterdetfast_0->myVars, multislaterdetfast_0->Optimizable, - multislaterdetfast_0->CI_Optimizable, cur); + if (nGroups == -1) + success = createMSDFast(multislaterdetfast_0->Dets, *multislaterdetfast_0->C2node, *multislaterdetfast_0->C, + *multislaterdetfast_0->CSFcoeff, *multislaterdetfast_0->DetsPerCSF, + *multislaterdetfast_0->CSFexpansion, multislaterdetfast_0->usingCSF, + *multislaterdetfast_0->myVars, multislaterdetfast_0->Optimizable, + multislaterdetfast_0->CI_Optimizable, cur, ab_format); + else + success = createMSDFast(multislaterdetfast_0->Dets, *multislaterdetfast_0->C2node, *multislaterdetfast_0->C, + *multislaterdetfast_0->CSFcoeff, *multislaterdetfast_0->DetsPerCSF, + *multislaterdetfast_0->CSFexpansion, multislaterdetfast_0->usingCSF, + *multislaterdetfast_0->myVars, multislaterdetfast_0->Optimizable, + multislaterdetfast_0->CI_Optimizable, cur, ab_format); // The primary purpose of this function is to create all the optimizable orbital rotation parameters. // But if orbital rotation parameters were supplied by the user it will also apply a unitary transformation @@ -595,7 +607,8 @@ bool SlaterDetBuilder::createMSDFast(std::vector& uniqueConfg_up, From 009b044af4d249068882b90f4d0d09ed6441fb2c Mon Sep 17 00:00:00 2001 From: camelto2 Date: Thu, 18 Feb 2021 10:50:00 -0700 Subject: [PATCH 07/16] remove redundant if statement --- .../Fermion/SlaterDetBuilder.cpp | 17 +++++------------ 1 file changed, 5 insertions(+), 12 deletions(-) diff --git a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp index 3c203a1c0e..bf8ee31f24 100644 --- a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp +++ b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp @@ -321,18 +321,11 @@ WaveFunctionComponent* SlaterDetBuilder::buildComponent(xmlNodePtr cur) bool ab_format = (nGroups == -1) ? true : false; multislaterdetfast_0->initialize(); - if (nGroups == -1) - success = createMSDFast(multislaterdetfast_0->Dets, *multislaterdetfast_0->C2node, *multislaterdetfast_0->C, - *multislaterdetfast_0->CSFcoeff, *multislaterdetfast_0->DetsPerCSF, - *multislaterdetfast_0->CSFexpansion, multislaterdetfast_0->usingCSF, - *multislaterdetfast_0->myVars, multislaterdetfast_0->Optimizable, - multislaterdetfast_0->CI_Optimizable, cur, ab_format); - else - success = createMSDFast(multislaterdetfast_0->Dets, *multislaterdetfast_0->C2node, *multislaterdetfast_0->C, - *multislaterdetfast_0->CSFcoeff, *multislaterdetfast_0->DetsPerCSF, - *multislaterdetfast_0->CSFexpansion, multislaterdetfast_0->usingCSF, - *multislaterdetfast_0->myVars, multislaterdetfast_0->Optimizable, - multislaterdetfast_0->CI_Optimizable, cur, ab_format); + success = createMSDFast(multislaterdetfast_0->Dets, *multislaterdetfast_0->C2node, *multislaterdetfast_0->C, + *multislaterdetfast_0->CSFcoeff, *multislaterdetfast_0->DetsPerCSF, + *multislaterdetfast_0->CSFexpansion, multislaterdetfast_0->usingCSF, + *multislaterdetfast_0->myVars, multislaterdetfast_0->Optimizable, + multislaterdetfast_0->CI_Optimizable, cur, ab_format); // The primary purpose of this function is to create all the optimizable orbital rotation parameters. // But if orbital rotation parameters were supplied by the user it will also apply a unitary transformation From 7f776880b47ac26e8f21e8458482833c89fc3f4b Mon Sep 17 00:00:00 2001 From: camelto2 Date: Thu, 18 Feb 2021 14:48:10 -0700 Subject: [PATCH 08/16] remove nGroups reading from input. Remove separate paths from createMSDFast, allow backwards compatibility input --- .../Fermion/SlaterDetBuilder.cpp | 149 ++++++++---------- .../Fermion/SlaterDetBuilder.h | 4 +- 2 files changed, 65 insertions(+), 88 deletions(-) diff --git a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp index bf8ee31f24..98843e31af 100644 --- a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp +++ b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp @@ -204,16 +204,20 @@ WaveFunctionComponent* SlaterDetBuilder::buildComponent(xmlNodePtr cur) if (multislaterdet_0 || multislaterdetfast_0) APP_ABORT("multideterminant is already instantiated."); - std::string spo_alpha_name; - std::string spo_beta_name; + int nGroups = targetPtcl.groups(); + std::vector spoNames(nGroups); + std::string fastAlg; - int nGroups = -1; OhmmsAttributeSet spoAttrib; - spoAttrib.add(spo_alpha_name, "spo_up"); - spoAttrib.add(spo_beta_name, "spo_dn"); + for (int grp = 0; grp < nGroups; grp++) + spoAttrib.add(spoNames[grp], "spo_" + std::to_string(grp)); + if (nGroups == 2) + { //for legacy + spoAttrib.add(spoNames[0], "spo_up"); + spoAttrib.add(spoNames[1], "spo_dn"); + } spoAttrib.add(fastAlg, "Fast", {"", "yes", "no"}, TagStatus::DEPRECATED); spoAttrib.add(msd_algorithm, "algorithm", {"", "precomputed_table_method", "table_method", "all_determinants"}); - spoAttrib.add(nGroups, "ngroups"); spoAttrib.put(cur); if (msd_algorithm.empty()) @@ -226,62 +230,23 @@ WaveFunctionComponent* SlaterDetBuilder::buildComponent(xmlNodePtr cur) msd_algorithm = "all_determinants"; } - //old format - SPOSetPtr spo_alpha; - SPOSetPtr spo_beta; - std::unique_ptr spo_alpha_clone, spo_beta_clone; //new format - std::vector spoNames; - std::vector spos; - std::vector> spo_clones; - if (nGroups == -1) - { - spo_alpha = sposet_builder_factory_.getSPOSet(spo_alpha_name); - spo_beta = sposet_builder_factory_.getSPOSet(spo_beta_name); - if (spo_alpha == nullptr) - { - app_error() << "In SlaterDetBuilder: SPOSet \"" << spo_alpha_name - << "\" is not found. Expected for MultiSlaterDeterminant.\n"; - abort(); - } - else - spo_alpha_clone.reset(spo_alpha->makeClone()); + std::vector spos(nGroups); + std::vector> spo_clones(nGroups); - if (spo_beta == nullptr) + for (int grp = 0; grp < nGroups; grp++) + { + spos[grp] = sposet_builder_factory_.getSPOSet(spoNames[grp]); + if (spos[grp] == nullptr) { - app_error() << "In SlaterDetBuilder: SPOSet \"" << spo_beta_name + app_error() << "In SlaterDetBuilder: SPOSet \"" << spoNames[grp] << "\" is not found. Expected for MultiSlaterDeterminant.\n"; abort(); } - else - spo_beta_clone.reset(spo_beta->makeClone()); - } - else - { - app_warning() << "Using new group method" << std::endl; - spoNames.resize(nGroups); - spos.resize(nGroups); - spo_clones.resize(nGroups); - OhmmsAttributeSet grpAttrib; - for (int grp = 0; grp < nGroups; grp++) - grpAttrib.add(spoNames[grp], "spo_" + std::to_string(grp)); - grpAttrib.put(cur); - app_log() << "GROUPS" << std::endl; - for (int grp = 0; grp < nGroups; grp++) - app_log() << spoNames[grp] << std::endl; - for (int grp = 0; grp < nGroups; grp++) - { - spos[grp] = sposet_builder_factory_.getSPOSet(spoNames[grp]); - if (spos[grp] == nullptr) - { - app_error() << "In SlaterDetBuilder: SPOSet \"" << spoNames[grp] - << "\" is not found. Expected for MultiSlaterDeterminant.\n"; - abort(); - } - spo_clones[grp].reset(spos[grp]->makeClone()); - } + spo_clones[grp].reset(spos[grp]->makeClone()); } + if (msd_algorithm == "precomputed_table_method" || msd_algorithm == "table_method") { app_summary() << " Using Bryan's table method." << std::endl; @@ -291,20 +256,10 @@ WaveFunctionComponent* SlaterDetBuilder::buildComponent(xmlNodePtr cur) } std::vector> dets; - if (nGroups == -1) + for (int grp = 0; grp < nGroups; grp++) { - app_log() << " Creating base determinant (up) for MSD expansion. \n"; - dets.emplace_back(std::make_unique(std::move(spo_alpha_clone), 0)); - app_log() << " Creating base determinant (down) for MSD expansion. \n"; - dets.emplace_back(std::make_unique(std::move(spo_beta_clone), 1)); - } - else - { - for (int grp = 0; grp < nGroups; grp++) - { - app_log() << " Creating base determinant (" << grp << ") for MSD expansion. \n"; - dets.emplace_back(std::make_unique(std::move(spo_clones[grp]), grp)); - } + app_log() << " Creating base determinant (" << grp << ") for MSD expansion. \n"; + dets.emplace_back(std::make_unique(std::move(spo_clones[grp]), grp)); } if (msd_algorithm == "precomputed_table_method") @@ -318,14 +273,12 @@ WaveFunctionComponent* SlaterDetBuilder::buildComponent(xmlNodePtr cur) multislaterdetfast_0 = new MultiSlaterDeterminantFast(targetPtcl, std::move(dets), false); } - bool ab_format = (nGroups == -1) ? true : false; - multislaterdetfast_0->initialize(); success = createMSDFast(multislaterdetfast_0->Dets, *multislaterdetfast_0->C2node, *multislaterdetfast_0->C, *multislaterdetfast_0->CSFcoeff, *multislaterdetfast_0->DetsPerCSF, *multislaterdetfast_0->CSFexpansion, multislaterdetfast_0->usingCSF, *multislaterdetfast_0->myVars, multislaterdetfast_0->Optimizable, - multislaterdetfast_0->CI_Optimizable, cur, ab_format); + multislaterdetfast_0->CI_Optimizable, cur); // The primary purpose of this function is to create all the optimizable orbital rotation parameters. // But if orbital rotation parameters were supplied by the user it will also apply a unitary transformation @@ -334,11 +287,16 @@ WaveFunctionComponent* SlaterDetBuilder::buildComponent(xmlNodePtr cur) } else { + if (nGroups != 2) + { + PRE.error("MSD using all_determinants algorithm requires two particle species."); + return nullptr; + } app_summary() << " Using a list of determinants for multi-deterimant expansion." << std::endl; auto spo_up = - std::make_unique(std::move(spo_alpha_clone), targetPtcl.first(0), targetPtcl.last(0)); + std::make_unique(std::move(spo_clones[0]), targetPtcl.first(0), targetPtcl.last(0)); auto spo_dn = - std::make_unique(std::move(spo_beta_clone), targetPtcl.first(1), targetPtcl.last(1)); + std::make_unique(std::move(spo_clones[1]), targetPtcl.first(1), targetPtcl.last(1)); if (UseBackflow) { app_summary() << " Using backflow transformation." << std::endl; @@ -600,8 +558,7 @@ bool SlaterDetBuilder::createMSDFast(std::vector tmp(ndets, N_int); temps.push_back(tmp); - hin.read(temps[grp], "CI_" + std::to_string(grp)); + std::string ci_str; + + if (!hin.readEntry(temps[grp], "CI_" + std::to_string(grp))) + { + //for backwards compatibility + if (grp == 0) + hin.read(temps[grp], "CI_Alpha"); + else if (grp == 1) + hin.read(temps[grp], "CI_Beta"); + else + APP_ABORT("Unknown HDF5 CI format"); + } } std::vector MyCIs(nGroups); diff --git a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.h b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.h index 8847e2e6f1..47f71a9507 100644 --- a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.h +++ b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.h @@ -88,8 +88,8 @@ class SlaterDetBuilder : public WaveFunctionComponentBuilder opt_variables_type& myVars, bool& Optimizable, bool& CI_Optimizable, - xmlNodePtr cur, - bool& ab_format); + xmlNodePtr cur); + bool readDetList(xmlNodePtr cur, std::vector& uniqueConfg_up, From ed5aec6b305b26e570ac15c02cb4aff6e1bb11e3 Mon Sep 17 00:00:00 2001 From: camelto2 Date: Mon, 22 Feb 2021 11:31:08 -0700 Subject: [PATCH 09/16] add unit test for arbitrary number of MSD species --- .../test_spo_collection_input_MSD_LCAO_h5.cpp | 107 +++++++++++++++++- 1 file changed, 103 insertions(+), 4 deletions(-) diff --git a/src/QMCWaveFunctions/tests/test_spo_collection_input_MSD_LCAO_h5.cpp b/src/QMCWaveFunctions/tests/test_spo_collection_input_MSD_LCAO_h5.cpp index a006fcfb47..ded9ae5650 100644 --- a/src/QMCWaveFunctions/tests/test_spo_collection_input_MSD_LCAO_h5.cpp +++ b/src/QMCWaveFunctions/tests/test_spo_collection_input_MSD_LCAO_h5.cpp @@ -26,7 +26,10 @@ using std::string; namespace qmcplusplus { -void test_LiH_msd_xml_input(const std::string& spo_xml_string, const std::string& check_sponame, int check_spo_size, int check_basisset_size) +void test_LiH_msd_xml_input(const std::string& spo_xml_string, + const std::string& check_sponame, + int check_spo_size, + int check_basisset_size) { Communicate* c; c = OHMMS::Controller; @@ -37,9 +40,9 @@ void test_LiH_msd_xml_input(const std::string& spo_xml_string, const std::string ParticleSet& elec_(*elec_uptr); ions_.setName("ion0"); - ions_.create({1,1}); - ions_.R[0] = {0.0, 0.0, 0.0}; - ions_.R[1] = {0.0, 0.0, 3.0139239693}; + ions_.create({1, 1}); + ions_.R[0] = {0.0, 0.0, 0.0}; + ions_.R[1] = {0.0, 0.0, 3.0139239693}; SpeciesSet& ispecies = ions_.getSpeciesSet(); int LiIdx = ispecies.addSpecies("Li"); int HIdx = ispecies.addSpecies("H"); @@ -198,4 +201,100 @@ TEST_CASE("SPO input spline from xml LiH_msd", "[wavefunction]") "; test_LiH_msd_xml_input(spo_xml_string3, "spo-up", 85, 105); } + +void test_LiH_msd_xml_input_with_positron(const std::string& spo_xml_string, + const std::string& check_sponame, + int check_spo_size, + int check_basisset_size) +{ + Communicate* c; + c = OHMMS::Controller; + + auto ions_uptr = std::make_unique(); + auto elec_uptr = std::make_unique(); + ParticleSet& ions_(*ions_uptr); + ParticleSet& elec_(*elec_uptr); + + ions_.setName("ion0"); + ions_.create({1, 1}); + ions_.R[0] = {0.0, 0.0, 0.0}; + ions_.R[1] = {0.0, 0.0, 3.0139239693}; + SpeciesSet& ispecies = ions_.getSpeciesSet(); + int LiIdx = ispecies.addSpecies("Li"); + int HIdx = ispecies.addSpecies("H"); + + elec_.setName("elec"); + elec_.create({2, 2, 1}); + elec_.R[0] = {0.5, 0.5, 0.5}; + elec_.R[1] = {0.1, 0.1, 1.1}; + elec_.R[2] = {-0.5, -0.5, -0.5}; + elec_.R[3] = {-0.1, -0.1, 1.5}; + elec_.R[4] = {0.0, -1.0, 2}; + + SpeciesSet& tspecies = elec_.getSpeciesSet(); + int upIdx = tspecies.addSpecies("u"); + int downIdx = tspecies.addSpecies("d"); + int positronIdx = tspecies.addSpecies("pos"); + int massIdx = tspecies.addAttribute("mass"); + int chargeIdx = tspecies.addAttribute("charge"); + tspecies(massIdx, upIdx) = 1.0; + tspecies(massIdx, downIdx) = 1.0; + tspecies(massIdx, positronIdx) = 1.0; + tspecies(chargeIdx, upIdx) = -1.0; + tspecies(chargeIdx, downIdx) = -1.0; + tspecies(chargeIdx, positronIdx) = 1.0; + + // Need 1 electron and 1 proton, somehow + //ParticleSet target = ParticleSet(); + ParticleSetPool ptcl = ParticleSetPool(c); + ptcl.addParticleSet(std::move(elec_uptr)); + ptcl.addParticleSet(std::move(ions_uptr)); + + Libxml2Document doc; + bool okay = doc.parseFromString(spo_xml_string); + REQUIRE(okay); + + xmlNodePtr ein_xml = doc.getRoot(); + + WaveFunctionFactory wf_factory("psi0", elec_, ptcl.getPool(), c); + wf_factory.put(ein_xml); + + SPOSet* spo_ptr(wf_factory.getSPOSet(check_sponame)); + REQUIRE(spo_ptr != nullptr); + REQUIRE(spo_ptr->getOrbitalSetSize() == check_spo_size); + REQUIRE(spo_ptr->getBasisSetSize() == check_basisset_size); +} + +TEST_CASE("SPO input spline from xml LiH_msd arbitrary species", "[wavefunction]") +{ + app_log() << "-------------------------------------------------------------" << std::endl; + app_log() << "LiH_msd with positron xml input style" << std::endl; + app_log() << "-------------------------------------------------------------" << std::endl; + const char* spo_xml_string1 = " \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ +"; + test_LiH_msd_xml_input_with_positron(spo_xml_string1, "spo-ps", 5, 105); +} } // namespace qmcplusplus From acf14bac7adc8a401ee0d694c2bca089140edc7a Mon Sep 17 00:00:00 2001 From: camelto2 Date: Mon, 22 Feb 2021 12:46:32 -0700 Subject: [PATCH 10/16] Remove need for overloaded readDetList and readDetListH5 --- .../Fermion/SlaterDetBuilder.cpp | 718 +----------------- .../Fermion/SlaterDetBuilder.h | 28 +- 2 files changed, 19 insertions(+), 727 deletions(-) diff --git a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp index 98843e31af..f560ce3b43 100644 --- a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp +++ b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp @@ -698,11 +698,13 @@ bool SlaterDetBuilder::createMSDFast(std::vector uniqueConfg_up, uniqueConfg_dn; + std::vector> uniqueConfgs(2); std::vector CItags; bool optimizeCI; - int nels_up = multiSD->nels_up; - int nels_dn = multiSD->nels_dn; + std::vector nels(2); + nels[0] = multiSD->nels_up; + nels[1] = multiSD->nels_dn; + std::vector> C2nodes(2); //Check id multideterminants are in HDF5 @@ -722,23 +724,24 @@ bool SlaterDetBuilder::createMSD(MultiSlaterDeterminant* multiSD, xmlNodePtr cur if (HDF5Path != "") { app_log() << "Found Multideterminants in H5 File" << std::endl; - success = readDetListH5(cur, uniqueConfg_up, uniqueConfg_dn, multiSD->C2node_up, multiSD->C2node_dn, CItags, - multiSD->C, optimizeCI, nels_up, nels_dn); + success = readDetListH5(cur, uniqueConfgs, C2nodes, CItags, multiSD->C, optimizeCI, nels); } else - success = readDetList(cur, uniqueConfg_up, uniqueConfg_dn, multiSD->C2node_up, multiSD->C2node_dn, CItags, - multiSD->C, optimizeCI, nels_up, nels_dn, multiSD->CSFcoeff, multiSD->DetsPerCSF, - multiSD->CSFexpansion, multiSD->usingCSF); + success = readDetList(cur, uniqueConfgs, C2nodes, CItags, multiSD->C, optimizeCI, nels, multiSD->CSFcoeff, + multiSD->DetsPerCSF, multiSD->CSFexpansion, multiSD->usingCSF); if (!success) return false; - multiSD->resize(uniqueConfg_up.size(), uniqueConfg_dn.size()); + + multiSD->C2node_up = C2nodes[0]; + multiSD->C2node_dn = C2nodes[1]; + multiSD->resize(uniqueConfgs[0].size(), uniqueConfgs[1].size()); { auto& spo = multiSD->spo_up; - spo->occup.resize(uniqueConfg_up.size(), multiSD->nels_up); - for (int i = 0; i < uniqueConfg_up.size(); i++) + spo->occup.resize(uniqueConfgs[0].size(), multiSD->nels_up); + for (int i = 0; i < uniqueConfgs[0].size(); i++) { int nq = 0; - ci_configuration& ci = uniqueConfg_up[i]; + ci_configuration& ci = uniqueConfgs[0][i]; for (int k = 0; k < ci.occup.size(); k++) { if (ci.occup[k]) @@ -761,11 +764,11 @@ bool SlaterDetBuilder::createMSD(MultiSlaterDeterminant* multiSD, xmlNodePtr cur } { auto& spo = multiSD->spo_dn; - spo->occup.resize(uniqueConfg_dn.size(), multiSD->nels_dn); - for (int i = 0; i < uniqueConfg_dn.size(); i++) + spo->occup.resize(uniqueConfgs[1].size(), multiSD->nels_dn); + for (int i = 0; i < uniqueConfgs[1].size(); i++) { int nq = 0; - ci_configuration& ci = uniqueConfg_dn[i]; + ci_configuration& ci = uniqueConfgs[1][i]; for (int k = 0; k < ci.occup.size(); k++) { if (ci.occup[k]) @@ -835,438 +838,6 @@ bool SlaterDetBuilder::createMSD(MultiSlaterDeterminant* multiSD, xmlNodePtr cur return success; } - -bool SlaterDetBuilder::readDetList(xmlNodePtr cur, - std::vector& uniqueConfg_up, - std::vector& uniqueConfg_dn, - std::vector& C2node_up, - std::vector& C2node_dn, - std::vector& CItags, - std::vector& coeff, - bool& optimizeCI, - int nels_up, - int nels_dn, - std::vector& CSFcoeff, - std::vector& DetsPerCSF, - std::vector& CSFexpansion, - bool& usingCSF) -{ - bool success = true; - uniqueConfg_up.clear(); - uniqueConfg_dn.clear(); - C2node_up.clear(); - C2node_dn.clear(); - CItags.clear(); - coeff.clear(); - CSFcoeff.clear(); - DetsPerCSF.clear(); - CSFexpansion.clear(); - std::vector confgList_up; - std::vector confgList_dn; - std::string optCI = "no"; - RealType cutoff = 0.0; - RealType zero_cutoff = 0.0; - OhmmsAttributeSet ciAttrib; - ciAttrib.add(optCI, "optimize"); - ciAttrib.add(optCI, "Optimize"); - ciAttrib.put(cur); - optimizeCI = (optCI == "yes"); - xmlNodePtr curRoot = cur, DetListNode = nullptr; - std::string cname, cname0; - cur = curRoot->children; - while (cur != NULL) //check the basis set - { - getNodeName(cname, cur); - if (cname == "detlist") - { - DetListNode = cur; - app_log() << "Found determinant list. \n"; - } - cur = cur->next; - } - size_t NCA = 0, NCB = 0; - size_t NEA = 0, NEB = 0; - size_t nstates = 0; - size_t ndets = 0; - size_t count = 0; - size_t cnt0 = 0; - std::string Dettype = "DETS"; - std::string CSFChoice = "qchem_coeff"; - OhmmsAttributeSet spoAttrib; - spoAttrib.add(NCA, "nca"); - spoAttrib.add(NCB, "ncb"); - spoAttrib.add(NEA, "nea"); - spoAttrib.add(NEB, "neb"); - spoAttrib.add(ndets, "size"); - spoAttrib.add(Dettype, "type"); - spoAttrib.add(nstates, "nstates"); - spoAttrib.add(cutoff, "cutoff"); - spoAttrib.add(zero_cutoff, "zero_cutoff"); - spoAttrib.add(zero_cutoff, "zerocutoff"); - spoAttrib.add(CSFChoice, "sortby"); - spoAttrib.put(DetListNode); - - if (ndets == 0) - { - APP_ABORT("size==0 in detlist is not allowed. Use slaterdeterminant in this case.\n"); - } - - if (Dettype == "DETS" || Dettype == "Determinants") - usingCSF = false; - else if (Dettype == "CSF") - usingCSF = true; - else - { - APP_ABORT("Only allowed type in detlist is DETS or CSF.\n"); - } - - if (zero_cutoff > 0) - app_log() << " Initializing CI coeffs less than " << zero_cutoff << " to zero." << std::endl; - - if (NEA == 0) - NEA = nels_up - NCA; - else if (NEA != nels_up - NCA) - throw std::runtime_error("nea is not equal to nels_up - nca"); - if (NEB == 0) - NEB = nels_dn - NCB; - else if (NEB != nels_dn - NCB) - throw std::runtime_error("neb is not equal to nels_dn - ncb"); - - // mmorales: a little messy for now, clean up later - cur = DetListNode->children; - ci_configuration dummyC_alpha; - ci_configuration dummyC_beta; - dummyC_alpha.occup.resize(NCA + nstates, false); - for (size_t i = 0; i < NCA + NEA; i++) - dummyC_alpha.occup[i] = true; - dummyC_beta.occup.resize(NCB + nstates, false); - for (size_t i = 0; i < NCB + NEB; i++) - dummyC_beta.occup[i] = true; - RealType sumsq_qc = 0.0; - RealType sumsq = 0.0; - //app_log() <<"alpha reference: \n" < cutoff)) || ((CSFChoice == "coeff") && (std::abs(ci) < cutoff))) - { - cur = cur->next; - cnt0++; - continue; - } - cnt0++; - if (std::abs(qc_ci) < zero_cutoff) -#ifdef QMC_COMPLEX - ci_real = 0.0; -#else - ci = 0.0; -#endif - CSFcoeff.push_back(ci); - sumsq_qc += qc_ci * qc_ci; - DetsPerCSF.push_back(0); - CItags.push_back(tag); - count++; - xmlNodePtr csf = cur->children; - while (csf != NULL) - { - getNodeName(cname0, csf); - if (cname0 == "det") - { - std::string alpha, beta, tag0; - RealType coef = 0.0; - OhmmsAttributeSet detAttrib; - detAttrib.add(tag0, "id"); - detAttrib.add(coef, "coeff"); - detAttrib.add(beta, "beta"); - detAttrib.add(alpha, "alpha"); - detAttrib.put(csf); - size_t nq = 0; - if (alpha.size() < nstates) - { - std::cerr << "alpha: " << alpha << std::endl; - APP_ABORT("Found incorrect alpha determinant label. size < nca+nstates"); - } - for (size_t i = 0; i < nstates; i++) - { - if (alpha[i] != '0' && alpha[i] != '1') - { - std::cerr << alpha << std::endl; - APP_ABORT("Found incorrect determinant label."); - } - if (alpha[i] == '1') - nq++; - } - if (nq != NEA) - { - std::cerr << "alpha: " << alpha << std::endl; - APP_ABORT("Found incorrect alpha determinant label. noccup != nca+nea"); - } - nq = 0; - if (beta.size() < nstates) - { - std::cerr << "beta: " << beta << std::endl; - APP_ABORT("Found incorrect beta determinant label. size < ncb+nstates"); - } - for (size_t i = 0; i < nstates; i++) - { - if (beta[i] != '0' && beta[i] != '1') - { - std::cerr << beta << std::endl; - APP_ABORT("Found incorrect determinant label."); - } - if (beta[i] == '1') - nq++; - } - if (nq != NEB) - { - std::cerr << "beta: " << beta << std::endl; - APP_ABORT("Found incorrect beta determinant label. noccup != ncb+neb"); - } - //app_log() <<" " << std::endl; - DetsPerCSF.back()++; - CSFexpansion.push_back(coef); - coeff.push_back(coef * ci); - confgList_up.push_back(dummyC_alpha); - for (size_t i = 0; i < NCA; i++) - confgList_up.back().occup[i] = true; - for (size_t i = NCA; i < NCA + nstates; i++) - confgList_up.back().occup[i] = (alpha[i - NCA] == '1'); - confgList_dn.push_back(dummyC_beta); - for (size_t i = 0; i < NCB; i++) - confgList_dn.back().occup[i] = true; - for (size_t i = NCB; i < NCB + nstates; i++) - confgList_dn.back().occup[i] = (beta[i - NCB] == '1'); - } // if(name=="det") - csf = csf->next; - } // csf loop - if (DetsPerCSF.back() == 0) - { - APP_ABORT("Found empty CSF (no det blocks)."); - } - } // if (name == "csf") - cur = cur->next; - } - if (cnt0 != ndets) - { - std::cerr << "count, ndets: " << cnt0 << " " << ndets << std::endl; - APP_ABORT("Problems reading determinant ci_configurations. Found a number of determinants inconsistent with xml " - "file size parameter.\n"); - } - //if(!usingCSF) - // if(confgList_up.size() != ndets || confgList_dn.size() != ndets || coeff.size() != ndets) { - // APP_ABORT("Problems reading determinant ci_configurations."); - // } - C2node_up.resize(coeff.size()); - C2node_dn.resize(coeff.size()); - app_log() << "Found " << coeff.size() << " terms in the MSD expansion.\n"; - RealType sumsq = 0.0; - for (size_t i = 0; i < coeff.size(); i++) - sumsq += std::abs(coeff[i] * coeff[i]); - app_log() << "Norm of ci vector (sum of ci^2): " << sumsq << std::endl; - app_log() << "Norm of qchem ci vector (sum of qchem_ci^2): " << sumsq_qc << std::endl; - for (size_t i = 0; i < confgList_up.size(); i++) - { - bool found = false; - size_t k = -1; - for (size_t j = 0; j < uniqueConfg_up.size(); j++) - { - if (confgList_up[i] == uniqueConfg_up[j]) - { - found = true; - k = j; - break; - } - } - if (found) - { - C2node_up[i] = k; - } - else - { - uniqueConfg_up.push_back(confgList_up[i]); - C2node_up[i] = uniqueConfg_up.size() - 1; - } - } - for (size_t i = 0; i < confgList_dn.size(); i++) - { - bool found = false; - size_t k = -1; - for (size_t j = 0; j < uniqueConfg_dn.size(); j++) - { - if (confgList_dn[i] == uniqueConfg_dn[j]) - { - found = true; - k = j; - break; - } - } - if (found) - { - C2node_dn[i] = k; - } - else - { - uniqueConfg_dn.push_back(confgList_dn[i]); - C2node_dn[i] = uniqueConfg_dn.size() - 1; - } - } - } - else - { - app_log() << "Reading CI expansion." << std::endl; - - int cntup = 0; - int cntdn = 0; - std::unordered_map MyMapUp; - std::unordered_map MyMapDn; - while (cur != NULL) //check the basis set - { - getNodeName(cname, cur); - if (cname == "configuration" || cname == "ci") - { - RealType qc_ci = 0.0; - std::string alpha, beta, tag; - OhmmsAttributeSet confAttrib; -#ifdef QMC_COMPLEX - RealType ci_real = 0.0, ci_imag = 0.0; - confAttrib.add(ci_real, "coeff_real"); - confAttrib.add(ci_imag, "coeff_imag"); -#else - RealType ci = 0.0; - confAttrib.add(ci, "coeff"); -#endif - confAttrib.add(qc_ci, "qchem_coeff"); - confAttrib.add(alpha, "alpha"); - confAttrib.add(beta, "beta"); - confAttrib.add(tag, "id"); - confAttrib.put(cur); - -#ifdef QMC_COMPLEX - ValueType ci(ci_real, ci_imag); -#endif - - //Will always loop through the whole determinant set as no assumption on the order of the determinant is made - if (std::abs(ci) < cutoff) - { - cur = cur->next; - continue; - } - - for (size_t i = 0; i < nstates; i++) - { - if (alpha[i] != '0' && alpha[i] != '1') - { - std::cerr << alpha << std::endl; - APP_ABORT("Found incorrect determinant label."); - } - if (beta[i] != '0' && beta[i] != '1') - { - std::cerr << beta << std::endl; - APP_ABORT("Found incorrect determinant label."); - } - } - - if (alpha.size() < nstates) - { - std::cerr << "alpha: " << alpha << std::endl; - APP_ABORT("Found incorrect alpha determinant label. size < nca+nstates"); - } - if (beta.size() < nstates) - { - std::cerr << "beta: " << beta << std::endl; - APP_ABORT("Found incorrect beta determinant label. size < ncb+nstates"); - } - - coeff.push_back(ci); - CItags.push_back(tag); - - std::unordered_map::const_iterator gotup = MyMapUp.find(alpha); - - - if (gotup == MyMapUp.end()) - { - uniqueConfg_up.push_back(dummyC_alpha); - uniqueConfg_up.back().add_occupation(alpha); - C2node_up.push_back(cntup); - MyMapUp.insert(std::pair(alpha, cntup)); - cntup++; - } - else - { - C2node_up.push_back(gotup->second); - } - - std::unordered_map::const_iterator gotdn = MyMapDn.find(beta); - - if (gotdn == MyMapDn.end()) - { - uniqueConfg_dn.push_back(dummyC_beta); - uniqueConfg_dn.back().add_occupation(beta); - C2node_dn.push_back(cntdn); - MyMapDn.insert(std::pair(beta, cntdn)); - cntdn++; - } - else - { - C2node_dn.push_back(gotdn->second); - } - - //if (qc_ci == 0.0) - // qc_ci = ci; - - cnt0++; - sumsq_qc += qc_ci * qc_ci; - sumsq += std::abs(ci * ci); - } - cur = cur->next; - } - - app_log() << "Found " << coeff.size() << " terms in the MSD expansion.\n"; - app_log() << "Norm of ci vector (sum of ci^2): " << sumsq << std::endl; - app_log() << "Norm of qchem ci vector (sum of qchem_ci^2): " << sumsq_qc << std::endl; - - } //usingCSF - - app_log() << "Found " << uniqueConfg_up.size() << " unique up determinants.\n"; - app_log() << "Found " << uniqueConfg_dn.size() << " unique down determinants.\n"; - - return success; -} - bool SlaterDetBuilder::readDetList(xmlNodePtr cur, std::vector>& uniqueConfgs, std::vector>& C2nodes, @@ -1642,259 +1213,6 @@ bool SlaterDetBuilder::readDetList(xmlNodePtr cur, return success; } - -bool SlaterDetBuilder::readDetListH5(xmlNodePtr cur, - std::vector& uniqueConfg_up, - std::vector& uniqueConfg_dn, - std::vector& C2node_up, - std::vector& C2node_dn, - std::vector& CItags, - std::vector& coeff, - bool& optimizeCI, - int nels_up, - int nels_dn) -{ - bool success = true; - int extlevel(0); - uniqueConfg_up.clear(); - uniqueConfg_dn.clear(); - C2node_up.clear(); - C2node_dn.clear(); - CItags.clear(); - coeff.clear(); - std::string CICoeffH5path(""); - std::vector confgList_up; - std::vector confgList_dn; - std::vector CIcoeff; - std::vector ConfigTag; - std::string optCI = "no"; - RealType cutoff = 0.0; - OhmmsAttributeSet ciAttrib; - ciAttrib.add(optCI, "optimize"); - ciAttrib.put(cur); - optimizeCI = (optCI == "yes"); - xmlNodePtr curRoot = cur, DetListNode = nullptr; - std::string cname, cname0, multidetH5path; - cur = curRoot->children; - while (cur != NULL) //check the basis set - { - getNodeName(cname, cur); - if (cname == "detlist") - { - DetListNode = cur; - app_log() << "Found determinant list. \n"; - } - cur = cur->next; - } - - app_log() << " H5 code path implicitly assumes NCA = NCB = 0, " - << "NEA = " << nels_up << " and NEB = " << nels_dn << std::endl; - size_t nstates = 0; - size_t ndets = 0; - size_t H5_ndets, H5_nstates; - /// 64 bit fixed width integer - const unsigned bit_kind = 64; - static_assert(bit_kind == sizeof(int64_t) * 8, "Must be 64 bit fixed width integer"); - /// the number of 64 bit integers which represent the binary string for occupation - int N_int; - std::string Dettype = "DETS"; - ValueType sumsq = 0.0; - OhmmsAttributeSet spoAttrib; - spoAttrib.add(ndets, "size"); - spoAttrib.add(Dettype, "type"); - spoAttrib.add(nstates, "nstates"); - spoAttrib.add(extlevel, "ext_level"); - spoAttrib.add(cutoff, "cutoff"); - spoAttrib.add(multidetH5path, "href"); - spoAttrib.add(CICoeffH5path, "opt_coeffs"); - spoAttrib.put(DetListNode); - if (ndets == 0) - { - APP_ABORT("size==0 in detlist is not allowed. Use slaterdeterminant in this case.\n"); - } - if (Dettype != "DETS" && Dettype != "Determinants") - APP_ABORT("Reading from HDF5 is only enabled for CI DETS. Must be accessed through (type=\"DETS\") or " - "(type=\"Determinants\") .\n"); - app_log() << "Reading CI expansion from HDF5:" << multidetH5path << std::endl; - - hdf_archive hin; - if (!hin.open(multidetH5path.c_str(), H5F_ACC_RDONLY)) - { - std::cerr << "Could not open H5 file" << std::endl; - abort(); - } - - if (!hin.push("MultiDet")) - { - std::cerr << "Could not open Multidet Group in H5 file" << std::endl; - abort(); - } - - hin.read(H5_ndets, "NbDet"); - if (ndets != H5_ndets) - { - std::cerr << "Number of determinants in H5 file (" << H5_ndets << ") different from number of dets in XML (" - << ndets << ")" << std::endl; - abort(); - } - - hin.read(H5_nstates, "nstate"); - if (nstates == 0) - nstates = H5_nstates; - else if (nstates != H5_nstates) - { - std::cerr << "Number of states/orbitals in H5 file (" << H5_nstates - << ") different from number of states/orbitals in XML (" << nstates << ")" << std::endl; - abort(); - } - - hin.read(N_int, "Nbits"); - CIcoeff.resize(ndets); - ConfigTag.resize(ndets); - - readCoeffs(hin, CIcoeff, ndets, extlevel); - - ///IF OPTIMIZED COEFFICIENTS ARE PRESENT IN opt_coeffs Path - ///THEY ARE READ FROM DIFFERENT HDF5 the replace the previous coeff - ///It is important to still read all old coeffs and only replace the optimized ones - ///in order to keep coherence with the cutoff on the number of determinants - ///REMEMBER!! FIRST COEFF IS FIXED. THEREFORE WE DO NOT REPLACE IT!!! - if (CICoeffH5path != "") - { - int OptCiSize = 0; - std::vector CIcoeffopt; - hdf_archive coeffin; - if (!coeffin.open(CICoeffH5path.c_str(), H5F_ACC_RDONLY)) - { - std::cerr << "Could not open H5 file containing Optimized Coefficients" << std::endl; - abort(); - } - - if (!coeffin.push("MultiDet")) - { - std::cerr << "Could not open Multidet Group in H5 file" << std::endl; - abort(); - } - coeffin.read(OptCiSize, "NbDet"); - CIcoeffopt.resize(OptCiSize); - - readCoeffs(coeffin, CIcoeffopt, ndets, extlevel); - - coeffin.close(); - - for (int i = 0; i < OptCiSize; i++) - CIcoeff[i + 1] = CIcoeffopt[i]; - - app_log() << "The first " << OptCiSize - << " Optimized coefficients were substituted to the original set of coefficients." << std::endl; - } - - Matrix tempAlpha(ndets, N_int); - hin.read(tempAlpha, "CI_Alpha"); - - Matrix tempBeta(ndets, N_int); - hin.read(tempBeta, "CI_Beta"); - - std::string MyCIAlpha, MyCIBeta; - MyCIAlpha.resize(nstates); - MyCIBeta.resize(nstates); - - ci_configuration dummyC_alpha; - ci_configuration dummyC_beta; - - dummyC_alpha.occup.resize(nstates, false); - dummyC_beta.occup.resize(nstates, false); - - for (size_t i = 0; i < nstates; i++) - dummyC_alpha.occup[i] = true; - - for (size_t i = 0; i < nstates; i++) - dummyC_beta.occup[i] = true; - - hin.close(); - app_log() << " Done reading " << ndets << " CIs from H5!" << std::endl; - - int cntup = 0; - int cntdn = 0; - std::unordered_map MyMapUp; - std::unordered_map MyMapDn; - - app_log() << " Sorting unique CIs" << std::endl; - ///This loop will find all unique Determinants in and store them "unsorted" in a new container uniqueConfg_up - ///and uniqueConfg_dn. The sorting is not done here - for (int ni = 0; ni < ndets; ni++) - { - if (std::abs(CIcoeff[ni]) < cutoff) - continue; - - int j = 0; - for (int k = 0; k < N_int; k++) - { - int64_t a = tempAlpha[ni][k]; - std::bitset a2 = a; - - auto b = tempBeta[ni][k]; - auto b2 = std::bitset(b); - - for (int i = 0; i < bit_kind; i++) - { - if (j < nstates) - { - MyCIAlpha[j] = a2[i] ? '1' : '0'; - MyCIBeta[j] = b2[i] ? '1' : '0'; - j++; - } - } - } - coeff.push_back(CIcoeff[ni]); - std::ostringstream h5tag; - h5tag << "CIcoeff_" << ni; - CItags.push_back(h5tag.str()); - - std::unordered_map::const_iterator gotup = MyMapUp.find(MyCIAlpha); - - if (gotup == MyMapUp.end()) - { - uniqueConfg_up.push_back(dummyC_alpha); - uniqueConfg_up.back().add_occupation(MyCIAlpha); - C2node_up.push_back(cntup); - MyMapUp.insert(std::pair(MyCIAlpha, cntup)); - cntup++; - } - else - { - C2node_up.push_back(gotup->second); - } - - std::unordered_map::const_iterator gotdn = MyMapDn.find(MyCIBeta); - - if (gotdn == MyMapDn.end()) - { - uniqueConfg_dn.push_back(dummyC_beta); - uniqueConfg_dn.back().add_occupation(MyCIBeta); - C2node_dn.push_back(cntdn); - MyMapDn.insert(std::pair(MyCIBeta, cntdn)); - cntdn++; - } - else - { - C2node_dn.push_back(gotdn->second); - } - - sumsq += CIcoeff[ni] * CIcoeff[ni]; - } - - app_log() << " Done Sorting unique CIs" << std::endl; - app_log() << "Found " << coeff.size() << " terms in the MSD expansion.\n"; - app_log() << "Norm of ci vector (sum of ci^2): " << sumsq << std::endl; - - app_log() << "Found " << uniqueConfg_up.size() << " unique up determinants.\n"; - app_log() << "Found " << uniqueConfg_dn.size() << " unique down determinants.\n"; - - - return success; -} - bool SlaterDetBuilder::readDetListH5(xmlNodePtr cur, std::vector>& uniqueConfgs, std::vector>& C2nodes, diff --git a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.h b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.h index 47f71a9507..ebda1eb284 100644 --- a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.h +++ b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.h @@ -89,22 +89,8 @@ class SlaterDetBuilder : public WaveFunctionComponentBuilder bool& Optimizable, bool& CI_Optimizable, xmlNodePtr cur); - - bool readDetList(xmlNodePtr cur, - std::vector& uniqueConfg_up, - std::vector& uniqueConfg_dn, - std::vector& C2node_up, - std::vector& C2node_dn, - std::vector& CItags, - std::vector& coeff, - bool& optimizeCI, - int nels_up, - int nels_dn, - std::vector& CSFcoeff, - std::vector& DetsPerCSF, - std::vector& CSFexpansion, - bool& usingCSF); + bool readDetList(xmlNodePtr cur, std::vector>& uniqueConfgs, std::vector>& C2nodes, @@ -117,18 +103,6 @@ class SlaterDetBuilder : public WaveFunctionComponentBuilder std::vector& CSFexpansion, bool& usingCSF); - - bool readDetListH5(xmlNodePtr cur, - std::vector& uniqueConfg_up, - std::vector& uniqueConfg_dn, - std::vector& C2node_up, - std::vector& C2node_dn, - std::vector& CItags, - std::vector& coeff, - bool& optimizeCI, - int nels_up, - int nels_dn); - bool readDetListH5(xmlNodePtr cur, std::vector>& uniqueConfgs, std::vector>& C2nodes, From 2463b1870d7dd29d04a82248e5fbffab297ea990 Mon Sep 17 00:00:00 2001 From: camelto2 Date: Mon, 22 Feb 2021 15:23:52 -0700 Subject: [PATCH 11/16] nGroups to const --- src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp index f560ce3b43..9223ba3421 100644 --- a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp +++ b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp @@ -204,7 +204,7 @@ WaveFunctionComponent* SlaterDetBuilder::buildComponent(xmlNodePtr cur) if (multislaterdet_0 || multislaterdetfast_0) APP_ABORT("multideterminant is already instantiated."); - int nGroups = targetPtcl.groups(); + const int nGroups = targetPtcl.groups(); std::vector spoNames(nGroups); std::string fastAlg; @@ -565,7 +565,7 @@ bool SlaterDetBuilder::createMSDFast(std::vector nptcls(nGroups); for (int grp = 0; grp < nGroups; grp++) @@ -852,7 +852,7 @@ bool SlaterDetBuilder::readDetList(xmlNodePtr cur, { bool success = true; - int nGroups = uniqueConfgs.size(); + const int nGroups = uniqueConfgs.size(); for (int grp = 0; grp < nGroups; grp++) { uniqueConfgs[grp].clear(); @@ -1223,7 +1223,7 @@ bool SlaterDetBuilder::readDetListH5(xmlNodePtr cur, { bool success = true; int extlevel(0); - int nGroups = uniqueConfgs.size(); + const int nGroups = uniqueConfgs.size(); for (int grp = 0; grp < nGroups; grp++) { uniqueConfgs[grp].clear(); From d0d528ec794560673c5e22966c3efcf067814115 Mon Sep 17 00:00:00 2001 From: camelto2 Date: Mon, 22 Feb 2021 16:02:45 -0700 Subject: [PATCH 12/16] change abort to barrier_and_abort --- src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp index 9223ba3421..7e74a172db 100644 --- a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp +++ b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp @@ -239,9 +239,9 @@ WaveFunctionComponent* SlaterDetBuilder::buildComponent(xmlNodePtr cur) spos[grp] = sposet_builder_factory_.getSPOSet(spoNames[grp]); if (spos[grp] == nullptr) { - app_error() << "In SlaterDetBuilder: SPOSet \"" << spoNames[grp] - << "\" is not found. Expected for MultiSlaterDeterminant.\n"; - abort(); + std::stringstream err_msg; + err_msg << "In SlaterDetBuilder: SPOSet \"" << spoNames[grp] << "\" is not found. Expected for MultiSlaterDeterminant." << std::endl; + myComm->barrier_and_abort(err_msg.str()); } spo_clones[grp].reset(spos[grp]->makeClone()); } From 5b9fcb4861da8862af0edfd6bd5e9d23dd8693e2 Mon Sep 17 00:00:00 2001 From: camelto2 Date: Tue, 23 Feb 2021 10:26:33 -0700 Subject: [PATCH 13/16] change spo_clones --- src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp index 7e74a172db..dfc9c48dcc 100644 --- a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp +++ b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp @@ -232,7 +232,7 @@ WaveFunctionComponent* SlaterDetBuilder::buildComponent(xmlNodePtr cur) //new format std::vector spos(nGroups); - std::vector> spo_clones(nGroups); + std::vector> spo_clones; for (int grp = 0; grp < nGroups; grp++) { @@ -243,7 +243,8 @@ WaveFunctionComponent* SlaterDetBuilder::buildComponent(xmlNodePtr cur) err_msg << "In SlaterDetBuilder: SPOSet \"" << spoNames[grp] << "\" is not found. Expected for MultiSlaterDeterminant." << std::endl; myComm->barrier_and_abort(err_msg.str()); } - spo_clones[grp].reset(spos[grp]->makeClone()); + std::unique_ptr tmp(spos[grp]->makeClone()); + spo_clones.emplace_back(std::move(tmp)); } From 8c0bf87a5a17cbe80785da8c5797588bc79c00a9 Mon Sep 17 00:00:00 2001 From: camelto2 Date: Tue, 23 Feb 2021 12:13:59 -0700 Subject: [PATCH 14/16] remove spos vector --- src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp index dfc9c48dcc..e9effacf1a 100644 --- a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp +++ b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp @@ -231,20 +231,18 @@ WaveFunctionComponent* SlaterDetBuilder::buildComponent(xmlNodePtr cur) } //new format - std::vector spos(nGroups); std::vector> spo_clones; for (int grp = 0; grp < nGroups; grp++) { - spos[grp] = sposet_builder_factory_.getSPOSet(spoNames[grp]); - if (spos[grp] == nullptr) + SPOSetPtr spo_tmp = sposet_builder_factory_.getSPOSet(spoNames[grp]); + if (spo_tmp == nullptr) { std::stringstream err_msg; err_msg << "In SlaterDetBuilder: SPOSet \"" << spoNames[grp] << "\" is not found. Expected for MultiSlaterDeterminant." << std::endl; myComm->barrier_and_abort(err_msg.str()); } - std::unique_ptr tmp(spos[grp]->makeClone()); - spo_clones.emplace_back(std::move(tmp)); + spo_clones.emplace_back(spo_tmp->makeClone()); } From 5ebface033177140eb301d1e7f3240f27f8897d6 Mon Sep 17 00:00:00 2001 From: camelto2 Date: Tue, 23 Feb 2021 14:57:54 -0700 Subject: [PATCH 15/16] fix unit test to work for real build --- .../test_spo_collection_input_MSD_LCAO_h5.cpp | 28 +++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/src/QMCWaveFunctions/tests/test_spo_collection_input_MSD_LCAO_h5.cpp b/src/QMCWaveFunctions/tests/test_spo_collection_input_MSD_LCAO_h5.cpp index ded9ae5650..ccc7090126 100644 --- a/src/QMCWaveFunctions/tests/test_spo_collection_input_MSD_LCAO_h5.cpp +++ b/src/QMCWaveFunctions/tests/test_spo_collection_input_MSD_LCAO_h5.cpp @@ -270,6 +270,7 @@ TEST_CASE("SPO input spline from xml LiH_msd arbitrary species", "[wavefunction] app_log() << "-------------------------------------------------------------" << std::endl; app_log() << "LiH_msd with positron xml input style" << std::endl; app_log() << "-------------------------------------------------------------" << std::endl; +#ifdef QMC_COMPLEX const char* spo_xml_string1 = " \ \ \ @@ -295,6 +296,33 @@ TEST_CASE("SPO input spline from xml LiH_msd arbitrary species", "[wavefunction] \ \ "; +#else + const char* spo_xml_string1 = " \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ +"; +#endif test_LiH_msd_xml_input_with_positron(spo_xml_string1, "spo-ps", 5, 105); } } // namespace qmcplusplus From 321b5961de8475d86456d86568bf65c2840d469c Mon Sep 17 00:00:00 2001 From: camelto2 Date: Tue, 23 Feb 2021 16:55:53 -0700 Subject: [PATCH 16/16] change unit test and CI coeff parsing in readDetList --- .../Fermion/SlaterDetBuilder.cpp | 42 +++++++++---------- .../test_spo_collection_input_MSD_LCAO_h5.cpp | 28 ------------- 2 files changed, 20 insertions(+), 50 deletions(-) diff --git a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp index e9effacf1a..d673d748eb 100644 --- a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp +++ b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp @@ -239,7 +239,8 @@ WaveFunctionComponent* SlaterDetBuilder::buildComponent(xmlNodePtr cur) if (spo_tmp == nullptr) { std::stringstream err_msg; - err_msg << "In SlaterDetBuilder: SPOSet \"" << spoNames[grp] << "\" is not found. Expected for MultiSlaterDeterminant." << std::endl; + err_msg << "In SlaterDetBuilder: SPOSet \"" << spoNames[grp] + << "\" is not found. Expected for MultiSlaterDeterminant." << std::endl; myComm->barrier_and_abort(err_msg.str()); } spo_clones.emplace_back(spo_tmp->makeClone()); @@ -961,28 +962,27 @@ bool SlaterDetBuilder::readDetList(xmlNodePtr cur, RealType exctLvl, qc_ci = 0.0; OhmmsAttributeSet confAttrib; std::string tag, OccString; -#ifdef QMC_COMPLEX RealType ci_real = 0.0, ci_imag = 0.0; + confAttrib.add(ci_real, "coeff"); confAttrib.add(ci_real, "coeff_real"); confAttrib.add(ci_imag, "coeff_imag"); -#else - RealType ci = 0.0; - confAttrib.add(ci, "coeff"); -#endif confAttrib.add(qc_ci, "qchem_coeff"); confAttrib.add(tag, "id"); confAttrib.add(OccString, "occ"); confAttrib.add(exctLvl, "exctLvl"); confAttrib.put(cur); + if (qc_ci == 0.0) -#ifdef QMC_COMPLEX qc_ci = ci_real; -#else - qc_ci = ci; -#endif + ValueType ci; #ifdef QMC_COMPLEX - ValueType ci(ci_real, ci_imag); + ci = ValueType(ci_real, ci_imag); +#else + if (ci_imag != RealType(0)) + myComm->barrier_and_abort( + "SlaterDetBuilder::readDetList. Build with QMC_COMPLEX if using complex CI expansion coefficients."); + ci = ci_real; #endif //Can discriminate based on any of 3 criterion if (((std::abs(qc_ci) < cutoff) && (CSFChoice == "qchem_coeff")) || @@ -994,11 +994,7 @@ bool SlaterDetBuilder::readDetList(xmlNodePtr cur, } cnt0++; if (std::abs(qc_ci) < zero_cutoff) -#ifdef QMC_COMPLEX - ci_real = 0.0; -#else - ci = 0.0; -#endif + ci = 0.0; CSFcoeff.push_back(ci); sumsq_qc += qc_ci * qc_ci; DetsPerCSF.push_back(0); @@ -1126,14 +1122,10 @@ bool SlaterDetBuilder::readDetList(xmlNodePtr cur, std::vector occs(nGroups); std::string tag; OhmmsAttributeSet confAttrib; -#ifdef QMC_COMPLEX RealType ci_real = 0.0, ci_imag = 0.0; + confAttrib.add(ci_real, "coeff"); confAttrib.add(ci_real, "coeff_real"); confAttrib.add(ci_imag, "coeff_imag"); -#else - RealType ci = 0.0; - confAttrib.add(ci, "coeff"); -#endif confAttrib.add(qc_ci, "qchem_coeff"); for (int grp = 0; grp < nGroups; grp++) confAttrib.add(occs[grp], "occ" + std::to_string(grp)); @@ -1145,8 +1137,14 @@ bool SlaterDetBuilder::readDetList(xmlNodePtr cur, confAttrib.add(tag, "id"); confAttrib.put(cur); + ValueType ci; #ifdef QMC_COMPLEX - ValueType ci(ci_real, ci_imag); + ci = ValueType(ci_real, ci_imag); +#else + if (ci_imag != RealType(0)) + myComm->barrier_and_abort( + "SlaterDetBuilder::readDetList. Build with QMC_COMPLEX if using complex CI expansion coefficients."); + ci = ci_real; #endif //Will always loop through the whole determinant set as no assumption on the order of the determinant is made diff --git a/src/QMCWaveFunctions/tests/test_spo_collection_input_MSD_LCAO_h5.cpp b/src/QMCWaveFunctions/tests/test_spo_collection_input_MSD_LCAO_h5.cpp index ccc7090126..bcca7bda57 100644 --- a/src/QMCWaveFunctions/tests/test_spo_collection_input_MSD_LCAO_h5.cpp +++ b/src/QMCWaveFunctions/tests/test_spo_collection_input_MSD_LCAO_h5.cpp @@ -270,33 +270,6 @@ TEST_CASE("SPO input spline from xml LiH_msd arbitrary species", "[wavefunction] app_log() << "-------------------------------------------------------------" << std::endl; app_log() << "LiH_msd with positron xml input style" << std::endl; app_log() << "-------------------------------------------------------------" << std::endl; -#ifdef QMC_COMPLEX - const char* spo_xml_string1 = " \ - \ - \ - \ - \ - \ - \ - \ - \ - \ - \ - \ - \ - \ - \ - \ - \ - \ - \ - \ - \ - \ - \ - \ -"; -#else const char* spo_xml_string1 = " \ \ \ @@ -322,7 +295,6 @@ TEST_CASE("SPO input spline from xml LiH_msd arbitrary species", "[wavefunction] \ \ "; -#endif test_LiH_msd_xml_input_with_positron(spo_xml_string1, "spo-ps", 5, 105); } } // namespace qmcplusplus