From cabbe466772edd823751f4ec47617bf158161eb4 Mon Sep 17 00:00:00 2001 From: Philip Fackler Date: Wed, 16 Aug 2023 12:08:51 -0400 Subject: [PATCH] PWOribitalSetT and PWBasisT Add FullRealType in SPOSet and RotatedSPOs Move generic definition after specialization add implicit implementations Fix some errors initial commit of templated PWOribitSetT that compiles cleanup templateitze PWBasis as well, as is dependancy remove inaccurate comment remove polluted commit --- src/QMCWaveFunctions/CMakeLists.txt | 4 +- src/QMCWaveFunctions/PlaneWave/PWBasisT.cpp | 197 ++++++++++ src/QMCWaveFunctions/PlaneWave/PWBasisT.h | 343 ++++++++++++++++++ .../PlaneWave/PWOrbitalSetT.cpp | 145 ++++++++ .../PlaneWave/PWOrbitalSetT.h | 130 +++++++ 5 files changed, 817 insertions(+), 2 deletions(-) create mode 100644 src/QMCWaveFunctions/PlaneWave/PWBasisT.cpp create mode 100644 src/QMCWaveFunctions/PlaneWave/PWBasisT.h create mode 100644 src/QMCWaveFunctions/PlaneWave/PWOrbitalSetT.cpp create mode 100644 src/QMCWaveFunctions/PlaneWave/PWOrbitalSetT.h diff --git a/src/QMCWaveFunctions/CMakeLists.txt b/src/QMCWaveFunctions/CMakeLists.txt index 6ad7f54cb5..1bbfc0520f 100644 --- a/src/QMCWaveFunctions/CMakeLists.txt +++ b/src/QMCWaveFunctions/CMakeLists.txt @@ -111,9 +111,9 @@ if(OHMMS_DIM MATCHES 3) endif(HAVE_EINSPLINE) # plane wave SPO - set(FERMION_SRCS ${FERMION_SRCS} PlaneWave/PWBasis.cpp PlaneWave/PWParameterSet.cpp PlaneWave/PWOrbitalBuilder.cpp) + set(FERMION_SRCS ${FERMION_SRCS} PlaneWave/PWBasis.cpp PlaneWave/PWBasisT.cpp PlaneWave/PWParameterSet.cpp PlaneWave/PWOrbitalBuilder.cpp) if(QMC_COMPLEX) - set(FERMION_SRCS ${FERMION_SRCS} PlaneWave/PWOrbitalSet.cpp) + set(FERMION_SRCS ${FERMION_SRCS} PlaneWave/PWOrbitalSet.cpp PlaneWave/PWOrbitalSetT.cpp) else() set(FERMION_SRCS ${FERMION_SRCS} PlaneWave/PWRealOrbitalSet.cpp PlaneWave/PWRealOrbitalSetT.cpp) endif(QMC_COMPLEX) diff --git a/src/QMCWaveFunctions/PlaneWave/PWBasisT.cpp b/src/QMCWaveFunctions/PlaneWave/PWBasisT.cpp new file mode 100644 index 0000000000..fe00655309 --- /dev/null +++ b/src/QMCWaveFunctions/PlaneWave/PWBasisT.cpp @@ -0,0 +1,197 @@ +////////////////////////////////////////////////////////////////////////////////////// +// 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 PWBasisT.cpp + * @brief Definition of member functions of Plane-wave basis set + */ +#include "PWBasisT.h" + +namespace qmcplusplus +{ +template +int PWBasisT::readbasis(hdf_archive& h5basisgroup, + RealType ecutoff, + const ParticleLayout& lat, + const std::string& pwname, + const std::string& pwmultname, + bool resizeContainer) +{ + ///make a local copy + Lattice = lat; + ecut = ecutoff; + app_log() << " PWBasisT::" << pwmultname << " is found " << std::endl; + h5basisgroup.read(gvecs, "/electrons/kpoint_0/gvectors"); + NumPlaneWaves = std::max(gvecs.size(), kplusgvecs_cart.size()); + if (NumPlaneWaves == 0) + { + app_error() << " PWBasisT::readbasis Basis is missing. Abort " << std::endl; + abort(); //FIX_ABORT + } + if (kplusgvecs_cart.empty()) + { + kplusgvecs_cart.resize(NumPlaneWaves); + for (int i = 0; i < NumPlaneWaves; i++) + kplusgvecs_cart[i] = Lattice.k_cart(gvecs[i]); + } + //app_log() << " Gx Gy Gz " << std::endl; + //for(int i=0; i(std::cout,"\n")); + return NumPlaneWaves; +} + +template +void PWBasisT::setTwistAngle(const PosType& tang) +{ + PosType dang = twist - tang; + bool sameTwist = dot(dang, dang) < std::numeric_limits::epsilon(); + if (maxmaxg && sameTwist) + return; + twist = tang; + reset(); +} + +template +void PWBasisT::reset() +{ + trimforecut(); + //logC.resize(3,2*maxmaxg+1); + Z.resize(NumPlaneWaves, 2 + DIM); + Zv.resize(NumPlaneWaves); + phi.resize(NumPlaneWaves); +} + +/** Remove basis elements if kinetic energy > ecut. + * + * Keep and indexmap so we know how to match coefficients on read. + */ +template +void PWBasisT::trimforecut() +{ + //Convert the twist angle to Cartesian coordinates. + twist_cart = Lattice.k_cart(twist); + inputmap.resize(NumPlaneWaves); + app_log() << " PWBasisT::TwistAngle (unit) =" << twist << std::endl; + app_log() << " PWBasisT::TwistAngle (cart) =" << twist_cart << std::endl; + app_log() << " PWBasisT::trimforecut NumPlaneWaves (before) =" << NumPlaneWaves << std::endl; + std::vector gvecCopy(gvecs); + std::vector gcartCopy(kplusgvecs_cart); + gvecs.clear(); + kplusgvecs_cart.clear(); + minusModKplusG2.reserve(NumPlaneWaves); + // RealType kcutoff2 = 2.0*ecut; //std::sqrt(2.0*ecut); + int ngIn = NumPlaneWaves; + for (int ig = 0, newig = 0; ig < ngIn; ig++) + { + //PosType tempvec = Lattice.k_cart(gvecCopy[ig]+twist); + PosType tempvec = gcartCopy[ig] + twist_cart; + RealType mod2 = dot(tempvec, tempvec); + + // Keep all the g-vectors + // The cutoff energy is not stored in the HDF file now. + // Is truncating the gvectors to a spherical shell necessary? + if (true) + { + gvecs.push_back(gvecCopy[ig]); + kplusgvecs_cart.push_back(tempvec); + minusModKplusG2.push_back(-mod2); + //Remember which position in the HDF5 file this came from...for coefficients + inputmap[ig] = newig++; + } +#if 0 + if(mod2<=kcutoff2) + { + gvecs.push_back(gvecCopy[ig]); + kplusgvecs_cart.push_back(tempvec); + minusModKplusG2.push_back(-mod2); + //Remember which position in the HDF5 file this came from...for coefficients + inputmap[ig] = newig++; + } + else + { + inputmap[ig] = -1; //Temporary value...need to know final NumPlaneWaves. + NumPlaneWaves--; + } +#endif + } +#if defined(PWBasisT_USE_RECURSIVE) + //Store the maximum number of translations, within ecut, of any reciprocal cell vector. + for (int ig = 0; ig < NumPlaneWaves; ig++) + for (int i = 0; i < OHMMS_DIM; i++) + if (std::abs(gvecs[ig][i]) > maxg[i]) + maxg[i] = std::abs(gvecs[ig][i]); + gvecs_shifted.resize(NumPlaneWaves); + for (int ig = 0; ig < NumPlaneWaves; ig++) + gvecs_shifted[ig] = gvecs[ig] + maxg; + maxmaxg = std::max(maxg[0], std::max(maxg[1], maxg[2])); + //changes the order???? ok + C.resize(3, 2 * maxmaxg + 2); +#else + maxmaxg = 1; +#endif + // //make a copy of input to gvecCopy + //// for(int ig=0, newig=0; ig 0) + //// negative.push_back(1); + //// else { //gx == 0, test gy + //// if(gvecCopy[ig][1] < 0) + //// negative.push_back(0); + //// else if(gvecCopy[ig][1] > 0) + //// negative.push_back(1); + //// else { //gx == gy == 0; test gz. If gz==0 also, take negative=1 (arbitrary) + //// if(gvecCopy[ig][2] < 0) + //// negative.push_back(0); + //// else + //// negative.push_back(1); + //// } + //// } + ////#endif + // } else { + // inputmap[ig] = -1; //Temporary value...need to know final NumPlaneWaves. + // NumPlaneWaves--; + // } + // } + //Finalize the basis. Fix temporary values of inputmap. + //for(int ig=0; igecut + app_log() << " NumPlaneWaves (after) =" << NumPlaneWaves << std::endl; +} +// template class PWBasisT; +// template class PWBasisT; +template class PWBasisT>; +template class PWBasisT>; +} // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/PlaneWave/PWBasisT.h b/src/QMCWaveFunctions/PlaneWave/PWBasisT.h new file mode 100644 index 0000000000..54592b9ba7 --- /dev/null +++ b/src/QMCWaveFunctions/PlaneWave/PWBasisT.h @@ -0,0 +1,343 @@ +////////////////////////////////////////////////////////////////////////////////////// +// 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 +// +// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign +////////////////////////////////////////////////////////////////////////////////////// + + +/** @file PWBasis.h + * @brief Declaration of Plane-wave basis set + */ +#ifndef QMCPLUSPLUS_PLANEWAVEBASIS_BLAS_H +#define QMCPLUSPLUS_PLANEWAVEBASIS_BLAS_H + +#include "Configuration.h" +#include "Particle/ParticleSet.h" +#include "Message/Communicate.h" +#include "type_traits/complex_help.hpp" +#include "CPU/e2iphi.h" +#include "hdf/hdf_archive.h" + +/** If defined, use recursive method to build the basis set for each position + * + * performance improvement is questionable: load vs sin/cos + */ +//#define PWBASIS_USE_RECURSIVE + +namespace qmcplusplus +{ +/** Plane-wave basis set + * + * Rewrite of PlaneWaveBasis to utilize blas II or III + * Support more general input tags + */ +template +class PWBasisT : public QMCTraits +{ +public: + using RealType = typename RealAlias_impl::value_type; + using ComplexType = T; + using PosType = TinyVector; + using IndexType = QMCTraits::IndexType; + using ParticleLayout = ParticleSet::ParticleLayout; + using GIndex_t = TinyVector; + +private: + ///max of maxg[i] + int maxmaxg; + //Need to store the maximum translation in each dimension to use recursive PW generation. + GIndex_t maxg; + //The PlaneWave data - keep all of these strictly private to prevent inconsistencies. + RealType ecut; + ///twist angle in reduced + PosType twist; + ///twist angle in cartesian + PosType twist_cart; //Twist angle in reduced and Cartesian. + + ///gvecs in reduced coordiates + std::vector gvecs; + ///Reduced coordinates with offset gvecs_shifted[][idim]=gvecs[][idim]+maxg[idim] + std::vector gvecs_shifted; + + std::vector minusModKplusG2; + std::vector kplusgvecs_cart; //Cartesian. + + Matrix C; + //Real wavefunctions here. Now the basis states are cos(Gr) or sin(Gr), not exp(iGr) + //We need a way of switching between them for G -> -G, otherwise the + //determinant will have multiple rows that are equal (to within a constant factor) + //of others, giving a zero determinant. For this, we build a vector (negative) which + //stores whether a vector is "+" or "-" (with some criterion, to be defined). We + //the switch from cos() to sin() based on the value of this input. + std::vector negative; + +public: + //enumeration for the value, laplacian, gradients and size + enum + { + PW_VALUE, + PW_LAP, + PW_GRADX, + PW_GRADY, + PW_GRADZ, + PW_MAXINDEX + }; + + Matrix Z; + + Vector Zv; + /* inputmap is used for a memory efficient way of + * + * importing the basis-set and coefficients when the desired energy cutoff may be + * lower than that represented by all data in the wavefunction input file. + * The steps taken are: + * - Read all basis data. + * - Create map. inputmap[i] = j; j is correct PW index, i is input coef index. + * For basis elements outside cutoff, inputmap[i] = gvecs.size(); + * - Coefficients are in same order as PWs in inputfile => simply file into + * storage matrix using the map as the input. All excess coefficients are + * put into [gvecs.size()] and not used. i.e. coefs need to be allocated 1 higher. + * Such an approach is not needed for Gamma-point only calculations because the + * basis is spherically ordered. However, when a twist-angle is used, the "sphere" + * of allowed planewaves is shifted. + */ + + Vector phi; + + std::vector inputmap; + + ///total number of basis functions + int NumPlaneWaves; + + ///local copy of Lattice + ParticleLayout Lattice; + + ///default constructor + PWBasisT() : maxmaxg(0), NumPlaneWaves(0) {} + + ///constructor + PWBasisT(const PosType& twistangle) : maxmaxg(0), twist(twistangle), NumPlaneWaves(0) {} + + ~PWBasisT() {} + + ///set the twist angle + void setTwistAngle(const PosType& tang); + + ///reset + void reset(); + + /** Read basisset from hdf5 file. Apply ecut. + * @param h5basisgroup h5 node where basis is located + * @param ecutoff cutoff energy + * @param lat CrystalLattice + * @param resizeContainer if true, resize internal storage. + * @return the number of plane waves + */ + int readbasis(hdf_archive& h5basisgroup, + RealType ecutoff, + const ParticleLayout& lat, + const std::string& pwname = "planewaves", + const std::string& pwmultname = "multipliers", + bool resizeContainer = true); + + /** Remove basis elements if kinetic energy > ecut. + * + * Keep and indexmap so we know how to match coefficients on read. + */ + void trimforecut(); + +#if defined(PWBASIS_USE_RECURSIVE) + /** Fill the recursion coefficients matrix. + * + * @todo Generalize to non-orthorohmbic cells + */ + inline void BuildRecursionCoefs(const PosType& pos) + { + PosType tau_red(Lattice.toUnit(pos)); +// RealType phi=TWOPI*tau_red[0]; +// RealType nphi=maxg0*phi; +// ComplexType ct0(std::cos(phi),std::sin(phi)); +// ComplexType t(std::cos(nphi),-std::sin(nphi)); +// C0[0]=t; +// for(int n=1; n<=2*maxg0; n++) C0[n] = (t *= ct0); +// +// phi=TWOPI*tau_red[1]; +// nphi=maxg1*phi; +// ct0=ComplexType(std::cos(phi),std::sin(phi)); +// t=ComplexType(std::cos(nphi),-std::sin(nphi)); +// C1[0]=t; +// for(int n=1; n<=2*maxg1; n++) C1[n] = (t *= ct0); +// +// phi=TWOPI*tau_red[2]; +// nphi=maxg2*phi; +// ct0=ComplexType(std::cos(phi),std::sin(phi)); +// t=ComplexType(std::cos(nphi),-std::sin(nphi)); +// C2[0]=t; +// for(int n=1; n<=2*maxg2; n++) C2[n] = (t *= ct0); +#pragma ivdep + for (int idim = 0; idim < 3; idim++) + { + int ng = maxg[idim]; + RealType phi = TWOPI * tau_red[idim]; + RealType nphi = ng * phi; + ComplexType Ctemp(std::cos(phi), std::sin(phi)); + ComplexType t(std::cos(nphi), -std::sin(nphi)); + ComplexType* restrict cp_ptr = C[idim]; + *cp_ptr++ = t; + for (int n = 1; n <= 2 * ng; n++) + { + *cp_ptr++ = (t *= Ctemp); + } + } + //Base version + //#pragma ivdep + // for(int idim=0; idim<3; idim++){ + // RealType phi=TWOPI*tau_red[idim]; + // ComplexType Ctemp(std::cos(phi),std::sin(phi)); + // int ng=maxg[idim]; + // ComplexType* restrict cp_ptr=C[idim]+ng; + // ComplexType* restrict cn_ptr=C[idim]+ng-1; + // *cp_ptr=1.0; + // for(int n=1; n<=ng; n++,cn_ptr--){ + // ComplexType t(Ctemp*(*cp_ptr++)); + // *cp_ptr = t; + // *cn_ptr = conj(t); + // } + // } + //Not valid for general supercell + // // Cartesian of twist for 1,1,1 (reduced coordinates) + // PosType G111(1.0,1.0,1.0); + // G111 = Lattice.k_cart(G111); + // + // //Precompute a small number of complex factors (PWs along b1,b2,b3 lines) + // //using a fast recursion algorithm + //#pragma ivdep + // for(int idim=0; idim<3; idim++){ + // //start the recursion with the 111 vector. + // RealType phi = pos[idim] * G111[idim]; + // register ComplexType Ctemp(std::cos(phi), std::sin(phi)); + // int ng=maxg[idim]; + // ComplexType* restrict cp_ptr=C[idim]+ng; + // ComplexType* restrict cn_ptr=C[idim]+ng-1; + // *cp_ptr=1.0; + // for(int n=1; n<=ng; n++,cn_ptr--){ + // ComplexType t(Ctemp*(*cp_ptr++)); + // *cp_ptr = t; + // *cn_ptr = conj(t); + // } + // } + } + + inline void evaluate(const PosType& pos) + { + BuildRecursionCoefs(pos); + RealType twistdotr = dot(twist_cart, pos); + ComplexType pw0(std::cos(twistdotr), std::sin(twistdotr)); + //Evaluate the planewaves for particle iat. + for (int ig = 0; ig < NumPlaneWaves; ig++) + { + //PW is initialized as exp(i*twist.r) so that the final basis evaluations are for (twist+G).r + ComplexType pw(pw0); //std::cos(twistdotr),std::sin(twistdotr)); + for (int idim = 0; idim < 3; idim++) + pw *= C(idim, gvecs_shifted[ig][idim]); + //pw *= C0[gvecs_shifted[ig][0]]; + //pw *= C1[gvecs_shifted[ig][1]]; + //pw *= C2[gvecs_shifted[ig][2]]; + Zv[ig] = pw; + } + } + /** Evaluate all planewaves and derivatives for the iat-th particle + * + * The basis functions are evaluated for particles iat: first <= iat < last + * Evaluate the plane-waves at current particle coordinates using a fast + * recursion algorithm. Order of Y,dY and d2Y is kept correct. + * These can be "dotted" with coefficients later to complete orbital evaluations. + */ + inline void evaluateAll(const ParticleSet& P, int iat) + { + const PosType& r(P.activeR(iat)); + BuildRecursionCoefs(r); + RealType twistdotr = dot(twist_cart, r); + ComplexType pw0(std::cos(twistdotr), std::sin(twistdotr)); + //Evaluate the planewaves and derivatives. + ComplexType* restrict zptr = Z.data(); + for (int ig = 0; ig < NumPlaneWaves; ig++, zptr += 5) + { + //PW is initialized as exp(i*twist.r) so that the final basis evaluations + //are for (twist+G).r + ComplexType pw(pw0); + // THE INDEX ORDER OF C DOESN'T LOOK TOO GOOD: this could be fixed + for (int idim = 0; idim < 3; idim++) + pw *= C(idim, gvecs_shifted[ig][idim]); + //pw *= C0[gvecs_shifted[ig][0]]; + //pw *= C1[gvecs_shifted[ig][1]]; + //pw *= C2[gvecs_shifted[ig][2]]; + zptr[0] = pw; + zptr[1] = minusModKplusG2[ig] * pw; + zptr[2] = kplusgvecs_cart[ig][0] * ComplexType(-pw.imag(), pw.real()); + zptr[3] = kplusgvecs_cart[ig][1] * ComplexType(-pw.imag(), pw.real()); + zptr[4] = kplusgvecs_cart[ig][2] * ComplexType(-pw.imag(), pw.real()); + } + } +#else + inline void evaluate(const PosType& pos) + { + //Evaluate the planewaves for particle iat. + for (int ig = 0; ig < NumPlaneWaves; ig++) + phi[ig] = dot(kplusgvecs_cart[ig], pos); + eval_e2iphi(NumPlaneWaves, phi.data(), Zv.data()); + } + inline void evaluateAll(const ParticleSet& P, int iat) + { + const PosType& r(P.activeR(iat)); + evaluate(r); + ComplexType* restrict zptr = Z.data(); + for (int ig = 0; ig < NumPlaneWaves; ig++, zptr += 5) + { + //PW is initialized as exp(i*twist.r) so that the final basis evaluations + //are for (twist+G).r + ComplexType& pw = Zv[ig]; + zptr[0] = pw; + zptr[1] = minusModKplusG2[ig] * pw; + zptr[2] = kplusgvecs_cart[ig][0] * ComplexType(-pw.imag(), pw.real()); + zptr[3] = kplusgvecs_cart[ig][1] * ComplexType(-pw.imag(), pw.real()); + zptr[4] = kplusgvecs_cart[ig][2] * ComplexType(-pw.imag(), pw.real()); + } + } +#endif + // /** Fill the recursion coefficients matrix. + // * + // * @todo Generalize to non-orthorohmbic cells + // */ + // void BuildRecursionCoefsByAdd(const PosType& pos) + // { + // // Cartesian of twist for 1,1,1 (reduced coordinates) + // PosType G111(1.0,1.0,1.0); + // G111 = Lattice.k_cart(G111); + // //PosType redP=P.Lattice.toUnit(P.R[iat]); + // //Precompute a small number of complex factors (PWs along b1,b2,b3 lines) + // for(int idim=0; idim<3; idim++){ + // //start the recursion with the 111 vector. + // RealType phi = pos[idim] * G111[idim]; + // int ng(maxg[idim]); + // RealType* restrict cp_ptr=logC[idim]+ng; + // RealType* restrict cn_ptr=logC[idim]+ng-1; + // *cp_ptr=0.0; + // //add INTEL vectorization + // for(int n=1; n<=ng; n++,cn_ptr--){ + // RealType t(phi+*cp_ptr++); + // *cp_ptr = t; + // *cn_ptr = -t; + // } + // } + // } +}; +} // namespace qmcplusplus +#endif diff --git a/src/QMCWaveFunctions/PlaneWave/PWOrbitalSetT.cpp b/src/QMCWaveFunctions/PlaneWave/PWOrbitalSetT.cpp new file mode 100644 index 0000000000..a3b1e135ec --- /dev/null +++ b/src/QMCWaveFunctions/PlaneWave/PWOrbitalSetT.cpp @@ -0,0 +1,145 @@ +////////////////////////////////////////////////////////////////////////////////////// +// 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 +////////////////////////////////////////////////////////////////////////////////////// + + +#include "Message/Communicate.h" +#include "PWOrbitalSetT.h" +#include "Numerics/MatrixOperators.h" + +namespace qmcplusplus +{ +template +PWOrbitalSetT::~PWOrbitalSetT() +{ + if (OwnBasisSet && myBasisSet) + delete myBasisSet; + if (!IsCloned && this->C != nullptr) + delete this->C; +} + +template +std::unique_ptr> PWOrbitalSetT::makeClone() const +{ + auto myclone = std::make_unique>(*this); + myclone->myBasisSet = new PWBasisT(*myBasisSet); + myclone->IsCloned = true; + return myclone; +} + +template +void PWOrbitalSetT::setOrbitalSetSize(int norbs) {} + +template +void PWOrbitalSetT::resize(PWBasisPtr bset, int nbands, bool cleanup) +{ + myBasisSet = bset; + this->OrbitalSetSize = nbands; + OwnBasisSet = cleanup; + BasisSetSize = myBasisSet->NumPlaneWaves; + this->C = new ValueMatrix(this->OrbitalSetSize, BasisSetSize); + this->Temp.resize(this->OrbitalSetSize, PW_MAXINDEX); + app_log() << " PWOrbitalSetT::resize OrbitalSetSize =" << this->OrbitalSetSize << " BasisSetSize = " << BasisSetSize + << std::endl; +} + +template +void PWOrbitalSetT::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) + (*(this->C))[jorb][inputmap[ig]] = coefs[ig]; + } +} + +template +void PWOrbitalSetT::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) + (*(this->C))[jorb][inputmap[ig]] = coefs[ig]; + } +} + +template +void PWOrbitalSetT::evaluateValue(const ParticleSet& P, int iat, ValueVector& psi) +{ + //Evaluate every orbital for particle iat. + //Evaluate the basis-set at these coordinates: + //myBasisSet->evaluate(P,iat); + myBasisSet->evaluate(P.activeR(iat)); + MatrixOperators::product(*(this->C), myBasisSet->Zv, psi); +} + +template +void PWOrbitalSetT::evaluateVGL(const ParticleSet& P, int iat, ValueVector& psi, GradVector& dpsi, ValueVector& d2psi) +{ + //Evaluate the orbitals and derivatives for particle iat only. + myBasisSet->evaluateAll(P, iat); + MatrixOperators::product(*(this->C), myBasisSet->Z, this->Temp); + const T* restrict tptr = this->Temp.data(); + for (int j = 0; j < this->OrbitalSetSize; j++, tptr += PW_MAXINDEX) + { + psi[j] = tptr[PW_VALUE]; + d2psi[j] = tptr[PW_LAP]; + dpsi[j] = GradType(tptr[PW_GRADX], tptr[PW_GRADY], tptr[PW_GRADZ]); + } +} + +template +void PWOrbitalSetT::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(*(this->C), myBasisSet->Z, this->Temp); + const T* restrict tptr = this->Temp.data(); + for (int j = 0; j < this->OrbitalSetSize; j++, tptr += PW_MAXINDEX) + { + logdet(i, j) = tptr[PW_VALUE]; + d2logdet(i, j) = tptr[PW_LAP]; + dlogdet(i, j) = GradType(tptr[PW_GRADX], tptr[PW_GRADY], tptr[PW_GRADZ]); + } + } +} + +// Class concrete types from T +// NOTE: This class only gets compiled if QMC_COMPLEX is defined, thus it is inherently complex +// template class PWOrbitalSetT; +// template class PWOrbitalSetT; +template class PWOrbitalSetT>; +template class PWOrbitalSetT>; +} // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/PlaneWave/PWOrbitalSetT.h b/src/QMCWaveFunctions/PlaneWave/PWOrbitalSetT.h new file mode 100644 index 0000000000..39d67f70b1 --- /dev/null +++ b/src/QMCWaveFunctions/PlaneWave/PWOrbitalSetT.h @@ -0,0 +1,130 @@ +////////////////////////////////////////////////////////////////////////////////////// +// 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 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 PWOrbitalSetT.h + * @brief Definition of member functions of Plane-wave basis set + */ +#ifndef QMCPLUSPLUS_PLANEWAVE_ORBITALSET_BLAS_H +#define QMCPLUSPLUS_PLANEWAVE_ORBITALSET_BLAS_H + +#include "QMCWaveFunctions/PlaneWave/PWBasisT.h" +#include "type_traits/complex_help.hpp" +#include "QMCWaveFunctions/SPOSetT.h" +#include "CPU/BLAS.hpp" + +namespace qmcplusplus +{ + +template +class PWOrbitalSetT : public SPOSetT +{ + +public: + using RealType = typename RealAlias_impl::value_type; + using ComplexType = T; + using PosType = QMCTraits::PosType; + using ValueVector = typename SPOSetT::ValueVector; + using GradVector = typename SPOSetT::GradVector; + using ValueMatrix = typename SPOSetT::ValueMatrix; + using GradMatrix = typename SPOSetT::GradMatrix; + using GradType = QMCTraits::GradType; + using IndexType = QMCTraits::IndexType; + + using BasisSet_t = PWBasisT; + using PWBasisPtr = PWBasisT*; + + /** 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 + */ + PWOrbitalSetT(const std::string& my_name) + : SPOSetT(my_name), OwnBasisSet(false), myBasisSet(nullptr), BasisSetSize(0), C(nullptr), IsCloned(false) + {} + + std::string getClassName() const override { return "PWOrbitalSetT"; } + + + /** delete BasisSet only it owns this + * + * Builder takes care of who owns what + */ + ~PWOrbitalSetT() 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); + + /** Builder class takes care of the assertion + */ + void addVector(const std::vector& coefs, int jorb); + 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 BLAS::dot(BasisSetSize, (*C)[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; + + /** boolean + * + * If true, this has to delete the BasisSet + */ + bool OwnBasisSet; + ///TwistAngle of this PWOrbitalSetT + PosType TwistAngle; + ///My basis set + PWBasisPtr myBasisSet; + ///number of basis + IndexType BasisSetSize; + /** pointer to matrix containing the coefficients + * + * makeClone makes a shallow copy and flag IsCloned + */ + ValueMatrix* C; + ///if true, do not clean up + bool IsCloned; + + /** temporary array to perform gemm operation */ + Matrix Temp; +}; +} // namespace qmcplusplus +#endif