From 17e55cce00939fb755ca0dc4e77a83b4bf7a72bd Mon Sep 17 00:00:00 2001 From: Steven Hahn Date: Fri, 4 Aug 2023 14:57:43 -0400 Subject: [PATCH] Add FakeSPOT class Move SpinorSet to a templated class Refactor FreeOrbital class Base typed aliases on SPOSet on OrbitalSetTraits Add FullRealType in SPOSet and RotatedSPOs Add this in templated meta class Add explicit function instantions for FreeOrbital Add templated class SHOSetT Signed-off-by: Steven Hahn Add PWRealOrbitalSetT template class Revert test_RotatedSPOs.cpp --- src/QMCWaveFunctions/CMakeLists.txt | 9 +- .../ElectronGas/FreeOrbitalT.cpp | 718 +++++++ .../ElectronGas/FreeOrbitalT.h | 88 + .../HarmonicOscillator/SHOSetT.cpp | 577 ++++++ .../HarmonicOscillator/SHOSetT.h | 158 ++ .../PlaneWave/PWRealOrbitalSetT.cpp | 165 ++ .../PlaneWave/PWRealOrbitalSetT.h | 143 ++ src/QMCWaveFunctions/RotatedSPOsT.cpp | 1690 +++++++++++++++++ src/QMCWaveFunctions/RotatedSPOsT.h | 420 ++++ src/QMCWaveFunctions/SPOSetT.cpp | 6 +- src/QMCWaveFunctions/SPOSetT.h | 6 +- src/QMCWaveFunctions/SpinorSetT.cpp | 586 ++++++ src/QMCWaveFunctions/SpinorSetT.h | 229 +++ src/QMCWaveFunctions/tests/CMakeLists.txt | 2 +- src/QMCWaveFunctions/tests/FakeSPOT.cpp | 160 ++ src/QMCWaveFunctions/tests/FakeSPOT.h | 62 + .../tests/test_RotatedSPOs.cpp | 63 + 17 files changed, 5073 insertions(+), 9 deletions(-) create mode 100644 src/QMCWaveFunctions/ElectronGas/FreeOrbitalT.cpp create mode 100644 src/QMCWaveFunctions/ElectronGas/FreeOrbitalT.h create mode 100644 src/QMCWaveFunctions/HarmonicOscillator/SHOSetT.cpp create mode 100644 src/QMCWaveFunctions/HarmonicOscillator/SHOSetT.h create mode 100644 src/QMCWaveFunctions/PlaneWave/PWRealOrbitalSetT.cpp create mode 100644 src/QMCWaveFunctions/PlaneWave/PWRealOrbitalSetT.h create mode 100644 src/QMCWaveFunctions/RotatedSPOsT.cpp create mode 100644 src/QMCWaveFunctions/RotatedSPOsT.h create mode 100644 src/QMCWaveFunctions/SpinorSetT.cpp create mode 100644 src/QMCWaveFunctions/SpinorSetT.h create mode 100644 src/QMCWaveFunctions/tests/FakeSPOT.cpp create mode 100644 src/QMCWaveFunctions/tests/FakeSPOT.h diff --git a/src/QMCWaveFunctions/CMakeLists.txt b/src/QMCWaveFunctions/CMakeLists.txt index 11da0d4cf9..959e7743ae 100644 --- a/src/QMCWaveFunctions/CMakeLists.txt +++ b/src/QMCWaveFunctions/CMakeLists.txt @@ -35,16 +35,17 @@ set(WFBASE_SRCS SPOSetT.cpp CompositeSPOSet.cpp HarmonicOscillator/SHOSet.cpp + HarmonicOscillator/SHOSetT.cpp HarmonicOscillator/SHOSetBuilder.cpp ExampleHeBuilder.cpp ExampleHeComponent.cpp) if(NOT QMC_COMPLEX) - set(WFBASE_SRCS ${WFBASE_SRCS} RotatedSPOs.cpp) + set(WFBASE_SRCS ${WFBASE_SRCS} RotatedSPOs.cpp RotatedSPOsT.cpp) endif(NOT QMC_COMPLEX) if(QMC_COMPLEX) - set(WFBASE_SRCS ${WFBASE_SRCS} SpinorSet.cpp) + set(WFBASE_SRCS ${WFBASE_SRCS} SpinorSet.cpp SpinorSetT.cpp) endif(QMC_COMPLEX) ######################## # build jastrows @@ -63,7 +64,7 @@ set(JASTROW_SRCS set(JASTROW_OMPTARGET_SRCS Jastrow/TwoBodyJastrow.cpp Jastrow/BsplineFunctor.cpp) -set(FERMION_SRCS ${FERMION_SRCS} ElectronGas/FreeOrbital.cpp ElectronGas/FreeOrbitalBuilder.cpp) +set(FERMION_SRCS ${FERMION_SRCS} ElectronGas/FreeOrbital.cpp ElectronGas/FreeOrbitalT.cpp ElectronGas/FreeOrbitalBuilder.cpp) # wavefunctions only availbale to 3-dim problems if(OHMMS_DIM MATCHES 3) @@ -114,7 +115,7 @@ if(OHMMS_DIM MATCHES 3) if(QMC_COMPLEX) set(FERMION_SRCS ${FERMION_SRCS} PlaneWave/PWOrbitalSet.cpp) else() - set(FERMION_SRCS ${FERMION_SRCS} PlaneWave/PWRealOrbitalSet.cpp) + set(FERMION_SRCS ${FERMION_SRCS} PlaneWave/PWRealOrbitalSet.cpp PlaneWave/PWRealOrbitalSetT.cpp) endif(QMC_COMPLEX) if(NOT QMC_COMPLEX) diff --git a/src/QMCWaveFunctions/ElectronGas/FreeOrbitalT.cpp b/src/QMCWaveFunctions/ElectronGas/FreeOrbitalT.cpp new file mode 100644 index 0000000000..81bc37cc79 --- /dev/null +++ b/src/QMCWaveFunctions/ElectronGas/FreeOrbitalT.cpp @@ -0,0 +1,718 @@ +////////////////////////////////////////////////////////////////////////////////////// +// This file is distributed under the University of Illinois/NCSA Open Source License. +// See LICENSE file in top directory for details. +// +// Copyright (c) 2022 QMCPACK developers. +// +// File developed by: Miguel Morales, moralessilva2@llnl.gov, Lawrence Livermore National Laboratory +// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign +// Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign +// Jaron T. Krogel, krogeljt@ornl.gov, Oak Ridge National Laboratory +// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory +// Yubo "Paul" Yang, yubo.paul.yang@gmail.com, CCQ @ Flatiron +// William F Godoy, godoywf@ornl.gov, Oak Ridge National Laboratory +// +// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign +////////////////////////////////////////////////////////////////////////////////////// + +#include "FreeOrbitalT.h" + +namespace qmcplusplus +{ + +template +FreeOrbitalT::FreeOrbitalT(const std::string& my_name, const std::vector& kpts_cart) : SPOSetT(my_name) +{} + +//Explicit template specialization +template<> +FreeOrbitalT::FreeOrbitalT(const std::string& my_name, const std::vector& kpts_cart) + : SPOSetT(my_name), + kvecs(kpts_cart), + mink(1), // treat k=0 as special case + maxk(kpts_cart.size()), + k2neg(maxk) +{ + this->OrbitalSetSize = 2 * maxk - 1; // k=0 has no (cos, sin) split, SPOSet member + for (int ik = 0; ik < maxk; ik++) + k2neg[ik] = -dot(kvecs[ik], kvecs[ik]); +} + +template<> +FreeOrbitalT::FreeOrbitalT(const std::string& my_name, const std::vector& kpts_cart) + : SPOSetT(my_name), + kvecs(kpts_cart), + mink(1), // treat k=0 as special case + maxk(kpts_cart.size()), + k2neg(maxk) +{ + this->OrbitalSetSize = 2 * maxk - 1; // k=0 has no (cos, sin) split, SPOSet member + for (int ik = 0; ik < maxk; ik++) + k2neg[ik] = -dot(kvecs[ik], kvecs[ik]); +} + +template<> +FreeOrbitalT>::FreeOrbitalT(const std::string& my_name, const std::vector& kpts_cart) + : SPOSetT>(my_name), + kvecs(kpts_cart), + mink(0), // treat k=0 as special case + maxk(kpts_cart.size()), + k2neg(maxk) +{ + this->OrbitalSetSize = maxk; // SPOSet member + for (int ik = 0; ik < maxk; ik++) + k2neg[ik] = -dot(kvecs[ik], kvecs[ik]); +} + +template<> +FreeOrbitalT>::FreeOrbitalT(const std::string& my_name, const std::vector& kpts_cart) + : SPOSetT>(my_name), + kvecs(kpts_cart), + mink(0), // treat k=0 as special case + maxk(kpts_cart.size()), + k2neg(maxk) +{ + this->OrbitalSetSize = maxk; // SPOSet member + for (int ik = 0; ik < maxk; ik++) + k2neg[ik] = -dot(kvecs[ik], kvecs[ik]); +} + + +template +void FreeOrbitalT::evaluateVGL(const ParticleSet& P, + int iat, + ValueVector& pvec, + GradVector& dpvec, + ValueVector& d2pvec) +{} + +template<> +void FreeOrbitalT::evaluateVGL(const ParticleSet& P, + int iat, + ValueVector& pvec, + GradVector& dpvec, + ValueVector& d2pvec) +{ + const PosType& r = P.activeR(iat); + RealType sinkr, coskr; + for (int ik = mink; ik < maxk; ik++) + { + sincos(dot(kvecs[ik], r), &sinkr, &coskr); + const int j2 = 2 * ik; + const int j1 = j2 - 1; + pvec[j1] = coskr; + pvec[j2] = sinkr; + dpvec[j1] = -sinkr * kvecs[ik]; + dpvec[j2] = coskr * kvecs[ik]; + d2pvec[j1] = k2neg[ik] * coskr; + d2pvec[j2] = k2neg[ik] * sinkr; + } + pvec[0] = 1.0; + dpvec[0] = 0.0; + d2pvec[0] = 0.0; +} + +template<> +void FreeOrbitalT::evaluateVGL(const ParticleSet& P, + int iat, + ValueVector& pvec, + GradVector& dpvec, + ValueVector& d2pvec) +{ + const PosType& r = P.activeR(iat); + RealType sinkr, coskr; + for (int ik = mink; ik < maxk; ik++) + { + sincos(dot(kvecs[ik], r), &sinkr, &coskr); + const int j2 = 2 * ik; + const int j1 = j2 - 1; + pvec[j1] = coskr; + pvec[j2] = sinkr; + dpvec[j1] = -sinkr * kvecs[ik]; + dpvec[j2] = coskr * kvecs[ik]; + d2pvec[j1] = k2neg[ik] * coskr; + d2pvec[j2] = k2neg[ik] * sinkr; + } + pvec[0] = 1.0; + dpvec[0] = 0.0; + d2pvec[0] = 0.0; +} + + +template<> +void FreeOrbitalT>::evaluateVGL(const ParticleSet& P, + int iat, + ValueVector& pvec, + GradVector& dpvec, + ValueVector& d2pvec) +{ + const PosType& r = P.activeR(iat); + RealType sinkr, coskr; + for (int ik = mink; ik < maxk; ik++) + { + sincos(dot(kvecs[ik], r), &sinkr, &coskr); + + pvec[ik] = ValueType(coskr, sinkr); + dpvec[ik] = ValueType(-sinkr, coskr) * kvecs[ik]; + d2pvec[ik] = ValueType(k2neg[ik] * coskr, k2neg[ik] * sinkr); + } +} + +template<> +void FreeOrbitalT>::evaluateVGL(const ParticleSet& P, + int iat, + ValueVector& pvec, + GradVector& dpvec, + ValueVector& d2pvec) +{ + const PosType& r = P.activeR(iat); + RealType sinkr, coskr; + for (int ik = mink; ik < maxk; ik++) + { + sincos(dot(kvecs[ik], r), &sinkr, &coskr); + + pvec[ik] = ValueType(coskr, sinkr); + dpvec[ik] = ValueType(-sinkr, coskr) * kvecs[ik]; + d2pvec[ik] = ValueType(k2neg[ik] * coskr, k2neg[ik] * sinkr); + } +} + + +template<> +void FreeOrbitalT::evaluateValue(const ParticleSet& P, int iat, ValueVector& pvec) +{ + const PosType& r = P.activeR(iat); + RealType sinkr, coskr; + for (int ik = mink; ik < maxk; ik++) + { + sincos(dot(kvecs[ik], r), &sinkr, &coskr); + const int j2 = 2 * ik; + const int j1 = j2 - 1; + pvec[j1] = coskr; + pvec[j2] = sinkr; + } + pvec[0] = 1.0; +} + +template<> +void FreeOrbitalT::evaluateValue(const ParticleSet& P, int iat, ValueVector& pvec) +{ + const PosType& r = P.activeR(iat); + RealType sinkr, coskr; + for (int ik = mink; ik < maxk; ik++) + { + sincos(dot(kvecs[ik], r), &sinkr, &coskr); + const int j2 = 2 * ik; + const int j1 = j2 - 1; + pvec[j1] = coskr; + pvec[j2] = sinkr; + } + pvec[0] = 1.0; +} + +template<> +void FreeOrbitalT>::evaluateValue(const ParticleSet& P, int iat, ValueVector& pvec) +{ + const PosType& r = P.activeR(iat); + RealType sinkr, coskr; + for (int ik = mink; ik < maxk; ik++) + { + sincos(dot(kvecs[ik], r), &sinkr, &coskr); + + pvec[ik] = std::complex(coskr, sinkr); + const int j2 = 2 * ik; + const int j1 = j2 - 1; + pvec[j1] = coskr; + pvec[j2] = sinkr; + } +} + +template<> +void FreeOrbitalT>::evaluateValue(const ParticleSet& P, int iat, ValueVector& pvec) +{ + const PosType& r = P.activeR(iat); + RealType sinkr, coskr; + for (int ik = mink; ik < maxk; ik++) + { + sincos(dot(kvecs[ik], r), &sinkr, &coskr); + + pvec[ik] = std::complex(coskr, sinkr); + const int j2 = 2 * ik; + const int j1 = j2 - 1; + pvec[j1] = coskr; + pvec[j2] = sinkr; + } +} + +template +void FreeOrbitalT::evaluate_notranspose(const ParticleSet& P, + int first, + int last, + ValueMatrix& phi, + GradMatrix& dphi, + HessMatrix& d2phi_mat) +{} + + +template<> +void FreeOrbitalT::evaluate_notranspose(const ParticleSet& P, + int first, + int last, + ValueMatrix& phi, + GradMatrix& dphi, + HessMatrix& d2phi_mat) +{ + RealType sinkr, coskr; + float phi_of_r; + for (int iat = first, i = 0; iat < last; iat++, i++) + { + ValueVector p(phi[i], this->OrbitalSetSize); + GradVector dp(dphi[i], this->OrbitalSetSize); + HessVector hess(d2phi_mat[i], this->OrbitalSetSize); + + const PosType& r = P.activeR(iat); + for (int ik = mink; ik < maxk; ik++) + { + sincos(dot(kvecs[ik], r), &sinkr, &coskr); + const int j2 = 2 * ik; + const int j1 = j2 - 1; + p[j1] = coskr; + p[j2] = sinkr; + dp[j1] = -sinkr * kvecs[ik]; + dp[j2] = coskr * kvecs[ik]; + for (int la = 0; la < OHMMS_DIM; la++) + { + hess[j1](la, la) = -coskr * (kvecs[ik])[la] * (kvecs[ik])[la]; + hess[j2](la, la) = -sinkr * (kvecs[ik])[la] * (kvecs[ik])[la]; + for (int lb = la + 1; lb < OHMMS_DIM; lb++) + { + hess[j1](la, lb) = -coskr * (kvecs[ik])[la] * (kvecs[ik])[lb]; + hess[j2](la, lb) = -sinkr * (kvecs[ik])[la] * (kvecs[ik])[lb]; + hess[j1](lb, la) = hess[j1](la, lb); + hess[j2](lb, la) = hess[j2](la, lb); + } + } + } + p[0] = 1.0; + dp[0] = 0.0; + hess[0] = 0.0; + } +} + +template<> +void FreeOrbitalT::evaluate_notranspose(const ParticleSet& P, + int first, + int last, + ValueMatrix& phi, + GradMatrix& dphi, + HessMatrix& d2phi_mat) +{ + RealType sinkr, coskr; + double phi_of_r; + for (int iat = first, i = 0; iat < last; iat++, i++) + { + ValueVector p(phi[i], this->OrbitalSetSize); + GradVector dp(dphi[i], this->OrbitalSetSize); + HessVector hess(d2phi_mat[i], this->OrbitalSetSize); + + const PosType& r = P.activeR(iat); + for (int ik = mink; ik < maxk; ik++) + { + sincos(dot(kvecs[ik], r), &sinkr, &coskr); + const int j2 = 2 * ik; + const int j1 = j2 - 1; + p[j1] = coskr; + p[j2] = sinkr; + dp[j1] = -sinkr * kvecs[ik]; + dp[j2] = coskr * kvecs[ik]; + for (int la = 0; la < OHMMS_DIM; la++) + { + hess[j1](la, la) = -coskr * (kvecs[ik])[la] * (kvecs[ik])[la]; + hess[j2](la, la) = -sinkr * (kvecs[ik])[la] * (kvecs[ik])[la]; + for (int lb = la + 1; lb < OHMMS_DIM; lb++) + { + hess[j1](la, lb) = -coskr * (kvecs[ik])[la] * (kvecs[ik])[lb]; + hess[j2](la, lb) = -sinkr * (kvecs[ik])[la] * (kvecs[ik])[lb]; + hess[j1](lb, la) = hess[j1](la, lb); + hess[j2](lb, la) = hess[j2](la, lb); + } + } + } + p[0] = 1.0; + dp[0] = 0.0; + hess[0] = 0.0; + } +} + + +template<> +void FreeOrbitalT>::evaluate_notranspose(const ParticleSet& P, + int first, + int last, + ValueMatrix& phi, + GradMatrix& dphi, + HessMatrix& d2phi_mat) +{ + RealType sinkr, coskr; + std::complex phi_of_r; + for (int iat = first, i = 0; iat < last; iat++, i++) + { + ValueVector p(phi[i], this->OrbitalSetSize); + GradVector dp(dphi[i], this->OrbitalSetSize); + HessVector hess(d2phi_mat[i], this->OrbitalSetSize); + + const PosType& r = P.activeR(iat); + for (int ik = mink; ik < maxk; ik++) + { + sincos(dot(kvecs[ik], r), &sinkr, &coskr); + + phi_of_r = std::complex(coskr, sinkr); + p[ik] = phi_of_r; + + dp[ik] = std::complex(-sinkr, coskr) * kvecs[ik]; + for (int la = 0; la < OHMMS_DIM; la++) + { + hess[ik](la, la) = -phi_of_r * (kvecs[ik])[la] * (kvecs[ik])[la]; + for (int lb = la + 1; lb < OHMMS_DIM; lb++) + { + hess[ik](la, lb) = -phi_of_r * (kvecs[ik])[la] * (kvecs[ik])[lb]; + hess[ik](lb, la) = hess[ik](la, lb); + } + } + } + } +} + +template<> +void FreeOrbitalT>::evaluate_notranspose(const ParticleSet& P, + int first, + int last, + ValueMatrix& phi, + GradMatrix& dphi, + HessMatrix& d2phi_mat) +{ + RealType sinkr, coskr; + std::complex phi_of_r; + for (int iat = first, i = 0; iat < last; iat++, i++) + { + ValueVector p(phi[i], this->OrbitalSetSize); + GradVector dp(dphi[i], this->OrbitalSetSize); + HessVector hess(d2phi_mat[i], this->OrbitalSetSize); + + const PosType& r = P.activeR(iat); + for (int ik = mink; ik < maxk; ik++) + { + sincos(dot(kvecs[ik], r), &sinkr, &coskr); + + phi_of_r = std::complex(coskr, sinkr); + p[ik] = phi_of_r; + + dp[ik] = std::complex(-sinkr, coskr) * kvecs[ik]; + for (int la = 0; la < OHMMS_DIM; la++) + { + hess[ik](la, la) = -phi_of_r * (kvecs[ik])[la] * (kvecs[ik])[la]; + for (int lb = la + 1; lb < OHMMS_DIM; lb++) + { + hess[ik](la, lb) = -phi_of_r * (kvecs[ik])[la] * (kvecs[ik])[lb]; + hess[ik](lb, la) = hess[ik](la, lb); + } + } + } + } +} + +template +void FreeOrbitalT::evaluate_notranspose(const ParticleSet& P, + int first, + int last, + ValueMatrix& phi, + GradMatrix& dphi, + HessMatrix& d2phi_mat, + GGGMatrix& d3phi_mat) +{} + +template<> +void FreeOrbitalT::evaluate_notranspose(const ParticleSet& P, + int first, + int last, + ValueMatrix& phi, + GradMatrix& dphi, + HessMatrix& d2phi_mat, + GGGMatrix& d3phi_mat) +{ + RealType sinkr, coskr; + ValueType phi_of_r; + for (int iat = first, i = 0; iat < last; iat++, i++) + { + ValueVector p(phi[i], OrbitalSetSize); + GradVector dp(dphi[i], OrbitalSetSize); + HessVector hess(d2phi_mat[i], OrbitalSetSize); + GGGVector ggg(d3phi_mat[i], OrbitalSetSize); + + const PosType& r = P.activeR(iat); + for (int ik = mink; ik < maxk; ik++) + { + sincos(dot(kvecs[ik], r), &sinkr, &coskr); + const int j2 = 2 * ik; + const int j1 = j2 - 1; + p[j1] = coskr; + p[j2] = sinkr; + dp[j1] = -sinkr * kvecs[ik]; + dp[j2] = coskr * kvecs[ik]; + for (int la = 0; la < OHMMS_DIM; la++) + { + hess[j1](la, la) = -coskr * (kvecs[ik])[la] * (kvecs[ik])[la]; + hess[j2](la, la) = -sinkr * (kvecs[ik])[la] * (kvecs[ik])[la]; + ggg[j1][la](la, la) = sinkr * (kvecs[ik])[la] * (kvecs[ik])[la] * (kvecs[ik])[la]; + ggg[j2][la](la, la) = -coskr * (kvecs[ik])[la] * (kvecs[ik])[la] * (kvecs[ik])[la]; + for (int lb = la + 1; lb < OHMMS_DIM; lb++) + { + hess[j1](la, lb) = -coskr * (kvecs[ik])[la] * (kvecs[ik])[lb]; + hess[j2](la, lb) = -sinkr * (kvecs[ik])[la] * (kvecs[ik])[lb]; + hess[j1](lb, la) = hess[j1](la, lb); + hess[j2](lb, la) = hess[j2](la, lb); + ggg[j1][la](lb, la) = sinkr * (kvecs[ik])[la] * (kvecs[ik])[lb] * (kvecs[ik])[la]; + ggg[j2][la](lb, la) = -coskr * (kvecs[ik])[la] * (kvecs[ik])[lb] * (kvecs[ik])[la]; + ggg[j1][la](la, lb) = ggg[j1][la](lb, la); + ggg[j2][la](la, lb) = ggg[j2][la](lb, la); + ggg[j1][lb](la, la) = ggg[j1][la](lb, la); + ggg[j2][lb](la, la) = ggg[j2][la](lb, la); + ggg[j1][la](lb, lb) = sinkr * (kvecs[ik])[la] * (kvecs[ik])[lb] * (kvecs[ik])[lb]; + ggg[j2][la](lb, lb) = -coskr * (kvecs[ik])[la] * (kvecs[ik])[lb] * (kvecs[ik])[lb]; + ggg[j1][lb](la, lb) = ggg[j1][la](lb, lb); + ggg[j2][lb](la, lb) = ggg[j2][la](lb, lb); + ggg[j1][lb](lb, la) = ggg[j1][la](lb, lb); + ggg[j2][lb](lb, la) = ggg[j2][la](lb, lb); + for (int lc = lb + 1; lc < OHMMS_DIM; lc++) + { + ggg[j1][la](lb, lc) = sinkr * (kvecs[ik])[la] * (kvecs[ik])[lb] * (kvecs[ik])[lc]; + ggg[j2][la](lb, lc) = -coskr * (kvecs[ik])[la] * (kvecs[ik])[lb] * (kvecs[ik])[lc]; + ggg[j1][la](lc, lb) = ggg[j1][la](lb, lc); + ggg[j2][la](lc, lb) = ggg[j2][la](lb, lc); + ggg[j1][lb](la, lc) = ggg[j1][la](lb, lc); + ggg[j2][lb](la, lc) = ggg[j2][la](lb, lc); + ggg[j1][lb](lc, la) = ggg[j1][la](lb, lc); + ggg[j2][lb](lc, la) = ggg[j2][la](lb, lc); + ggg[j1][lc](la, lb) = ggg[j1][la](lb, lc); + ggg[j2][lc](la, lb) = ggg[j2][la](lb, lc); + ggg[j1][lc](lb, la) = ggg[j1][la](lb, lc); + ggg[j2][lc](lb, la) = ggg[j2][la](lb, lc); + } + } + } + } + + p[0] = 1.0; + dp[0] = 0.0; + hess[0] = 0.0; + ggg[0] = 0.0; + } +} + +template<> +void FreeOrbitalT::evaluate_notranspose(const ParticleSet& P, + int first, + int last, + ValueMatrix& phi, + GradMatrix& dphi, + HessMatrix& d2phi_mat, + GGGMatrix& d3phi_mat) +{ + RealType sinkr, coskr; + ValueType phi_of_r; + for (int iat = first, i = 0; iat < last; iat++, i++) + { + ValueVector p(phi[i], OrbitalSetSize); + GradVector dp(dphi[i], OrbitalSetSize); + HessVector hess(d2phi_mat[i], OrbitalSetSize); + GGGVector ggg(d3phi_mat[i], OrbitalSetSize); + + const PosType& r = P.activeR(iat); + for (int ik = mink; ik < maxk; ik++) + { + sincos(dot(kvecs[ik], r), &sinkr, &coskr); + const int j2 = 2 * ik; + const int j1 = j2 - 1; + p[j1] = coskr; + p[j2] = sinkr; + dp[j1] = -sinkr * kvecs[ik]; + dp[j2] = coskr * kvecs[ik]; + for (int la = 0; la < OHMMS_DIM; la++) + { + hess[j1](la, la) = -coskr * (kvecs[ik])[la] * (kvecs[ik])[la]; + hess[j2](la, la) = -sinkr * (kvecs[ik])[la] * (kvecs[ik])[la]; + ggg[j1][la](la, la) = sinkr * (kvecs[ik])[la] * (kvecs[ik])[la] * (kvecs[ik])[la]; + ggg[j2][la](la, la) = -coskr * (kvecs[ik])[la] * (kvecs[ik])[la] * (kvecs[ik])[la]; + for (int lb = la + 1; lb < OHMMS_DIM; lb++) + { + hess[j1](la, lb) = -coskr * (kvecs[ik])[la] * (kvecs[ik])[lb]; + hess[j2](la, lb) = -sinkr * (kvecs[ik])[la] * (kvecs[ik])[lb]; + hess[j1](lb, la) = hess[j1](la, lb); + hess[j2](lb, la) = hess[j2](la, lb); + ggg[j1][la](lb, la) = sinkr * (kvecs[ik])[la] * (kvecs[ik])[lb] * (kvecs[ik])[la]; + ggg[j2][la](lb, la) = -coskr * (kvecs[ik])[la] * (kvecs[ik])[lb] * (kvecs[ik])[la]; + ggg[j1][la](la, lb) = ggg[j1][la](lb, la); + ggg[j2][la](la, lb) = ggg[j2][la](lb, la); + ggg[j1][lb](la, la) = ggg[j1][la](lb, la); + ggg[j2][lb](la, la) = ggg[j2][la](lb, la); + ggg[j1][la](lb, lb) = sinkr * (kvecs[ik])[la] * (kvecs[ik])[lb] * (kvecs[ik])[lb]; + ggg[j2][la](lb, lb) = -coskr * (kvecs[ik])[la] * (kvecs[ik])[lb] * (kvecs[ik])[lb]; + ggg[j1][lb](la, lb) = ggg[j1][la](lb, lb); + ggg[j2][lb](la, lb) = ggg[j2][la](lb, lb); + ggg[j1][lb](lb, la) = ggg[j1][la](lb, lb); + ggg[j2][lb](lb, la) = ggg[j2][la](lb, lb); + for (int lc = lb + 1; lc < OHMMS_DIM; lc++) + { + ggg[j1][la](lb, lc) = sinkr * (kvecs[ik])[la] * (kvecs[ik])[lb] * (kvecs[ik])[lc]; + ggg[j2][la](lb, lc) = -coskr * (kvecs[ik])[la] * (kvecs[ik])[lb] * (kvecs[ik])[lc]; + ggg[j1][la](lc, lb) = ggg[j1][la](lb, lc); + ggg[j2][la](lc, lb) = ggg[j2][la](lb, lc); + ggg[j1][lb](la, lc) = ggg[j1][la](lb, lc); + ggg[j2][lb](la, lc) = ggg[j2][la](lb, lc); + ggg[j1][lb](lc, la) = ggg[j1][la](lb, lc); + ggg[j2][lb](lc, la) = ggg[j2][la](lb, lc); + ggg[j1][lc](la, lb) = ggg[j1][la](lb, lc); + ggg[j2][lc](la, lb) = ggg[j2][la](lb, lc); + ggg[j1][lc](lb, la) = ggg[j1][la](lb, lc); + ggg[j2][lc](lb, la) = ggg[j2][la](lb, lc); + } + } + } + } + + p[0] = 1.0; + dp[0] = 0.0; + hess[0] = 0.0; + ggg[0] = 0.0; + } +} + +template<> +void FreeOrbitalT>::evaluate_notranspose(const ParticleSet& P, + int first, + int last, + ValueMatrix& phi, + GradMatrix& dphi, + HessMatrix& d2phi_mat, + GGGMatrix& d3phi_mat) +{ + RealType sinkr, coskr; + ValueType phi_of_r; + for (int iat = first, i = 0; iat < last; iat++, i++) + { + ValueVector p(phi[i], OrbitalSetSize); + GradVector dp(dphi[i], OrbitalSetSize); + HessVector hess(d2phi_mat[i], OrbitalSetSize); + GGGVector ggg(d3phi_mat[i], OrbitalSetSize); + + const PosType& r = P.activeR(iat); + for (int ik = mink; ik < maxk; ik++) + { + sincos(dot(kvecs[ik], r), &sinkr, &coskr); + const ValueType compi(0, 1); + phi_of_r = ValueType(coskr, sinkr); + p[ik] = phi_of_r; + dp[ik] = compi * phi_of_r * kvecs[ik]; + for (int la = 0; la < OHMMS_DIM; la++) + { + hess[ik](la, la) = -phi_of_r * (kvecs[ik])[la] * (kvecs[ik])[la]; + for (int lb = la + 1; lb < OHMMS_DIM; lb++) + { + hess[ik](la, lb) = -phi_of_r * (kvecs[ik])[la] * (kvecs[ik])[lb]; + hess[ik](lb, la) = hess[ik](la, lb); + } + } + for (int la = 0; la < OHMMS_DIM; la++) + { + ggg[ik][la] = compi * (kvecs[ik])[la] * hess[ik]; + } + } + } +} + +template<> +void FreeOrbitalT>::evaluate_notranspose(const ParticleSet& P, + int first, + int last, + ValueMatrix& phi, + GradMatrix& dphi, + HessMatrix& d2phi_mat, + GGGMatrix& d3phi_mat) +{ + RealType sinkr, coskr; + ValueType phi_of_r; + for (int iat = first, i = 0; iat < last; iat++, i++) + { + ValueVector p(phi[i], OrbitalSetSize); + GradVector dp(dphi[i], OrbitalSetSize); + HessVector hess(d2phi_mat[i], OrbitalSetSize); + GGGVector ggg(d3phi_mat[i], OrbitalSetSize); + + const PosType& r = P.activeR(iat); + for (int ik = mink; ik < maxk; ik++) + { + sincos(dot(kvecs[ik], r), &sinkr, &coskr); + const ValueType compi(0, 1); + phi_of_r = ValueType(coskr, sinkr); + p[ik] = phi_of_r; + dp[ik] = compi * phi_of_r * kvecs[ik]; + for (int la = 0; la < OHMMS_DIM; la++) + { + hess[ik](la, la) = -phi_of_r * (kvecs[ik])[la] * (kvecs[ik])[la]; + for (int lb = la + 1; lb < OHMMS_DIM; lb++) + { + hess[ik](la, lb) = -phi_of_r * (kvecs[ik])[la] * (kvecs[ik])[lb]; + hess[ik](lb, la) = hess[ik](la, lb); + } + } + for (int la = 0; la < OHMMS_DIM; la++) + { + ggg[ik][la] = compi * (kvecs[ik])[la] * hess[ik]; + } + } + } +} + +// generic implementation + +template +FreeOrbitalT::~FreeOrbitalT() +{} + +template +void FreeOrbitalT::evaluate_notranspose(const ParticleSet& P, + int first, + int last, + ValueMatrix& phi, + GradMatrix& dphi, + ValueMatrix& d2phi) +{ + for (int iat = first, i = 0; iat < last; iat++, i++) + { + ValueVector p(phi[i], this->OrbitalSetSize); + GradVector dp(dphi[i], this->OrbitalSetSize); + ValueVector d2p(d2phi[i], this->OrbitalSetSize); + evaluateVGL(P, iat, p, dp, d2p); + } +} + + +template +void FreeOrbitalT::report(const std::string& pad) const +{ + app_log() << pad << "FreeOrbital report" << std::endl; + for (int ik = 0; ik < kvecs.size(); ik++) + { + app_log() << pad << ik << " " << kvecs[ik] << std::endl; + } + app_log() << pad << "end FreeOrbital report" << std::endl; + app_log().flush(); +} + +template class FreeOrbitalT; +template class FreeOrbitalT; +template class FreeOrbitalT>; +template class FreeOrbitalT>; + + +} // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/ElectronGas/FreeOrbitalT.h b/src/QMCWaveFunctions/ElectronGas/FreeOrbitalT.h new file mode 100644 index 0000000000..c73ab26a2a --- /dev/null +++ b/src/QMCWaveFunctions/ElectronGas/FreeOrbitalT.h @@ -0,0 +1,88 @@ +////////////////////////////////////////////////////////////////////////////////////// +// This file is distributed under the University of Illinois/NCSA Open Source License. +// See LICENSE file in top directory for details. +// +// Copyright (c) 2022 QMCPACK developers. +// +// File developed by: Miguel Morales, moralessilva2@llnl.gov, Lawrence Livermore National Laboratory +// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign +// Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign +// Jaron T. Krogel, krogeljt@ornl.gov, Oak Ridge National Laboratory +// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory +// Yubo "Paul" Yang, yubo.paul.yang@gmail.com, CCQ @ Flatiron +// William F Godoy, godoywf@ornl.gov, Oak Ridge National Laboratory +// +// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign +////////////////////////////////////////////////////////////////////////////////////// + +#ifndef QMCPLUSPLUS_FREE_ORBITAL +#define QMCPLUSPLUS_FREE_ORBITAL + +#include "QMCWaveFunctions/SPOSetT.h" + +namespace qmcplusplus +{ +template +class FreeOrbitalT : public SPOSetT +{ +public: + using ValueVector = typename SPOSetT::ValueVector; + using GradVector = typename SPOSetT::GradVector; + using HessVector = typename SPOSetT::HessVector; + using ValueMatrix = typename SPOSetT::ValueMatrix; + using GradMatrix = typename SPOSetT::GradMatrix; + using HessMatrix = typename SPOSetT::HessMatrix; + using GGGMatrix = typename SPOSetT::GGGMatrix; + using RealType = typename SPOSetT::RealType; + using PosType = typename SPOSetT::PosType; + using ValueType = typename SPOSetT::ValueType; + + FreeOrbitalT(const std::string& my_name, const std::vector& kpts_cart); + ~FreeOrbitalT(); + + inline std::string getClassName() const final { return "FreeOrbital"; } + + // phi[i][j] is phi_j(r_i), i.e. electron i in orbital j + // i \in [first, last) + void evaluate_notranspose(const ParticleSet& P, + int first, + int last, + ValueMatrix& phi, + GradMatrix& dphi, + ValueMatrix& d2phi) final; + + // plug r_i into all orbitals + void evaluateVGL(const ParticleSet& P, int i, ValueVector& pvec, GradVector& dpvec, ValueVector& d2pvec) final; + void evaluateValue(const ParticleSet& P, int iat, ValueVector& pvec) final; + + // hessian matrix is needed by backflow + void evaluate_notranspose(const ParticleSet& P, + int first, + int last, + ValueMatrix& phi, + GradMatrix& dphi, + HessMatrix& d2phi_mat) final; + + // derivative of hessian is needed to optimize backflow + void evaluate_notranspose(const ParticleSet& P, + int first, + int last, + ValueMatrix& phi, + GradMatrix& dphi, + HessMatrix& d2phi_mat, + GGGMatrix& d3phi_mat) override; + + void report(const std::string& pad) const override; + // ---- begin required overrides + std::unique_ptr> makeClone() const final { return std::make_unique>(*this); } + void setOrbitalSetSize(int norbs) final { throw std::runtime_error("not implemented"); } + // required overrides end ---- +private: + const std::vector kvecs; // kvecs vectors + const int mink; // minimum k index + const int maxk; // maximum number of kvecs vectors + std::vector k2neg; // minus kvecs^2 +}; + +} // namespace qmcplusplus +#endif diff --git a/src/QMCWaveFunctions/HarmonicOscillator/SHOSetT.cpp b/src/QMCWaveFunctions/HarmonicOscillator/SHOSetT.cpp new file mode 100644 index 0000000000..76a606151d --- /dev/null +++ b/src/QMCWaveFunctions/HarmonicOscillator/SHOSetT.cpp @@ -0,0 +1,577 @@ +////////////////////////////////////////////////////////////////////////////////////// +// This file is distributed under the University of Illinois/NCSA Open Source License. +// See LICENSE file in top directory for details. +// +// Copyright (c) 2016 Jeongnim Kim and QMCPACK developers. +// +// File developed by: Jaron T. Krogel, krogeljt@ornl.gov, Oak Ridge National Laboratory +// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory +// +// File created by: Jaron T. Krogel, krogeljt@ornl.gov, Oak Ridge National Laboratory +////////////////////////////////////////////////////////////////////////////////////// + + +#include "SHOSetT.h" +#include "Utilities/string_utils.h" + +namespace qmcplusplus +{ +template +SHOSetT::SHOSetT(const std::string& my_name, RealType l, PosType c, const std::vector& sho_states) + : SPOSetT(my_name), length(l), center(c) +{ + state_info.resize(sho_states.size()); + for (int s = 0; s < sho_states.size(); ++s) + state_info[s] = *sho_states[s]; + initialize(); +} + +template +void SHOSetT::initialize() +{ + using std::sqrt; + + this->OrbitalSetSize = state_info.size(); + + qn_max = -1; + for (int s = 0; s < state_info.size(); ++s) + for (int d = 0; d < QMCTraits::DIM; ++d) + qn_max[d] = std::max(qn_max[d], state_info[s].quantum_number[d]); + qn_max += 1; + + nmax = -1; + for (int d = 0; d < QMCTraits::DIM; ++d) + nmax = std::max(nmax, qn_max[d]); + + prefactors.resize(nmax); + hermite.resize(QMCTraits::DIM, nmax); + bvalues.resize(QMCTraits::DIM, nmax); + + if (nmax > 0) + { + prefactors[0] = 1.0 / (sqrt(sqrt(M_PI) * length)); + for (int n = 1; n < nmax; ++n) + prefactors[n] = prefactors[n - 1] / sqrt(2. * n); + } +} + +template +SHOSetT::~SHOSetT() = default; + +template +std::unique_ptr> SHOSetT::makeClone() const { return std::make_unique>(*this); } + +template +void SHOSetT::report(const std::string& pad) const +{ + app_log() << pad << "SHOSet report" << std::endl; + app_log() << pad << " length = " << length << std::endl; + app_log() << pad << " center = " << center << std::endl; + app_log() << pad << " nmax = " << nmax << std::endl; + app_log() << pad << " qn_max = " << qn_max << std::endl; + app_log() << pad << " # states = " << state_info.size() << std::endl; + app_log() << pad << " states" << std::endl; + for (int s = 0; s < state_info.size(); ++s) + state_info[s].sho_report(pad + " " + int2string(s) + " "); + app_log() << pad << "end SHOSet report" << std::endl; + app_log().flush(); +} + +template +void SHOSetT::evaluateValue(const ParticleSet& P, int iat, ValueVector& psi) +{ + const PosType& r(P.activeR(iat)); + ValueVector p(&psi[0], this->size()); + evaluate_v(r, p); +} + +template +void SHOSetT::evaluateVGL(const ParticleSet& P, int iat, ValueVector& psi, GradVector& dpsi, ValueVector& d2psi) +{ + const PosType& r(P.activeR(iat)); + ValueVector p(&psi[0], this->size()); + GradVector dp(&dpsi[0], this->size()); + ValueVector d2p(&d2psi[0], this->size()); + evaluate_vgl(r, p, dp, d2p); +} + +template +void SHOSetT::evaluate_notranspose(const ParticleSet& P, + int first, + int last, + ValueMatrix& logdet, + GradMatrix& dlogdet, + ValueMatrix& d2logdet) +{ + for (int iat = first, i = 0; iat < last; ++iat, ++i) + { + ValueVector p(logdet[i], this->size()); + GradVector dp(dlogdet[i], this->size()); + ValueVector d2p(d2logdet[i], this->size()); + evaluate_vgl(P.R[iat], p, dp, d2p); + } +} + +template +void SHOSetT::evaluate_v(PosType r, ValueVector& psi) +{ + PosType x = (r - center) / length; + evaluate_hermite(x); + evaluate_d0(x, psi); +} + +template +void SHOSetT::evaluate_vgl(PosType r, ValueVector& psi, GradVector& dpsi, ValueVector& d2psi) +{ + PosType x = (r - center) / length; + evaluate_hermite(x); + evaluate_d0(x, psi); + evaluate_d1(x, psi, dpsi); + evaluate_d2(x, psi, d2psi); +} + +template +void SHOSetT::evaluate_hermite(const PosType& xpos) +{ + for (int d = 0; d < QMCTraits::DIM; ++d) + { + int nh = qn_max[d]; + if (nh > 0) + { + RealType x = xpos[d]; + hermite(d, 0) = 1.0; + RealType Hnm2 = 0.0; + RealType Hnm1 = 1.0; + for (int n = 1; n < nh; ++n) + { + RealType Hn = 2 * (x * Hnm1 - (n - 1) * Hnm2); + hermite(d, n) = Hn; + Hnm2 = Hnm1; + Hnm1 = Hn; + } + } + } +} + +template +void SHOSetT::evaluate_d0(const PosType& xpos, ValueVector& psi) +{ + using std::exp; + for (int d = 0; d < QMCTraits::DIM; ++d) + { + RealType x = xpos[d]; + RealType g = exp(-.5 * x * x); + for (int n = 0; n < qn_max[d]; ++n) + { + bvalues(d, n) = prefactors[n] * g * hermite(d, n); + } + } + for (int s = 0; s < state_info.size(); ++s) + { + const SHOState& state = state_info[s]; + RealType phi = 1.0; + for (int d = 0; d < QMCTraits::DIM; ++d) + phi *= bvalues(d, state.quantum_number[d]); + psi[s] = phi; + } +} + +template +void SHOSetT::evaluate_d1(const PosType& xpos, ValueVector& psi, GradVector& dpsi) +{ + RealType ol = 1.0 / length; + for (int d = 0; d < QMCTraits::DIM; ++d) + { + RealType x = xpos[d]; + RealType Hnm1 = 0.0; + for (int n = 0; n < qn_max[d]; ++n) + { + RealType Hn = hermite(d, n); + bvalues(d, n) = (-x + 2 * n * Hnm1 / Hn) * ol; + Hnm1 = Hn; + } + } + for (int s = 0; s < state_info.size(); ++s) + { + const SHOState& state = state_info[s]; + TinyVector dphi; + for (int d = 0; d < QMCTraits::DIM; ++d) + dphi[d] = bvalues(d, state.quantum_number[d]); + dphi *= psi[s]; + dpsi[s] = dphi; + } +} + +template +void SHOSetT::evaluate_d2(const PosType& xpos, ValueVector& psi, ValueVector& d2psi) +{ + RealType ol2 = 1.0 / (length * length); + for (int d = 0; d < QMCTraits::DIM; ++d) + { + RealType x = xpos[d]; + RealType x2 = x * x; + for (int n = 0; n < qn_max[d]; ++n) + { + bvalues(d, n) = (-1.0 + x2 - 2 * n) * ol2; + } + } + for (int s = 0; s < state_info.size(); ++s) + { + const SHOState& state = state_info[s]; + T d2phi = 0.0; + for (int d = 0; d < QMCTraits::DIM; ++d) + d2phi += bvalues(d, state.quantum_number[d]); + d2phi *= psi[s]; + d2psi[s] = d2phi; + } +} + +template +void SHOSetT::evaluate_check(PosType r, ValueVector& psi, GradVector& dpsi, ValueVector& d2psi) +{ + using std::exp; + using std::sqrt; + + evaluate_vgl(r, psi, dpsi, d2psi); + + const int N = 6; + RealType H[N], dH[N], d2H[N], pre[N]; + RealType p[N], dp[N], d2p[N]; + + pre[0] = 1.0 / (sqrt(sqrt(M_PI) * length)); + for (int n = 1; n < N; ++n) + pre[n] = pre[n - 1] / sqrt(2. * n); + + for (int d = 0; d < QMCTraits::DIM; ++d) + { + RealType x = (r[d] - center[d]) / length; + RealType x2 = x * x, x3 = x * x * x, x4 = x * x * x * x, x5 = x * x * x * x * x; + H[0] = 1; + dH[0] = 0; + d2H[0] = 0; + H[1] = 2 * x; + dH[1] = 2; + d2H[1] = 0; + H[2] = 4 * x2 - 2; + dH[2] = 8 * x; + d2H[2] = 8; + H[3] = 8 * x3 - 12 * x; + dH[3] = 24 * x2 - 12; + d2H[3] = 48 * x; + H[4] = 16 * x4 - 48 * x2 + 12; + dH[4] = 64 * x3 - 96 * x; + d2H[4] = 192 * x2 - 96; + H[5] = 32 * x5 - 160 * x3 + 120 * x; + dH[5] = 160 * x4 - 480 * x2 + 120; + d2H[5] = 640 * x3 - 960 * x; + RealType g = exp(-x2 / 2); + for (int n = 0; n < N; ++n) + { + p[n] = pre[n] * g * H[n]; + dp[n] = pre[n] * g * (-x * H[n] + dH[n]); + d2p[n] = pre[n] * g * ((x2 - 1) * H[n] - 2 * x * dH[n] + d2H[n]); + } + app_log() << "eval check dim = " << d << " x = " << x << std::endl; + app_log() << " hermite check" << std::endl; + for (int n = 0; n < qn_max[d]; ++n) + { + app_log() << " " << n << " " << H[n] << std::endl; + app_log() << " " << n << " " << hermite(d, n) << std::endl; + } + app_log() << " phi d0 check" << std::endl; + for (int n = 0; n < qn_max[d]; ++n) + { + app_log() << " " << n << " " << p[n] << std::endl; + app_log() << " " << n << " " << d0_values(d, n) << std::endl; + } + app_log() << " phi d1 check" << std::endl; + for (int n = 0; n < qn_max[d]; ++n) + { + app_log() << " " << n << " " << dp[n] / p[n] << std::endl; + app_log() << " " << n << " " << d1_values(d, n) << std::endl; + } + app_log() << " phi d2 check" << std::endl; + for (int n = 0; n < qn_max[d]; ++n) + { + app_log() << " " << n << " " << d2p[n] / p[n] << std::endl; + app_log() << " " << n << " " << d2_values(d, n) << std::endl; + } + } +} + +template +void SHOSetT::test_derivatives() +{ + int n = 3; + PosType c = 5.123; + PosType L = 1.0; + PosType drg = L / n; + PosType dr = L / 1000; + int nphi = state_info.size(); + + PosType o2dr, odr2; + + ValueVector vpsi, vpsitmp; + GradVector vdpsi, vdpsin; + ValueVector vd2psi, vd2psin; + + + vpsi.resize(nphi); + vdpsi.resize(nphi); + vd2psi.resize(nphi); + + vpsitmp.resize(nphi); + vdpsin.resize(nphi); + vd2psin.resize(nphi); + + + ValueVector psi(&vpsi[0], this->size()); + GradVector dpsi(&vdpsi[0], this->size()); + ValueVector d2psi(&vd2psi[0], this->size()); + + ValueVector psitmp(&vpsitmp[0], this->size()); + GradVector dpsin(&vdpsin[0], this->size()); + ValueVector d2psin(&vd2psin[0], this->size()); + + + app_log() << " loading dr" << std::endl; + + RealType odr2sum = 0.0; + for (int d = 0; d < QMCTraits::DIM; ++d) + { + RealType odr = 1.0 / dr[d]; + o2dr[d] = .5 * odr; + odr2[d] = odr * odr; + odr2sum += odr2[d]; + } + + app_log() << "SHOSet::test_derivatives" << std::endl; + + const SimulationCell simulation_cell; + ParticleSet Ps(simulation_cell); + + int p = 0; + PosType r, rtmp; + for (int i = 0; i < n; ++i) + { + r[0] = c[0] + i * drg[0]; + for (int j = 0; j < n; ++j) + { + r[1] = c[1] + j * drg[1]; + for (int k = 0; k < n; ++k) + { + r[2] = c[2] + k * drg[2]; + + evaluate_vgl(r, psi, dpsi, d2psi); + + for (int m = 0; m < nphi; ++m) + d2psin[m] = -2 * odr2sum * psi[m]; + for (int d = 0; d < QMCTraits::DIM; ++d) + { + rtmp = r; + rtmp[d] += dr[d]; + evaluate_v(rtmp, psitmp); + for (int m = 0; m < nphi; ++m) + { + T phi = psitmp[m]; + dpsin[m][d] = phi * o2dr[d]; + d2psin[m] += phi * odr2[d]; + } + rtmp = r; + rtmp[d] -= dr[d]; + evaluate_v(rtmp, psitmp); + for (int m = 0; m < nphi; ++m) + { + T phi = psitmp[m]; + dpsin[m][d] -= phi * o2dr[d]; + d2psin[m] += phi * odr2[d]; + } + } + + RealType dphi_diff = 0.0; + RealType d2phi_diff = 0.0; + for (int m = 0; m < nphi; ++m) + for (int d = 0; d < QMCTraits::DIM; ++d) + dphi_diff = std::max(dphi_diff, std::abs(dpsi[m][d] - dpsin[m][d]) / std::abs(dpsin[m][d])); + for (int m = 0; m < nphi; ++m) + d2phi_diff = std::max(d2phi_diff, std::abs(d2psi[m] - d2psin[m]) / std::abs(d2psin[m])); + app_log() << " " << p << " " << dphi_diff << " " << d2phi_diff << std::endl; + app_log() << " derivatives" << std::endl; + for (int m = 0; m < nphi; ++m) + { + std::string qn = ""; + for (int d = 0; d < QMCTraits::DIM; ++d) + qn += int2string(state_info[m].quantum_number[d]) + " "; + app_log() << " " << qn; + for (int d = 0; d < QMCTraits::DIM; ++d) + app_log() << real(dpsi[m][d]) << " "; + app_log() << std::endl; + app_log() << " " << qn; + for (int d = 0; d < QMCTraits::DIM; ++d) + app_log() << real(dpsin[m][d]) << " "; + app_log() << std::endl; + } + app_log() << " laplacians" << std::endl; + PosType x = r / length; + for (int m = 0; m < nphi; ++m) + { + std::string qn = ""; + for (int d = 0; d < QMCTraits::DIM; ++d) + qn += int2string(state_info[m].quantum_number[d]) + " "; + app_log() << " " << qn << real(d2psi[m] / psi[m]) << std::endl; + app_log() << " " << qn << real(d2psin[m] / psi[m]) << std::endl; + } + p++; + } + } + } + + app_log() << "end SHOSet::test_derivatives" << std::endl; +} + +template +void SHOSetT::test_overlap() +{ + app_log() << "SHOSet::test_overlap" << std::endl; + + + //linear + int d = 0; + + app_log() << " length = " << length << std::endl; + app_log() << " prefactors" << std::endl; + for (int n = 0; n < qn_max[d]; ++n) + app_log() << " " << n << " " << prefactors[n] << std::endl; + + app_log() << " 1d overlap" << std::endl; + + ValueVector vpsi; + vpsi.resize(this->size()); + ValueVector psi(&vpsi[0], this->size()); + + double xmax = 4.0; + double dx = .1; + double dr = length * dx; + + int nphi = qn_max[d]; + Array omat; + omat.resize(nphi, nphi); + for (int i = 0; i < nphi; ++i) + for (int j = 0; j < nphi; ++j) + omat(i, j) = 0.0; + + PosType xp = 0.0; + for (double x = -xmax; x < xmax; x += dx) + { + xp[d] = x; + evaluate_hermite(xp); + evaluate_d0(xp, psi); + + for (int i = 0; i < nphi; ++i) + for (int j = 0; j < nphi; ++j) + omat(i, j) += bvalues(d, i) * bvalues(d, j) * dr; + } + + for (int i = 0; i < nphi; ++i) + { + app_log() << std::endl; + for (int j = 0; j < nphi; ++j) + app_log() << omat(i, j) << " "; + } + app_log() << std::endl; + + + //volumetric + app_log() << " 3d overlap" << std::endl; + double dV = dr * dr * dr; + nphi = this->size(); + omat.resize(nphi, nphi); + for (int i = 0; i < nphi; ++i) + for (int j = 0; j < nphi; ++j) + omat(i, j) = 0.0; + for (double x = -xmax; x < xmax; x += dx) + for (double y = -xmax; y < xmax; y += dx) + for (double z = -xmax; z < xmax; z += dx) + { + xp[0] = x; + xp[1] = y; + xp[2] = z; + evaluate_hermite(xp); + evaluate_d0(xp, psi); + + for (int i = 0; i < nphi; ++i) + for (int j = 0; j < nphi; ++j) + omat(i, j) += std::abs(psi[i] * psi[j]) * dV; + } + for (int i = 0; i < nphi; ++i) + { + app_log() << std::endl; + for (int j = 0; j < nphi; ++j) + app_log() << omat(i, j) << " "; + } + app_log() << std::endl; + + + app_log() << "end SHOSet::test_overlap" << std::endl; +} + +template +void SHOSetT::evaluateThirdDeriv(const ParticleSet& P, int first, int last, GGGMatrix& grad_grad_grad_logdet) +{ + not_implemented("evaluateThirdDeriv(P,first,last,dddlogdet)"); +} + +template +void SHOSetT::evaluate_notranspose(const ParticleSet& P, + int first, + int last, + ValueMatrix& logdet, + GradMatrix& dlogdet, + HessMatrix& grad_grad_logdet) +{ + not_implemented("evaluate_notranspose(P,first,last,logdet,dlogdet,ddlogdet)"); +} + +template +void SHOSetT::evaluate_notranspose(const ParticleSet& P, + int first, + int last, + ValueMatrix& logdet, + GradMatrix& dlogdet, + HessMatrix& grad_grad_logdet, + GGGMatrix& grad_grad_grad_logdet) +{ + not_implemented("evaluate_notranspose(P,first,last,logdet,dlogdet,ddlogdet,dddlogdet)"); +} + +template +void SHOSetT::evaluateGradSource(const ParticleSet& P, + int first, + int last, + const ParticleSet& source, + int iat_src, + GradMatrix& gradphi) +{ + not_implemented("evaluateGradSource(P,first,last,source,iat,dphi)"); +} + +template +void SHOSetT::evaluateGradSource(const ParticleSet& P, + int first, + int last, + const ParticleSet& source, + int iat_src, + GradMatrix& grad_phi, + HessMatrix& grad_grad_phi, + GradMatrix& grad_lapl_phi) +{ + not_implemented("evaluateGradSource(P,first,last,source,iat,dphi,ddphi,dd2phi)"); +} + +// Class concrete types from ValueType +template class SHOSetT; +template class SHOSetT; +template class SHOSetT>; +template class SHOSetT>; + +} // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/HarmonicOscillator/SHOSetT.h b/src/QMCWaveFunctions/HarmonicOscillator/SHOSetT.h new file mode 100644 index 0000000000..bd4870a63c --- /dev/null +++ b/src/QMCWaveFunctions/HarmonicOscillator/SHOSetT.h @@ -0,0 +1,158 @@ +////////////////////////////////////////////////////////////////////////////////////// +// This file is distributed under the University of Illinois/NCSA Open Source License. +// See LICENSE file in top directory for details. +// +// Copyright (c) 2016 Jeongnim Kim and QMCPACK developers. +// +// File developed by: Jaron T. Krogel, krogeljt@ornl.gov, Oak Ridge National Laboratory +// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory +// +// File created by: Jaron T. Krogel, krogeljt@ornl.gov, Oak Ridge National Laboratory +////////////////////////////////////////////////////////////////////////////////////// + + +#ifndef QMCPLUSPLUS_SHOSET_H +#define QMCPLUSPLUS_SHOSET_H + +#include "QMCWaveFunctions/SPOSetT.h" +#include "QMCWaveFunctions/SPOInfo.h" + +namespace qmcplusplus +{ +struct SHOState : public SPOInfo +{ + TinyVector quantum_number; + + SHOState() + { + quantum_number = -1; + energy = 0.0; + } + + ~SHOState() override {} + + inline void set(TinyVector qn, RealType e) + { + quantum_number = qn; + energy = e; + } + + inline void sho_report(const std::string& pad = "") const + { + app_log() << pad << "qn=" << quantum_number << " e=" << energy << std::endl; + } +}; + +template +class SHOSetT : public SPOSetT +{ +public: + using GradVector = typename SPOSetT::GradVector; + using ValueVector = typename SPOSetT::ValueVector; + using ValueMatrix = typename SPOSetT::ValueMatrix; + using GradMatrix = typename SPOSetT::GradMatrix; + using value_type = typename ValueMatrix::value_type; + using grad_type = typename GradMatrix::value_type; + using RealType = typename SPOSetT::RealType; + using PosType = TinyVector; + using HessType = typename OrbitalSetTraits::HessType; + using HessMatrix = typename OrbitalSetTraits::HessMatrix; + using GGGType = TinyVector; + using GGGVector = Vector; + using GGGMatrix = Matrix; + + RealType length; + PosType center; + + int nmax; + TinyVector qn_max; + std::vector state_info; + std::vector prefactors; + Array hermite; + Array bvalues; + Array d0_values; + Array d1_values; + Array d2_values; + + //construction/destruction + SHOSetT(const std::string& my_name, RealType l, PosType c, const std::vector& sho_states); + + ~SHOSetT() override; + + std::string getClassName() const override { return "SHOSet"; } + + void initialize(); + + //SPOSet interface methods + std::unique_ptr> makeClone() const override; + + void evaluateValue(const ParticleSet& P, int iat, ValueVector& psi) override; + + void evaluateVGL(const ParticleSet& P, int iat, ValueVector& psi, GradVector& dpsi, ValueVector& d2psi) override; + + void evaluate_notranspose(const ParticleSet& P, + int first, + int last, + ValueMatrix& logdet, + GradMatrix& dlogdet, + ValueMatrix& d2logdet) override; + + + //local functions + void evaluate_v(PosType r, ValueVector& psi); + void evaluate_vgl(PosType r, ValueVector& psi, GradVector& dpsi, ValueVector& d2psi); + void evaluate_hermite(const PosType& xpos); + void evaluate_d0(const PosType& xpos, ValueVector& psi); + void evaluate_d1(const PosType& xpos, ValueVector& psi, GradVector& dpsi); + void evaluate_d2(const PosType& xpos, ValueVector& psi, ValueVector& d2psi); + void report(const std::string& pad = "") const override; + void test_derivatives(); + void test_overlap(); + void evaluate_check(PosType r, ValueVector& psi, GradVector& dpsi, ValueVector& d2psi); + + //empty methods + /// number of orbitals is determined only by initial request + inline void setOrbitalSetSize(int norbs) override {} + + ///unimplemented functions call this to abort + inline void not_implemented(const std::string& method) + { + APP_ABORT("SHOSet::" + method + " has not been implemented."); + } + + + //methods to be implemented in the future (possibly) + void evaluateThirdDeriv(const ParticleSet& P, int first, int last, GGGMatrix& dddlogdet) override; + void evaluate_notranspose(const ParticleSet& P, + int first, + int last, + ValueMatrix& logdet, + GradMatrix& dlogdet, + HessMatrix& ddlogdet) override; + void evaluate_notranspose(const ParticleSet& P, + int first, + int last, + ValueMatrix& logdet, + GradMatrix& dlogdet, + HessMatrix& ddlogdet, + GGGMatrix& dddlogdet) override; + void evaluateGradSource(const ParticleSet& P, + int first, + int last, + const ParticleSet& source, + int iat_src, + GradMatrix& gradphi) override; + void evaluateGradSource(const ParticleSet& P, + int first, + int last, + const ParticleSet& source, + int iat_src, + GradMatrix& dphi, + HessMatrix& ddphi, + GradMatrix& dlapl_phi) override; +}; + +} // namespace qmcplusplus + + +#endif diff --git a/src/QMCWaveFunctions/PlaneWave/PWRealOrbitalSetT.cpp b/src/QMCWaveFunctions/PlaneWave/PWRealOrbitalSetT.cpp new file mode 100644 index 0000000000..3286624090 --- /dev/null +++ b/src/QMCWaveFunctions/PlaneWave/PWRealOrbitalSetT.cpp @@ -0,0 +1,165 @@ +////////////////////////////////////////////////////////////////////////////////////// +// This file is distributed under the University of Illinois/NCSA Open Source License. +// See LICENSE file in top directory for details. +// +// Copyright (c) 2016 Jeongnim Kim and QMCPACK developers. +// +// File developed by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign +// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign +// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory +// Mark Dewing, markdewing@gmail.com, University of Illinois at Urbana-Champaign +// +// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign +////////////////////////////////////////////////////////////////////////////////////// + + +/** @file PWRealOrbitalSetT.cpp + * @brief declaration of the member functions of PWRealOrbitalSetT + * + * Not the most optimized method to use wavefunctions in a plane-wave basis. + */ +#include "Message/Communicate.h" +#include "PWRealOrbitalSetT.h" +#include "Numerics/MatrixOperators.h" +#include "type_traits/ConvertToReal.h" + +namespace qmcplusplus +{ +template +PWRealOrbitalSetT::~PWRealOrbitalSetT() +{ + if (OwnBasisSet && myBasisSet) + delete myBasisSet; +} + +template +std::unique_ptr> PWRealOrbitalSetT::makeClone() const +{ + auto myclone = std::make_unique>(*this); + myclone->myBasisSet = new PWBasis(*(this->myBasisSet)); + return myclone; +} + +template +void PWRealOrbitalSetT::setOrbitalSetSize(int norbs) +{} + +template +void PWRealOrbitalSetT::resize(PWBasisPtr bset, int nbands, bool cleanup) +{ + myBasisSet = bset; + this->OrbitalSetSize = nbands; + OwnBasisSet = cleanup; + BasisSetSize = myBasisSet->NumPlaneWaves; + CC.resize(this->OrbitalSetSize, BasisSetSize); + Temp.resize(this->OrbitalSetSize, PW_MAXINDEX); + tempPsi.resize(this->OrbitalSetSize); + app_log() << " PWRealOrbitalSetT::resize OrbitalSetSize =" << this->OrbitalSetSize + << " BasisSetSize = " << BasisSetSize << std::endl; +} + +template +void PWRealOrbitalSetT::addVector(const std::vector& coefs, int jorb) +{ + int ng = myBasisSet->inputmap.size(); + if (ng != coefs.size()) + { + app_error() << " Input G map does not match the basis size of wave functions " << std::endl; + OHMMS::Controller->abort(); + } + //drop G points for the given TwistAngle + const std::vector& inputmap(myBasisSet->inputmap); + for (int ig = 0; ig < ng; ig++) + { + if (inputmap[ig] > -1) + CC[jorb][inputmap[ig]] = coefs[ig]; + } +} + +template +void PWRealOrbitalSetT::addVector(const std::vector& coefs, int jorb) +{ + int ng = myBasisSet->inputmap.size(); + if (ng != coefs.size()) + { + app_error() << " Input G map does not match the basis size of wave functions " << std::endl; + OHMMS::Controller->abort(); + } + //drop G points for the given TwistAngle + const std::vector& inputmap(myBasisSet->inputmap); + for (int ig = 0; ig < ng; ig++) + { + if (inputmap[ig] > -1) + CC[jorb][inputmap[ig]] = coefs[ig]; + } +} + +template +void PWRealOrbitalSetT::evaluateValue(const ParticleSet& P, int iat, ValueVector& psi) +{ + myBasisSet->evaluate(P.activeR(iat)); + MatrixOperators::product(CC, myBasisSet->Zv, tempPsi); + for (int j = 0; j < this->OrbitalSetSize; j++) + psi[j] = tempPsi[j].real(); +} + +template +void PWRealOrbitalSetT::evaluateVGL(const ParticleSet& P, + int iat, + ValueVector& psi, + GradVector& dpsi, + ValueVector& d2psi) +{ + myBasisSet->evaluateAll(P, iat); + MatrixOperators::product(CC, myBasisSet->Z, Temp); + const ComplexType* restrict tptr = Temp.data(); + for (int j = 0; j < this->OrbitalSetSize; j++, tptr += PW_MAXINDEX) + { + psi[j] = tptr[PW_VALUE].real(); + d2psi[j] = tptr[PW_LAP].real(); +#if OHMMS_DIM == 3 + dpsi[j] = GradType(tptr[PW_GRADX].real(), tptr[PW_GRADY].real(), tptr[PW_GRADZ].real()); +#elif OHMMS_DIM == 2 + dpsi[j] = GradType(tptr[PW_GRADX].real(), tptr[PW_GRADY].real()); +#elif OHMMS_DIM == 1 + dpsi[j] = GradType(tptr[PW_GRADX].real()); +#else +#error "Only physical dimensions 1/2/3 are supported." +#endif + } +} + +template +void PWRealOrbitalSetT::evaluate_notranspose(const ParticleSet& P, + int first, + int last, + ValueMatrix& logdet, + GradMatrix& dlogdet, + ValueMatrix& d2logdet) +{ + for (int iat = first, i = 0; iat < last; iat++, i++) + { + myBasisSet->evaluateAll(P, iat); + MatrixOperators::product(CC, myBasisSet->Z, Temp); + const ComplexType* restrict tptr = Temp.data(); + for (int j = 0; j < this->OrbitalSetSize; j++, tptr += PW_MAXINDEX) + { + convertToReal(tptr[PW_VALUE], logdet(i, j)); + convertToReal(tptr[PW_LAP], d2logdet(i, j)); +#if OHMMS_DIM == 3 + dlogdet(i, j) = GradType(tptr[PW_GRADX].real(), tptr[PW_GRADY].real(), tptr[PW_GRADZ].real()); +#elif OHMMS_DIM == 2 + dlogdet(i, j) = GradType(tptr[PW_GRADX].real(), tptr[PW_GRADY].real()); +#elif OHMMS_DIM == 1 + dlogdet(i, j) = GradType(tptr[PW_GRADX].real()); +#else +#error "Only physical dimensions 1/2/3 are supported." +#endif + } + } +} + +template class SPOSetT; +template class SPOSetT; + +} // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/PlaneWave/PWRealOrbitalSetT.h b/src/QMCWaveFunctions/PlaneWave/PWRealOrbitalSetT.h new file mode 100644 index 0000000000..29e484f3ff --- /dev/null +++ b/src/QMCWaveFunctions/PlaneWave/PWRealOrbitalSetT.h @@ -0,0 +1,143 @@ +////////////////////////////////////////////////////////////////////////////////////// +// This file is distributed under the University of Illinois/NCSA Open Source License. +// See LICENSE file in top directory for details. +// +// Copyright (c) 2016 Jeongnim Kim and QMCPACK developers. +// +// File developed by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign +// Miguel Morales, moralessilva2@llnl.gov, Lawrence Livermore National Laboratory +// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign +// Mark Dewing, markdewing@gmail.com, University of Illinois at Urbana-Champaign +// +// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign +////////////////////////////////////////////////////////////////////////////////////// + + +/** @file PWRealOrbitalSetT.h + * @brief Define PWRealOrbitalSetT derived from SPOSetT + * + * This is a specialized single-particle orbital set for real trial + * wavefunctions and enabled with QMC_COMPLEX=0 + */ +#ifndef QMCPLUSPLUS_PLANEWAVE_REALORBITALSET_BLAS_H +#define QMCPLUSPLUS_PLANEWAVE_REALORBITALSET_BLAS_H + +#include "QMCWaveFunctions/PlaneWave/PWBasis.h" +#include "QMCWaveFunctions/SPOSetT.h" +#include "CPU/BLAS.hpp" + +namespace qmcplusplus +{ +template +class PWRealOrbitalSetT : public SPOSetT +{ +public: + using BasisSet_t = PWBasis; + using PWBasisPtr = PWBasis*; + + using IndexType = typename SPOSetT::IndexType; + using RealType = typename SPOSetT::RealType; + using ComplexType = typename SPOSetT::ComplexType; + using ValueVector = typename SPOSetT::ValueVector; + using ValueMatrix = typename SPOSetT::ValueMatrix; + using GradVector = typename SPOSetT::GradVector; + using GradMatrix = typename SPOSetT::GradMatrix; + using HessMatrix = typename SPOSetT::HessMatrix; + using PosType = typename SPOSetT::PosType; + + /** inherit the enum of BasisSet_t */ + enum + { + PW_VALUE = BasisSet_t::PW_VALUE, + PW_LAP = BasisSet_t::PW_LAP, + PW_GRADX = BasisSet_t::PW_GRADX, + PW_GRADY = BasisSet_t::PW_GRADY, + PW_GRADZ = BasisSet_t::PW_GRADZ, + PW_MAXINDEX = BasisSet_t::PW_MAXINDEX + }; + + /** default constructor + */ + PWRealOrbitalSetT(const std::string& my_name) + : SPOSetT(my_name), OwnBasisSet(false), myBasisSet(nullptr), BasisSetSize(0) + {} + + std::string getClassName() const override { return "PWRealOrbitalSetT"; } + + /** delete BasisSet only it owns this + * + * Builder takes care of who owns what + */ + ~PWRealOrbitalSetT() override; + + std::unique_ptr> makeClone() const override; + + /** resize the orbital base + * @param bset PWBasis + * @param nbands number of bands + * @param cleaup if true, owns PWBasis. Will clean up. + */ + void resize(PWBasisPtr bset, int nbands, bool cleanup = false); + + /** add eigenstate for jorb-th orbital + * @param coefs real input data + * @param jorb orbital index + */ + void addVector(const std::vector& coefs, int jorb); + + /** add eigenstate for jorb-th orbital + * @param coefs complex input data + * @param jorb orbital index + */ + void addVector(const std::vector& coefs, int jorb); + + void setOrbitalSetSize(int norbs) override; + + inline T evaluate(int ib, const PosType& pos) + { + myBasisSet->evaluate(pos); + return real(BLAS::dot(BasisSetSize, CC[ib], myBasisSet->Zv.data())); + } + + void evaluateValue(const ParticleSet& P, int iat, ValueVector& psi) override; + + void evaluateVGL(const ParticleSet& P, int iat, ValueVector& psi, GradVector& dpsi, ValueVector& d2psi) override; + + void evaluate_notranspose(const ParticleSet& P, + int first, + int last, + ValueMatrix& logdet, + GradMatrix& dlogdet, + ValueMatrix& d2logdet) override; + + void evaluate_notranspose(const ParticleSet& P, + int first, + int last, + ValueMatrix& logdet, + GradMatrix& dlogdet, + HessMatrix& grad_grad_logdet) override + { + APP_ABORT("Need specialization of evaluate_notranspose() for grad_grad_logdet. \n"); + } + + + /** boolean + * + * If true, this has to delete the BasisSet + */ + bool OwnBasisSet; + ///TwistAngle of this PWRealOrbitalSet + PosType TwistAngle; + ///My basis set + PWBasisPtr myBasisSet; + ///number of basis + IndexType BasisSetSize; + ///Plane-wave coefficients of complex: (iband,g-vector) + Matrix CC; + /// temporary array to perform gemm operation + Matrix Temp; + ///temporary complex vector before assigning to a real psi + Vector tempPsi; +}; +} // namespace qmcplusplus +#endif diff --git a/src/QMCWaveFunctions/RotatedSPOsT.cpp b/src/QMCWaveFunctions/RotatedSPOsT.cpp new file mode 100644 index 0000000000..5a992ebce8 --- /dev/null +++ b/src/QMCWaveFunctions/RotatedSPOsT.cpp @@ -0,0 +1,1690 @@ +////////////////////////////////////////////////////////////////////////////////////// +//// This file is distributed under the University of Illinois/NCSA Open Source +/// License. / See LICENSE file in top directory for details. +//// +//// Copyright (c) QMCPACK developers. +//// +//// File developed by: Sergio D. Pineda Flores, +/// sergio_pinedaflores@berkeley.edu, University of California, Berkeley / Eric +/// Neuscamman, eneuscamman@berkeley.edu, University of California, Berkeley / +/// Ye Luo, yeluo@anl.gov, Argonne National Laboratory +//// +//// File created by: Sergio D. Pineda Flores, sergio_pinedaflores@berkeley.edu, +/// University of California, Berkeley +//////////////////////////////////////////////////////////////////////////////////////// +#include "RotatedSPOsT.h" + +#include "CPU/BLAS.hpp" +#include "Numerics/DeterminantOperators.h" +#include "Numerics/MatrixOperators.h" +#include "io/hdf/hdf_archive.h" + +namespace qmcplusplus +{ +template +RotatedSPOsT::RotatedSPOsT(const std::string& my_name, std::unique_ptr>&& spos) + : SPOSetT(my_name), OptimizableObject(my_name), Phi(std::move(spos)), nel_major_(0), params_supplied(false) +{ + this->OrbitalSetSize = Phi->getOrbitalSetSize(); +} + +template +RotatedSPOsT::~RotatedSPOsT() +{} + +template +void RotatedSPOsT::setRotationParameters(const std::vector& param_list) +{ + params = param_list; + params_supplied = true; +} + +template +void RotatedSPOsT::createRotationIndices(int nel, int nmo, RotationIndices& rot_indices) +{ + for (int i = 0; i < nel; i++) + for (int j = nel; j < nmo; j++) + rot_indices.emplace_back(i, j); +} + +template +void RotatedSPOsT::createRotationIndicesFull(int nel, int nmo, RotationIndices& rot_indices) +{ + rot_indices.reserve(nmo * (nmo - 1) / 2); + + // start with core-active rotations - put them at the beginning of the list + // so it matches the other list of rotation indices + for (int i = 0; i < nel; i++) + for (int j = nel; j < nmo; j++) + rot_indices.emplace_back(i, j); + + // Add core-core rotations - put them at the end of the list + for (int i = 0; i < nel; i++) + for (int j = i + 1; j < nel; j++) + rot_indices.emplace_back(i, j); + + // Add active-active rotations - put them at the end of the list + for (int i = nel; i < nmo; i++) + for (int j = i + 1; j < nmo; j++) + rot_indices.emplace_back(i, j); +} + +template +void RotatedSPOsT::constructAntiSymmetricMatrix(const RotationIndices& rot_indices, + const std::vector& param, + ValueMatrix& rot_mat) +{ + assert(rot_indices.size() == param.size()); + // Assumes rot_mat is of the correct size + + rot_mat = 0.0; + + for (int i = 0; i < rot_indices.size(); i++) + { + const int p = rot_indices[i].first; + const int q = rot_indices[i].second; + const RealType x = param[i]; + + rot_mat[q][p] = x; + rot_mat[p][q] = -x; + } +} + +template +void RotatedSPOsT::extractParamsFromAntiSymmetricMatrix(const RotationIndices& rot_indices, + const ValueMatrix& rot_mat, + std::vector& param) +{ + assert(rot_indices.size() == param.size()); + // Assumes rot_mat is of the correct size + + for (int i = 0; i < rot_indices.size(); i++) + { + const int p = rot_indices[i].first; + const int q = rot_indices[i].second; + param[i] = rot_mat[q][p]; + } +} + +template +void RotatedSPOsT::resetParametersExclusive(const opt_variables_type& active) +{ + std::vector delta_param(m_act_rot_inds.size()); + + size_t psize = m_act_rot_inds.size(); + + if (use_global_rot_) + { + psize = m_full_rot_inds.size(); + assert(psize >= m_act_rot_inds.size()); + } + + std::vector old_param(psize); + std::vector new_param(psize); + + for (int i = 0; i < m_act_rot_inds.size(); i++) + { + int loc = this->myVars.where(i); + delta_param[i] = active[loc] - this->myVars[i]; + this->myVars[i] = active[loc]; + } + + if (use_global_rot_) + { + for (int i = 0; i < m_full_rot_inds.size(); i++) + old_param[i] = myVarsFull[i]; + + applyDeltaRotation(delta_param, old_param, new_param); + + // Save the the params + for (int i = 0; i < m_full_rot_inds.size(); i++) + myVarsFull[i] = new_param[i]; + } + else + { + apply_rotation(delta_param, false); + + // Save the parameters in the history list + history_params_.push_back(delta_param); + } +} + +template +void RotatedSPOsT::writeVariationalParameters(hdf_archive& hout) +{ + hout.push("RotatedSPOsT"); + if (use_global_rot_) + { + hout.push("rotation_global"); + std::string rot_global_name = std::string("rotation_global_") + SPOSetT::getName(); + + int nparam_full = myVarsFull.size(); + std::vector full_params(nparam_full); + for (int i = 0; i < nparam_full; i++) + full_params[i] = myVarsFull[i]; + + hout.write(full_params, rot_global_name); + hout.pop(); + } + else + { + hout.push("rotation_history"); + size_t rows = history_params_.size(); + size_t cols = 0; + if (rows > 0) + cols = history_params_[0].size(); + + Matrix tmp(rows, cols); + for (size_t i = 0; i < rows; i++) + for (size_t j = 0; j < cols; j++) + tmp(i, j) = history_params_[i][j]; + + std::string rot_hist_name = std::string("rotation_history_") + SPOSetT::getName(); + hout.write(tmp, rot_hist_name); + hout.pop(); + } + + // Save myVars in order to restore object state exactly + // The values aren't meaningful, but they need to match those saved in + // VariableSet + hout.push("rotation_params"); + std::string rot_params_name = std::string("rotation_params_") + SPOSetT::getName(); + + int nparam = this->myVars.size(); + std::vector params(nparam); + for (int i = 0; i < nparam; i++) + params[i] = this->myVars[i]; + + hout.write(params, rot_params_name); + hout.pop(); + + hout.pop(); +} + +template +void RotatedSPOsT::readVariationalParameters(hdf_archive& hin) +{ + hin.push("RotatedSPOsT", false); + + bool grp_hist_exists = hin.is_group("rotation_history"); + bool grp_global_exists = hin.is_group("rotation_global"); + if (!grp_hist_exists && !grp_global_exists) + app_warning() << "Rotation parameters not found in VP file"; + + if (grp_global_exists) + { + hin.push("rotation_global", false); + std::string rot_global_name = std::string("rotation_global_") + SPOSetT::getName(); + + std::vector sizes(1); + if (!hin.getShape(rot_global_name, sizes)) + throw std::runtime_error("Failed to read rotation_global in VP file"); + + int nparam_full_actual = sizes[0]; + int nparam_full = myVarsFull.size(); + + if (nparam_full != nparam_full_actual) + { + std::ostringstream tmp_err; + tmp_err << "Expected number of full rotation parameters (" << nparam_full << ") does not match number in file (" + << nparam_full_actual << ")"; + throw std::runtime_error(tmp_err.str()); + } + std::vector full_params(nparam_full); + hin.read(full_params, rot_global_name); + for (int i = 0; i < nparam_full; i++) + myVarsFull[i] = full_params[i]; + + hin.pop(); + + applyFullRotation(full_params, true); + } + else if (grp_hist_exists) + { + hin.push("rotation_history", false); + std::string rot_hist_name = std::string("rotation_history_") + SPOSetT::getName(); + std::vector sizes(2); + if (!hin.getShape(rot_hist_name, sizes)) + throw std::runtime_error("Failed to read rotation history in VP file"); + + int rows = sizes[0]; + int cols = sizes[1]; + history_params_.resize(rows); + Matrix tmp(rows, cols); + hin.read(tmp, rot_hist_name); + for (size_t i = 0; i < rows; i++) + { + history_params_[i].resize(cols); + for (size_t j = 0; j < cols; j++) + history_params_[i][j] = tmp(i, j); + } + + hin.pop(); + + applyRotationHistory(); + } + + hin.push("rotation_params", false); + std::string rot_param_name = std::string("rotation_params_") + SPOSetT::getName(); + + std::vector sizes(1); + if (!hin.getShape(rot_param_name, sizes)) + throw std::runtime_error("Failed to read rotation_params in VP file"); + + int nparam_actual = sizes[0]; + int nparam = this->myVars.size(); + if (nparam != nparam_actual) + { + std::ostringstream tmp_err; + tmp_err << "Expected number of rotation parameters (" << nparam << ") does not match number in file (" + << nparam_actual << ")"; + throw std::runtime_error(tmp_err.str()); + } + + std::vector params(nparam); + hin.read(params, rot_param_name); + for (int i = 0; i < nparam; i++) + this->myVars[i] = params[i]; + + hin.pop(); + + hin.pop(); +} + +template +void RotatedSPOsT::buildOptVariables(const size_t nel) +{ +#if !defined(QMC_COMPLEX) + /* Only rebuild optimized variables if more after-rotation orbitals are + * needed Consider ROHF, there is only one set of SPO for both spin up and + * down Nup > Ndown. nel_major_ will be set Nup. + * + * Use the size of myVars as a flag to avoid building the rotation + * parameters again when a clone is made (the DiracDeterminant constructor + * calls buildOptVariables) + */ + if (nel > nel_major_ && this->myVars.size() == 0) + { + nel_major_ = nel; + + const size_t nmo = Phi->getOrbitalSetSize(); + + // create active rotation parameter indices + RotationIndices created_m_act_rot_inds; + + RotationIndices created_full_rot_inds; + if (use_global_rot_) + createRotationIndicesFull(nel, nmo, created_full_rot_inds); + + createRotationIndices(nel, nmo, created_m_act_rot_inds); + + buildOptVariables(created_m_act_rot_inds, created_full_rot_inds); + } +#endif +} + +template +void RotatedSPOsT::buildOptVariables(const RotationIndices& rotations, const RotationIndices& full_rotations) +{ +#if !defined(QMC_COMPLEX) + const size_t nmo = Phi->getOrbitalSetSize(); + + // create active rotations + m_act_rot_inds = rotations; + + if (use_global_rot_) + m_full_rot_inds = full_rotations; + + if (use_global_rot_) + app_log() << "Orbital rotation using global rotation" << std::endl; + else + app_log() << "Orbital rotation using history" << std::endl; + + // This will add the orbital rotation parameters to myVars + // and will also read in initial parameter values supplied in input file + int p, q; + int nparams_active = m_act_rot_inds.size(); + + app_log() << "nparams_active: " << nparams_active << " params2.size(): " << params.size() << std::endl; + if (params_supplied) + if (nparams_active != params.size()) + throw std::runtime_error("The number of supplied orbital rotation parameters does not " + "match number prdouced by the slater " + "expansion. \n"); + + this->myVars.clear(); + for (int i = 0; i < nparams_active; i++) + { + p = m_act_rot_inds[i].first; + q = m_act_rot_inds[i].second; + std::stringstream sstr; + sstr << SPOSetT::getName() << "_orb_rot_" << (p < 10 ? "0" : "") << (p < 100 ? "0" : "") << (p < 1000 ? "0" : "") + << p << "_" << (q < 10 ? "0" : "") << (q < 100 ? "0" : "") << (q < 1000 ? "0" : "") << q; + + // If the user input parameters, use those. Otherwise, initialize the + // parameters to zero + if (params_supplied) + { + this->myVars.insert(sstr.str(), params[i]); + } + else + { + this->myVars.insert(sstr.str(), 0.0); + } + } + + if (use_global_rot_) + { + myVarsFull.clear(); + for (int i = 0; i < m_full_rot_inds.size(); i++) + { + p = m_full_rot_inds[i].first; + q = m_full_rot_inds[i].second; + std::stringstream sstr; + sstr << SPOSetT::getName() << "_orb_rot_" << (p < 10 ? "0" : "") << (p < 100 ? "0" : "") + << (p < 1000 ? "0" : "") << p << "_" << (q < 10 ? "0" : "") << (q < 100 ? "0" : "") << (q < 1000 ? "0" : "") + << q; + + if (params_supplied && i < m_act_rot_inds.size()) + myVarsFull.insert(sstr.str(), params[i]); + else + myVarsFull.insert(sstr.str(), 0.0); + } + } + + // Printing the parameters + if (true) + { + app_log() << std::string(16, ' ') << "Parameter name" << std::string(15, ' ') << "Value\n"; + this->myVars.print(app_log()); + } + + if (params_supplied) + { + std::vector param(m_act_rot_inds.size()); + for (int i = 0; i < m_act_rot_inds.size(); i++) + param[i] = this->myVars[i]; + apply_rotation(param, false); + } +#endif +} + +template +void RotatedSPOsT::apply_rotation(const std::vector& param, bool use_stored_copy) +{ + assert(param.size() == m_act_rot_inds.size()); + + const size_t nmo = Phi->getOrbitalSetSize(); + ValueMatrix rot_mat(nmo, nmo); + + constructAntiSymmetricMatrix(m_act_rot_inds, param, rot_mat); + + /* + rot_mat is now an anti-hermitian matrix. Now we convert + it into a unitary matrix via rot_mat = exp(-rot_mat). + Finally, apply unitary matrix to orbs. + */ + exponentiate_antisym_matrix(rot_mat); + Phi->applyRotation(rot_mat, use_stored_copy); +} + +template +void RotatedSPOsT::applyDeltaRotation(const std::vector& delta_param, + const std::vector& old_param, + std::vector& new_param) +{ + const size_t nmo = Phi->getOrbitalSetSize(); + ValueMatrix new_rot_mat(nmo, nmo); + constructDeltaRotation(delta_param, old_param, m_act_rot_inds, m_full_rot_inds, new_param, new_rot_mat); + + Phi->applyRotation(new_rot_mat, true); +} + +template +void RotatedSPOsT::constructDeltaRotation(const std::vector& delta_param, + const std::vector& old_param, + const RotationIndices& act_rot_inds, + const RotationIndices& full_rot_inds, + std::vector& new_param, + ValueMatrix& new_rot_mat) +{ + assert(delta_param.size() == act_rot_inds.size()); + assert(old_param.size() == full_rot_inds.size()); + assert(new_param.size() == full_rot_inds.size()); + + const size_t nmo = new_rot_mat.rows(); + assert(new_rot_mat.rows() == new_rot_mat.cols()); + + ValueMatrix old_rot_mat(nmo, nmo); + + constructAntiSymmetricMatrix(full_rot_inds, old_param, old_rot_mat); + exponentiate_antisym_matrix(old_rot_mat); + + ValueMatrix delta_rot_mat(nmo, nmo); + + constructAntiSymmetricMatrix(act_rot_inds, delta_param, delta_rot_mat); + exponentiate_antisym_matrix(delta_rot_mat); + + // Apply delta rotation to old rotation. + BLAS::gemm('N', 'N', nmo, nmo, nmo, 1.0, delta_rot_mat.data(), nmo, old_rot_mat.data(), nmo, 0.0, new_rot_mat.data(), + nmo); + + ValueMatrix log_rot_mat(nmo, nmo); + log_antisym_matrix(new_rot_mat, log_rot_mat); + extractParamsFromAntiSymmetricMatrix(full_rot_inds, log_rot_mat, new_param); +} + +template +void RotatedSPOsT::applyFullRotation(const std::vector& full_param, bool use_stored_copy) +{ + assert(full_param.size() == m_full_rot_inds.size()); + + const size_t nmo = Phi->getOrbitalSetSize(); + ValueMatrix rot_mat(nmo, nmo); + rot_mat = T(0); + + constructAntiSymmetricMatrix(m_full_rot_inds, full_param, rot_mat); + + /* + rot_mat is now an anti-hermitian matrix. Now we convert + it into a unitary matrix via rot_mat = exp(-rot_mat). + Finally, apply unitary matrix to orbs. + */ + exponentiate_antisym_matrix(rot_mat); + Phi->applyRotation(rot_mat, use_stored_copy); +} + +template +void RotatedSPOsT::applyRotationHistory() +{ + for (auto delta_param : history_params_) + { + apply_rotation(delta_param, false); + } +} + +// compute exponential of a real, antisymmetric matrix by diagonalizing and +// exponentiating eigenvalues +template +void RotatedSPOsT::exponentiate_antisym_matrix(ValueMatrix& mat) +{ + const int n = mat.rows(); + std::vector> mat_h(n * n, 0); + std::vector eval(n, 0); + std::vector> work(2 * n, 0); + std::vector rwork(3 * n, 0); + std::vector> mat_d(n * n, 0); + std::vector> mat_t(n * n, 0); + // exponentiating e^X = e^iY (Y hermitian) + // i(-iX) = X, so -iX is hermitian + // diagonalize -iX = UDU^T, exponentiate e^iD, and return U e^iD U^T + // construct hermitian analogue of mat by multiplying by -i + for (int i = 0; i < n; ++i) + { + for (int j = i; j < n; ++j) + { + mat_h[i + n * j] = std::complex(0, -1.0 * mat[j][i]); + mat_h[j + n * i] = std::complex(0, 1.0 * mat[j][i]); + } + } + // diagonalize the matrix + char JOBZ('V'); + char UPLO('U'); + int N(n); + int LDA(n); + int LWORK(2 * n); + int info = 0; + LAPACK::heev(JOBZ, UPLO, N, &mat_h.at(0), LDA, &eval.at(0), &work.at(0), LWORK, &rwork.at(0), info); + if (info != 0) + { + std::ostringstream msg; + msg << "heev failed with info = " << info << " in RotatedSPOsT::exponentiate_antisym_matrix"; + throw std::runtime_error(msg.str()); + } + // iterate through diagonal matrix, exponentiate terms + for (int i = 0; i < n; ++i) + { + for (int j = 0; j < n; ++j) + { + mat_d[i + j * n] = (i == j) ? std::exp(std::complex(0.0, eval[i])) : std::complex(0.0, 0.0); + } + } + // perform matrix multiplication + // assume row major + BLAS::gemm('N', 'C', n, n, n, std::complex(1.0, 0), &mat_d.at(0), n, &mat_h.at(0), n, + std::complex(0.0, 0.0), &mat_t.at(0), n); + BLAS::gemm('N', 'N', n, n, n, std::complex(1.0, 0), &mat_h.at(0), n, &mat_t.at(0), n, + std::complex(0.0, 0.0), &mat_d.at(0), n); + for (int i = 0; i < n; ++i) + for (int j = 0; j < n; ++j) + { + if (mat_d[i + n * j].imag() > 1e-12) + { + app_log() << "warning: large imaginary value in orbital " + "rotation matrix: (i,j) = (" + << i << "," << j << "), im = " << mat_d[i + n * j].imag() << std::endl; + } + mat[j][i] = mat_d[i + n * j].real(); + } +} + +template +void RotatedSPOsT::log_antisym_matrix(const ValueMatrix& mat, ValueMatrix& output) +{ + const int n = mat.rows(); + std::vector mat_h(n * n, 0); + std::vector eval_r(n, 0); + std::vector eval_i(n, 0); + std::vector mat_l(n * n, 0); + std::vector work(4 * n, 0); + + std::vector> mat_cd(n * n, 0); + std::vector> mat_cl(n * n, 0); + std::vector> mat_ch(n * n, 0); + + for (int i = 0; i < n; ++i) + for (int j = 0; j < n; ++j) + mat_h[i + n * j] = mat[i][j]; + + // diagonalize the matrix + char JOBL('V'); + char JOBR('N'); + int N(n); + int LDA(n); + int LWORK(4 * n); + int info = 0; + LAPACK::geev(&JOBL, &JOBR, &N, &mat_h.at(0), &LDA, &eval_r.at(0), &eval_i.at(0), &mat_l.at(0), &LDA, nullptr, &LDA, + &work.at(0), &LWORK, &info); + if (info != 0) + { + std::ostringstream msg; + msg << "heev failed with info = " << info << " in RotatedSPOsT::log_antisym_matrix"; + throw std::runtime_error(msg.str()); + } + + // iterate through diagonal matrix, take log + for (int i = 0; i < n; ++i) + { + for (int j = 0; j < n; ++j) + { + auto tmp = (i == j) ? std::log(std::complex(eval_r[i], eval_i[i])) : std::complex(0.0, 0.0); + mat_cd[i + j * n] = tmp; + + if (eval_i[j] > 0.0) + { + mat_cl[i + j * n] = std::complex(mat_l[i + j * n], mat_l[i + (j + 1) * n]); + mat_cl[i + (j + 1) * n] = std::complex(mat_l[i + j * n], -mat_l[i + (j + 1) * n]); + } + else if (!(eval_i[j] < 0.0)) + { + mat_cl[i + j * n] = std::complex(mat_l[i + j * n], 0.0); + } + } + } + + RealType one(1.0); + RealType zero(0.0); + BLAS::gemm('N', 'N', n, n, n, one, &mat_cl.at(0), n, &mat_cd.at(0), n, zero, &mat_ch.at(0), n); + BLAS::gemm('N', 'C', n, n, n, one, &mat_ch.at(0), n, &mat_cl.at(0), n, zero, &mat_cd.at(0), n); + + for (int i = 0; i < n; ++i) + for (int j = 0; j < n; ++j) + { + if (mat_cd[i + n * j].imag() > 1e-12) + { + app_log() << "warning: large imaginary value in antisymmetric " + "matrix: (i,j) = (" + << i << "," << j << "), im = " << mat_cd[i + n * j].imag() << std::endl; + } + output[i][j] = mat_cd[i + n * j].real(); + } +} + +template +void RotatedSPOsT::evaluateDerivRatios(const VirtualParticleSet& VP, + const opt_variables_type& optvars, + ValueVector& psi, + const ValueVector& psiinv, + std::vector& ratios, + Matrix& dratios, + int FirstIndex, + int LastIndex) +{ + Phi->evaluateDetRatios(VP, psi, psiinv, ratios); + + const size_t nel = LastIndex - FirstIndex; + const size_t nmo = Phi->getOrbitalSetSize(); + + psiM_inv.resize(nel, nel); + psiM_all.resize(nel, nmo); + dpsiM_all.resize(nel, nmo); + d2psiM_all.resize(nel, nmo); + + psiM_inv = 0; + psiM_all = 0; + dpsiM_all = 0; + d2psiM_all = 0; + + const ParticleSet& P = VP.getRefPS(); + int iel = VP.refPtcl; + + Phi->evaluate_notranspose(P, FirstIndex, LastIndex, psiM_all, dpsiM_all, d2psiM_all); + + for (int i = 0; i < nel; i++) + for (int j = 0; j < nel; j++) + psiM_inv(i, j) = psiM_all(i, j); + + Invert(psiM_inv.data(), nel, nel); + + const T* const A(psiM_all.data()); + const T* const Ainv(psiM_inv.data()); + ValueMatrix T_orig; + T_orig.resize(nel, nmo); + + BLAS::gemm('N', 'N', nmo, nel, nel, T(1.0), A, nmo, Ainv, nel, T(0.0), T_orig.data(), nmo); + + ValueMatrix T_mat; + T_mat.resize(nel, nmo); + + ValueVector tmp_psi; + tmp_psi.resize(nmo); + + for (int iat = 0; iat < VP.getTotalNum(); iat++) + { + Phi->evaluateValue(VP, iat, tmp_psi); + + for (int j = 0; j < nmo; j++) + psiM_all(iel - FirstIndex, j) = tmp_psi[j]; + + for (int i = 0; i < nel; i++) + for (int j = 0; j < nel; j++) + psiM_inv(i, j) = psiM_all(i, j); + + Invert(psiM_inv.data(), nel, nel); + + const T* const A(psiM_all.data()); + const T* const Ainv(psiM_inv.data()); + + // The matrix A is rectangular. Ainv is the inverse of the square part + // of the matrix. The multiply of Ainv and the square part of A is just + // the identity. This multiply could be reduced to Ainv and the + // non-square part of A. + BLAS::gemm('N', 'N', nmo, nel, nel, T(1.0), A, nmo, Ainv, nel, T(0.0), T_mat.data(), nmo); + + for (int i = 0; i < m_act_rot_inds.size(); i++) + { + int kk = this->myVars.where(i); + if (kk >= 0) + { + const int p = m_act_rot_inds.at(i).first; + const int q = m_act_rot_inds.at(i).second; + dratios(iat, kk) = T_mat(p, q) - T_orig(p, q); // dratio size is (nknot, num_vars) + } + } + } +} + +template +void RotatedSPOsT::evaluateDerivativesWF(ParticleSet& P, + const opt_variables_type& optvars, + Vector& dlogpsi, + int FirstIndex, + int LastIndex) +{ + const size_t nel = LastIndex - FirstIndex; + const size_t nmo = Phi->getOrbitalSetSize(); + + //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~PART1 + + psiM_inv.resize(nel, nel); + psiM_all.resize(nel, nmo); + dpsiM_all.resize(nel, nmo); + d2psiM_all.resize(nel, nmo); + + psiM_inv = 0; + psiM_all = 0; + dpsiM_all = 0; + d2psiM_all = 0; + + Phi->evaluate_notranspose(P, FirstIndex, LastIndex, psiM_all, dpsiM_all, d2psiM_all); + + for (int i = 0; i < nel; i++) + for (int j = 0; j < nel; j++) + psiM_inv(i, j) = psiM_all(i, j); + + Invert(psiM_inv.data(), nel, nel); + + //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~PART2 + const T* const A(psiM_all.data()); + const T* const Ainv(psiM_inv.data()); + ValueMatrix T_mat; + T_mat.resize(nel, nmo); + + BLAS::gemm('N', 'N', nmo, nel, nel, T(1.0), A, nmo, Ainv, nel, T(0.0), T_mat.data(), nmo); + + for (int i = 0; i < m_act_rot_inds.size(); i++) + { + int kk = this->myVars.where(i); + if (kk >= 0) + { + const int p = m_act_rot_inds.at(i).first; + const int q = m_act_rot_inds.at(i).second; + dlogpsi[kk] = T_mat(p, q); + } + } +} + +template +void RotatedSPOsT::evaluateDerivatives(ParticleSet& P, + const opt_variables_type& optvars, + Vector& dlogpsi, + Vector& dhpsioverpsi, + const int& FirstIndex, + const int& LastIndex) +{ + const size_t nel = LastIndex - FirstIndex; + const size_t nmo = Phi->getOrbitalSetSize(); + + //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~PART1 + myG_temp.resize(nel); + myG_J.resize(nel); + myL_temp.resize(nel); + myL_J.resize(nel); + + myG_temp = 0; + myG_J = 0; + myL_temp = 0; + myL_J = 0; + + Bbar.resize(nel, nmo); + psiM_inv.resize(nel, nel); + psiM_all.resize(nel, nmo); + dpsiM_all.resize(nel, nmo); + d2psiM_all.resize(nel, nmo); + + Bbar = 0; + psiM_inv = 0; + psiM_all = 0; + dpsiM_all = 0; + d2psiM_all = 0; + + Phi->evaluate_notranspose(P, FirstIndex, LastIndex, psiM_all, dpsiM_all, d2psiM_all); + + for (int i = 0; i < nel; i++) + for (int j = 0; j < nel; j++) + psiM_inv(i, j) = psiM_all(i, j); + + Invert(psiM_inv.data(), nel, nel); + + // current value of Gradient and Laplacian + // gradient components + for (int a = 0; a < nel; a++) + for (int i = 0; i < nel; i++) + for (int k = 0; k < 3; k++) + myG_temp[a][k] += psiM_inv(i, a) * dpsiM_all(a, i)[k]; + // laplacian components + for (int a = 0; a < nel; a++) + { + for (int i = 0; i < nel; i++) + myL_temp[a] += psiM_inv(i, a) * d2psiM_all(a, i); + } + + // calculation of myG_J which will be used to represent + // \frac{\nabla\psi_{J}}{\psi_{J}} calculation of myL_J will be used to + // represent \frac{\nabla^2\psi_{J}}{\psi_{J}} IMPORTANT NOTE: The value of + // P.L holds \nabla^2 ln[\psi] but we need \frac{\nabla^2 \psi}{\psi} and + // this is what myL_J will hold + for (int a = 0, iat = FirstIndex; a < nel; a++, iat++) + { + myG_J[a] = (P.G[iat] - myG_temp[a]); + myL_J[a] = (P.L[iat] + dot(P.G[iat], P.G[iat]) - myL_temp[a]); + } + // possibly replace wit BLAS calls + for (int i = 0; i < nel; i++) + for (int j = 0; j < nmo; j++) + Bbar(i, j) = d2psiM_all(i, j) + 2 * dot(myG_J[i], dpsiM_all(i, j)) + myL_J[i] * psiM_all(i, j); + + //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~PART2 + const T* const A(psiM_all.data()); + const T* const Ainv(psiM_inv.data()); + const T* const B(Bbar.data()); + ValueMatrix T_mat; + ValueMatrix Y1; + ValueMatrix Y2; + ValueMatrix Y3; + ValueMatrix Y4; + T_mat.resize(nel, nmo); + Y1.resize(nel, nel); + Y2.resize(nel, nmo); + Y3.resize(nel, nmo); + Y4.resize(nel, nmo); + + BLAS::gemm('N', 'N', nmo, nel, nel, T(1.0), A, nmo, Ainv, nel, T(0.0), T_mat.data(), nmo); + BLAS::gemm('N', 'N', nel, nel, nel, T(1.0), B, nmo, Ainv, nel, T(0.0), Y1.data(), nel); + BLAS::gemm('N', 'N', nmo, nel, nel, T(1.0), T_mat.data(), nmo, Y1.data(), nel, T(0.0), Y2.data(), nmo); + BLAS::gemm('N', 'N', nmo, nel, nel, T(1.0), B, nmo, Ainv, nel, T(0.0), Y3.data(), nmo); + + // possibly replace with BLAS call + Y4 = Y3 - Y2; + + for (int i = 0; i < m_act_rot_inds.size(); i++) + { + int kk = this->myVars.where(i); + if (kk >= 0) + { + const int p = m_act_rot_inds.at(i).first; + const int q = m_act_rot_inds.at(i).second; + dlogpsi[kk] += T_mat(p, q); + dhpsioverpsi[kk] += T(-0.5) * Y4(p, q); + } + } +} + +template +void RotatedSPOsT::evaluateDerivatives(ParticleSet& P, + const opt_variables_type& optvars, + Vector& dlogpsi, + Vector& dhpsioverpsi, + const T& psiCurrent, + const std::vector& Coeff, + const std::vector& C2node_up, + const std::vector& C2node_dn, + const ValueVector& detValues_up, + const ValueVector& detValues_dn, + const GradMatrix& grads_up, + const GradMatrix& grads_dn, + const ValueMatrix& lapls_up, + const ValueMatrix& lapls_dn, + const ValueMatrix& M_up, + const ValueMatrix& M_dn, + const ValueMatrix& Minv_up, + const ValueMatrix& Minv_dn, + const GradMatrix& B_grad, + const ValueMatrix& B_lapl, + const std::vector& detData_up, + const size_t N1, + const size_t N2, + const size_t NP1, + const size_t NP2, + const std::vector>& lookup_tbl) +{ + bool recalculate(false); + for (int k = 0; k < this->myVars.size(); ++k) + { + int kk = this->myVars.where(k); + if (kk < 0) + continue; + if (optvars.recompute(kk)) + recalculate = true; + } + if (recalculate) + { + ParticleSet::ParticleGradient myG_temp, myG_J; + ParticleSet::ParticleLaplacian myL_temp, myL_J; + const int NP = P.getTotalNum(); + myG_temp.resize(NP); + myG_temp = 0.0; + myL_temp.resize(NP); + myL_temp = 0.0; + myG_J.resize(NP); + myG_J = 0.0; + myL_J.resize(NP); + myL_J = 0.0; + const size_t nmo = Phi->getOrbitalSetSize(); + const size_t nel = P.last(0) - P.first(0); + + const T* restrict C_p = Coeff.data(); + for (int i = 0; i < Coeff.size(); i++) + { + const size_t upC = C2node_up[i]; + const size_t dnC = C2node_dn[i]; + const T tmp1 = C_p[i] * detValues_dn[dnC]; + const T tmp2 = C_p[i] * detValues_up[upC]; + for (size_t k = 0, j = N1; k < NP1; k++, j++) + { + myG_temp[j] += tmp1 * grads_up(upC, k); + myL_temp[j] += tmp1 * lapls_up(upC, k); + } + for (size_t k = 0, j = N2; k < NP2; k++, j++) + { + myG_temp[j] += tmp2 * grads_dn(dnC, k); + myL_temp[j] += tmp2 * lapls_dn(dnC, k); + } + } + + myG_temp *= (1 / psiCurrent); + myL_temp *= (1 / psiCurrent); + + // calculation of myG_J which will be used to represent + // \frac{\nabla\psi_{J}}{\psi_{J}} calculation of myL_J will be used to + // represent \frac{\nabla^2\psi_{J}}{\psi_{J}} IMPORTANT NOTE: The + // value of P.L holds \nabla^2 ln[\psi] but we need \frac{\nabla^2 + // \psi}{\psi} and this is what myL_J will hold + for (int iat = 0; iat < (myL_temp.size()); iat++) + { + myG_J[iat] = (P.G[iat] - myG_temp[iat]); + myL_J[iat] = (P.L[iat] + dot(P.G[iat], P.G[iat]) - myL_temp[iat]); + } + + table_method_eval(dlogpsi, dhpsioverpsi, myL_J, myG_J, nel, nmo, psiCurrent, Coeff, C2node_up, C2node_dn, + detValues_up, detValues_dn, grads_up, grads_dn, lapls_up, lapls_dn, M_up, M_dn, Minv_up, Minv_dn, + B_grad, B_lapl, detData_up, N1, N2, NP1, NP2, lookup_tbl); + } +} + +template +void RotatedSPOsT::evaluateDerivativesWF(ParticleSet& P, + const opt_variables_type& optvars, + Vector& dlogpsi, + const FullRealType& psiCurrent, + const std::vector& Coeff, + const std::vector& C2node_up, + const std::vector& C2node_dn, + const ValueVector& detValues_up, + const ValueVector& detValues_dn, + const ValueMatrix& M_up, + const ValueMatrix& M_dn, + const ValueMatrix& Minv_up, + const ValueMatrix& Minv_dn, + const std::vector& detData_up, + const std::vector>& lookup_tbl) +{ + bool recalculate(false); + for (int k = 0; k < this->myVars.size(); ++k) + { + int kk = this->myVars.where(k); + if (kk < 0) + continue; + if (optvars.recompute(kk)) + recalculate = true; + } + if (recalculate) + { + const size_t nmo = Phi->getOrbitalSetSize(); + const size_t nel = P.last(0) - P.first(0); + + table_method_evalWF(dlogpsi, nel, nmo, psiCurrent, Coeff, C2node_up, C2node_dn, detValues_up, detValues_dn, M_up, + M_dn, Minv_up, Minv_dn, detData_up, lookup_tbl); + } +} + +template +void RotatedSPOsT::table_method_eval(Vector& dlogpsi, + Vector& dhpsioverpsi, + const ParticleSet::ParticleLaplacian& myL_J, + const ParticleSet::ParticleGradient& myG_J, + const size_t nel, + const size_t nmo, + const T& psiCurrent, + const std::vector& Coeff, + const std::vector& C2node_up, + const std::vector& C2node_dn, + const ValueVector& detValues_up, + const ValueVector& detValues_dn, + const GradMatrix& grads_up, + const GradMatrix& grads_dn, + const ValueMatrix& lapls_up, + const ValueMatrix& lapls_dn, + const ValueMatrix& M_up, + const ValueMatrix& M_dn, + const ValueMatrix& Minv_up, + const ValueMatrix& Minv_dn, + const GradMatrix& B_grad, + const ValueMatrix& B_lapl, + const std::vector& detData_up, + const size_t N1, + const size_t N2, + const size_t NP1, + const size_t NP2, + const std::vector>& lookup_tbl) +/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +GUIDE TO THE MATICES BEING BUILT +---------------------------------------------- +The idea here is that there is a loop over all unique determinants. For each +determiant the table method is employed to calculate the contributions to the +parameter derivatives (dhpsioverpsi/dlogpsi) + + loop through unquie determinants + loop through parameters + evaluate contributaion to dlogpsi and dhpsioverpsi +\noindent + + BLAS GUIDE for matrix multiplication of [ alpha * A.B + beta * C = C ] + Matrix A is of dimensions a1,a2 and Matrix B is b1,b2 in which a2=b1 + The BLAS command is as follows... + + BLAS::gemm('N','N', b2, a1, a2 ,alpha, B, b2, A, a2, beta, C, b2); + +Below is a human readable format for the matrix multiplications performed +below... + +This notation is inspired by http://dx.doi.org/10.1063/1.4948778 +\newline +\hfill\break +$ + A_{i,j}=\phi_j(r_{i}) \\ + T = A^{-1} \widetilde{A} \\ + B_{i,j} =\nabla^2 \phi_{j}(r_i) + \frac{\nabla_{i}J}{J} \cdot \nabla +\phi_{j}(r_{i}) + \frac{\nabla^2_i J}{J} \phi_{j}(r_{i}) \\ + \hat{O_{I}} = \hat{O}D_{I} \\ + D_{I}=det(A_{I}) \newline + \psi_{MS} = \sum_{I=0} C_{I} D_{I\uparrow}D_{I\downarrow} \\ + \Psi_{total} = \psi_{J}\psi_{MS} \\ + \alpha_{I} = P^{T}_{I}TQ_{I} \\ + M_{I} = P^{T}_{I} \widetilde{M} Q_{I} = P^{T}_{I} (A^{-1}\widetilde{B} - +A^{-1} B A^{-1}\widetilde{A} )Q_{I} \\ +$ +\newline +There are three constants I use in the expressions for dhpsioverpsi and dlogpsi +\newline +\hfill\break +$ + const0 = C_{0}*det(A_{0\downarrow})+\sum_{I=1} C_{I}*det(A_{I\downarrow})* +det(\alpha_{I\uparrow}) \\ + const1 = C_{0}*\hat{O} det(A_{0\downarrow})+\sum_{I=1} +C_{I}*\hat{O}det(A_{I\downarrow})* det(\alpha_{I\uparrow}) \\ + const2 = \sum_{I=1} C_{I}*det(A_{I\downarrow})* +Tr[\alpha_{I}^{-1}M_{I}]*det(\alpha_{I}) \\ +$ +\newline +Below is a translation of the shorthand I use to represent matrices independent +of ``excitation matrix". \newline \hfill\break +$ + Y_{1} = A^{-1}B \\ + Y_{2} = A^{-1}BA^{-1}\widetilde{A} \\ + Y_{3} = A^{-1}\widetilde{B} \\ + Y_{4} = \widetilde{M} = (A^{-1}\widetilde{B} - A^{-1} B A^{-1}\widetilde{A} +)\\ +$ +\newline +Below is a translation of the shorthand I use to represent matrices dependent on +``excitation" with respect to the reference Matrix and sums of matrices. Above +this line I have represented these excitation matrices with a subscript ``I" but +from this point on The subscript will be omitted and it is clear that whenever a +matrix depends on $P^{T}_I$ and $Q_{I}$ that this is an excitation matrix. The +reference matrix is always $A_{0}$ and is always the Hartree Fock Matrix. +\newline +\hfill\break +$ + Y_{5} = TQ \\ + Y_{6} = (P^{T}TQ)^{-1} = \alpha_{I}^{-1}\\ + Y_{7} = \alpha_{I}^{-1} P^{T} \\ + Y_{11} = \widetilde{M}Q \\ + Y_{23} = P^{T}\widetilde{M}Q \\ + Y_{24} = \alpha_{I}^{-1}P^{T}\widetilde{M}Q \\ + Y_{25} = \alpha_{I}^{-1}P^{T}\widetilde{M}Q\alpha_{I}^{-1} \\ + Y_{26} = \alpha_{I}^{-1}P^{T}\widetilde{M}Q\alpha_{I}^{-1}P^{T}\\ +$ +\newline +So far you will notice that I have not included up or down arrows to specify +what spin the matrices are of. This is because we are calculating the derivative +of all up or all down spin orbital rotation parameters at a time. If we are +finding the up spin derivatives then any term that is down spin will be +constant. The following assumes that we are taking up-spin MO rotation parameter +derivatives. Of course the down spin expression can be retrieved by swapping the +up and down arrows. I have dubbed any expression with lowercase p prefix as a +"precursor" to an expression actually used... \newline \hfill\break +$ + \dot{C_{I}} = C_{I}*det(A_{I\downarrow})\\ + \ddot{C_{I}} = C_{I}*\hat{O}det(A_{I\downarrow}) \\ + pK1 = \sum_{I=1} \dot{C_{I}} det(\alpha_{I}) Tr[\alpha_{I}^{-1}M_{I}] +(Q\alpha_{I}^{-1}P^{T}) \\ + pK2 = \sum_{I=1} \dot{C_{I}} det(\alpha_{I}) (Q\alpha_{I}^{-1}P^{T}) \\ + pK3 = \sum_{I=1} \ddot{C_{I}} det(\alpha_{I}) (Q\alpha_{I}^{-1}P^{T}) \\ + pK4 = \sum_{I=1} \dot{C_{I}} det(A_{I}) (Q\alpha_{I}^{-1}P^{T}) \\ + pK5 = \sum_{I=1} \dot{C_{I}} det(\alpha_{I}) (Q\alpha_{I}^{-1} M_{I} +\alpha_{I}^{-1}P^{T}) \\ +$ +\newline +Now these p matrices will be used to make various expressions via BLAS commands. +\newline +\hfill\break +$ + K1T = const0^{-1}*pK1.T =const0^{-1} \sum_{I=1} \dot{C_{I}} det(\alpha_{I}) +Tr[\alpha_{I}^{-1}M_{I}] (Q\alpha_{I}^{-1}P^{T}T) \\ + TK1T = T.K1T = const0^{-1} \sum_{I=1} \dot{C_{I}} det(\alpha_{I}) +Tr[\alpha_{I}^{-1}M_{I}] (TQ\alpha_{I}^{-1}P^{T}T)\\ \\ + K2AiB = const0^{-1} \sum_{I=1} \dot{C_{I}} det(\alpha_{I}) +(Q\alpha_{I}^{-1}P^{T}A^{-1}\widetilde{B})\\ + TK2AiB = T.K2AiB = const0^{-1} \sum_{I=1} \dot{C_{I}} det(\alpha_{I}) +(TQ\alpha_{I}^{-1}P^{T}A^{-1}\widetilde{B})\\ + K2XA = const0^{-1} \sum_{I=1} \dot{C_{I}} det(\alpha_{I}) +(Q\alpha_{I}^{-1}P^{T}X\widetilde{A})\\ + TK2XA = T.K2XA = const0^{-1} \sum_{I=1} \dot{C_{I}} det(\alpha_{I}) +(TQ\alpha_{I}^{-1}P^{T}X\widetilde{A})\\ \\ + K2T = \frac{const1}{const0^{2}} \sum_{I=1} \dot{C_{I}} det(\alpha_{I}) +(Q\alpha_{I}^{-1}P^{T}T) \\ + TK2T = T.K2T =\frac{const1}{const0^{2}} \sum_{I=1} \dot{C_{I}} +det(\alpha_{I}) (TQ\alpha_{I}^{-1}P^{T}T) \\ + MK2T = \frac{const0}{const1} Y_{4}.K2T= const0^{-1} \sum_{I=1} \dot{C_{I}} +det(\alpha_{I}) (\widetilde{M}Q\alpha_{I}^{-1}P^{T}T)\\ \\ + K3T = const0^{-1} \sum_{I=1} \ddot{C_{I}} det(\alpha_{I}) +(Q\alpha_{I}^{-1}P^{T}T) \\ + TK3T = T.K3T = const0^{-1} \sum_{I=1} \ddot{C_{I}} det(\alpha_{I}) +(TQ\alpha_{I}^{-1}P^{T}T)\\ \\ + K4T = \sum_{I=1} \dot{C_{I}} det(A_{I}) (Q\alpha_{I}^{-1}P^{T}T) \\ + TK4T = T.K4T = \sum_{I=1} \dot{C_{I}} det(A_{I}) (TQ\alpha_{I}^{-1}P^{T}T) +\\ \\ + K5T = const0^{-1} \sum_{I=1} \dot{C_{I}} det(\alpha_{I}) (Q\alpha_{I}^{-1} +M_{I} \alpha_{I}^{-1}P^{T} T) \\ + TK5T = T.K5T = \sum_{I=1} \dot{C_{I}} det(\alpha_{I}) (T Q\alpha_{I}^{-1} +M_{I} \alpha_{I}^{-1}P^{T} T) \\ +$ +\newline +Now with all these matrices and constants the expressions of dhpsioverpsi and +dlogpsi can be created. + + + + +In addition I will be using a special generalization of the kinetic operator +which I will denote as O. Our Slater matrix with the special O operator applied +to each element will be called B_bar + +$ +``Bbar"_{i,j} =\nabla^2 \phi_{j}(r_i) + \frac{\nabla_{i}J}{J} \cdot \nabla +\phi_{j}(r_{i}) + \frac{\nabla^2_i J}{J} \phi_{j}(r_{i}) +$ +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ +{ + ValueMatrix Table; + ValueMatrix Bbar; + ValueMatrix Y1, Y2, Y3, Y4, Y5, Y6, Y7, Y11, Y23, Y24, Y25, Y26; + ValueMatrix pK1, K1T, TK1T, pK2, K2AiB, TK2AiB, K2XA, TK2XA, K2T, TK2T, MK2T, pK3, K3T, TK3T, pK5, K5T, TK5T; + + Table.resize(nel, nmo); + + Bbar.resize(nel, nmo); + + Y1.resize(nel, nel); + Y2.resize(nel, nmo); + Y3.resize(nel, nmo); + Y4.resize(nel, nmo); + + pK1.resize(nmo, nel); + K1T.resize(nmo, nmo); + TK1T.resize(nel, nmo); + + pK2.resize(nmo, nel); + K2AiB.resize(nmo, nmo); + TK2AiB.resize(nel, nmo); + K2XA.resize(nmo, nmo); + TK2XA.resize(nel, nmo); + K2T.resize(nmo, nmo); + TK2T.resize(nel, nmo); + MK2T.resize(nel, nmo); + + pK3.resize(nmo, nel); + K3T.resize(nmo, nmo); + TK3T.resize(nel, nmo); + + pK5.resize(nmo, nel); + K5T.resize(nmo, nmo); + TK5T.resize(nel, nmo); + + const int parameters_size(m_act_rot_inds.size()); + const int parameter_start_index(0); + + const size_t num_unique_up_dets(detValues_up.size()); + const size_t num_unique_dn_dets(detValues_dn.size()); + + const T* restrict cptr = Coeff.data(); + const size_t nc = Coeff.size(); + const size_t* restrict upC(C2node_up.data()); + const size_t* restrict dnC(C2node_dn.data()); + // B_grad holds the gradient operator + // B_lapl holds the laplacian operator + // B_bar will hold our special O operator + + const int offset1(N1); + const int offset2(N2); + const int NPother(NP2); + + T* T_(Table.data()); + + // possibly replace wit BLAS calls + for (int i = 0; i < nel; i++) + for (int j = 0; j < nmo; j++) + Bbar(i, j) = B_lapl(i, j) + 2 * dot(myG_J[i + offset1], B_grad(i, j)) + myL_J[i + offset1] * M_up(i, j); + + const T* restrict B(Bbar.data()); + const T* restrict A(M_up.data()); + const T* restrict Ainv(Minv_up.data()); + // IMPORTANT NOTE: THE Dets[0]->psiMinv OBJECT DOES NOT HOLD THE INVERSE IF + // THE MULTIDIRACDETERMINANTBASE ONLY CONTAINS ONE ELECTRON. NEED A FIX FOR + // THIS CASE + // The T matrix should be calculated and stored for use + // T = A^{-1} \widetilde A + // REMINDER: that the ValueMatrix "matrix" stores data in a row major order + // and that BLAS commands assume column major + BLAS::gemm('N', 'N', nmo, nel, nel, T(1.0), A, nmo, Ainv, nel, RealType(0.0), T_, nmo); + + BLAS::gemm('N', 'N', nel, nel, nel, T(1.0), B, nmo, Ainv, nel, RealType(0.0), Y1.data(), nel); + BLAS::gemm('N', 'N', nmo, nel, nel, T(1.0), T_, nmo, Y1.data(), nel, RealType(0.0), Y2.data(), nmo); + BLAS::gemm('N', 'N', nmo, nel, nel, T(1.0), B, nmo, Ainv, nel, RealType(0.0), Y3.data(), nmo); + + // possibly replace with BLAS call + Y4 = Y3 - Y2; + + // Need to create the constants: (Oi, const0, const1, const2)to take + // advantage of minimal BLAS commands; Oi is the special operator applied to + // the slater matrix "A subscript i" from the total CI expansion \hat{O_{i}} + //= \hat{O}D_{i} with D_{i}=det(A_{i}) and Multi-Slater component defined as + //\sum_{i=0} C_{i} D_{i\uparrow}D_{i\downarrow} + std::vector Oi(num_unique_dn_dets); + + for (int index = 0; index < num_unique_dn_dets; index++) + for (int iat = 0; iat < NPother; iat++) + Oi[index] += lapls_dn(index, iat) + 2 * dot(grads_dn(index, iat), myG_J[offset2 + iat]) + + myL_J[offset2 + iat] * detValues_dn[index]; + + // const0 = C_{0}*det(A_{0\downarrow})+\sum_{i=1} + // C_{i}*det(A_{i\downarrow})* det(\alpha_{i\uparrow}) const1 = + // C_{0}*\hat{O} det(A_{0\downarrow})+\sum_{i=1} + // C_{i}*\hat{O}det(A_{i\downarrow})* det(\alpha_{i\uparrow}) const2 = + // \sum_{i=1} C_{i}*det(A_{i\downarrow})* + // Tr[\alpha_{i}^{-1}M_{i}]*det(\alpha_{i}) + RealType const0(0.0), const1(0.0), const2(0.0); + for (size_t i = 0; i < nc; ++i) + { + const RealType c = cptr[i]; + const size_t up = upC[i]; + const size_t down = dnC[i]; + + const0 += c * detValues_dn[down] * (detValues_up[up] / detValues_up[0]); + const1 += c * Oi[down] * (detValues_up[up] / detValues_up[0]); + } + + std::fill(pK1.begin(), pK1.end(), 0.0); + std::fill(pK2.begin(), pK2.end(), 0.0); + std::fill(pK3.begin(), pK3.end(), 0.0); + std::fill(pK5.begin(), pK5.end(), 0.0); + + // Now we are going to loop through all unique determinants. + // The few lines above are for the reference matrix contribution. + // Although I start the loop below from index 0, the loop only performs + // actions when the index is >= 1 the detData object contains all the + // information about the P^T and Q matrices (projection matrices) needed in + // the table method + const int* restrict data_it = detData_up.data(); + for (int index = 0, datum = 0; index < num_unique_up_dets; index++) + { + const int k = data_it[datum]; + + if (k == 0) + { + datum += 3 * k + 1; + } + + else + { + // Number of rows and cols of P^T + const int prows = k; + const int pcols = nel; + // Number of rows and cols of Q + const int qrows = nmo; + const int qcols = k; + + Y5.resize(nel, k); + Y6.resize(k, k); + + // Any matrix multiplication of P^T or Q is simply a projection + // Explicit matrix multiplication can be avoided; instead column or + // row copying can be done BlAS::copy(size of col/row being copied, + // Matrix pointer + place to begin copying, + // storage spacing (number of elements btw next row/col + // element), Pointer to resultant matrix + place to begin + // pasting, storage spacing of resultant matrix) + // For example the next 4 lines is the matrix multiplication of T*Q + // = Y5 + std::fill(Y5.begin(), Y5.end(), 0.0); + for (int i = 0; i < k; i++) + { + BLAS::copy(nel, T_ + data_it[datum + 1 + k + i], nmo, Y5.data() + i, k); + } + + std::fill(Y6.begin(), Y6.end(), 0.0); + for (int i = 0; i < k; i++) + { + BLAS::copy(k, Y5.data() + (data_it[datum + 1 + i]) * k, 1, (Y6.data() + i * k), 1); + } + + Vector WS; + Vector Piv; + WS.resize(k); + Piv.resize(k); + std::complex logdet = 0.0; + InvertWithLog(Y6.data(), k, k, WS.data(), Piv.data(), logdet); + + Y11.resize(nel, k); + Y23.resize(k, k); + Y24.resize(k, k); + Y25.resize(k, k); + Y26.resize(k, nel); + + std::fill(Y11.begin(), Y11.end(), 0.0); + for (int i = 0; i < k; i++) + { + BLAS::copy(nel, Y4.data() + (data_it[datum + 1 + k + i]), nmo, Y11.data() + i, k); + } + + std::fill(Y23.begin(), Y23.end(), 0.0); + for (int i = 0; i < k; i++) + { + BLAS::copy(k, Y11.data() + (data_it[datum + 1 + i]) * k, 1, (Y23.data() + i * k), 1); + } + + BLAS::gemm('N', 'N', k, k, k, RealType(1.0), Y23.data(), k, Y6.data(), k, RealType(0.0), Y24.data(), k); + BLAS::gemm('N', 'N', k, k, k, RealType(1.0), Y6.data(), k, Y24.data(), k, RealType(0.0), Y25.data(), k); + + Y26.resize(k, nel); + + std::fill(Y26.begin(), Y26.end(), 0.0); + for (int i = 0; i < k; i++) + { + BLAS::copy(k, Y25.data() + i, k, Y26.data() + (data_it[datum + 1 + i]), nel); + } + + Y7.resize(k, nel); + + std::fill(Y7.begin(), Y7.end(), 0.0); + for (int i = 0; i < k; i++) + { + BLAS::copy(k, Y6.data() + i, k, Y7.data() + (data_it[datum + 1 + i]), nel); + } + + // c_Tr_AlphaI_MI is a constant contributing to constant const2 + // c_Tr_AlphaI_MI = Tr[\alpha_{I}^{-1}(P^{T}\widetilde{M} Q)] + RealType c_Tr_AlphaI_MI = 0.0; + for (int i = 0; i < k; i++) + { + c_Tr_AlphaI_MI += Y24(i, i); + } + + for (int p = 0; p < lookup_tbl[index].size(); p++) + { + // el_p is the element position that contains information about + // the CI coefficient, and det up/dn values associated with the + // current unique determinant + const int el_p(lookup_tbl[index][p]); + const RealType c = cptr[el_p]; + const size_t up = upC[el_p]; + const size_t down = dnC[el_p]; + + const RealType alpha_1(c * detValues_dn[down] * detValues_up[up] / detValues_up[0] * c_Tr_AlphaI_MI); + const RealType alpha_2(c * detValues_dn[down] * detValues_up[up] / detValues_up[0]); + const RealType alpha_3(c * Oi[down] * detValues_up[up] / detValues_up[0]); + + const2 += alpha_1; + + for (int i = 0; i < k; i++) + { + BLAS::axpy(nel, alpha_1, Y7.data() + i * nel, 1, pK1.data() + (data_it[datum + 1 + k + i]) * nel, 1); + BLAS::axpy(nel, alpha_2, Y7.data() + i * nel, 1, pK2.data() + (data_it[datum + 1 + k + i]) * nel, 1); + BLAS::axpy(nel, alpha_3, Y7.data() + i * nel, 1, pK3.data() + (data_it[datum + 1 + k + i]) * nel, 1); + BLAS::axpy(nel, alpha_2, Y26.data() + i * nel, 1, pK5.data() + (data_it[datum + 1 + k + i]) * nel, 1); + } + } + datum += 3 * k + 1; + } + } + + BLAS::gemm('N', 'N', nmo, nmo, nel, 1.0 / const0, T_, nmo, pK1.data(), nel, RealType(0.0), K1T.data(), nmo); + BLAS::gemm('N', 'N', nmo, nel, nmo, RealType(1.0), K1T.data(), nmo, T_, nmo, RealType(0.0), TK1T.data(), nmo); + + BLAS::gemm('N', 'N', nmo, nmo, nel, 1.0 / const0, Y3.data(), nmo, pK2.data(), nel, RealType(0.0), K2AiB.data(), nmo); + BLAS::gemm('N', 'N', nmo, nel, nmo, RealType(1.0), K2AiB.data(), nmo, T_, nmo, RealType(0.0), TK2AiB.data(), nmo); + BLAS::gemm('N', 'N', nmo, nmo, nel, 1.0 / const0, Y2.data(), nmo, pK2.data(), nel, RealType(0.0), K2XA.data(), nmo); + BLAS::gemm('N', 'N', nmo, nel, nmo, RealType(1.0), K2XA.data(), nmo, T_, nmo, RealType(0.0), TK2XA.data(), nmo); + + BLAS::gemm('N', 'N', nmo, nmo, nel, const1 / (const0 * const0), T_, nmo, pK2.data(), nel, RealType(0.0), K2T.data(), + nmo); + BLAS::gemm('N', 'N', nmo, nel, nmo, RealType(1.0), K2T.data(), nmo, T_, nmo, RealType(0.0), TK2T.data(), nmo); + BLAS::gemm('N', 'N', nmo, nel, nmo, const0 / const1, K2T.data(), nmo, Y4.data(), nmo, RealType(0.0), MK2T.data(), + nmo); + + BLAS::gemm('N', 'N', nmo, nmo, nel, 1.0 / const0, T_, nmo, pK3.data(), nel, RealType(0.0), K3T.data(), nmo); + BLAS::gemm('N', 'N', nmo, nel, nmo, RealType(1.0), K3T.data(), nmo, T_, nmo, RealType(0.0), TK3T.data(), nmo); + + BLAS::gemm('N', 'N', nmo, nmo, nel, 1.0 / const0, T_, nmo, pK5.data(), nel, RealType(0.0), K5T.data(), nmo); + BLAS::gemm('N', 'N', nmo, nel, nmo, RealType(1.0), K5T.data(), nmo, T_, nmo, RealType(0.0), TK5T.data(), nmo); + + for (int mu = 0, k = parameter_start_index; k < (parameter_start_index + parameters_size); k++, mu++) + { + int kk = this->myVars.where(k); + if (kk >= 0) + { + const int i(m_act_rot_inds[mu].first), j(m_act_rot_inds[mu].second); + if (i <= nel - 1 && j > nel - 1) + { + dhpsioverpsi[kk] += + T(-0.5 * Y4(i, j) - + 0.5 * + (-K5T(i, j) + K5T(j, i) + TK5T(i, j) + K2AiB(i, j) - K2AiB(j, i) - TK2AiB(i, j) - K2XA(i, j) + + K2XA(j, i) + TK2XA(i, j) - MK2T(i, j) + K1T(i, j) - K1T(j, i) - TK1T(i, j) - + const2 / const1 * K2T(i, j) + const2 / const1 * K2T(j, i) + const2 / const1 * TK2T(i, j) + + K3T(i, j) - K3T(j, i) - TK3T(i, j) - K2T(i, j) + K2T(j, i) + TK2T(i, j))); + } + else if (i <= nel - 1 && j <= nel - 1) + { + dhpsioverpsi[kk] += + T(-0.5 * (Y4(i, j) - Y4(j, i)) - + 0.5 * + (-K5T(i, j) + K5T(j, i) + TK5T(i, j) - TK5T(j, i) + K2AiB(i, j) - K2AiB(j, i) - TK2AiB(i, j) + + TK2AiB(j, i) - K2XA(i, j) + K2XA(j, i) + TK2XA(i, j) - TK2XA(j, i) - MK2T(i, j) + MK2T(j, i) + + K1T(i, j) - K1T(j, i) - TK1T(i, j) + TK1T(j, i) - const2 / const1 * K2T(i, j) + + const2 / const1 * K2T(j, i) + const2 / const1 * TK2T(i, j) - const2 / const1 * TK2T(j, i) + + K3T(i, j) - K3T(j, i) - TK3T(i, j) + TK3T(j, i) - K2T(i, j) + K2T(j, i) + TK2T(i, j) - TK2T(j, i))); + } + else + { + dhpsioverpsi[kk] += T(-0.5 * + (-K5T(i, j) + K5T(j, i) + K2AiB(i, j) - K2AiB(j, i) - K2XA(i, j) + K2XA(j, i) + + + K1T(i, j) - K1T(j, i) - const2 / const1 * K2T(i, j) + const2 / const1 * K2T(j, i) + + K3T(i, j) - K3T(j, i) - K2T(i, j) + K2T(j, i))); + } + } + } +} + +template +void RotatedSPOsT::table_method_evalWF(Vector& dlogpsi, + const size_t nel, + const size_t nmo, + const T& psiCurrent, + const std::vector& Coeff, + const std::vector& C2node_up, + const std::vector& C2node_dn, + const ValueVector& detValues_up, + const ValueVector& detValues_dn, + const ValueMatrix& M_up, + const ValueMatrix& M_dn, + const ValueMatrix& Minv_up, + const ValueMatrix& Minv_dn, + const std::vector& detData_up, + const std::vector>& lookup_tbl) +{ + ValueMatrix Table; + ValueMatrix Y5, Y6, Y7; + ValueMatrix pK4, K4T, TK4T; + + Table.resize(nel, nmo); + + Bbar.resize(nel, nmo); + + pK4.resize(nmo, nel); + K4T.resize(nmo, nmo); + TK4T.resize(nel, nmo); + + const int parameters_size(m_act_rot_inds.size()); + const int parameter_start_index(0); + + const size_t num_unique_up_dets(detValues_up.size()); + const size_t num_unique_dn_dets(detValues_dn.size()); + + const T* restrict cptr = Coeff.data(); + const size_t nc = Coeff.size(); + const size_t* restrict upC(C2node_up.data()); + const size_t* restrict dnC(C2node_dn.data()); + + T* T_(Table.data()); + + const T* restrict A(M_up.data()); + const T* restrict Ainv(Minv_up.data()); + // IMPORTANT NOTE: THE Dets[0]->psiMinv OBJECT DOES NOT HOLD THE INVERSE IF + // THE MULTIDIRACDETERMINANTBASE ONLY CONTAINS ONE ELECTRON. NEED A FIX FOR + // THIS CASE + // The T matrix should be calculated and stored for use + // T = A^{-1} \widetilde A + // REMINDER: that the ValueMatrix "matrix" stores data in a row major order + // and that BLAS commands assume column major + BLAS::gemm('N', 'N', nmo, nel, nel, RealType(1.0), A, nmo, Ainv, nel, RealType(0.0), T_, nmo); + + // const0 = C_{0}*det(A_{0\downarrow})+\sum_{i=1} + // C_{i}*det(A_{i\downarrow})* det(\alpha_{i\uparrow}) + RealType const0(0.0), const1(0.0), const2(0.0); + for (size_t i = 0; i < nc; ++i) + { + const RealType c = cptr[i]; + const size_t up = upC[i]; + const size_t down = dnC[i]; + + const0 += c * detValues_dn[down] * (detValues_up[up] / detValues_up[0]); + } + + std::fill(pK4.begin(), pK4.end(), 0.0); + + // Now we are going to loop through all unique determinants. + // The few lines above are for the reference matrix contribution. + // Although I start the loop below from index 0, the loop only performs + // actions when the index is >= 1 the detData object contains all the + // information about the P^T and Q matrices (projection matrices) needed in + // the table method + const int* restrict data_it = detData_up.data(); + for (int index = 0, datum = 0; index < num_unique_up_dets; index++) + { + const int k = data_it[datum]; + + if (k == 0) + { + datum += 3 * k + 1; + } + + else + { + // Number of rows and cols of P^T + const int prows = k; + const int pcols = nel; + // Number of rows and cols of Q + const int qrows = nmo; + const int qcols = k; + + Y5.resize(nel, k); + Y6.resize(k, k); + + // Any matrix multiplication of P^T or Q is simply a projection + // Explicit matrix multiplication can be avoided; instead column or + // row copying can be done BlAS::copy(size of col/row being copied, + // Matrix pointer + place to begin copying, + // storage spacing (number of elements btw next row/col + // element), Pointer to resultant matrix + place to begin + // pasting, storage spacing of resultant matrix) + // For example the next 4 lines is the matrix multiplication of T*Q + // = Y5 + std::fill(Y5.begin(), Y5.end(), 0.0); + for (int i = 0; i < k; i++) + { + BLAS::copy(nel, T_ + data_it[datum + 1 + k + i], nmo, Y5.data() + i, k); + } + + std::fill(Y6.begin(), Y6.end(), 0.0); + for (int i = 0; i < k; i++) + { + BLAS::copy(k, Y5.data() + (data_it[datum + 1 + i]) * k, 1, (Y6.data() + i * k), 1); + } + + Vector WS; + Vector Piv; + WS.resize(k); + Piv.resize(k); + std::complex logdet = 0.0; + InvertWithLog(Y6.data(), k, k, WS.data(), Piv.data(), logdet); + + Y7.resize(k, nel); + + std::fill(Y7.begin(), Y7.end(), 0.0); + for (int i = 0; i < k; i++) + { + BLAS::copy(k, Y6.data() + i, k, Y7.data() + (data_it[datum + 1 + i]), nel); + } + + for (int p = 0; p < lookup_tbl[index].size(); p++) + { + // el_p is the element position that contains information about + // the CI coefficient, and det up/dn values associated with the + // current unique determinant + const int el_p(lookup_tbl[index][p]); + const RealType c = cptr[el_p]; + const size_t up = upC[el_p]; + const size_t down = dnC[el_p]; + + const RealType alpha_4(c * detValues_dn[down] * detValues_up[up] * (1 / psiCurrent)); + + for (int i = 0; i < k; i++) + { + BLAS::axpy(nel, alpha_4, Y7.data() + i * nel, 1, pK4.data() + (data_it[datum + 1 + k + i]) * nel, 1); + } + } + datum += 3 * k + 1; + } + } + + BLAS::gemm('N', 'N', nmo, nmo, nel, RealType(1.0), T_, nmo, pK4.data(), nel, RealType(0.0), K4T.data(), nmo); + BLAS::gemm('N', 'N', nmo, nel, nmo, RealType(1.0), K4T.data(), nmo, T_, nmo, RealType(0.0), TK4T.data(), nmo); + + for (int mu = 0, k = parameter_start_index; k < (parameter_start_index + parameters_size); k++, mu++) + { + int kk = this->myVars.where(k); + if (kk >= 0) + { + const int i(m_act_rot_inds[mu].first), j(m_act_rot_inds[mu].second); + if (i <= nel - 1 && j > nel - 1) + { + dlogpsi[kk] += + T(detValues_up[0] * (Table(i, j)) * const0 * (1 / psiCurrent) + (K4T(i, j) - K4T(j, i) - TK4T(i, j))); + } + else if (i <= nel - 1 && j <= nel - 1) + { + dlogpsi[kk] += T(detValues_up[0] * (Table(i, j) - Table(j, i)) * const0 * (1 / psiCurrent) + + (K4T(i, j) - TK4T(i, j) - K4T(j, i) + TK4T(j, i))); + } + else + { + dlogpsi[kk] += T((K4T(i, j) - K4T(j, i))); + } + } + } +} + +template +std::unique_ptr> RotatedSPOsT::makeClone() const +{ + auto myclone = std::make_unique(SPOSetT::getName(), std::unique_ptr>(Phi->makeClone())); + + myclone->params = this->params; + myclone->params_supplied = this->params_supplied; + myclone->m_act_rot_inds = this->m_act_rot_inds; + myclone->m_full_rot_inds = this->m_full_rot_inds; + myclone->myVars = this->myVars; + myclone->myVarsFull = this->myVarsFull; + myclone->history_params_ = this->history_params_; + myclone->use_global_rot_ = this->use_global_rot_; + return myclone; +} + +// Class concrete types from ValueType +template class RotatedSPOsT; +template class RotatedSPOsT; + +} // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/RotatedSPOsT.h b/src/QMCWaveFunctions/RotatedSPOsT.h new file mode 100644 index 0000000000..3273681455 --- /dev/null +++ b/src/QMCWaveFunctions/RotatedSPOsT.h @@ -0,0 +1,420 @@ +////////////////////////////////////////////////////////////////////////////////////// +//// This file is distributed under the University of Illinois/NCSA Open Source +/// License. / See LICENSE file in top directory for details. +//// +//// Copyright (c) QMCPACK developers. +//// +//// File developed by: Sergio D. Pineda Flores, +/// sergio_pinedaflores@berkeley.edu, University of California, Berkeley / Eric +/// Neuscamman, eneuscamman@berkeley.edu, University of California, Berkeley / +/// Ye Luo, yeluo@anl.gov, Argonne National Laboratory +//// +//// File created by: Sergio D. Pineda Flores, sergio_pinedaflores@berkeley.edu, +/// University of California, Berkeley +//////////////////////////////////////////////////////////////////////////////////////// +#ifndef QMCPLUSPLUS_ROTATEDSPOST_H +#define QMCPLUSPLUS_ROTATEDSPOST_H + +#include "QMCWaveFunctions/SPOSetT.h" + +namespace qmcplusplus +{ +template +class RotatedSPOsT; +namespace testing +{ +opt_variables_type& getMyVarsFull(RotatedSPOsT& rot); +opt_variables_type& getMyVarsFull(RotatedSPOsT& rot); +std::vector>& getHistoryParams(RotatedSPOsT& rot); +std::vector>& getHistoryParams(RotatedSPOsT& rot); +} // namespace testing + +template +class RotatedSPOsT : public SPOSetT, public OptimizableObject +{ +public: + using IndexType = typename SPOSetT::IndexType; + using RealType = typename SPOSetT::RealType; + using FullRealType = typename SPOSetT::FullRealType; + using ValueVector = typename SPOSetT::ValueVector; + using ValueMatrix = typename SPOSetT::ValueMatrix; + using GradVector = typename SPOSetT::GradVector; + using GradMatrix = typename SPOSetT::GradMatrix; + using HessVector = typename SPOSetT::HessVector; + using HessMatrix = typename SPOSetT::HessMatrix; + using GGGVector = typename SPOSetT::GGGVector; + using GGGMatrix = typename SPOSetT::GGGMatrix; + + // constructor + RotatedSPOsT(const std::string& my_name, std::unique_ptr>&& spos); + // destructor + ~RotatedSPOsT() override; + + std::string getClassName() const override { return "RotatedSPOsT"; } + bool isOptimizable() const override { return true; } + bool isOMPoffload() const override { return Phi->isOMPoffload(); } + bool hasIonDerivs() const override { return Phi->hasIonDerivs(); } + + // Vector of rotation matrix indices + using RotationIndices = std::vector>; + + // Active orbital rotation parameter indices + RotationIndices m_act_rot_inds; + + // Full set of rotation values for global rotation + RotationIndices m_full_rot_inds; + + // Construct a list of the matrix indices for non-zero rotation parameters. + // (The structure for a sparse representation of the matrix) + // Only core->active rotations are created. + static void createRotationIndices(int nel, int nmo, RotationIndices& rot_indices); + + // Construct a list for all the matrix indices, including core->active, + // core->core and active->active + static void createRotationIndicesFull(int nel, int nmo, RotationIndices& rot_indices); + + // Fill in antisymmetric matrix from the list of rotation parameter indices + // and a list of parameter values. + // This function assumes rot_mat is properly sized upon input and is set to + // zero. + static void constructAntiSymmetricMatrix(const RotationIndices& rot_indices, + const std::vector& param, + ValueMatrix& rot_mat); + + // Extract the list of rotation parameters from the entries in an + // antisymmetric matrix This function expects rot_indices and param are the + // same length. + static void extractParamsFromAntiSymmetricMatrix(const RotationIndices& rot_indices, + const ValueMatrix& rot_mat, + std::vector& param); + + // function to perform orbital rotations + void apply_rotation(const std::vector& param, bool use_stored_copy); + + // For global rotation, inputs are the old parameters and the delta + // parameters. The corresponding rotation matrices are constructed, + // multiplied together, and the new parameters extracted. The new rotation + // is applied to the underlying SPO coefficients + void applyDeltaRotation(const std::vector& delta_param, + const std::vector& old_param, + std::vector& new_param); + + // Perform the construction of matrices and extraction of parameters for a + // delta rotation. Split out and made static for testing. + static void constructDeltaRotation(const std::vector& delta_param, + const std::vector& old_param, + const RotationIndices& act_rot_inds, + const RotationIndices& full_rot_inds, + std::vector& new_param, + ValueMatrix& new_rot_mat); + + // When initializing the rotation from VP files + // This function applies the rotation history + void applyRotationHistory(); + + // This function applies the global rotation (similar to apply_rotation, but + // for the full set of rotation parameters) + void applyFullRotation(const std::vector& full_param, bool use_stored_copy); + + // Compute matrix exponential of an antisymmetric matrix (result is rotation + // matrix) + static void exponentiate_antisym_matrix(ValueMatrix& mat); + + // Compute matrix log of rotation matrix to produce antisymmetric matrix + static void log_antisym_matrix(const ValueMatrix& mat, ValueMatrix& output); + + // A particular SPOSet used for Orbitals + std::unique_ptr> Phi; + + /// Set the rotation parameters (usually from input file) + void setRotationParameters(const std::vector& param_list); + + /// the number of electrons of the majority spin + size_t nel_major_; + + std::unique_ptr> makeClone() const override; + + // myG_temp (myL_temp) is the Gradient (Laplacian) value of of the + // Determinant part of the wfn myG_J is the Gradient of the all other parts + // of the wavefunction (typically just the Jastrow). + // It represents \frac{\nabla\psi_{J}}{\psi_{J}} + // myL_J will be used to represent \frac{\nabla^2\psi_{J}}{\psi_{J}} . The + // Laplacian portion IMPORTANT NOTE: The value of P.L holds \nabla^2 + // ln[\psi] but we need \frac{\nabla^2 \psi}{\psi} and this is what myL_J + // will hold + ParticleSet::ParticleGradient myG_temp, myG_J; + ParticleSet::ParticleLaplacian myL_temp, myL_J; + + ValueMatrix Bbar; + ValueMatrix psiM_inv; + ValueMatrix psiM_all; + GradMatrix dpsiM_all; + ValueMatrix d2psiM_all; + + // Single Slater creation + void buildOptVariables(size_t nel); + + // For the MSD case rotations must be created in MultiSlaterDetTableMethod + // class + void buildOptVariables(const RotationIndices& rotations, const RotationIndices& full_rotations); + + void evaluateDerivatives(ParticleSet& P, + const opt_variables_type& optvars, + Vector& dlogpsi, + Vector& dhpsioverpsi, + const int& FirstIndex, + const int& LastIndex) override; + + void evaluateDerivativesWF(ParticleSet& P, + const opt_variables_type& optvars, + Vector& dlogpsi, + int FirstIndex, + int LastIndex) override; + + void evaluateDerivatives(ParticleSet& P, + const opt_variables_type& optvars, + Vector& dlogpsi, + Vector& dhpsioverpsi, + const T& psiCurrent, + const std::vector& Coeff, + const std::vector& C2node_up, + const std::vector& C2node_dn, + const ValueVector& detValues_up, + const ValueVector& detValues_dn, + const GradMatrix& grads_up, + const GradMatrix& grads_dn, + const ValueMatrix& lapls_up, + const ValueMatrix& lapls_dn, + const ValueMatrix& M_up, + const ValueMatrix& M_dn, + const ValueMatrix& Minv_up, + const ValueMatrix& Minv_dn, + const GradMatrix& B_grad, + const ValueMatrix& B_lapl, + const std::vector& detData_up, + const size_t N1, + const size_t N2, + const size_t NP1, + const size_t NP2, + const std::vector>& lookup_tbl) override; + + void evaluateDerivativesWF(ParticleSet& P, + const opt_variables_type& optvars, + Vector& dlogpsi, + const FullRealType& psiCurrent, + const std::vector& Coeff, + const std::vector& C2node_up, + const std::vector& C2node_dn, + const ValueVector& detValues_up, + const ValueVector& detValues_dn, + const ValueMatrix& M_up, + const ValueMatrix& M_dn, + const ValueMatrix& Minv_up, + const ValueMatrix& Minv_dn, + const std::vector& detData_up, + const std::vector>& lookup_tbl) override; + + // helper function to evaluatederivative; evaluate orbital rotation + // parameter derivative using table method + void table_method_eval(Vector& dlogpsi, + Vector& dhpsioverpsi, + const ParticleSet::ParticleLaplacian& myL_J, + const ParticleSet::ParticleGradient& myG_J, + const size_t nel, + const size_t nmo, + const T& psiCurrent, + const std::vector& Coeff, + const std::vector& C2node_up, + const std::vector& C2node_dn, + const ValueVector& detValues_up, + const ValueVector& detValues_dn, + const GradMatrix& grads_up, + const GradMatrix& grads_dn, + const ValueMatrix& lapls_up, + const ValueMatrix& lapls_dn, + const ValueMatrix& M_up, + const ValueMatrix& M_dn, + const ValueMatrix& Minv_up, + const ValueMatrix& Minv_dn, + const GradMatrix& B_grad, + const ValueMatrix& B_lapl, + const std::vector& detData_up, + const size_t N1, + const size_t N2, + const size_t NP1, + const size_t NP2, + const std::vector>& lookup_tbl); + + void table_method_evalWF(Vector& dlogpsi, + const size_t nel, + const size_t nmo, + const T& psiCurrent, + const std::vector& Coeff, + const std::vector& C2node_up, + const std::vector& C2node_dn, + const ValueVector& detValues_up, + const ValueVector& detValues_dn, + const ValueMatrix& M_up, + const ValueMatrix& M_dn, + const ValueMatrix& Minv_up, + const ValueMatrix& Minv_dn, + const std::vector& detData_up, + const std::vector>& lookup_tbl); + + void extractOptimizableObjectRefs(UniqueOptObjRefs& opt_obj_refs) override { opt_obj_refs.push_back(*this); } + + void checkInVariablesExclusive(opt_variables_type& active) override + { + if (this->myVars.size()) + active.insertFrom(this->myVars); + } + + void checkOutVariables(const opt_variables_type& active) override { this->myVars.getIndex(active); } + + /// reset + void resetParametersExclusive(const opt_variables_type& active) override; + + void writeVariationalParameters(hdf_archive& hout) override; + + void readVariationalParameters(hdf_archive& hin) override; + + //********************************************************************************* + // the following functions simply call Phi's corresponding functions + void setOrbitalSetSize(int norbs) override { Phi->setOrbitalSetSize(norbs); } + + void checkObject() const override { Phi->checkObject(); } + + void evaluateValue(const ParticleSet& P, int iat, ValueVector& psi) override + { + assert(psi.size() <= this->OrbitalSetSize); + Phi->evaluateValue(P, iat, psi); + } + + void evaluateVGL(const ParticleSet& P, int iat, ValueVector& psi, GradVector& dpsi, ValueVector& d2psi) override + { + assert(psi.size() <= this->OrbitalSetSize); + Phi->evaluateVGL(P, iat, psi, dpsi, d2psi); + } + + void evaluateDetRatios(const VirtualParticleSet& VP, + ValueVector& psi, + const ValueVector& psiinv, + std::vector& ratios) override + { + Phi->evaluateDetRatios(VP, psi, psiinv, ratios); + } + + void evaluateDerivRatios(const VirtualParticleSet& VP, + const opt_variables_type& optvars, + ValueVector& psi, + const ValueVector& psiinv, + std::vector& ratios, + Matrix& dratios, + int FirstIndex, + int LastIndex) override; + + void evaluateVGH(const ParticleSet& P, + int iat, + ValueVector& psi, + GradVector& dpsi, + HessVector& grad_grad_psi) override + { + assert(psi.size() <= this->OrbitalSetSize); + Phi->evaluateVGH(P, iat, psi, dpsi, grad_grad_psi); + } + + void evaluateVGHGH(const ParticleSet& P, + int iat, + ValueVector& psi, + GradVector& dpsi, + HessVector& grad_grad_psi, + GGGVector& grad_grad_grad_psi) override + { + Phi->evaluateVGHGH(P, iat, psi, dpsi, grad_grad_psi, grad_grad_grad_psi); + } + + void evaluate_notranspose(const ParticleSet& P, + int first, + int last, + ValueMatrix& logdet, + GradMatrix& dlogdet, + ValueMatrix& d2logdet) override + { + Phi->evaluate_notranspose(P, first, last, logdet, dlogdet, d2logdet); + } + + void evaluate_notranspose(const ParticleSet& P, + int first, + int last, + ValueMatrix& logdet, + GradMatrix& dlogdet, + HessMatrix& grad_grad_logdet) override + { + Phi->evaluate_notranspose(P, first, last, logdet, dlogdet, grad_grad_logdet); + } + + void evaluate_notranspose(const ParticleSet& P, + int first, + int last, + ValueMatrix& logdet, + GradMatrix& dlogdet, + HessMatrix& grad_grad_logdet, + GGGMatrix& grad_grad_grad_logdet) override + { + Phi->evaluate_notranspose(P, first, last, logdet, dlogdet, grad_grad_logdet, grad_grad_grad_logdet); + } + + void evaluateGradSource(const ParticleSet& P, + int first, + int last, + const ParticleSet& source, + int iat_src, + GradMatrix& grad_phi) override + { + Phi->evaluateGradSource(P, first, last, source, iat_src, grad_phi); + } + + void evaluateGradSource(const ParticleSet& P, + int first, + int last, + const ParticleSet& source, + int iat_src, + GradMatrix& grad_phi, + HessMatrix& grad_grad_phi, + GradMatrix& grad_lapl_phi) override + { + Phi->evaluateGradSource(P, first, last, source, iat_src, grad_phi, grad_grad_phi, grad_lapl_phi); + } + + // void evaluateThirdDeriv(const ParticleSet& P, int first, int last, + // GGGMatrix& grad_grad_grad_logdet) {Phi->evaluateThridDeriv(P, first, + // last, grad_grad_grad_logdet); } + + /// Use history list (false) or global rotation (true) + void set_use_global_rotation(bool use_global_rotation) { use_global_rot_ = use_global_rotation; } + +private: + /// true if SPO parameters (orbital rotation parameters) have been supplied + /// by input + bool params_supplied; + /// list of supplied orbital rotation parameters + std::vector params; + + /// Full set of rotation matrix parameters for use in global rotation method + opt_variables_type myVarsFull; + + /// List of previously applied parameters + std::vector> history_params_; + + /// Use global rotation or history list + bool use_global_rot_ = true; + + friend opt_variables_type& testing::getMyVarsFull(RotatedSPOsT& rot); + friend opt_variables_type& testing::getMyVarsFull(RotatedSPOsT& rot); + friend std::vector>& testing::getHistoryParams(RotatedSPOsT& rot); + friend std::vector>& testing::getHistoryParams(RotatedSPOsT& rot); +}; + +} // namespace qmcplusplus + +#endif diff --git a/src/QMCWaveFunctions/SPOSetT.cpp b/src/QMCWaveFunctions/SPOSetT.cpp index e61f4cace1..34c76bad82 100644 --- a/src/QMCWaveFunctions/SPOSetT.cpp +++ b/src/QMCWaveFunctions/SPOSetT.cpp @@ -134,7 +134,7 @@ void SPOSetT::mw_evaluateVGLandDetRatioGrads(const RefVectorWithLeader& grads) const { assert(this == &spo_list.getLeader()); - assert(phi_vgl_v.size(0) == DIM_VGL); + assert(phi_vgl_v.size(0) == QMCTraits::DIM_VGL); assert(phi_vgl_v.size(1) == spo_list.size()); const size_t nw = spo_list.size(); const size_t norb_requested = phi_vgl_v.size(2); @@ -360,7 +360,7 @@ template void SPOSetT::evaluateDerivativesWF(ParticleSet& P, const opt_variables_type& optvars, Vector& dlogpsi, - const T& psiCurrent, + const typename QTFull::ValueType& psiCurrent, const std::vector& Coeff, const std::vector& C2node_up, const std::vector& C2node_dn, @@ -435,4 +435,4 @@ template class SPOSetT; template class SPOSetT>; template class SPOSetT>; -} // namespace qmcplusplus \ No newline at end of file +} // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/SPOSetT.h b/src/QMCWaveFunctions/SPOSetT.h index 95643fd3c5..ddc14c6593 100644 --- a/src/QMCWaveFunctions/SPOSetT.h +++ b/src/QMCWaveFunctions/SPOSetT.h @@ -65,6 +65,10 @@ class SPOSetT : public QMCTraits using SPOMap = std::map>>; using OffloadMWVGLArray = Array>; // [VGL, walker, Orbs] using OffloadMWVArray = Array>; // [walker, Orbs] + using PosType = typename OrbitalSetTraits::PosType; + using RealType = typename OrbitalSetTraits::RealType; + using ValueType = typename OrbitalSetTraits::ValueType; + using FullRealType = typename OrbitalSetTraits::RealType; template using OffloadMatrix = Matrix>; @@ -176,7 +180,7 @@ class SPOSetT : public QMCTraits virtual void evaluateDerivativesWF(ParticleSet& P, const opt_variables_type& optvars, Vector& dlogpsi, - const T& psiCurrent, + const typename QTFull::ValueType& psiCurrent, const std::vector& Coeff, const std::vector& C2node_up, const std::vector& C2node_dn, diff --git a/src/QMCWaveFunctions/SpinorSetT.cpp b/src/QMCWaveFunctions/SpinorSetT.cpp new file mode 100644 index 0000000000..64d7d3d6b1 --- /dev/null +++ b/src/QMCWaveFunctions/SpinorSetT.cpp @@ -0,0 +1,586 @@ +////////////////////////////////////////////////////////////////////////////////////// +// This file is distributed under the University of Illinois/NCSA Open Source License. +// See LICENSE file in top directory for details. +// +// Copyright (c) 2022 QMCPACK developers +// +// File developed by: Raymond Clay III, rclay@sandia.gov, Sandia National Laboratories +// Cody A. Melton, cmelton@sandia.gov, Sandia National Laboratories +// +// File created by: Raymond Clay III, rclay@sandia.gov, Sandia National Laboratories +////////////////////////////////////////////////////////////////////////////////////// + +#include "SpinorSetT.h" +#include "Utilities/ResourceCollection.h" +#include "Platforms/OMPTarget/OMPTargetMath.hpp" + +namespace qmcplusplus +{ +template +struct SpinorSetT::SpinorSetMultiWalkerResource : public Resource +{ + SpinorSetMultiWalkerResource() : Resource("SpinorSet") {} + SpinorSetMultiWalkerResource(const SpinorSetMultiWalkerResource&) : SpinorSetMultiWalkerResource() {} + std::unique_ptr makeClone() const override { return std::make_unique(*this); } + OffloadMWVGLArray up_phi_vgl_v, dn_phi_vgl_v; + std::vector up_ratios, dn_ratios; + std::vector up_grads, dn_grads; + std::vector spins; +}; + +template +SpinorSetT::SpinorSetT(const std::string& my_name) : SPOSetT(my_name), spo_up(nullptr), spo_dn(nullptr) +{} + +template +SpinorSetT::~SpinorSetT() = default; + +template +void SpinorSetT::set_spos(std::unique_ptr>&& up, std::unique_ptr>&& dn) +{ + //Sanity check for input SPO's. They need to be the same size or + IndexType spo_size_up = up->getOrbitalSetSize(); + IndexType spo_size_down = dn->getOrbitalSetSize(); + + if (spo_size_up != spo_size_down) + throw std::runtime_error("SpinorSet::set_spos(...): up and down SPO components have different sizes."); + + setOrbitalSetSize(spo_size_up); + + spo_up = std::move(up); + spo_dn = std::move(dn); + + psi_work_up.resize(this->OrbitalSetSize); + psi_work_down.resize(this->OrbitalSetSize); + + dpsi_work_up.resize(this->OrbitalSetSize); + dpsi_work_down.resize(this->OrbitalSetSize); + + d2psi_work_up.resize(this->OrbitalSetSize); + d2psi_work_down.resize(this->OrbitalSetSize); +} + +template +void SpinorSetT::setOrbitalSetSize(int norbs) +{ + this->OrbitalSetSize = norbs; +}; + +template +void SpinorSetT::evaluateValue(const ParticleSet& P, int iat, ValueVector& psi) +{ + psi_work_up = 0.0; + psi_work_down = 0.0; + + spo_up->evaluateValue(P, iat, psi_work_up); + spo_dn->evaluateValue(P, iat, psi_work_down); + + ParticleSet::Scalar_t s = P.activeSpin(iat); + + RealType coss(0.0), sins(0.0); + + coss = std::cos(s); + sins = std::sin(s); + + //This is only supported in the complex build, so T is some complex number depending on the precision. + T eis(coss, sins); + T emis(coss, -sins); + + psi = eis * psi_work_up + emis * psi_work_down; +} + +template +void SpinorSetT::evaluateVGL(const ParticleSet& P, int iat, ValueVector& psi, GradVector& dpsi, ValueVector& d2psi) +{ + psi_work_up = 0.0; + psi_work_down = 0.0; + dpsi_work_up = 0.0; + dpsi_work_down = 0.0; + d2psi_work_up = 0.0; + d2psi_work_down = 0.0; + + spo_up->evaluateVGL(P, iat, psi_work_up, dpsi_work_up, d2psi_work_up); + spo_dn->evaluateVGL(P, iat, psi_work_down, dpsi_work_down, d2psi_work_down); + + ParticleSet::Scalar_t s = P.activeSpin(iat); + + RealType coss(0.0), sins(0.0); + + coss = std::cos(s); + sins = std::sin(s); + + T eis(coss, sins); + T emis(coss, -sins); + + psi = eis * psi_work_up + emis * psi_work_down; + dpsi = eis * dpsi_work_up + emis * dpsi_work_down; + d2psi = eis * d2psi_work_up + emis * d2psi_work_down; +} + +template +void SpinorSetT::evaluateVGL_spin(const ParticleSet& P, + int iat, + ValueVector& psi, + GradVector& dpsi, + ValueVector& d2psi, + ValueVector& dspin) +{ + psi_work_up = 0.0; + psi_work_down = 0.0; + dpsi_work_up = 0.0; + dpsi_work_down = 0.0; + d2psi_work_up = 0.0; + d2psi_work_down = 0.0; + + spo_up->evaluateVGL(P, iat, psi_work_up, dpsi_work_up, d2psi_work_up); + spo_dn->evaluateVGL(P, iat, psi_work_down, dpsi_work_down, d2psi_work_down); + + ParticleSet::Scalar_t s = P.activeSpin(iat); + + RealType coss(0.0), sins(0.0); + + coss = std::cos(s); + sins = std::sin(s); + + T eis(coss, sins); + T emis(coss, -sins); + T eye(0, 1.0); + + psi = eis * psi_work_up + emis * psi_work_down; + dpsi = eis * dpsi_work_up + emis * dpsi_work_down; + d2psi = eis * d2psi_work_up + emis * d2psi_work_down; + dspin = eye * (eis * psi_work_up - emis * psi_work_down); +} + +template +void SpinorSetT::mw_evaluateVGLWithSpin(const RefVectorWithLeader>& spo_list, + const RefVectorWithLeader& P_list, + int iat, + const RefVector& psi_v_list, + const RefVector& dpsi_v_list, + const RefVector& d2psi_v_list, + OffloadMatrix& mw_dspin) const +{ + auto& spo_leader = spo_list.template getCastedLeader>(); + auto& P_leader = P_list.getLeader(); + assert(this == &spo_leader); + + IndexType nw = spo_list.size(); + auto [up_spo_list, dn_spo_list] = extractSpinComponentRefList(spo_list); + auto& up_spo_leader = up_spo_list.getLeader(); + auto& dn_spo_leader = dn_spo_list.getLeader(); + + RefVector up_psi_v_list, dn_psi_v_list; + RefVector up_dpsi_v_list, dn_dpsi_v_list; + RefVector up_d2psi_v_list, dn_d2psi_v_list; + for (int iw = 0; iw < nw; iw++) + { + auto& spo = spo_list.template getCastedElement>(iw); + up_psi_v_list.push_back(spo.psi_work_up); + dn_psi_v_list.push_back(spo.psi_work_down); + up_dpsi_v_list.push_back(spo.dpsi_work_up); + dn_dpsi_v_list.push_back(spo.dpsi_work_down); + up_d2psi_v_list.push_back(spo.d2psi_work_up); + dn_d2psi_v_list.push_back(spo.d2psi_work_down); + } + + up_spo_leader.mw_evaluateVGL(up_spo_list, P_list, iat, up_psi_v_list, up_dpsi_v_list, up_d2psi_v_list); + dn_spo_leader.mw_evaluateVGL(dn_spo_list, P_list, iat, dn_psi_v_list, dn_dpsi_v_list, dn_d2psi_v_list); + + for (int iw = 0; iw < nw; iw++) + { + ParticleSet::Scalar_t s = P_list[iw].activeSpin(iat); + RealType coss = std::cos(s); + RealType sins = std::sin(s); + + T eis(coss, sins); + T emis(coss, -sins); + T eye(0, 1.0); + + psi_v_list[iw].get() = eis * up_psi_v_list[iw].get() + emis * dn_psi_v_list[iw].get(); + dpsi_v_list[iw].get() = eis * up_dpsi_v_list[iw].get() + emis * dn_dpsi_v_list[iw].get(); + d2psi_v_list[iw].get() = eis * up_d2psi_v_list[iw].get() + emis * dn_d2psi_v_list[iw].get(); + for (int iorb = 0; iorb < this->OrbitalSetSize; iorb++) + mw_dspin(iw, iorb) = eye * (eis * (up_psi_v_list[iw].get())[iorb] - emis * (dn_psi_v_list[iw].get())[iorb]); + } + //Data above is all on host, but since mw_dspin is DualMatrix we need to sync the host and device + mw_dspin.updateTo(); +} + +template +void SpinorSetT::mw_evaluateVGLandDetRatioGradsWithSpin(const RefVectorWithLeader>& spo_list, + const RefVectorWithLeader& P_list, + int iat, + const std::vector& invRow_ptr_list, + OffloadMWVGLArray& phi_vgl_v, + std::vector& ratios, + std::vector& grads, + std::vector& spingrads) const +{ + auto& spo_leader = spo_list.template getCastedLeader>(); + auto& P_leader = P_list.getLeader(); + assert(this == &spo_leader); + assert(phi_vgl_v.size(0) == DIM_VGL); + assert(phi_vgl_v.size(1) == spo_list.size()); + const size_t nw = spo_list.size(); + const size_t norb_requested = phi_vgl_v.size(2); + + auto& mw_res = spo_leader.mw_res_handle_.getResource(); + auto& up_phi_vgl_v = mw_res.up_phi_vgl_v; + auto& dn_phi_vgl_v = mw_res.dn_phi_vgl_v; + auto& up_ratios = mw_res.up_ratios; + auto& dn_ratios = mw_res.dn_ratios; + auto& up_grads = mw_res.up_grads; + auto& dn_grads = mw_res.dn_grads; + auto& spins = mw_res.spins; + + up_phi_vgl_v.resize(QMCTraits::DIM_VGL, nw, norb_requested); + dn_phi_vgl_v.resize(QMCTraits::DIM_VGL, nw, norb_requested); + up_ratios.resize(nw); + dn_ratios.resize(nw); + up_grads.resize(nw); + dn_grads.resize(nw); + spins.resize(nw); + + auto [up_spo_list, dn_spo_list] = extractSpinComponentRefList(spo_list); + auto& up_spo_leader = up_spo_list.getLeader(); + auto& dn_spo_leader = dn_spo_list.getLeader(); + + up_spo_leader.mw_evaluateVGLandDetRatioGrads(up_spo_list, P_list, iat, invRow_ptr_list, up_phi_vgl_v, up_ratios, + up_grads); + dn_spo_leader.mw_evaluateVGLandDetRatioGrads(dn_spo_list, P_list, iat, invRow_ptr_list, dn_phi_vgl_v, dn_ratios, + dn_grads); + for (int iw = 0; iw < nw; iw++) + { + ParticleSet::Scalar_t s = P_list[iw].activeSpin(iat); + spins[iw] = s; + RealType coss = std::cos(s); + RealType sins = std::sin(s); + + T eis(coss, sins); + T emis(coss, -sins); + T eye(0, 1.0); + + ratios[iw] = eis * up_ratios[iw] + emis * dn_ratios[iw]; + grads[iw] = (eis * up_grads[iw] * up_ratios[iw] + emis * dn_grads[iw] * dn_ratios[iw]) / ratios[iw]; + spingrads[iw] = eye * (eis * up_ratios[iw] - emis * dn_ratios[iw]) / ratios[iw]; + } + + auto* spins_ptr = spins.data(); + //This data lives on the device + auto* phi_vgl_ptr = phi_vgl_v.data(); + auto* up_phi_vgl_ptr = up_phi_vgl_v.data(); + auto* dn_phi_vgl_ptr = dn_phi_vgl_v.data(); + PRAGMA_OFFLOAD("omp target teams distribute map(to:spins_ptr[0:nw])") + for (int iw = 0; iw < nw; iw++) + { + RealType c, s; + omptarget::sincos(spins_ptr[iw], &s, &c); + T eis(c, s), emis(c, -s); + PRAGMA_OFFLOAD("omp parallel for collapse(2)") + for (int idim = 0; idim < QMCTraits::DIM_VGL; idim++) + for (int iorb = 0; iorb < norb_requested; iorb++) + { + auto offset = idim * nw * norb_requested + iw * norb_requested + iorb; + phi_vgl_ptr[offset] = eis * up_phi_vgl_ptr[offset] + emis * dn_phi_vgl_ptr[offset]; + } + } +} + +template +void SpinorSetT::evaluate_notranspose(const ParticleSet& P, + int first, + int last, + ValueMatrix& logdet, + GradMatrix& dlogdet, + ValueMatrix& d2logdet) +{ + IndexType nelec = P.getTotalNum(); + + logpsi_work_up.resize(nelec, this->OrbitalSetSize); + logpsi_work_down.resize(nelec, this->OrbitalSetSize); + + dlogpsi_work_up.resize(nelec, this->OrbitalSetSize); + dlogpsi_work_down.resize(nelec, this->OrbitalSetSize); + + d2logpsi_work_up.resize(nelec, this->OrbitalSetSize); + d2logpsi_work_down.resize(nelec, this->OrbitalSetSize); + + spo_up->evaluate_notranspose(P, first, last, logpsi_work_up, dlogpsi_work_up, d2logpsi_work_up); + spo_dn->evaluate_notranspose(P, first, last, logpsi_work_down, dlogpsi_work_down, d2logpsi_work_down); + + + for (int iat = 0; iat < nelec; iat++) + { + ParticleSet::Scalar_t s = P.activeSpin(iat); + + RealType coss(0.0), sins(0.0); + + coss = std::cos(s); + sins = std::sin(s); + + T eis(coss, sins); + T emis(coss, -sins); + + for (int no = 0; no < this->OrbitalSetSize; no++) + { + logdet(iat, no) = eis * logpsi_work_up(iat, no) + emis * logpsi_work_down(iat, no); + dlogdet(iat, no) = eis * dlogpsi_work_up(iat, no) + emis * dlogpsi_work_down(iat, no); + d2logdet(iat, no) = eis * d2logpsi_work_up(iat, no) + emis * d2logpsi_work_down(iat, no); + } + } +} + +template +void SpinorSetT::mw_evaluate_notranspose(const RefVectorWithLeader>& spo_list, + const RefVectorWithLeader& P_list, + int first, + int last, + const RefVector& logdet_list, + const RefVector& dlogdet_list, + const RefVector& d2logdet_list) const +{ + auto& spo_leader = spo_list.template getCastedLeader>(); + auto& P_leader = P_list.getLeader(); + assert(this == &spo_leader); + + IndexType nw = spo_list.size(); + IndexType nelec = P_leader.getTotalNum(); + + auto [up_spo_list, dn_spo_list] = extractSpinComponentRefList(spo_list); + auto& up_spo_leader = up_spo_list.getLeader(); + auto& dn_spo_leader = dn_spo_list.getLeader(); + + std::vector mw_up_logdet, mw_dn_logdet; + std::vector mw_up_dlogdet, mw_dn_dlogdet; + std::vector mw_up_d2logdet, mw_dn_d2logdet; + mw_up_logdet.reserve(nw); + mw_dn_logdet.reserve(nw); + mw_up_dlogdet.reserve(nw); + mw_dn_dlogdet.reserve(nw); + mw_up_d2logdet.reserve(nw); + mw_dn_d2logdet.reserve(nw); + + RefVector up_logdet_list, dn_logdet_list; + RefVector up_dlogdet_list, dn_dlogdet_list; + RefVector up_d2logdet_list, dn_d2logdet_list; + up_logdet_list.reserve(nw); + dn_logdet_list.reserve(nw); + up_dlogdet_list.reserve(nw); + dn_dlogdet_list.reserve(nw); + up_d2logdet_list.reserve(nw); + dn_d2logdet_list.reserve(nw); + + ValueMatrix tmp_val_mat(nelec, this->OrbitalSetSize); + GradMatrix tmp_grad_mat(nelec, this->OrbitalSetSize); + for (int iw = 0; iw < nw; iw++) + { + mw_up_logdet.emplace_back(tmp_val_mat); + up_logdet_list.emplace_back(mw_up_logdet.back()); + mw_dn_logdet.emplace_back(tmp_val_mat); + dn_logdet_list.emplace_back(mw_dn_logdet.back()); + + mw_up_dlogdet.emplace_back(tmp_grad_mat); + up_dlogdet_list.emplace_back(mw_up_dlogdet.back()); + mw_dn_dlogdet.emplace_back(tmp_grad_mat); + dn_dlogdet_list.emplace_back(mw_dn_dlogdet.back()); + + mw_up_d2logdet.emplace_back(tmp_val_mat); + up_d2logdet_list.emplace_back(mw_up_d2logdet.back()); + mw_dn_d2logdet.emplace_back(tmp_val_mat); + dn_d2logdet_list.emplace_back(mw_dn_d2logdet.back()); + } + + up_spo_leader.mw_evaluate_notranspose(up_spo_list, P_list, first, last, up_logdet_list, up_dlogdet_list, + up_d2logdet_list); + dn_spo_leader.mw_evaluate_notranspose(dn_spo_list, P_list, first, last, dn_logdet_list, dn_dlogdet_list, + dn_d2logdet_list); + + for (int iw = 0; iw < nw; iw++) + for (int iat = 0; iat < nelec; iat++) + { + ParticleSet::Scalar_t s = P_list[iw].activeSpin(iat); + RealType coss = std::cos(s); + RealType sins = std::sin(s); + T eis(coss, sins); + T emis(coss, -sins); + + for (int no = 0; no < this->OrbitalSetSize; no++) + { + logdet_list[iw].get()(iat, no) = + eis * up_logdet_list[iw].get()(iat, no) + emis * dn_logdet_list[iw].get()(iat, no); + dlogdet_list[iw].get()(iat, no) = + eis * up_dlogdet_list[iw].get()(iat, no) + emis * dn_dlogdet_list[iw].get()(iat, no); + d2logdet_list[iw].get()(iat, no) = + eis * up_d2logdet_list[iw].get()(iat, no) + emis * dn_d2logdet_list[iw].get()(iat, no); + } + } +} + +template +void SpinorSetT::evaluate_notranspose_spin(const ParticleSet& P, + int first, + int last, + ValueMatrix& logdet, + GradMatrix& dlogdet, + ValueMatrix& d2logdet, + ValueMatrix& dspinlogdet) +{ + IndexType nelec = P.getTotalNum(); + + logpsi_work_up.resize(nelec, this->OrbitalSetSize); + logpsi_work_down.resize(nelec, this->OrbitalSetSize); + + dlogpsi_work_up.resize(nelec, this->OrbitalSetSize); + dlogpsi_work_down.resize(nelec, this->OrbitalSetSize); + + d2logpsi_work_up.resize(nelec, this->OrbitalSetSize); + d2logpsi_work_down.resize(nelec, this->OrbitalSetSize); + + spo_up->evaluate_notranspose(P, first, last, logpsi_work_up, dlogpsi_work_up, d2logpsi_work_up); + spo_dn->evaluate_notranspose(P, first, last, logpsi_work_down, dlogpsi_work_down, d2logpsi_work_down); + + + for (int iat = 0; iat < nelec; iat++) + { + ParticleSet::Scalar_t s = P.activeSpin(iat); + + RealType coss(0.0), sins(0.0); + + coss = std::cos(s); + sins = std::sin(s); + + T eis(coss, sins); + T emis(coss, -sins); + T eye(0, 1.0); + + for (int no = 0; no < this->OrbitalSetSize; no++) + { + logdet(iat, no) = eis * logpsi_work_up(iat, no) + emis * logpsi_work_down(iat, no); + dlogdet(iat, no) = eis * dlogpsi_work_up(iat, no) + emis * dlogpsi_work_down(iat, no); + d2logdet(iat, no) = eis * d2logpsi_work_up(iat, no) + emis * d2logpsi_work_down(iat, no); + dspinlogdet(iat, no) = eye * (eis * logpsi_work_up(iat, no) - emis * logpsi_work_down(iat, no)); + } + } +} + +template +void SpinorSetT::evaluate_spin(const ParticleSet& P, int iat, ValueVector& psi, ValueVector& dpsi) +{ + psi_work_up = 0.0; + psi_work_down = 0.0; + + spo_up->evaluateValue(P, iat, psi_work_up); + spo_dn->evaluateValue(P, iat, psi_work_down); + + ParticleSet::Scalar_t s = P.activeSpin(iat); + + RealType coss(0.0), sins(0.0); + + coss = std::cos(s); + sins = std::sin(s); + + T eis(coss, sins); + T emis(coss, -sins); + T eye(0, 1.0); + + psi = eis * psi_work_up + emis * psi_work_down; + dpsi = eye * (eis * psi_work_up - emis * psi_work_down); +} + +template +void SpinorSetT::evaluateGradSource(const ParticleSet& P, + int first, + int last, + const ParticleSet& source, + int iat_src, + GradMatrix& gradphi) +{ + IndexType nelec = P.getTotalNum(); + + GradMatrix gradphi_up(nelec, this->OrbitalSetSize); + GradMatrix gradphi_dn(nelec, this->OrbitalSetSize); + spo_up->evaluateGradSource(P, first, last, source, iat_src, gradphi_up); + spo_dn->evaluateGradSource(P, first, last, source, iat_src, gradphi_dn); + + for (int iat = 0; iat < nelec; iat++) + { + ParticleSet::Scalar_t s = P.activeSpin(iat); + RealType coss = std::cos(s); + RealType sins = std::sin(s); + T eis(coss, sins); + T emis(coss, -sins); + for (int imo = 0; imo < this->OrbitalSetSize; imo++) + gradphi(iat, imo) = gradphi_up(iat, imo) * eis + gradphi_dn(iat, imo) * emis; + } +} + +template +std::unique_ptr> SpinorSetT::makeClone() const +{ + auto myclone = std::make_unique>(this->my_name_); + std::unique_ptr> cloneup(spo_up->makeClone()); + std::unique_ptr> clonedn(spo_dn->makeClone()); + myclone->set_spos(std::move(cloneup), std::move(clonedn)); + return myclone; +} + +template +void SpinorSetT::createResource(ResourceCollection& collection) const +{ + spo_up->createResource(collection); + spo_dn->createResource(collection); + auto index = collection.addResource(std::make_unique()); +} + +template +void SpinorSetT::acquireResource(ResourceCollection& collection, + const RefVectorWithLeader>& spo_list) const +{ + auto [up_spo_list, dn_spo_list] = extractSpinComponentRefList(spo_list); + auto& spo_leader = spo_list.template getCastedLeader>(); + auto& up_spo_leader = up_spo_list.getLeader(); + auto& dn_spo_leader = dn_spo_list.getLeader(); + up_spo_leader.acquireResource(collection, up_spo_list); + dn_spo_leader.acquireResource(collection, dn_spo_list); + spo_leader.mw_res_handle_ = collection.lendResource(); +} + +template +void SpinorSetT::releaseResource(ResourceCollection& collection, + const RefVectorWithLeader>& spo_list) const +{ + auto [up_spo_list, dn_spo_list] = extractSpinComponentRefList(spo_list); + auto& spo_leader = spo_list.template getCastedLeader>(); + auto& up_spo_leader = up_spo_list.getLeader(); + auto& dn_spo_leader = dn_spo_list.getLeader(); + up_spo_leader.releaseResource(collection, up_spo_list); + dn_spo_leader.releaseResource(collection, dn_spo_list); + collection.takebackResource(spo_leader.mw_res_handle_); +} + +template +std::pair>, RefVectorWithLeader>> SpinorSetT::extractSpinComponentRefList( + const RefVectorWithLeader>& spo_list) const +{ + SpinorSetT& spo_leader = spo_list.template getCastedLeader>(); + IndexType nw = spo_list.size(); + SPOSetT& up_spo_leader = *(spo_leader.spo_up); + SPOSetT& dn_spo_leader = *(spo_leader.spo_dn); + RefVectorWithLeader> up_spo_list(up_spo_leader); + RefVectorWithLeader> dn_spo_list(dn_spo_leader); + up_spo_list.reserve(nw); + dn_spo_list.reserve(nw); + for (int iw = 0; iw < nw; iw++) + { + SpinorSetT& spinor = spo_list.template getCastedElement>(iw); + up_spo_list.emplace_back(*(spinor.spo_up)); + dn_spo_list.emplace_back(*(spinor.spo_dn)); + } + return std::make_pair(up_spo_list, dn_spo_list); +} + +template class SpinorSetT>; +template class SpinorSetT>; + +} // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/SpinorSetT.h b/src/QMCWaveFunctions/SpinorSetT.h new file mode 100644 index 0000000000..bc59e610aa --- /dev/null +++ b/src/QMCWaveFunctions/SpinorSetT.h @@ -0,0 +1,229 @@ +////////////////////////////////////////////////////////////////////////////////////// +// This file is distributed under the University of Illinois/NCSA Open Source License. +// See LICENSE file in top directory for details. +// +// Copyright (c) 2022 QMCPACK developers +// +// File developed by: Raymond Clay III, rclay@sandia.gov, Sandia National Laboratories +// Cody A. Melton, cmelton@sandia.gov, Sandia National Laboratories +// +// File created by: Raymond Clay III, rclay@sandia.gov, Sandia National Laboratories +////////////////////////////////////////////////////////////////////////////////////// + +#ifndef QMCPLUSPLUS_SPINORSET_H +#define QMCPLUSPLUS_SPINORSET_H + +#include "QMCWaveFunctions/SPOSetT.h" +#include "ResourceHandle.h" + +namespace qmcplusplus +{ +/** Class for Melton & Mitas style Spinors. + * + */ +template +class SpinorSetT : public SPOSetT +{ +public: + using ValueMatrix = typename SPOSetT::ValueMatrix; + using ValueVector = typename SPOSetT::ValueVector; + using GradMatrix = typename SPOSetT::GradMatrix; + using GradType = typename SPOSetT::GradType; + using GradVector = typename SPOSetT::GradVector; + using OffloadMWVGLArray = Array>; // [VGL, walker, Orbs] + //using OffloadMWVGLArray = typename SPOSetT::template OffloadMWCGLArray; + template + using OffloadMatrix = typename SPOSetT::template OffloadMatrix
; + using RealType = typename SPOSetT::RealType; + using IndexType = OHMMS_INDEXTYPE; + + /** constructor */ + SpinorSetT(const std::string& my_name); + ~SpinorSetT() override; + + std::string getClassName() const override { return "SpinorSet"; } + bool isOptimizable() const override { return spo_up->isOptimizable() || spo_dn->isOptimizable(); } + bool isOMPoffload() const override { return spo_up->isOMPoffload() || spo_dn->isOMPoffload(); } + bool hasIonDerivs() const override { return spo_up->hasIonDerivs() || spo_dn->hasIonDerivs(); } + + //This class is initialized by separately building the up and down channels of the spinor set and + //then registering them. + void set_spos(std::unique_ptr>&& up, std::unique_ptr>&& dn); + + /** set the OrbitalSetSize + * @param norbs number of single-particle orbitals + */ + void setOrbitalSetSize(int norbs) override; + + /** evaluate the values of this spinor set + * @param P current ParticleSet + * @param iat active particle + * @param psi values of the SPO + */ + void evaluateValue(const ParticleSet& P, int iat, ValueVector& psi) override; + + /** evaluate the values, gradients and laplacians of this single-particle orbital set + * @param P current ParticleSet + * @param iat active particle + * @param psi values of the SPO + * @param dpsi gradients of the SPO + * @param d2psi laplacians of the SPO + */ + void evaluateVGL(const ParticleSet& P, int iat, ValueVector& psi, GradVector& dpsi, ValueVector& d2psi) override; + + /** evaluate the values, gradients and laplacians of this single-particle orbital set + * @param P current ParticleSet + * @param iat active particle + * @param psi values of the SPO + * @param dpsi gradients of the SPO + * @param d2psi laplacians of the SPO + * @param dspin spin gradient of the SPO + */ + void evaluateVGL_spin(const ParticleSet& P, + int iat, + ValueVector& psi, + GradVector& dpsi, + ValueVector& d2psi, + ValueVector& dspin) override; + + /** evaluate the values, gradients and laplacians and spin gradient of this single-particle orbital sets of multiple walkers + * @param spo_list the list of SPOSet pointers in a walker batch + * @param P_list the list of ParticleSet pointers in a walker batch + * @param iat active particle + * @param psi_v_list the list of value vector pointers in a walker batch + * @param dpsi_v_list the list of gradient vector pointers in a walker batch + * @param d2psi_v_list the list of laplacian vector pointers in a walker batch + * @param mw_dspin dual matrix of spin gradients. nw x num_orbitals + */ + void mw_evaluateVGLWithSpin(const RefVectorWithLeader>& spo_list, + const RefVectorWithLeader& P_list, + int iat, + const RefVector& psi_v_list, + const RefVector& dpsi_v_list, + const RefVector& d2psi_v_list, + OffloadMatrix& mw_dspin) const override; + + + /** evaluate the values, gradients and laplacians of this single-particle orbital sets and determinant ratio + * and grads of multiple walkers. Device data of phi_vgl_v must be up-to-date upon return. + * Includes spin gradients + * @param spo_list the list of SPOSet pointers in a walker batch + * @param P_list the list of ParticleSet pointers in a walker batch + * @param iat active particle + * @param phi_vgl_v orbital values, gradients and laplacians of all the walkers + * @param ratios, ratios of all walkers + * @param grads, spatial gradients of all walkers + * @param spingrads, spin gradients of all walkers + */ + void mw_evaluateVGLandDetRatioGradsWithSpin(const RefVectorWithLeader>& spo_list, + const RefVectorWithLeader& P_list, + int iat, + const std::vector& invRow_ptr_list, + OffloadMWVGLArray& phi_vgl_v, + std::vector& ratios, + std::vector& grads, + std::vector& spingrads) const override; + + /** evaluate the values, gradients and laplacians of this single-particle orbital for [first,last) particles + * @param P current ParticleSet + * @param first starting index of the particles + * @param last ending index of the particles + * @param logdet determinant matrix to be inverted + * @param dlogdet gradients + * @param d2logdet laplacians + * + */ + void evaluate_notranspose(const ParticleSet& P, + int first, + int last, + ValueMatrix& logdet, + GradMatrix& dlogdet, + ValueMatrix& d2logdet) override; + + void mw_evaluate_notranspose(const RefVectorWithLeader>& spo_list, + const RefVectorWithLeader& P_list, + int first, + int last, + const RefVector& logdet_list, + const RefVector& dlogdet_list, + const RefVector& d2logdet_list) const override; + + void evaluate_notranspose_spin(const ParticleSet& P, + int first, + int last, + ValueMatrix& logdet, + GradMatrix& dlogdet, + ValueMatrix& d2logdet, + ValueMatrix& dspinlogdet) override; + /** Evaluate the values, spin gradients, and spin laplacians of single particle spinors corresponding to electron iat. + * @param P current particle set. + * @param iat electron index. + * @param spinor values. + * @param spin gradient values. d/ds phi(r,s). + * + */ + void evaluate_spin(const ParticleSet& P, int iat, ValueVector& psi, ValueVector& dpsi) override; + + /** evaluate the gradients of this single-particle orbital + * for [first,last) target particles with respect to the given source particle + * @param P current ParticleSet + * @param first starting index of the particles + * @param last ending index of the particles + * @param iat_src source particle index + * @param gradphi gradients + * + */ + virtual void evaluateGradSource(const ParticleSet& P, + int first, + int last, + const ParticleSet& source, + int iat_src, + GradMatrix& gradphi) override; + + std::unique_ptr> makeClone() const override; + + void createResource(ResourceCollection& collection) const override; + + void acquireResource(ResourceCollection& collection, const RefVectorWithLeader>& spo_list) const override; + + void releaseResource(ResourceCollection& collection, const RefVectorWithLeader>& spo_list) const override; + + /// check if the multi walker resource is owned. For testing only. + bool isResourceOwned() const { return bool(mw_res_handle_); } + +private: + struct SpinorSetMultiWalkerResource; + ResourceHandle mw_res_handle_; + + std::pair>, RefVectorWithLeader>> extractSpinComponentRefList( + const RefVectorWithLeader>& spo_list) const; + + //Sposet for the up and down channels of our spinors. + std::unique_ptr> spo_up; + std::unique_ptr> spo_dn; + + //temporary arrays for holding the values of the up and down channels respectively. + ValueVector psi_work_up; + ValueVector psi_work_down; + + //temporary arrays for holding the gradients of the up and down channels respectively. + GradVector dpsi_work_up; + GradVector dpsi_work_down; + + //temporary arrays for holding the laplacians of the up and down channels respectively. + ValueVector d2psi_work_up; + ValueVector d2psi_work_down; + + //Same as above, but these are the full matrices containing all spinor/particle combinations. + ValueMatrix logpsi_work_up; + ValueMatrix logpsi_work_down; + + GradMatrix dlogpsi_work_up; + GradMatrix dlogpsi_work_down; + + ValueMatrix d2logpsi_work_up; + ValueMatrix d2logpsi_work_down; +}; + +} // namespace qmcplusplus +#endif diff --git a/src/QMCWaveFunctions/tests/CMakeLists.txt b/src/QMCWaveFunctions/tests/CMakeLists.txt index ec066f8735..5c9484da4c 100644 --- a/src/QMCWaveFunctions/tests/CMakeLists.txt +++ b/src/QMCWaveFunctions/tests/CMakeLists.txt @@ -139,7 +139,7 @@ set(DETERMINANT_SRC test_ci_configuration.cpp test_multi_slater_determinant.cpp) -add_library(sposets_for_testing FakeSPO.cpp ConstantSPOSet.cpp) +add_library(sposets_for_testing FakeSPOT.cpp FakeSPO.cpp ConstantSPOSet.cpp) target_include_directories(sposets_for_testing PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}") target_link_libraries(sposets_for_testing PUBLIC qmcwfs) diff --git a/src/QMCWaveFunctions/tests/FakeSPOT.cpp b/src/QMCWaveFunctions/tests/FakeSPOT.cpp new file mode 100644 index 0000000000..fcf1637682 --- /dev/null +++ b/src/QMCWaveFunctions/tests/FakeSPOT.cpp @@ -0,0 +1,160 @@ +////////////////////////////////////////////////////////////////////////////////////// +// 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: Mark Dewing, mdewing@anl.gov, Argonne National Laboratory +// +// File created by: Mark Dewing, mdewing@anl.gov, Argonne National Laboratory +////////////////////////////////////////////////////////////////////////////////////// + +#include "FakeSPOT.h" + +namespace qmcplusplus +{ +template +FakeSPOT::FakeSPOT() : SPOSetT("one_FakeSPO") +{ + a.resize(3, 3); + + a(0, 0) = 2.3; + a(0, 1) = 4.5; + a(0, 2) = 2.6; + a(1, 0) = 0.5; + a(1, 1) = 8.5; + a(1, 2) = 3.3; + a(2, 0) = 1.8; + a(2, 1) = 4.4; + a(2, 2) = 4.9; + + v.resize(3); + v[0] = 1.9; + v[1] = 2.0; + v[2] = 3.1; + + + a2.resize(4, 4); + a2(0, 0) = 2.3; + a2(0, 1) = 4.5; + a2(0, 2) = 2.6; + a2(0, 3) = 1.2; + a2(1, 0) = 0.5; + a2(1, 1) = 8.5; + a2(1, 2) = 3.3; + a2(1, 3) = 0.3; + a2(2, 0) = 1.8; + a2(2, 1) = 4.4; + a2(2, 2) = 4.9; + a2(2, 3) = 2.8; + a2(3, 0) = 0.8; + a2(3, 1) = 4.1; + a2(3, 2) = 3.2; + a2(3, 3) = 1.1; + + v2.resize(4, 4); + + v2(0, 0) = 3.2; + v2(0, 1) = 0.5; + v2(0, 2) = 5.9; + v2(0, 3) = 3.7; + v2(1, 0) = 0.3; + v2(1, 1) = 1.4; + v2(1, 2) = 3.9; + v2(1, 3) = 8.2; + v2(2, 0) = 3.3; + v2(2, 1) = 5.4; + v2(2, 2) = 4.9; + v2(2, 3) = 2.2; + v2(3, 1) = 5.4; + v2(3, 2) = 4.9; + v2(3, 3) = 2.2; + + gv.resize(4); + gv[0] = GradType(1.0, 0.0, 0.1); + gv[1] = GradType(1.0, 2.0, 0.1); + gv[2] = GradType(2.0, 1.0, 0.1); + gv[3] = GradType(0.4, 0.3, 0.1); +} +template +std::unique_ptr> FakeSPOT::makeClone() const +{ + return std::make_unique>(*this); +} + +template +void FakeSPOT::setOrbitalSetSize(int norbs) +{ + this->OrbitalSetSize = norbs; +} + +template +void FakeSPOT::evaluateValue(const ParticleSet& P, int iat, ValueVector& psi) +{ + if (iat < 0) + for (int i = 0; i < psi.size(); i++) + psi[i] = 1.2 * i - i * i; + else if (this->OrbitalSetSize == 3) + for (int i = 0; i < 3; i++) + psi[i] = a(iat, i); + else if (this->OrbitalSetSize == 4) + for (int i = 0; i < 4; i++) + psi[i] = a2(iat, i); +} + +template +void FakeSPOT::evaluateVGL(const ParticleSet& P, int iat, ValueVector& psi, GradVector& dpsi, ValueVector& d2psi) +{ + if (this->OrbitalSetSize == 3) + { + for (int i = 0; i < 3; i++) + { + psi[i] = v[i]; + dpsi[i] = gv[i]; + } + } + else if (this->OrbitalSetSize == 4) + { + for (int i = 0; i < 4; i++) + { + psi[i] = v2(iat, i); + dpsi[i] = gv[i]; + } + } +} + +template +void FakeSPOT::evaluate_notranspose(const ParticleSet& P, + int first, + int last, + ValueMatrix& logdet, + GradMatrix& dlogdet, + ValueMatrix& d2logdet) +{ + if (this->OrbitalSetSize == 3) + { + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + { + logdet(j, i) = a(i, j); + dlogdet[i][j] = gv[j] + GradType(i); + } + } + else if (this->OrbitalSetSize == 4) + { + for (int i = 0; i < 4; i++) + for (int j = 0; j < 4; j++) + { + logdet(j, i) = a2(i, j); + dlogdet[i][j] = gv[j] + GradType(i); + } + } +} + +// Class concrete types from ValueType +template class FakeSPOT; +template class FakeSPOT; +template class FakeSPOT>; +template class FakeSPOT>; + +} // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/tests/FakeSPOT.h b/src/QMCWaveFunctions/tests/FakeSPOT.h new file mode 100644 index 0000000000..ee452842f5 --- /dev/null +++ b/src/QMCWaveFunctions/tests/FakeSPOT.h @@ -0,0 +1,62 @@ +////////////////////////////////////////////////////////////////////////////////////// +// 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: Mark Dewing, mdewing@anl.gov, Argonne National Laboratory +// +// File created by: Mark Dewing, mdewing@anl.gov, Argonne National Laboratory +////////////////////////////////////////////////////////////////////////////////////// + + +#ifndef QMCPLUSPLUS_FAKESPOT_H +#define QMCPLUSPLUS_FAKESPOT_H + +#include "QMCWaveFunctions/SPOSetT.h" + +namespace qmcplusplus +{ +template +class FakeSPOT : public SPOSetT +{ +public: + Matrix a; + Matrix a2; + Vector v; + Matrix v2; + + using ValueVector = typename SPOSetT::ValueVector; + using ValueMatrix = typename SPOSetT::ValueMatrix; + using GradVector = typename SPOSetT::GradVector; + using GradMatrix = typename SPOSetT::GradMatrix; + using GradType = typename SPOSetT::GradType; + + typename SPOSetT::GradVector gv; + + FakeSPOT(); + + ~FakeSPOT() override = default; + + std::string getClassName() const override { return "FakeSPO"; } + + std::unique_ptr> makeClone() const override; + + virtual void report() {} + + void setOrbitalSetSize(int norbs) override; + + void evaluateValue(const ParticleSet& P, int iat, ValueVector& psi) override; + + void evaluateVGL(const ParticleSet& P, int iat, ValueVector& psi, GradVector& dpsi, ValueVector& d2psi) override; + + void evaluate_notranspose(const ParticleSet& P, + int first, + int last, + ValueMatrix& logdet, + GradMatrix& dlogdet, + ValueMatrix& d2logdet) override; +}; + +} // namespace qmcplusplus +#endif diff --git a/src/QMCWaveFunctions/tests/test_RotatedSPOs.cpp b/src/QMCWaveFunctions/tests/test_RotatedSPOs.cpp index c162b9985c..b126cfec86 100644 --- a/src/QMCWaveFunctions/tests/test_RotatedSPOs.cpp +++ b/src/QMCWaveFunctions/tests/test_RotatedSPOs.cpp @@ -20,6 +20,7 @@ #include "QMCWaveFunctions/WaveFunctionComponent.h" #include "QMCWaveFunctions/EinsplineSetBuilder.h" #include "QMCWaveFunctions/RotatedSPOs.h" +#include "QMCWaveFunctions/RotatedSPOsT.h" #include "QMCWaveFunctions/SPOSetT.h" #include "checkMatrix.hpp" #include "FakeSPO.h" @@ -651,7 +652,17 @@ opt_variables_type& getMyVars(SPOSetT& rot) { return rot.myVars; } opt_variables_type& getMyVars(SPOSetT>& rot) { return rot.myVars; } opt_variables_type& getMyVars(SPOSetT>& rot) { return rot.myVars; } opt_variables_type& getMyVarsFull(RotatedSPOs& rot) { return rot.myVarsFull; } +opt_variables_type& getMyVarsFull(RotatedSPOsT& rot) { return rot.myVarsFull; } +opt_variables_type& getMyVarsFull(RotatedSPOsT& rot) { return rot.myVarsFull; } std::vector>& getHistoryParams(RotatedSPOs& rot) { return rot.history_params_; } +std::vector>& getHistoryParams(RotatedSPOsT& rot) +{ + return rot.history_params_; +} +std::vector>& getHistoryParams(RotatedSPOsT& rot) +{ + return rot.history_params_; +} } // namespace testing // Test using global rotation @@ -706,6 +717,58 @@ TEST_CASE("RotatedSPOs read and write parameters", "[wavefunction]") CHECK(full_var[5] == Approx(0.0)); } +// Test using global rotation +TEMPLATE_TEST_CASE("RotatedSPOs read and write parameters", "[wavefunction][template]", double, float) +{ + auto fake_spo = std::make_unique>(); + fake_spo->setOrbitalSetSize(4); + RotatedSPOsT rot("fake_rot", std::move(fake_spo)); + int nel = 2; + rot.buildOptVariables(nel); + + optimize::VariableSet vs; + rot.checkInVariablesExclusive(vs); + vs[0] = 0.1; + vs[1] = 0.15; + vs[2] = 0.2; + vs[3] = 0.25; + rot.resetParametersExclusive(vs); + + { + hdf_archive hout; + vs.writeToHDF("rot_vp.h5", hout); + + rot.writeVariationalParameters(hout); + } + + auto fake_spo2 = std::make_unique>(); + fake_spo2->setOrbitalSetSize(4); + + RotatedSPOsT rot2("fake_rot", std::move(fake_spo2)); + rot2.buildOptVariables(nel); + + optimize::VariableSet vs2; + rot2.checkInVariablesExclusive(vs2); + + hdf_archive hin; + vs2.readFromHDF("rot_vp.h5", hin); + rot2.readVariationalParameters(hin); + + opt_variables_type& var = testing::getMyVars(rot2); + CHECK(var[0] == Approx(vs[0])); + CHECK(var[1] == Approx(vs[1])); + CHECK(var[2] == Approx(vs[2])); + CHECK(var[3] == Approx(vs[3])); + + opt_variables_type& full_var = testing::getMyVarsFull(rot2); + CHECK(full_var[0] == Approx(vs[0])); + CHECK(full_var[1] == Approx(vs[1])); + CHECK(full_var[2] == Approx(vs[2])); + CHECK(full_var[3] == Approx(vs[3])); + CHECK(full_var[4] == Approx(0.0)); + CHECK(full_var[5] == Approx(0.0)); +} + // Test using history list. TEST_CASE("RotatedSPOs read and write parameters history", "[wavefunction]") {