From c05b561baf7894d0a740511007f679f3960e1662 Mon Sep 17 00:00:00 2001 From: William F Godoy Date: Wed, 16 Aug 2023 10:21:53 -0400 Subject: [PATCH] Refactor FreeOrbital class Base typed aliases on SPOSet on OrbitalSetTraits --- src/QMCWaveFunctions/CMakeLists.txt | 2 +- .../ElectronGas/FreeOrbitalT.cpp | 683 ++++++++++++++++++ .../ElectronGas/FreeOrbitalT.h | 88 +++ src/QMCWaveFunctions/SPOSetT.h | 3 + 4 files changed, 775 insertions(+), 1 deletion(-) create mode 100644 src/QMCWaveFunctions/ElectronGas/FreeOrbitalT.cpp create mode 100644 src/QMCWaveFunctions/ElectronGas/FreeOrbitalT.h diff --git a/src/QMCWaveFunctions/CMakeLists.txt b/src/QMCWaveFunctions/CMakeLists.txt index e43ce7e02a2..f92a3f2a149 100644 --- a/src/QMCWaveFunctions/CMakeLists.txt +++ b/src/QMCWaveFunctions/CMakeLists.txt @@ -63,7 +63,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) diff --git a/src/QMCWaveFunctions/ElectronGas/FreeOrbitalT.cpp b/src/QMCWaveFunctions/ElectronGas/FreeOrbitalT.cpp new file mode 100644 index 00000000000..fe9e9eff058 --- /dev/null +++ b/src/QMCWaveFunctions/ElectronGas/FreeOrbitalT.cpp @@ -0,0 +1,683 @@ +////////////////////////////////////////////////////////////////////////////////////// +// 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 +{ + +//Constructors need to be fully specialized +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) +{ + 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) +{ + 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) +{ + 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) +{ + OrbitalSetSize = maxk; // SPOSet member + for (int ik = 0; ik < maxk; ik++) + k2neg[ik] = -dot(kvecs[ik], kvecs[ik]); +} + + +template +FreeOrbitalT::~FreeOrbitalT() +{} + +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, + 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::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) +{ + 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]; + } + } + } +} + +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 00000000000..c73ab26a2af --- /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/SPOSetT.h b/src/QMCWaveFunctions/SPOSetT.h index 95643fd3c54..2b7e2c7ca0b 100644 --- a/src/QMCWaveFunctions/SPOSetT.h +++ b/src/QMCWaveFunctions/SPOSetT.h @@ -65,6 +65,9 @@ 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; template using OffloadMatrix = Matrix>;