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