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