1 ////////////////////////////////////////////////////////////////////////////////////// 2 // This file is distributed under the University of Illinois/NCSA Open Source License. 3 // See LICENSE file in top directory for details. 4 // 5 // Copyright (c) 2016 Jeongnim Kim and QMCPACK developers. 6 // 7 // File developed by: Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign 8 // Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign 9 // Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign 10 // Jaron T. Krogel, krogeljt@ornl.gov, Oak Ridge National Laboratory 11 // Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory 12 // 13 // File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign 14 ////////////////////////////////////////////////////////////////////////////////////// 15 16 17 #ifndef QMCPLUSPLUS_NONLOCAL_ECPOTENTIAL_COMPONENT_H 18 #define QMCPLUSPLUS_NONLOCAL_ECPOTENTIAL_COMPONENT_H 19 #include "QMCHamiltonians/OperatorBase.h" 20 #include <ResourceCollection.h> 21 #include <TrialWaveFunction.h> 22 #include "Numerics/OneDimGridBase.h" 23 #include "Numerics/OneDimGridFunctor.h" 24 #include "Numerics/OneDimLinearSpline.h" 25 #include "Numerics/OneDimCubicSpline.h" 26 #include "NLPPJob.h" 27 28 namespace qmcplusplus 29 { 30 /** Contains a set of radial grid potentials around a center. 31 */ 32 class NonLocalECPComponent : public QMCTraits 33 { 34 private: 35 typedef std::vector<PosType> SpherGridType; 36 typedef OneDimGridBase<RealType> GridType; 37 typedef OneDimCubicSpline<RealType> RadialPotentialType; 38 39 ///Non Local part: angular momentum, potential and grid 40 int lmax; 41 ///the number of non-local channels 42 int nchannel; 43 ///the number of nknot 44 int nknot; 45 ///Maximum cutoff the non-local pseudopotential 46 RealType Rmax; 47 ///Angular momentum map 48 aligned_vector<int> angpp_m; 49 ///Weight of the angular momentum 50 aligned_vector<RealType> wgt_angpp_m; 51 /// Lfactor1[l]=(2*l+1)/(l+1) 52 RealType Lfactor1[8]; 53 /// Lfactor1[l]=(l)/(l+1) 54 RealType Lfactor2[8]; 55 ///Non-Local part of the pseudo-potential 56 std::vector<RadialPotentialType*> nlpp_m; 57 ///fixed Spherical Grid for species 58 SpherGridType sgridxyz_m; 59 ///randomized spherical grid 60 SpherGridType rrotsgrid_m; 61 ///weight of the spherical grid 62 std::vector<RealType> sgridweight_m; 63 ///Working arrays 64 std::vector<ValueType> wvec, Amat, dAmat; 65 66 //Position delta for virtual moves. 67 std::vector<PosType> deltaV; 68 //Array for P_l[cos(theta)]. 69 std::vector<RealType> lpol; 70 //Array for P'_l[cos(theta)] 71 std::vector<RealType> dlpol; 72 //Array for v_l(r). 73 std::vector<ValueType> vrad; 74 //Array for (2l+1)*v'_l(r)/r. 75 std::vector<RealType> dvrad; 76 //$\Psi(...q...)/\Psi(...r...)$ for all quadrature points q. 77 std::vector<ValueType> psiratio; 78 //$\nabla \Psi(...q...)/\Psi(...r...)$ for all quadrature points q. 79 // $\nabla$ is w.r.t. the electron coordinates involved in the quadrature. 80 std::vector<PosType> gradpsiratio; 81 //This stores gradient of v(r): 82 std::vector<PosType> vgrad; 83 //This stores the gradient of the cos(theta) term in force expression. 84 std::vector<PosType> cosgrad; 85 //This stores grad psi/psi - dot(u,grad psi) 86 std::vector<PosType> wfngrad; 87 //This stores potential contribution per knot: 88 std::vector<RealType> knot_pots; 89 90 /// scratch spaces used by evaluateValueAndDerivatives 91 Matrix<ValueType> dratio; 92 std::vector<ValueType> dlogpsi_vp; 93 94 // For Pulay correction to the force 95 std::vector<RealType> WarpNorm; 96 ParticleSet::ParticleGradient_t dG; 97 ParticleSet::ParticleLaplacian_t dL; 98 /// First index is knot, second is electron 99 Matrix<PosType> Gnew; 100 ///The gradient of the wave function w.r.t. the ion position 101 ParticleSet::ParticleGradient_t Gion; 102 103 ///virtual particle set: delayed initialization 104 VirtualParticleSet* VP; 105 106 /// build QP position deltas from the reference electron using internally stored random grid points 107 void buildQuadraturePointDeltaPositions(RealType r, const PosType& dr, std::vector<PosType>& deltaV) const; 108 109 /** finalize the calculation of $\frac{V\Psi_T}{\Psi_T}$ 110 */ 111 RealType calculateProjector(RealType r, const PosType& dr); 112 113 public: 114 NonLocalECPComponent(); 115 116 ///destructor 117 ~NonLocalECPComponent(); 118 119 NonLocalECPComponent* makeClone(const ParticleSet& qp); 120 121 ///add a new Non Local component 122 void add(int l, RadialPotentialType* pp); 123 124 ///add knots to the spherical grid addknot(const PosType & xyz,RealType weight)125 void addknot(const PosType& xyz, RealType weight) 126 { 127 sgridxyz_m.push_back(xyz); 128 sgridweight_m.push_back(weight); 129 } 130 131 void resize_warrays(int n, int m, int l); 132 133 void randomize_grid(RandomGenerator_t& myRNG); 134 template<typename T> 135 void randomize_grid(std::vector<T>& sphere, RandomGenerator_t& myRNG); 136 137 /** contribute local non-local move data 138 * @param iel reference electron id. 139 * @param Txy nonlocal move data. 140 */ 141 void contributeTxy(int iel, std::vector<NonLocalData>& Txy) const; 142 143 /** @brief Evaluate the nonlocal pp contribution via randomized quadrature grid 144 * to total energy from ion "iat" and electron "iel". 145 * 146 * @param W electron particle set. 147 * @param iat index of ion. 148 * @param Psi trial wave function object 149 * @param iel index of electron 150 * @param r the distance between ion iat and electron iel. 151 * @param dr displacement from ion iat to electron iel. 152 * @param use_DLA if ture, use determinant localization approximation (DLA). 153 * 154 * @return RealType Contribution to $\frac{V\Psi_T}{\Psi_T}$ from ion iat and electron iel. 155 */ 156 RealType evaluateOne(ParticleSet& W, 157 int iat, 158 TrialWaveFunction& Psi, 159 int iel, 160 RealType r, 161 const PosType& dr, 162 bool use_DLA); 163 164 /** @brief Evaluate the nonlocal pp contribution via randomized quadrature grid 165 * to total energy from ion "iat" and electron "iel" for a batch of walkers. 166 * 167 * @param ecp_component_list a list of ECP components 168 * @param p_list a list of electron particle set. 169 * @param psi_list a list of trial wave function object 170 * @param joblist a list of ion-electron pairs 171 * @param pairpots a list of contribution to $\frac{V\Psi_T}{\Psi_T}$ from ion iat and electron iel. 172 * @param use_DLA if ture, use determinant localization approximation (DLA). 173 * 174 * Note: ecp_component_list allows including different NLPP component for different walkers. 175 * electrons in iel_list must be of the same group (spin) 176 */ 177 static void mw_evaluateOne(const RefVectorWithLeader<NonLocalECPComponent>& ecp_component_list, 178 const RefVectorWithLeader<ParticleSet>& p_list, 179 const RefVectorWithLeader<TrialWaveFunction>& psi_list, 180 const RefVector<const NLPPJob<RealType>>& joblist, 181 std::vector<RealType>& pairpots, 182 ResourceCollection& collection, 183 bool use_DLA); 184 185 /** @brief Evaluate the nonlocal pp contribution via randomized quadrature grid 186 * to total energy from ion "iat" and electron "iel". 187 * 188 * @param W electron particle set. 189 * @param iat index of ion. 190 * @param Psi trial wave function object 191 * @param iel index of electron 192 * @param r the distance between ion iat and electron iel. 193 * @param dr displacement from ion iat to electron iel. 194 * @param force_iat 3d vector for Hellman-Feynman contribution. This gets modified. 195 * 196 * @return RealType Contribution to $\frac{V\Psi_T}{\Psi_T}$ from ion iat and electron iel. 197 */ 198 RealType evaluateOneWithForces(ParticleSet& W, 199 int iat, 200 TrialWaveFunction& Psi, 201 int iel, 202 RealType r, 203 const PosType& dr, 204 PosType& force_iat); 205 206 /** @brief Evaluate the nonlocal pp energy, Hellman-Feynman force, and "Pulay" force contribution 207 * via randomized quadrature grid from ion "iat" and electron "iel". 208 * 209 * @param W electron particle set. 210 * @param ions ion particle set. 211 * @param iat index of ion. 212 * @param Psi trial wave function object 213 * @param iel index of electron 214 * @param r the distance between ion iat and electron iel. 215 * @param dr displacement from ion iat to electron iel. 216 * @param force_iat 3d vector for Hellman-Feynman contribution. This gets modified. 217 * @param pulay_terms Nion x 3 object, holding a contribution for each ionic gradient from \Psi_T. 218 * 219 * @return RealType Contribution to $\frac{V\Psi_T}{\Psi_T}$ from ion iat and electron iel. 220 */ 221 RealType evaluateOneWithForces(ParticleSet& W, 222 ParticleSet& ions, 223 int iat, 224 TrialWaveFunction& Psi, 225 int iel, 226 RealType r, 227 const PosType& dr, 228 PosType& force_iat, 229 ParticleSet::ParticlePos_t& pulay_terms); 230 231 // This function needs to be updated to SoA. myTableIndex is introduced temporarily. 232 RealType evaluateValueAndDerivatives(ParticleSet& P, 233 int iat, 234 TrialWaveFunction& psi, 235 int iel, 236 RealType r, 237 const PosType& dr, 238 const opt_variables_type& optvars, 239 const std::vector<ValueType>& dlogpsi, 240 std::vector<ValueType>& dhpsioverpsi); 241 242 void print(std::ostream& os); 243 244 void initVirtualParticle(const ParticleSet& qp); 245 setRmax(int rmax)246 inline void setRmax(int rmax) { Rmax = rmax; } getRmax()247 inline RealType getRmax() const { return Rmax; } getNknot()248 inline int getNknot() const { return nknot; } setLmax(int Lmax)249 inline void setLmax(int Lmax) { lmax = Lmax; } getLmax()250 inline int getLmax() const { return lmax; } getVP()251 const VirtualParticleSet* getVP() const { return VP; }; 252 253 // copy sgridxyz_m to rrotsgrid_m without rotation. For testing only. 254 friend void copyGridUnrotatedForTest(NonLocalECPComponent& nlpp); 255 256 friend struct ECPComponentBuilder; 257 // a lazy temporal solution 258 friend class NonLocalECPotential_CUDA; 259 }; //end of RadialPotentialSet 260 261 } // namespace qmcplusplus 262 #endif 263