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 // Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
9 // Jeremy McMinnis, jmcminis@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_BAREKINETICENERGY_H
18 #define QMCPLUSPLUS_BAREKINETICENERGY_H
19 #include "Particle/ParticleSet.h"
20 #include "Particle/WalkerSetRef.h"
21 #include "QMCHamiltonians/OperatorBase.h"
22 #include "ParticleBase/ParticleAttribOps.h"
23 #include "QMCWaveFunctions/TrialWaveFunction.h"
24 #include "QMCDrivers/WalkerProperties.h"
25 #ifdef QMC_CUDA
26 #include "Particle/MCWalkerConfiguration.h"
27 #endif
28 #include "type_traits/scalar_traits.h"
29
30 namespace qmcplusplus
31 {
32 using WP = WalkerProperties::Indexes;
33
34 /** compute real(laplacian)
35 */
36 template<typename T, unsigned D>
laplacian(const TinyVector<T,D> & g,T l)37 inline T laplacian(const TinyVector<T, D>& g, T l)
38 {
39 return dot(g, g) + l;
40 }
41
42 /** specialization of laplacian with complex g & l
43 */
44 template<typename T, unsigned D>
laplacian(const TinyVector<std::complex<T>,D> & g,const std::complex<T> & l)45 inline T laplacian(const TinyVector<std::complex<T>, D>& g, const std::complex<T>& l)
46 {
47 return l.real() + OTCDot<T, T, D>::apply(g, g);
48 }
49
50 /** Convenience function to compute \f$\Re( \nabla^2_i \partial \Psi_T/\Psi_T)\f$
51 * @param g OHMMS_DIM dimensional vector for \f$\nabla_i \ln \Psi_T\f$ .
52 * @param l A number, representing \f$\nabla^2_i \ln \Psi_T\f$ .
53 * @param gg OHMMS_DIM dimensional vector containing \f$\nabla_i \partial \ln \Psi_T\f$ .
54 * @param gl A number, representing \f$\nabla^2_i \partial \ln \Psi_T\f$
55 * @param ideriv A number, representing \f$\partial \ln \Psi_T\f$
56 *
57 * @return A number corresponding to \f$\Re( \nabla^2_i \partial \Psi_T/\Psi_T)\f$
58 */
59
60 template<typename T, unsigned D>
dlaplacian(const TinyVector<T,D> & g,const T l,const TinyVector<T,D> & gg,const T gl,const T ideriv)61 inline T dlaplacian(const TinyVector<T, D>& g, const T l, const TinyVector<T, D>& gg, const T gl, const T ideriv)
62 {
63 return gl + l * ideriv + 2.0 * dot(g, gg) + dot(g, g) * ideriv;
64 }
65
66 template<typename T, unsigned D>
dlaplacian(const TinyVector<std::complex<T>,D> & g,const std::complex<T> l,const TinyVector<std::complex<T>,D> & gg,const std::complex<T> gl,const std::complex<T> ideriv)67 inline T dlaplacian(const TinyVector<std::complex<T>, D>& g,
68 const std::complex<T> l,
69 const TinyVector<std::complex<T>, D>& gg,
70 const std::complex<T> gl,
71 const std::complex<T> ideriv)
72 {
73 std::complex<T> l_times_ideriv = l * ideriv;
74 TinyVector<std::complex<T>, D> g_times_ideriv = g * ideriv;
75
76 return gl.real() + l_times_ideriv.real() + 2.0 * OTCDot<T, T, D>::apply(g, gg) +
77 OTCDot<T, T, D>::apply(g, g_times_ideriv);
78 }
79
80
81 /** @ingroup hamiltonian
82 @brief Evaluate the kinetic energy with a single mass
83
84 *The unit of the mass is AU, i.e., the electron mass \f$ m_e = 1 \f$.
85 * To evaluate the Bare Kinetic part of the local energy
86 \f$E_L({\bf R}) = \Psi^{-1}({\bf R})\hat{H}\Psi({\bf R}),\f$
87 it is useful to use the following trick
88 \f{eqnarray*}
89 \nabla^2\Psi({\bf R}) &=& \nabla^2(\exp(\ln \Psi({\bf R})))\\
90 &=&\nabla\cdot(\nabla\exp(\ln \Psi({\bf R}))) \\
91 &=&\nabla\cdot(\nabla\ln \Psi({\bf R}))\exp(\ln \Psi({\bf R}))\\
92 -\frac{1}{2}\frac{\nabla^2\Psi({\bf R})}{\Psi({\bf R})} &=&
93 -\frac{1}{2}\nabla^2\ln \Psi({\bf R})
94 -\frac{1}{2}(\nabla\ln \Psi({\bf R}))^2
95 \f}
96 */
97
98
int2string(const int & i)99 inline std::string int2string(const int& i)
100 {
101 std::stringstream ss;
102 ss << i;
103 return ss.str();
104 }
105
106
107 template<typename T>
108 struct BareKineticEnergy : public OperatorBase
109 {
110 ///true, if all the species have the same mass
111 bool SameMass;
112 ///mass of the particle
113 T M;
114 ///\f$ 1/(2 m^*) \f$
115 T OneOver2M;
116 ///MinusOver2M[i] = \f$ -1/2m[i]\f$ for the ith species
117 std::vector<T> MinusOver2M;
118
119 ParticleSet::ParticleGradient_t Gtmp;
120 ParticleSet::ParticleLaplacian_t Ltmp;
121
122 ///single particle trace samples
123 bool streaming_kinetic;
124 bool streaming_kinetic_comp;
125 bool streaming_momentum;
126
127 #if !defined(REMOVE_TRACEMANAGER)
128 Array<TraceReal, 1>* T_sample;
129 Array<TraceComp, 1>* T_sample_comp;
130 Array<TraceComp, 2>* p_sample;
131 #endif
132 ParticleSet& Ps;
133
134 /** constructor
135 *
136 * Kinetic operators need to be re-evaluated during optimization.
137 */
SameMassBareKineticEnergy138 BareKineticEnergy(RealType m = 1.0) : SameMass(true), M(m), OneOver2M(0.5 / m)
139 {
140 set_energy_domain(kinetic);
141 set_quantum_domain(quantum);
142 streaming_kinetic = false;
143 streaming_kinetic_comp = false;
144 streaming_momentum = false;
145 UpdateMode.set(OPTIMIZABLE, 1);
146 }
147
148 /** constructor with particleset
149 * @param target particleset
150 *
151 * Store mass per species and use SameMass to choose the methods.
152 * if SameMass, probably faster and easy to vectorize but no impact on the performance.
153 */
BareKineticEnergyBareKineticEnergy154 BareKineticEnergy(ParticleSet& p) : Ps(p)
155 {
156 set_energy_domain(kinetic);
157 one_body_quantum_domain(p);
158 streaming_kinetic = false;
159 streaming_kinetic_comp = false;
160 streaming_momentum = false;
161 UpdateMode.set(OPTIMIZABLE, 1);
162 SpeciesSet& tspecies(p.getSpeciesSet());
163 MinusOver2M.resize(tspecies.size());
164 int massind = tspecies.addAttribute("mass");
165 SameMass = true;
166 M = tspecies(massind, 0);
167 OneOver2M = 0.5 / M;
168 for (int i = 0; i < tspecies.size(); ++i)
169 {
170 SameMass &= (std::abs(tspecies(massind, i) - M) < 1e-6);
171 MinusOver2M[i] = -1.0 / (2.0 * tspecies(massind, i));
172 }
173 }
174 ///destructor
~BareKineticEnergyBareKineticEnergy175 ~BareKineticEnergy() {}
176
resetTargetParticleSetBareKineticEnergy177 void resetTargetParticleSet(ParticleSet& P) {}
178
179
180 #if !defined(REMOVE_TRACEMANAGER)
contribute_particle_quantitiesBareKineticEnergy181 virtual void contribute_particle_quantities()
182 {
183 request.contribute_array(myName);
184 request.contribute_array(myName + "_complex");
185 request.contribute_array("momentum");
186 }
187
checkout_particle_quantitiesBareKineticEnergy188 virtual void checkout_particle_quantities(TraceManager& tm)
189 {
190 streaming_particles = request.streaming_array(myName) || request.streaming_array(myName + "_complex") ||
191 request.streaming_array("momentum");
192 if (streaming_particles)
193 {
194 T_sample = tm.checkout_real<1>(myName, Ps);
195 T_sample_comp = tm.checkout_complex<1>(myName + "_complex", Ps);
196 p_sample = tm.checkout_complex<2>("momentum", Ps, DIM);
197 }
198 }
199
delete_particle_quantitiesBareKineticEnergy200 virtual void delete_particle_quantities()
201 {
202 if (streaming_particles)
203 {
204 delete T_sample;
205 delete T_sample_comp;
206 delete p_sample;
207 }
208 }
209 #endif
210
211
evaluateBareKineticEnergy212 inline Return_t evaluate(ParticleSet& P)
213 {
214 #if !defined(REMOVE_TRACEMANAGER)
215 if (streaming_particles)
216 {
217 Value = evaluate_sp(P);
218 }
219 else
220 #endif
221 if (SameMass)
222 {
223 //app_log() << "Here" << std::endl;
224 #ifdef QMC_COMPLEX
225 Value = std::real( CplxDot(P.G, P.G) + CplxSum(P.L) );
226 Value *= -OneOver2M;
227 #else
228 Value = Dot(P.G, P.G) + Sum(P.L);
229 Value *= -OneOver2M;
230 #endif
231 }
232 else
233 {
234 Value = 0.0;
235 for (int i = 0; i < MinusOver2M.size(); ++i)
236 {
237 T x = 0.0;
238 for (int j = P.first(i); j < P.last(i); ++j)
239 x += laplacian(P.G[j], P.L[j]);
240 Value += x * MinusOver2M[i];
241 }
242 }
243 return Value;
244 }
245
246 /**@brief Function to compute the value, direct ionic gradient terms, and pulay terms for the local kinetic energy.
247 *
248 * This general function represents the OperatorBase interface for computing. For an operator \hat{O}, this
249 * function will return \frac{\hat{O}\Psi_T}{\Psi_T}, \frac{\partial(\hat{O})\Psi_T}{\Psi_T}, and
250 * \frac{\hat{O}\partial\Psi_T}{\Psi_T} - \frac{\hat{O}\Psi_T}{\Psi_T}\frac{\partial \Psi_T}{\Psi_T}. These are
251 * referred to as Value, HF term, and pulay term respectively.
252 *
253 * @param P electron particle set.
254 * @param ions ion particle set
255 * @param psi Trial wave function object.
256 * @param hf_terms 3Nion dimensional object. All direct force terms, or ionic gradient of operator itself.
257 * Contribution of this operator is ADDED onto hf_terms.
258 * @param pulay_terms The terms coming from ionic gradients of trial wavefunction. Contribution of this operator is
259 * ADDED onto pulay_terms.
260 * @return Value of kinetic energy operator at electron/ion positions given by P and ions. The force contributions from
261 * this operator are added into hf_terms and pulay_terms.
262 */
evaluateWithIonDerivsBareKineticEnergy263 inline Return_t evaluateWithIonDerivs(ParticleSet& P,
264 ParticleSet& ions,
265 TrialWaveFunction& psi,
266 ParticleSet::ParticlePos_t& hf_terms,
267 ParticleSet::ParticlePos_t& pulay_terms)
268 {
269 typedef ParticleSet::ParticlePos_t ParticlePos_t;
270 typedef ParticleSet::ParticleGradient_t ParticleGradient_t;
271 typedef ParticleSet::ParticleLaplacian_t ParticleLaplacian_t;
272
273 int Nions = ions.getTotalNum();
274 int Nelec = P.getTotalNum();
275
276 //These are intermediate arrays for potentially complex math.
277 ParticleLaplacian_t term2_(Nelec);
278 ParticleGradient_t term4_(Nelec);
279
280 //Potentially complex temporary array for \partial \psi/\psi and \nabla^2 \partial \psi / \psi
281 ParticleGradient_t iongradpsi_(Nions), pulaytmp_(Nions);
282 //temporary arrays that will be explicitly real.
283 ParticlePos_t pulaytmpreal_(Nions), iongradpsireal_(Nions);
284
285
286 TinyVector<ParticleGradient_t, OHMMS_DIM> iongrad_grad_;
287 TinyVector<ParticleLaplacian_t, OHMMS_DIM> iongrad_lapl_;
288
289 for (int iondim = 0; iondim < OHMMS_DIM; iondim++)
290 {
291 iongrad_grad_[iondim].resize(Nelec);
292 iongrad_lapl_[iondim].resize(Nelec);
293 }
294
295 iongradpsi_ = 0;
296 iongradpsireal_ = 0;
297 pulaytmpreal_ = 0;
298 pulaytmp_ = 0;
299
300 RealType logpsi_ = psi.evaluateLog(P);
301 for (int iat = 0; iat < Nions; iat++)
302 {
303 //reset the iongrad_X containers.
304 for (int iondim = 0; iondim < OHMMS_DIM; iondim++)
305 {
306 iongrad_grad_[iondim] = 0;
307 iongrad_lapl_[iondim] = 0;
308 }
309 iongradpsi_[iat] = psi.evalGradSource(P, ions, iat, iongrad_grad_, iongrad_lapl_);
310 //conversion from potentially complex to definitely real.
311 convert(iongradpsi_[iat], iongradpsireal_[iat]);
312 if (SameMass)
313 {
314 for (int iondim = 0; iondim < OHMMS_DIM; iondim++)
315 {
316 //These term[24]_ variables exist because I want to do complex math first, and then take the real part at the
317 //end. Sum() and Dot() perform the needed functions and spit out the real part at the end.
318 term2_ = P.L * iongradpsi_[iat][iondim];
319 term4_ = P.G * iongradpsi_[iat][iondim];
320 pulaytmp_[iat][iondim] = -OneOver2M *
321 (Sum(iongrad_lapl_[iondim]) + Sum(term2_) + 2.0 * Dot(iongrad_grad_[iondim], P.G) + Dot(P.G, term4_));
322 }
323 }
324 else
325 {
326 for (int iondim = 0; iondim < OHMMS_DIM; iondim++)
327 {
328 for (int g = 0; g < MinusOver2M.size(); g++)
329 {
330 for (int iel = P.first(g); iel < P.last(g); iel++)
331 {
332 pulaytmp_[iat][iondim] += MinusOver2M[g] *
333 (dlaplacian(P.G[iel], P.L[iel], iongrad_grad_[iondim][iel], iongrad_lapl_[iondim][iel],
334 iongradpsi_[iat][iondim]));
335 }
336 }
337 }
338 }
339 //convert to real.
340 convert(pulaytmp_[iat], pulaytmpreal_[iat]);
341 }
342
343 if (SameMass)
344 {
345 Value = Dot(P.G, P.G) + Sum(P.L);
346 Value *= -OneOver2M;
347 }
348 else
349 {
350 Value = 0.0;
351 for (int i = 0; i < MinusOver2M.size(); ++i)
352 {
353 T x = 0.0;
354 for (int j = P.first(i); j < P.last(i); ++j)
355 x += laplacian(P.G[j], P.L[j]);
356 Value += x * MinusOver2M[i];
357 }
358 }
359 pulaytmpreal_ -= Value * iongradpsireal_;
360
361
362 pulay_terms += pulaytmpreal_;
363 return Value;
364 }
365
366
367 #if !defined(REMOVE_TRACEMANAGER)
evaluate_spBareKineticEnergy368 inline Return_t evaluate_sp(ParticleSet& P)
369 {
370 Array<RealType, 1>& T_samp = *T_sample;
371 Array<std::complex<RealType>, 1>& T_samp_comp = *T_sample_comp;
372 Array<std::complex<RealType>, 2>& p_samp = *p_sample;
373 std::complex<RealType> t1 = 0.0;
374 const RealType clambda(-OneOver2M);
375 Value = 0.0;
376 if (SameMass)
377 {
378 for (int i = 0; i < P.getTotalNum(); i++)
379 {
380 t1 = P.L[i] + dot(P.G[i], P.G[i]);
381 t1 = clambda * t1;
382 T_samp(i) = real(t1);
383 T_samp_comp(i) = t1;
384 for (int d = 0; d < DIM; ++d)
385 p_samp(i, d) = P.G[i][d];
386 Value += real(t1);
387 }
388 }
389 else
390 {
391 for (int s = 0; s < MinusOver2M.size(); ++s)
392 {
393 T mlambda = MinusOver2M[s];
394 for (int i = P.first(s); i < P.last(s); ++i)
395 {
396 //t1 = mlambda*( P.L[i] + dot(P.G[i],P.G[i]) );
397 t1 = P.L[i] + dot(P.G[i], P.G[i]);
398 t1 *= mlambda;
399 T_samp(i) = real(t1);
400 T_samp_comp(i) = t1;
401 for (int d = 0; d < DIM; ++d)
402 p_samp(i, d) = P.G[i][d];
403 Value += real(t1);
404 }
405 }
406 }
407 #if defined(TRACE_CHECK)
408 RealType Vnow = Value;
409 RealType Vsum = T_samp.sum();
410 RealType Vold = evaluate_orig(P);
411 if (std::abs(Vsum - Vnow) > TraceManager::trace_tol)
412 {
413 app_log() << "accumtest: BareKineticEnergy::evaluate()" << std::endl;
414 app_log() << "accumtest: tot:" << Vnow << std::endl;
415 app_log() << "accumtest: sum:" << Vsum << std::endl;
416 APP_ABORT("Trace check failed");
417 }
418 if (std::abs(Vold - Vnow) > TraceManager::trace_tol)
419 {
420 app_log() << "versiontest: BareKineticEnergy::evaluate()" << std::endl;
421 app_log() << "versiontest: orig:" << Vold << std::endl;
422 app_log() << "versiontest: mod:" << Vnow << std::endl;
423 APP_ABORT("Trace check failed");
424 }
425 #endif
426 return Value;
427 }
428
429 #endif
430
evaluate_origBareKineticEnergy431 inline Return_t evaluate_orig(ParticleSet& P)
432 {
433 if (SameMass)
434 {
435 Value = Dot(P.G, P.G) + Sum(P.L);
436 Value *= -OneOver2M;
437 }
438 else
439 {
440 Value = 0.0;
441 for (int i = 0; i < MinusOver2M.size(); ++i)
442 {
443 T x = 0.0;
444 for (int j = P.first(i); j < P.last(i); ++j)
445 x += laplacian(P.G[j], P.L[j]);
446 Value += x * MinusOver2M[i];
447 }
448 }
449 return Value;
450 }
451
452 /** implements the virtual function.
453 *
454 * Nothing is done but should check the mass
455 */
putBareKineticEnergy456 bool put(xmlNodePtr) { return true; }
457
getBareKineticEnergy458 bool get(std::ostream& os) const
459 {
460 os << "Kinetic energy";
461 return true;
462 }
463
makeCloneBareKineticEnergy464 OperatorBase* makeClone(ParticleSet& qp, TrialWaveFunction& psi) { return new BareKineticEnergy(*this); }
465
466 #ifdef QMC_CUDA
467 ////////////////////////////////
468 // Vectorized version for GPU //
469 ////////////////////////////////
470 // Nothing is done on GPU here, just copy into vector
addEnergyBareKineticEnergy471 void addEnergy(MCWalkerConfiguration& W, std::vector<RealType>& LocalEnergy)
472 {
473 std::vector<Walker_t*>& walkers = W.WalkerList;
474 for (int iw = 0; iw < walkers.size(); iw++)
475 {
476 Walker_t& w = *(walkers[iw]);
477 RealType KE = -OneOver2M * (Dot(w.G, w.G) + Sum(w.L));
478 w.getPropertyBase()[WP::NUMPROPERTIES + myIndex] = KE;
479 LocalEnergy[iw] += KE;
480 }
481 }
482 #endif
483
484 //Not used anymore
485 //void evaluate(WalkerSetRef& W, ValueVectorType& LE) {
486 // for(int iw=0; iw< W.walkers(); iw++) {
487 // RealType ke = 0.0;
488 // for(int iat=0; iat< W.particles(); iat++) {
489 // ke += dot(W.G(iw,iat),W.G(iw,iat)) + W.L(iw,iat);
490 // }
491 // LE[iw] -= M*ke;
492 // }
493 //}
494 };
495 } // namespace qmcplusplus
496 #endif
497