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 "LatticeDeviationEstimator.h"
13 #include "OhmmsData/AttributeSet.h"
14
15 namespace qmcplusplus
16 {
LatticeDeviationEstimator(ParticleSet & P,ParticleSet & sP,const std::string & tgroup_in,const std::string & sgroup_in)17 LatticeDeviationEstimator::LatticeDeviationEstimator(ParticleSet& P,
18 ParticleSet& sP,
19 const std::string& tgroup_in,
20 const std::string& sgroup_in)
21 : tspecies(P.getSpeciesSet()),
22 sspecies(sP.getSpeciesSet()),
23 tpset(P),
24 spset(sP),
25 tgroup(tgroup_in),
26 sgroup(sgroup_in),
27 hdf5_out(false),
28 per_xyz(false),
29 myTableID_(P.addTable(sP))
30 {
31 // calculate number of source particles to use as lattice sites
32 int src_species_id = sspecies.findSpecies(sgroup);
33 num_sites = spset.last(src_species_id) - spset.first(src_species_id);
34 int tar_species_id = tspecies.findSpecies(tgroup);
35 int num_tars = tpset.last(tar_species_id) - tpset.first(tar_species_id);
36 if (num_tars != num_sites)
37 {
38 app_log() << "number of target particles = " << num_tars << std::endl;
39 app_log() << "number of source particles = " << num_sites << std::endl;
40 APP_ABORT("nsource != ntargets in LatticeDeviationEstimator");
41 }
42 }
43
put(xmlNodePtr cur)44 bool LatticeDeviationEstimator::put(xmlNodePtr cur)
45 {
46 input_xml = cur;
47 std::string hdf5_flag = "no";
48 std::string xyz_flag = "no";
49 OhmmsAttributeSet attrib;
50 attrib.add(hdf5_flag, "hdf5");
51 attrib.add(xyz_flag, "per_xyz");
52 attrib.put(cur);
53
54 if (hdf5_flag == "yes")
55 {
56 hdf5_out = true;
57 }
58 else if (hdf5_flag == "no")
59 {
60 hdf5_out = false;
61 }
62 else
63 {
64 APP_ABORT("LatticeDeviationEstimator::put() - Please choose \"yes\" or \"no\" for hdf5 flag");
65 } // end if hdf5_flag
66 if (hdf5_out)
67 {
68 // change the default behavior of addValue() called by addObservables()
69 // YY: does this still matter if I have overwritten addObservables()?
70 UpdateMode.set(COLLECTABLE, 1);
71 }
72
73 if (xyz_flag == "yes")
74 {
75 per_xyz = true;
76 xyz2.resize(OHMMS_DIM);
77 }
78 else if (xyz_flag == "no")
79 {
80 per_xyz = false;
81 }
82 else
83 {
84 APP_ABORT("LatticeDeviationEstimator::put() - Please choose \"yes\" or \"no\" for per_xyz flag");
85 } // end if xyz_flag
86
87 return true;
88 }
89
get(std::ostream & os) const90 bool LatticeDeviationEstimator::get(std::ostream& os) const
91 { // class description
92 os << "LatticeDeviationEstimator: " << myName << "lattice = " << spset.getName();
93 return true;
94 }
95
evaluate(ParticleSet & P)96 LatticeDeviationEstimator::Return_t LatticeDeviationEstimator::evaluate(ParticleSet& P)
97 { // calculate <r^2> averaged over lattice sites
98 Value = 0.0;
99 std::fill(xyz2.begin(), xyz2.end(), 0.0);
100
101 RealType wgt = tWalker->Weight;
102 const auto& d_table = P.getDistTable(myTableID_);
103
104 // temp variables
105 RealType r, r2;
106 PosType dr;
107
108 int nsite(0); // site index
109 int cur_jat(-1); // target particle index
110 for (int iat = 0; iat < spset.getTotalNum(); iat++)
111 { // for each desired source particle
112 if (sspecies.speciesName[spset.GroupID[iat]] == sgroup)
113 { // find desired species
114 for (int jat = cur_jat + 1; jat < tpset.getTotalNum(); jat++)
115 { // find corresponding (!!!! assume next) target particle
116 if (tspecies.speciesName[tpset.GroupID[jat]] == tgroup)
117 {
118 // distance between particle iat in source pset, and jat in target pset
119 r = d_table.getDistRow(jat)[iat];
120 r2 = r * r;
121 Value += r2;
122
123 if (hdf5_out & !per_xyz)
124 { // store deviration for each lattice site if h5 file is available
125 P.Collectables[h5_index + nsite] = wgt * r2;
126 }
127
128 if (per_xyz)
129 {
130 dr = d_table.getDisplRow(jat)[iat];
131 for (int idir = 0; idir < OHMMS_DIM; idir++)
132 {
133 RealType dir2 = dr[idir] * dr[idir];
134 xyz2[idir] += dir2;
135 if (hdf5_out)
136 {
137 P.Collectables[h5_index + nsite * OHMMS_DIM + idir] = wgt * dir2;
138 }
139 }
140 }
141
142 cur_jat = jat;
143 break;
144 }
145 }
146 nsite += 1; // count the number of sites, for checking only
147 } // found desired species (source particle)
148 }
149
150 if (nsite != num_sites)
151 {
152 app_log() << "num_sites = " << num_sites << " nsites = " << nsite << std::endl;
153 APP_ABORT("Mismatch in LatticeDeivationEstimator.");
154 }
155
156 // average per site
157 Value /= num_sites;
158 if (per_xyz)
159 {
160 for (int idir = 0; idir < OHMMS_DIM; idir++)
161 {
162 xyz2[idir] /= num_sites;
163 }
164 }
165
166 return Value;
167 }
168
addObservables(PropertySetType & plist,BufferType & collectables)169 void LatticeDeviationEstimator::addObservables(PropertySetType& plist, BufferType& collectables)
170 {
171 // get myIndex for scalar.dat
172 if (per_xyz)
173 {
174 myIndex = plist.size();
175 for (int idir = 0; idir < OHMMS_DIM; idir++)
176 {
177 std::stringstream ss;
178 ss << idir;
179 plist.add(myName + "_dir" + ss.str());
180 }
181 }
182 else
183 {
184 myIndex = plist.add(myName); // same as OperatorBase::addObservables
185 }
186
187 // get h5_index for stat.h5
188 if (hdf5_out)
189 {
190 h5_index = collectables.size();
191 std::vector<RealType> tmp;
192 if (per_xyz)
193 {
194 tmp.resize(num_sites * OHMMS_DIM);
195 }
196 else
197 {
198 tmp.resize(num_sites);
199 }
200 collectables.add(tmp.begin(), tmp.end());
201 }
202 }
203
setObservables(PropertySetType & plist)204 void LatticeDeviationEstimator::setObservables(PropertySetType& plist)
205 {
206 if (per_xyz)
207 {
208 copy(xyz2.begin(), xyz2.end(), plist.begin() + myIndex);
209 }
210 else
211 {
212 plist[myIndex] = Value; // default behavior
213 }
214 }
215
resetTargetParticleSet(ParticleSet & P)216 void LatticeDeviationEstimator::resetTargetParticleSet(ParticleSet& P) {}
217
makeClone(ParticleSet & qp,TrialWaveFunction & psi)218 OperatorBase* LatticeDeviationEstimator::makeClone(ParticleSet& qp, TrialWaveFunction& psi)
219 {
220 // default constructor does not work with threads
221 //LatticeDeviationEstimator* myclone = new LatticeDeviationEstimator(*this);
222 LatticeDeviationEstimator* myclone = new LatticeDeviationEstimator(qp, spset, tgroup, sgroup);
223 myclone->put(input_xml);
224
225 return myclone;
226 }
227
registerCollectables(std::vector<observable_helper * > & h5desc,hid_t gid) const228 void LatticeDeviationEstimator::registerCollectables(std::vector<observable_helper*>& h5desc, hid_t gid) const
229 {
230 if (hdf5_out)
231 {
232 // one scalar per lattice site (i.e. source particle)
233 std::vector<int> ndim(1, num_sites);
234 if (per_xyz)
235 {
236 // one scalar per lattice site per dimension
237 ndim[0] = num_sites * OHMMS_DIM;
238 }
239
240 // open hdf5 entry and resize
241 observable_helper* h5o = new observable_helper(myName);
242 h5o->set_dimensions(ndim, h5_index);
243 h5o->open(gid);
244
245 // add to h5 file
246 h5desc.push_back(h5o);
247 }
248 }
249
250
251 } // namespace qmcplusplus
252