QMCPACK
test_J1OrbitalSoA.cpp
Go to the documentation of this file.
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) 2021 QMCPACK developers.
6 //
7 // File developed by: Shiv Upadhyay, shivnupadhyay@gmail.com, University of Pittsburgh
8 //
9 // File created by: Shiv Upadhyay, shivnupadhyay@gmail.com, University of Pittsburgh
10 //////////////////////////////////////////////////////////////////////////////////////
11 
12 #include "catch.hpp"
13 
14 #include "Particle/ParticleSetPool.h"
15 #include "QMCWaveFunctions/Jastrow/J1OrbitalSoA.h"
16 #include "QMCWaveFunctions/Jastrow/RadialJastrowBuilder.h"
17 #include "QMCWaveFunctions/WaveFunctionFactory.h"
18 #include "Utilities/RuntimeOptions.h"
19 
20 namespace qmcplusplus
21 {
22 TEST_CASE("J1 evaluate derivatives Jastrow", "[wavefunction]")
23 {
24  Communicate* c = OHMMS::Controller;
25 
26  ParticleSetPool ptcl = ParticleSetPool(c);
27  auto ions_uptr = std::make_unique<ParticleSet>(ptcl.getSimulationCell());
28  auto elec_uptr = std::make_unique<ParticleSet>(ptcl.getSimulationCell());
29  ParticleSet& ions_(*ions_uptr);
30  ParticleSet& elec_(*elec_uptr);
31 
32  ions_.setName("ion0");
33  ptcl.addParticleSet(std::move(ions_uptr));
34  ions_.create({1});
35  ions_.R[0] = {0.0, 0.0, 0.0};
36  SpeciesSet& ispecies = ions_.getSpeciesSet();
37  int HIdx = ispecies.addSpecies("H");
38  int ichargeIdx = ispecies.addAttribute("charge");
39  ispecies(ichargeIdx, HIdx) = 1.0;
40 
41  elec_.setName("e");
42  ptcl.addParticleSet(std::move(elec_uptr));
43  elec_.create({1, 1});
44  elec_.R[0] = {0.5, 0.5, 0.5};
45  elec_.R[1] = {-0.5, -0.5, -0.5};
46 
47  SpeciesSet& tspecies = elec_.getSpeciesSet();
48  int upIdx = tspecies.addSpecies("u");
49  int downIdx = tspecies.addSpecies("d");
50  int massIdx = tspecies.addAttribute("mass");
51  int chargeIdx = tspecies.addAttribute("charge");
52  tspecies(massIdx, upIdx) = 1.0;
53  tspecies(massIdx, downIdx) = 1.0;
54  tspecies(chargeIdx, upIdx) = -1.0;
55  tspecies(massIdx, downIdx) = -1.0;
56  // Necessary to set mass
57  elec_.resetGroups();
58 
59  ions_.update();
60  elec_.addTable(elec_);
61  elec_.addTable(ions_);
62  elec_.update();
63 
64  ions_.get(app_log());
65  elec_.get(app_log());
66 
67  const char* jasxml = R"(<wavefunction name="psi0" target="e">
68  <jastrow name="J1" type="One-Body" function="Bspline" print="yes" source="ion0">
69  <correlation elementType="H" cusp="0.0" size="2" rcut="5.0">
70  <coefficients id="J1H" type="Array"> 0.5 0.1 </coefficients>
71  </correlation>
72  </jastrow>
73 </wavefunction>
74 )";
75  Libxml2Document doc;
76  bool okay = doc.parseFromString(jasxml);
77  REQUIRE(okay);
78 
79  xmlNodePtr jas1 = doc.getRoot();
80 
81  // update all distance tables
82  elec_.update();
83  RuntimeOptions runtime_options;
84  WaveFunctionFactory wf_factory(elec_, ptcl.getPool(), c);
85  auto twf_ptr = wf_factory.buildTWF(jas1, runtime_options);
86  auto& twf(*twf_ptr);
87  twf.setMassTerm(elec_);
88  twf.evaluateLog(elec_);
89  twf.prepareGroup(elec_, 0);
90 
91  auto& twf_component_list = twf.getOrbitals();
92 
93  opt_variables_type active;
94  twf.checkInVariables(active);
95  active.removeInactive();
96  int nparam = active.size_of_active();
97  REQUIRE(nparam == 2);
98 
99  using ValueType = QMCTraits::ValueType;
100  Vector<ValueType> dlogpsi(nparam);
101  Vector<ValueType> dhpsioverpsi(nparam);
102  //twf.evaluateDerivatives(elec_, active, dlogpsi, dhpsioverpsi);
103  twf_component_list[0]->evaluateDerivatives(elec_, active, dlogpsi, dhpsioverpsi);
104 
105  // Numbers not validated
106  std::vector<ValueType> expected_dlogpsi = {-0.9336294487, -1.0196051794};
107  std::vector<ValueType> expected_dhpsioverpsi = {-1.1596433096, 0.7595492539};
108  for (int i = 0; i < nparam; i++)
109  {
110  CHECK(dlogpsi[i] == ValueApprox(expected_dlogpsi[i]));
111  CHECK(dhpsioverpsi[i] == ValueApprox(expected_dhpsioverpsi[i]));
112  }
113 }
114 
115 TEST_CASE("J1 evaluate derivatives Jastrow with two species", "[wavefunction]")
116 {
117  Communicate* c = OHMMS::Controller;
118  ParticleSetPool ptcl = ParticleSetPool(c);
119  auto ions_uptr = std::make_unique<ParticleSet>(ptcl.getSimulationCell());
120  auto elec_uptr = std::make_unique<ParticleSet>(ptcl.getSimulationCell());
121  ParticleSet& ions_(*ions_uptr);
122  ParticleSet& elec_(*elec_uptr);
123 
124  ions_.setName("ion0");
125  ptcl.addParticleSet(std::move(ions_uptr));
126  ions_.create({1, 1});
127  ions_.R[0] = {0.0, 0.0, 1.0};
128  ions_.R[1] = {0.0, 0.0, 0.0};
129  SpeciesSet& ispecies = ions_.getSpeciesSet();
130  int OIdx = ispecies.addSpecies("O");
131  int HIdx = ispecies.addSpecies("H");
132  int imassIdx = ispecies.addAttribute("mass");
133  int ichargeIdx = ispecies.addAttribute("charge");
134  ispecies(ichargeIdx, HIdx) = 1.0;
135  ispecies(ichargeIdx, OIdx) = 8.0;
136  ispecies(imassIdx, HIdx) = 1.0;
137  ispecies(imassIdx, OIdx) = 16.0;
138 
139  elec_.setName("e");
140  ptcl.addParticleSet(std::move(elec_uptr));
141  elec_.create({1, 1});
142  elec_.R[0] = {0.5, 0.5, 0.5};
143  elec_.R[1] = {-0.5, -0.5, -0.5};
144 
145  SpeciesSet& tspecies = elec_.getSpeciesSet();
146  int upIdx = tspecies.addSpecies("u");
147  int downIdx = tspecies.addSpecies("d");
148  int massIdx = tspecies.addAttribute("mass");
149  int chargeIdx = tspecies.addAttribute("charge");
150  tspecies(massIdx, upIdx) = 1.0;
151  tspecies(massIdx, downIdx) = 1.0;
152  tspecies(chargeIdx, upIdx) = -1.0;
153  tspecies(massIdx, downIdx) = -1.0;
154  // Necessary to set mass
155  elec_.resetGroups();
156 
157  ions_.update();
158  elec_.addTable(elec_);
159  elec_.addTable(ions_);
160  elec_.update();
161 
162  ions_.get(app_log());
163  elec_.get(app_log());
164 
165  const char* jasxml = R"(<wavefunction name="psi0" target="e">
166  <jastrow name="J1" type="One-Body" function="Bspline" print="yes" source="ion0">
167  <correlation elementType="H" cusp="0.0" size="2" rcut="5.0">
168  <coefficients id="J1H" type="Array"> 0.5 0.1 </coefficients>
169  </correlation>
170  <correlation elementType="O" cusp="0.0" size="2" rcut="5.0">
171  <coefficients id="J1O" type="Array"> 0.2 0.1 </coefficients>
172  </correlation>
173  </jastrow>
174 </wavefunction>
175 )";
176  Libxml2Document doc;
177  bool okay = doc.parseFromString(jasxml);
178  REQUIRE(okay);
179 
180  xmlNodePtr jas1 = doc.getRoot();
181 
182  // update all distance tables
183  elec_.update();
184  RuntimeOptions runtime_options;
185  WaveFunctionFactory wf_factory(elec_, ptcl.getPool(), c);
186  auto twf_ptr = wf_factory.buildTWF(jas1, runtime_options);
187  auto& twf(*twf_ptr);
188  twf.setMassTerm(elec_);
189  twf.evaluateLog(elec_);
190  twf.prepareGroup(elec_, 0);
191 
192  auto& twf_component_list = twf.getOrbitals();
193 
194  opt_variables_type active;
195  twf.checkInVariables(active);
196  active.removeInactive();
197  int nparam = active.size_of_active();
198  REQUIRE(nparam == 4);
199 
200  using ValueType = QMCTraits::ValueType;
201  Vector<ValueType> dlogpsi(nparam);
202  Vector<ValueType> dhpsioverpsi(nparam);
203  //twf.evaluateDerivatives(elec_, active, dlogpsi, dhpsioverpsi);
204  twf_component_list[0]->evaluateDerivatives(elec_, active, dlogpsi, dhpsioverpsi);
205 
206  // Numbers not validated independently
207  std::vector<ValueType> expected_dlogpsi = {-0.6360001724, -1.1764442146, -0.9336294487, -1.0196051794};
208  std::vector<ValueType> expected_dhpsioverpsi = {-0.6225838942, 0.0099980417, -1.1853702074, 0.7798000176};
209  for (int i = 0; i < nparam; i++)
210  {
211  CHECK(dlogpsi[i] == ValueApprox(expected_dlogpsi[i]));
212  CHECK(dhpsioverpsi[i] == ValueApprox(expected_dhpsioverpsi[i]));
213  }
214 }
215 
216 TEST_CASE("J1 evaluate derivatives Jastrow with two species one without Jastrow", "[wavefunction]")
217 {
218  Communicate* c = OHMMS::Controller;
219  ParticleSetPool ptcl = ParticleSetPool(c);
220  auto ions_uptr = std::make_unique<ParticleSet>(ptcl.getSimulationCell());
221  auto elec_uptr = std::make_unique<ParticleSet>(ptcl.getSimulationCell());
222  ParticleSet& ions_(*ions_uptr);
223  ParticleSet& elec_(*elec_uptr);
224 
225  ions_.setName("ion0");
226  ptcl.addParticleSet(std::move(ions_uptr));
227  ions_.create({1, 1});
228  ions_.R[0] = {0.0, 0.0, 1.0};
229  ions_.R[1] = {0.0, 0.0, 0.0};
230  SpeciesSet& ispecies = ions_.getSpeciesSet();
231  int OIdx = ispecies.addSpecies("O");
232  int HIdx = ispecies.addSpecies("H");
233  int imassIdx = ispecies.addAttribute("mass");
234  int ichargeIdx = ispecies.addAttribute("charge");
235  ispecies(ichargeIdx, HIdx) = 1.0;
236  ispecies(ichargeIdx, OIdx) = 8.0;
237  ispecies(imassIdx, HIdx) = 1.0;
238  ispecies(imassIdx, OIdx) = 16.0;
239 
240  elec_.setName("e");
241  ptcl.addParticleSet(std::move(elec_uptr));
242  elec_.create({1, 1});
243  elec_.R[0] = {0.5, 0.5, 0.5};
244  elec_.R[1] = {-0.5, -0.5, -0.5};
245 
246  SpeciesSet& tspecies = elec_.getSpeciesSet();
247  int upIdx = tspecies.addSpecies("u");
248  int downIdx = tspecies.addSpecies("d");
249  int massIdx = tspecies.addAttribute("mass");
250  int chargeIdx = tspecies.addAttribute("charge");
251  tspecies(massIdx, upIdx) = 1.0;
252  tspecies(massIdx, downIdx) = 1.0;
253  tspecies(chargeIdx, upIdx) = -1.0;
254  tspecies(massIdx, downIdx) = -1.0;
255  // Necessary to set mass
256  elec_.resetGroups();
257 
258  ions_.update();
259  elec_.addTable(elec_);
260  elec_.addTable(ions_);
261  elec_.update();
262 
263  ions_.get(app_log());
264  elec_.get(app_log());
265 
266  const char* jasxml = R"(<wavefunction name="psi0" target="e">
267  <jastrow name="J1" type="One-Body" function="Bspline" print="yes" source="ion0">
268  <correlation elementType="H" cusp="0.0" size="2" rcut="5.0">
269  <coefficients id="J1H" type="Array"> 0.5 0.1 </coefficients>
270  </correlation>
271  </jastrow>
272 </wavefunction>
273 )";
274  Libxml2Document doc;
275  bool okay = doc.parseFromString(jasxml);
276  REQUIRE(okay);
277 
278  xmlNodePtr jas1 = doc.getRoot();
279 
280  // update all distance tables
281  elec_.update();
282  RuntimeOptions runtime_options;
283  WaveFunctionFactory wf_factory(elec_, ptcl.getPool(), c);
284  auto twf_ptr = wf_factory.buildTWF(jas1, runtime_options);
285  auto& twf(*twf_ptr);
286  twf.setMassTerm(elec_);
287  twf.evaluateLog(elec_);
288  twf.prepareGroup(elec_, 0);
289 
290  auto& twf_component_list = twf.getOrbitals();
291 
292  opt_variables_type active;
293  twf.checkInVariables(active);
294  active.removeInactive();
295  int nparam = active.size_of_active();
296  REQUIRE(nparam == 2);
297 
298  using ValueType = QMCTraits::ValueType;
299  Vector<ValueType> dlogpsi(nparam);
300  Vector<ValueType> dhpsioverpsi(nparam);
301  //twf.evaluateDerivatives(elec_, active, dlogpsi, dhpsioverpsi);
302  twf_component_list[0]->evaluateDerivatives(elec_, active, dlogpsi, dhpsioverpsi);
303 
304  // Numbers not validated independently
305  std::vector<ValueType> expected_dlogpsi = {-0.9336294487, -1.0196051794};
306  std::vector<ValueType> expected_dhpsioverpsi = {-1.1596433096, 0.7595492539};
307  for (int i = 0; i < nparam; i++)
308  {
309  CHECK(dlogpsi[i] == ValueApprox(expected_dlogpsi[i]));
310  CHECK(dhpsioverpsi[i] == ValueApprox(expected_dhpsioverpsi[i]));
311  }
312 }
313 } // namespace qmcplusplus
qmcplusplus::ParticleSet::setName
void setName(const std::string &aname)
Definition: ParticleSet.h:237
Libxml2Document
class that handles xmlDoc
Definition: Libxml2Doc.h:76
qmcplusplus::SpeciesSet::addSpecies
int addSpecies(const std::string &aname)
When a name species does not exist, add a new species.
Definition: SpeciesSet.cpp:33
qmcplusplus
helper functions for EinsplineSetBuilder
Definition: Configuration.h:43
qmcplusplus::app_log
std::ostream & app_log()
Definition: OutputManager.h:65
qmcplusplus::TEST_CASE
TEST_CASE("complex_helper", "[type_traits]")
Definition: test_complex_helper.cpp:38
qmcplusplus::ParticleSetPool::addParticleSet
void addParticleSet(std::unique_ptr< ParticleSet > &&p)
add a ParticleSet* to the pool with its ownership transferred ParticleSet built outside the ParticleS...
Definition: ParticleSetPool.cpp:67
Libxml2Document::getRoot
xmlNodePtr getRoot()
Definition: Libxml2Doc.h:88
qmcplusplus::ParticleSet::update
void update(bool skipSK=false)
update the internal data
Definition: ParticleSet.cpp:349
OHMMS::Controller
Communicate * Controller
Global Communicator for a process.
Definition: Communicate.cpp:35
qmcplusplus::SpeciesSet::addAttribute
int addAttribute(const std::string &aname)
for a new attribute, allocate the data, !More often used to get the index of a species ...
Definition: SpeciesSet.cpp:45
Communicate
Wrapping information on parallelism.
Definition: Communicate.h:68
qmcplusplus::ParticleSet::addTable
int addTable(const ParticleSet &psrc, DTModes modes=DTModes::ALL_OFF)
add a distance table
Definition: ParticleSet.cpp:307
qmcplusplus::ParticleSet
Specialized paritlce class for atomistic simulations.
Definition: ParticleSet.h:55
qmcplusplus::QMCTraits::ValueType
QTBase::ValueType ValueType
Definition: Configuration.h:60
qmcplusplus::WaveFunctionFactory::buildTWF
std::unique_ptr< TrialWaveFunction > buildTWF(xmlNodePtr cur, const RuntimeOptions &runtime_options)
read from xmlNode
Definition: WaveFunctionFactory.cpp:44
qmcplusplus::REQUIRE
REQUIRE(std::filesystem::exists(filename))
qmcplusplus::ParticleSetPool
Manage a collection of ParticleSet objects.
Definition: ParticleSetPool.h:34
qmcplusplus::WaveFunctionFactory
Factory class to build a many-body wavefunction.
Definition: WaveFunctionFactory.h:31
optimize::VariableSet
class to handle a set of variables that can be modified during optimizations
Definition: VariableSet.h:49
qmcplusplus::ParticleSet::R
ParticlePos R
Position.
Definition: ParticleSet.h:79
optimize::VariableSet::size_of_active
int size_of_active() const
return the number of active variables
Definition: VariableSet.h:78
qmcplusplus::ParticleSet::getSpeciesSet
SpeciesSet & getSpeciesSet()
retrun the SpeciesSet of this particle set
Definition: ParticleSet.h:231
qmcplusplus::ParticleSet::create
void create(const std::vector< int > &agroup)
create grouped particles
Definition: ParticleSet.cpp:125
Libxml2Document::parseFromString
bool parseFromString(const std::string_view data)
Definition: Libxml2Doc.cpp:204
qmcplusplus::CHECK
CHECK(log_values[0]==ComplexApprox(std::complex< double >{ 5.603777579195571, -6.1586603331188225 }))
qmcplusplus::ValueType
LatticeGaussianProduct::ValueType ValueType
Definition: LatticeGaussianProduct.cpp:20
WaveFunctionFactory.h
Declaration of a WaveFunctionFactory.
qmcplusplus::SpeciesSet
Custom container for set of attributes for a set of species.
Definition: SpeciesSet.h:33
qmcplusplus::ParticleSetPool::getPool
const PoolType & getPool() const
get the Pool object
Definition: ParticleSetPool.h:92
optimize::VariableSet::removeInactive
void removeInactive()
remove inactive variables and trim the internal data
Definition: VariableSet.cpp:135
qmcplusplus::ParticleSet::get
bool get(std::ostream &os) const override
dummy. For satisfying OhmmsElementBase.
Definition: ParticleSet.cpp:300
qmcplusplus::ParticleSetPool::getSimulationCell
const auto & getSimulationCell() const
get simulation cell
Definition: ParticleSetPool.h:95
ParticleSetPool.h
Declaration of ParticleSetPool.