diff --git a/src/QMCWaveFunctions/CMakeLists.txt b/src/QMCWaveFunctions/CMakeLists.txt index 40e8e5a17d..9d6eace715 100644 --- a/src/QMCWaveFunctions/CMakeLists.txt +++ b/src/QMCWaveFunctions/CMakeLists.txt @@ -87,7 +87,7 @@ if(OHMMS_DIM MATCHES 3) LCAO/AOBasisBuilder.cpp LCAO/SoaLocalizedBasisSet.cpp) if(QMC_COMPLEX) - set(FERMION_SRCS ${FERMION_SRCS} LCAO/LCAOSpinorBuilder.cpp) + set(FERMION_SRCS ${FERMION_SRCS} LCAO/LCAOSpinorBuilder.cpp LCAO/LCAOSpinorBuilder.cpp) else(QMC_COMPLEX) #LCAO cusp correction is not ready for complex set(FERMION_SRCS ${FERMION_SRCS} diff --git a/src/QMCWaveFunctions/LCAO/LCAOSpinorBuilderT.cpp b/src/QMCWaveFunctions/LCAO/LCAOSpinorBuilderT.cpp new file mode 100644 index 0000000000..4e5a3fd2b0 --- /dev/null +++ b/src/QMCWaveFunctions/LCAO/LCAOSpinorBuilderT.cpp @@ -0,0 +1,207 @@ +////////////////////////////////////////////////////////////////////////////////////// +// This file is distributed under the University of Illinois/NCSA Open Source License. +// See LICENSE file in top directory for details. +// +// Copyright (c) 2020 QMCPACK developers. +// +// File developed by: Cody A. Melton, cmelton@sandia.gov, Sandia National Laboratories +// +// File created by: Cody A. Melton, cmelton@sandia.gov, Sandia National Laboratories +////////////////////////////////////////////////////////////////////////////////////// + +#include "LCAOSpinorBuilder.h" +#include "QMCWaveFunctions/SpinorSet.h" +#include "OhmmsData/AttributeSet.h" +#include "Utilities/ProgressReportEngine.h" +#include "hdf/hdf_archive.h" +#include "Message/CommOperators.h" + +namespace qmcplusplus +{ +template +LCAOSpinorBuilderT::LCAOSpinorBuilderT(ParticleSet& els, ParticleSet& ions, Communicate* comm, xmlNodePtr cur) + : LCAOrbitalBuilder(els, ions, comm, cur) +{ + ClassName = "LCAOSpinorBuilder"; + + if (h5_path == "") + myComm->barrier_and_abort("LCAOSpinorBuilder only works with href"); +} + +template +std::unique_ptr> LCAOSpinorBuilderT::createSPOSetFromXML(xmlNodePtr cur) +{ + ReportEngine PRE(ClassName, "createSPO(xmlNodePtr)"); + std::string spo_name(""), optimize("no"); + std::string basisset_name("LCAOBSet"); + OhmmsAttributeSet spoAttrib; + spoAttrib.add(spo_name, "name"); + spoAttrib.add(optimize, "optimize"); + spoAttrib.add(basisset_name, "basisset"); + spoAttrib.put(cur); + + BasisSet_t* myBasisSet = nullptr; + if (basisset_map_.find(basisset_name) == basisset_map_.end()) + myComm->barrier_and_abort("basisset \"" + basisset_name + "\" cannot be found\n"); + else + myBasisSet = basisset_map_[basisset_name].get(); + + if (optimize == "yes") + app_log() << " SPOSet " << spo_name << " is optimizable\n"; + + std::unique_ptr upspo = + std::make_unique(spo_name + "_up", std::unique_ptr(myBasisSet->makeClone())); + std::unique_ptr dnspo = + std::make_unique(spo_name + "_dn", std::unique_ptr(myBasisSet->makeClone())); + + loadMO(*upspo, *dnspo, cur); + + //create spinor and register up/dn + auto spinor_set = std::make_unique(spo_name); + spinor_set->set_spos(std::move(upspo), std::move(dnspo)); + return spinor_set; +} + +template +bool LCAOSpinorBuilderT::loadMO(LCAOrbitalSetT& up, LCAOrbitalSetT& dn, xmlNodePtr cur) +{ + bool PBC = false; + int norb = up.getBasisSetSize(); + std::string debugc("no"); + OhmmsAttributeSet aAttrib; + aAttrib.add(norb, "size"); + aAttrib.add(debugc, "debug"); + aAttrib.put(cur); + + up.setOrbitalSetSize(norb); + dn.setOrbitalSetSize(norb); + + xmlNodePtr occ_ptr = nullptr; + cur = cur->xmlChildrenNode; + while (cur != nullptr) + { + std::string cname((const char*)(cur->name)); + if (cname == "occupation") + { + occ_ptr = cur; + } + cur = cur->next; + } + + hdf_archive hin(myComm); + if (myComm->rank() == 0) + { + if (!hin.open(h5_path, H5F_ACC_RDONLY)) + myComm->barrier_and_abort("LCAOSpinorBuilder::loadMO missing or incorrect path to H5 file."); + hin.push("PBC"); + PBC = false; + hin.read(PBC, "PBC"); + hin.close(); + } + myComm->bcast(PBC); + if (PBC) + myComm->barrier_and_abort("LCAOSpinorBuilder::loadMO lcao spinors not implemented in PBC"); + + bool success = putFromH5(up, dn, occ_ptr); + + + if (debugc == "yes") + { + app_log() << "UP: Single-particle orbital coefficients dims=" << up.C->rows() << " x " << up.C->cols() + << std::endl; + app_log() << *up.C << std::endl; + app_log() << "DN: Single-particle orbital coefficients dims=" << dn.C->rows() << " x " << dn.C->cols() + << std::endl; + app_log() << *dn.C << std::endl; + } + return success; +} + +template +bool LCAOSpinorBuilderT::putFromH5(LCAOrbitalSetT& up, LCAOrbitalSetT& dn, xmlNodePtr occ_ptr) +{ + if (up.getBasisSetSize() == 0 || dn.getBasisSetSize() == 0) + { + myComm->barrier_and_abort("LCASpinorBuilder::loadMO detected ZERO BasisSetSize"); + return false; + } + + bool success = true; + hdf_archive hin(myComm); + if (myComm->rank() == 0) + { + istd::string setname = "/Super_Twist/eigenset_0"; + readRealMatrixFromH5(hin, setname, upReal); + setname += "_imag"; + readRealMatrixFromH5(hin, setname, upImag); + + af(!hin.open(h5_path, H5F_ACC_RDONLY)) + myComm->barrier_and_abort("LCAOSpinorBuilder::putFromH5 missing or incorrect path to H5 file"); + + Matrix upReal; + Matrix upImag; + ssert(upReal.rows() == upImag.rows()); + assert(upReal.cols() == upImag.cols()); + + Matrix upTemp(upReal.rows(), upReal.cols()); + for (int i = 0; i < upTemp.rows(); i++) + { + for (int j = 0; j < upTemp.cols(); j++) + { + upTemp[i][j] = ValueType(upReal[i][j], upImag[i][j]); + } + } + + Matrix dnReal; + Matrix dnImag; + setname = "/Super_Twist/eigenset_1"; + readRealMatrixFromH5(hin, setname, dnReal); + setname += "_imag"; + readRealMatrixFromH5(hin, setname, dnImag); + + assert(dnReal.rows() == dnImag.rows()); + assert(dnReal.cols() == dnImag.cols()); + + Matrix dnTemp(dnReal.rows(), dnReal.cols()); + for (int i = 0; i < dnTemp.rows(); i++) + { + for (int j = 0; j < dnTemp.cols(); j++) + { + dnTemp[i][j] = ValueType(dnReal[i][j], dnImag[i][j]); + } + } + + assert(upReal.rows() == dnReal.rows()); + assert(upReal.cols() == dnReal.cols()); + + Occ.resize(upReal.rows()); + success = putOccupation(up, occ_ptr); + + int norbs = up.getOrbitalSetSize(); + + int n = 0, i = 0; + while (i < norbs) + { + if (Occ[n] > 0.0) + { + std::copy(upTemp[n], upTemp[n + 1], (*up.C)[i]); + std::copy(dnTemp[n], dnTemp[n + 1], (*dn.C)[i]); + i++; + } + n++; + } + + hin.close(); + } + +#ifdef HAVE_MPI + myComm->comm.broadcast_n(up.C->data(), up.C->size()); + myComm->comm.broadcast_n(dn.C->data(), dn.C->size()); +#endif + + return success; +} + +template class LCAOSpinorBuilderT>; +template class LCAOSpinorBuilderT>; +} // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/LCAO/LCAOSpinorBuilderT.h b/src/QMCWaveFunctions/LCAO/LCAOSpinorBuilderT.h new file mode 100644 index 0000000000..62b40b43b1 --- /dev/null +++ b/src/QMCWaveFunctions/LCAO/LCAOSpinorBuilderT.h @@ -0,0 +1,64 @@ +////////////////////////////////////////////////////////////////////////////////////// +// This file is distributed under the University of Illinois/NCSA Open Source License. +// See LICENSE file in top directory for details. +// +// Copyright (c) 2020 QMCPACK developers. +// +// File developed by: Cody A. Melton, cmelton@sandia.gov, Sandia National Laboratories +// +// File created by: Cody A. Melton, cmelton@sandia.gov, Sandia National Laboratories +////////////////////////////////////////////////////////////////////////////////////// + + +#ifndef QMCPLUSPLUS_SOA_LCAO_SPINOR_BUILDER_H +#define QMCPLUSPLUS_SOA_LCAO_SPINOR_BUILDER_H + +#include "QMCWaveFunctions/LCAO/LCAOrbitalBuilderT.h" + +namespace qmcplusplus +{ +/** @file LCAOSpinorBuidler.h + * + * Derives from LCAOrbitalBuilder.h. Overrides createSPOSetFromXML method to read up and + * down channel from HDF5 and construct SpinorSet + * + */ +template +class LCAOSpinorBuilderT : public LCAOrbitalBuilderT +{ +public: + /** constructor + * \param els reference to the electrons + * \param ions reference to the ions + * + * Derives from LCAOrbitalBuilder, but will require an h5_path to be set + */ + LCAOSpinorBuilderT(ParticleSet& els, ParticleSet& ions, Communicate* comm, xmlNodePtr cur); + + /** creates and returns SpinorSet + * + * Creates an up and down LCAOrbitalSet + * calls LCAOSpinorBuilder::loadMO to build up and down from the H5 file + * registers up and down into a SpinorSet and returns + */ + std::unique_ptr> createSPOSetFromXML(xmlNodePtr cur) override; + +private: + /** load the up and down MO sets + * + * checks to make sure not PBC and initialize the Occ vector. + * call putFromH5 to parse the up and down MO coefficients + */ + bool loadMO(LCAOrbitalSetT& up, LCAOrbitalSetT& dn, xmlNodePtr cur); + + /** parse h5 file for spinor info + * + * assumes the h5 file has KPTS_0/eigenset_0(_imag) for the real/imag part of up component of spinor + * assumes the h5 file as KPTS_0/eigenset_1(_imag) for the real/imag part of dn component of spinor + * reads the various coefficient matricies and broadcast + * after this, we have up/dn LCAOrbitalSet that can be registered to the SpinorSet + */ + bool putFromH5(LCAOrbitalSetT& up, LCAOrbitalSetT& dn, xmlNodePtr); +}; +} // namespace qmcplusplus +#endif diff --git a/src/QMCWaveFunctions/LCAO/LCAOrbitalBuilderT.cpp b/src/QMCWaveFunctions/LCAO/LCAOrbitalBuilderT.cpp index 4a9de5f847..4e1e4f6bd1 100644 --- a/src/QMCWaveFunctions/LCAO/LCAOrbitalBuilderT.cpp +++ b/src/QMCWaveFunctions/LCAO/LCAOrbitalBuilderT.cpp @@ -563,7 +563,7 @@ LCAOrbitalBuilderT::createSPOSetFromXML(xmlNodePtr cur) #ifdef HAVE_MPI for (int orb_idx = 0; orb_idx < orbital_set_size; orb_idx++) for (int center_idx = 0; center_idx < num_centers; center_idx++) - broadcastCuspInfo( + CuspCorrectionConstructionT::broadcastCuspInfo( info(center_idx, orb_idx), *this->myComm, 0); #endif } diff --git a/src/QMCWaveFunctions/SPOSetBuilderFactoryT.cpp b/src/QMCWaveFunctions/SPOSetBuilderFactoryT.cpp index 284eaa7efe..b98952f779 100644 --- a/src/QMCWaveFunctions/SPOSetBuilderFactoryT.cpp +++ b/src/QMCWaveFunctions/SPOSetBuilderFactoryT.cpp @@ -30,7 +30,7 @@ #if defined(QMC_COMPLEX) #include "QMCWaveFunctions/EinsplineSpinorSetBuilder.h" -#include "QMCWaveFunctions/LCAO/LCAOSpinorBuilder.h" +#include "QMCWaveFunctions/LCAO/LCAOSpinorBuilderT.h" #endif #if defined(HAVE_EINSPLINE) @@ -154,7 +154,7 @@ SPOSetBuilderFactoryT::createSPOSetBuilder(xmlNodePtr rootNode) ions = pit->second.get(); if (targetPtcl.isSpinor()) #ifdef QMC_COMPLEX - bb = std::make_unique( + bb = std::make_unique>( targetPtcl, *ions, myComm, rootNode); #else PRE.error("Use of lcao spinors requires QMC_COMPLEX=1. Rebuild " @@ -253,9 +253,11 @@ SPOSetBuilderFactoryT::addSPOSet(std::unique_ptr> spo) template std::string SPOSetBuilderFactoryT::basisset_tag = "basisset"; -template class SPOSetBuilderFactoryT; -template class SPOSetBuilderFactoryT; +#ifdef QMC_COMPLEX template class SPOSetBuilderFactoryT>; template class SPOSetBuilderFactoryT>; - +#else +template class SPOSetBuilderFactoryT; +template class SPOSetBuilderFactoryT; +#endif } // namespace qmcplusplus