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: Yubo Yang, paul.young.0414@gmail.com, University of Illinois at Urbana-Champaign
8 //
9 // File created by: Yubo Yang, paul.young.0414@gmail.com, University of Illinois at Urbana-Champaign
10 //////////////////////////////////////////////////////////////////////////////////////
11
12 #include "SpeciesKineticEnergy.h"
13 #include "QMCHamiltonians/BareKineticEnergy.h" // laplaician() defined here
14 #include "OhmmsData/AttributeSet.h"
15
16 namespace qmcplusplus
17 {
SpeciesKineticEnergy(ParticleSet & P)18 SpeciesKineticEnergy::SpeciesKineticEnergy(ParticleSet& P) : tpset(P), hdf5_out(false)
19 {
20 SpeciesSet& tspecies(P.getSpeciesSet());
21 int massind = tspecies.getAttribute("mass");
22
23 num_species = tspecies.size();
24 app_log() << "SpeciesKineticEnergy is tracking " << num_species << " species." << std::endl;
25
26 species_kinetic.resize(num_species);
27 vec_minus_over_2m.resize(num_species);
28 species_names.resize(num_species);
29 for (int ispec = 0; ispec < num_species; ispec++)
30 {
31 species_names[ispec] = tspecies.speciesName[ispec];
32 RealType mass = tspecies(massind, ispec);
33 vec_minus_over_2m[ispec] = -1. / (2. * mass);
34 app_log() << " " << species_names[ispec] << " mass = " << mass << std::endl;
35 }
36 };
37
put(xmlNodePtr cur)38 bool SpeciesKineticEnergy::put(xmlNodePtr cur)
39 {
40 // read hdf5="yes" or "no"
41 OhmmsAttributeSet attrib;
42 std::string hdf5_flag = "no";
43 attrib.add(hdf5_flag, "hdf5");
44 attrib.put(cur);
45
46 if (hdf5_flag == "yes")
47 hdf5_out = true;
48 else if (hdf5_flag == "no")
49 hdf5_out = false;
50 else
51 APP_ABORT("SpeciesKineticEnergy::put() - Please choose \"yes\" or \"no\" for hdf5 flag");
52 // end if
53
54 if (hdf5_out)
55 { // add this estimator to stat.h5
56 UpdateMode.set(COLLECTABLE, 1);
57 }
58 return true;
59 }
60
get(std::ostream & os) const61 bool SpeciesKineticEnergy::get(std::ostream& os) const
62 { // class description
63 os << "SpeciesKineticEnergy: " << myName;
64 return true;
65 }
66
addObservables(PropertySetType & plist,BufferType & collectables)67 void SpeciesKineticEnergy::addObservables(PropertySetType& plist, BufferType& collectables)
68 {
69 myIndex = plist.size();
70 for (int ispec = 0; ispec < num_species; ispec++)
71 { // make columns named "$myName_u", "$myName_d" etc.
72 plist.add(myName + "_" + species_names[ispec]);
73 }
74 if (hdf5_out)
75 { // make room in h5 file and save its index
76 h5_index = collectables.size();
77 std::vector<RealType> tmp(num_species);
78 collectables.add(tmp.begin(), tmp.end());
79 }
80 }
81
registerCollectables(std::vector<observable_helper * > & h5desc,hid_t gid) const82 void SpeciesKineticEnergy::registerCollectables(std::vector<observable_helper*>& h5desc, hid_t gid) const
83 {
84 if (hdf5_out)
85 {
86 std::vector<int> ndim(1, num_species);
87 observable_helper* h5o = new observable_helper(myName);
88 h5o->set_dimensions(ndim, h5_index);
89 h5o->open(gid);
90 h5desc.push_back(h5o);
91 }
92 }
93
setObservables(PropertySetType & plist)94 void SpeciesKineticEnergy::setObservables(PropertySetType& plist)
95 { // slots in plist must be allocated by addObservables() first
96 copy(species_kinetic.begin(), species_kinetic.end(), plist.begin() + myIndex);
97 }
98
evaluate(ParticleSet & P)99 SpeciesKineticEnergy::Return_t SpeciesKineticEnergy::evaluate(ParticleSet& P)
100 {
101 std::fill(species_kinetic.begin(), species_kinetic.end(), 0.0);
102 RealType wgt = tWalker->Weight;
103
104 for (int iat = 0; iat < P.getTotalNum(); iat++)
105 {
106 int ispec = P.GroupID[iat];
107 RealType my_kinetic = vec_minus_over_2m[ispec] * laplacian(P.G[iat], P.L[iat]);
108 if (hdf5_out)
109 {
110 P.Collectables[h5_index + ispec] += my_kinetic * wgt;
111 }
112 species_kinetic[ispec] += my_kinetic;
113 }
114
115 Value = 0.0; // Value is no longer used
116 return Value;
117 }
118
makeClone(ParticleSet & qp,TrialWaveFunction & psi)119 OperatorBase* SpeciesKineticEnergy::makeClone(ParticleSet& qp, TrialWaveFunction& psi)
120 {
121 return new SpeciesKineticEnergy(*this);
122 }
123
124 } // namespace qmcplusplus
125