QMCPACK
test_einset_NiO_a16.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) 2016 Jeongnim Kim and QMCPACK developers.
6 //
7 // File developed by: Mark Dewing, markdewing@gmail.com, University of Illinois at Urbana-Champaign
8 //
9 // File created by: Mark Dewing, markdewing@gmail.com, University of Illinois at Urbana-Champaign
10 //////////////////////////////////////////////////////////////////////////////////////
11 
12 
13 #include "catch.hpp"
14 
15 #include "OhmmsData/Libxml2Doc.h"
16 #include "OhmmsPETE/OhmmsMatrix.h"
17 #include "Particle/ParticleSet.h"
18 #include "Particle/ParticleSetPool.h"
19 #include "QMCWaveFunctions/WaveFunctionComponent.h"
20 #include "BsplineFactory/EinsplineSetBuilder.h"
21 #include "BsplineFactory/EinsplineSpinorSetBuilder.h"
22 #include <ResourceCollection.h>
23 
24 #include <stdio.h>
25 #include <string>
26 #include <limits>
27 
28 using std::string;
29 
30 namespace qmcplusplus
31 {
32 
33 TEST_CASE("Einspline SPO from HDF NiO a16 97 electrons", "[wavefunction]")
34 {
35  Communicate* c = OHMMS::Controller;
36 
37  ParticleSet::ParticleLayout lattice;
38  lattice.R = {3.94055, 3.94055, 7.8811, 3.94055, 3.94055, -7.8811, -7.8811, 7.8811, 0};
39 
40  ParticleSetPool ptcl = ParticleSetPool(c);
41  ptcl.setSimulationCell(lattice);
42  auto ions_uptr = std::make_unique<ParticleSet>(ptcl.getSimulationCell());
43  auto elec_uptr = std::make_unique<ParticleSet>(ptcl.getSimulationCell());
44  ParticleSet& ions_(*ions_uptr);
45  ParticleSet& elec_(*elec_uptr);
46 
47  ions_.setName("ion");
48  ptcl.addParticleSet(std::move(ions_uptr));
49  ions_.create({8, 8});
50  ions_.R[0] = {0.75, 0.25, 0};
51  ions_.R[1] = {0.75, 0.25, 0.5};
52  ions_.R[2] = {0.75, 0.75, 0.25};
53  ions_.R[3] = {0.75, 0.75, 0.75};
54  ions_.R[4] = {0.25, 0.75, 0};
55  ions_.R[5] = {0.25, 0.25, 0.25};
56  ions_.R[6] = {0.25, 0.75, 0.5};
57  ions_.R[7] = {0.25, 0.25, 0.75};
58  ions_.R[8] = {0, 0, 0};
59  ions_.R[9] = {0, 0, 0.5};
60  ions_.R[10] = {0, 0.5, 0.25};
61  ions_.R[11] = {0, 0.5, 0.75};
62  ions_.R[12] = {0.5, 0.5, 0};
63  ions_.R[13] = {0.5, 0, 0.25};
64  ions_.R[14] = {0.5, 0.5, 0.5};
65  ions_.R[15] = {0.5, 0, 0.75};
66 
67  // convert to cartesian
68  ions_.R.setUnit(PosUnit::Lattice);
69  ions_.convert2Cart(ions_.R);
70 
71  // printout coordinates
72  SpeciesSet& ispecies = ions_.getSpeciesSet();
73  int Oidx = ispecies.addSpecies("O");
74  int Niidx = ispecies.addSpecies("Ni");
75  ions_.print(app_log());
76 
77  elec_.setName("elec");
78  ptcl.addParticleSet(std::move(elec_uptr));
79  elec_.create({97});
80  elec_.R[0] = {0.0, 0.0, 0.0};
81  elec_.R[1] = {0.0, 1.0, 0.0};
82  elec_.R[2] = {0.0, 1.1, 0.0};
83  elec_.R[3] = {0.0, 1.2, 0.0};
84  elec_.R[4] = {0.0, 1.3, 0.0};
85 
86  SpeciesSet& tspecies = elec_.getSpeciesSet();
87  int upIdx = tspecies.addSpecies("u");
88  int chargeIdx = tspecies.addAttribute("charge");
89  tspecies(chargeIdx, upIdx) = -1;
90 
91  //diamondC_2x1x1
92  const char* particles = R"(<tmp>
93 <determinantset type="einspline" href="NiO-fcc-supertwist111-supershift000-S4.h5" tilematrix="1 0 0 0 1 1 0 2 -2" twistnum="0" source="ion" meshfactor="1.0" precision="float" size="97" gpu="omptarget"/>
94 </tmp>
95 )";
96 
97  Libxml2Document doc;
98  bool okay = doc.parseFromString(particles);
99  REQUIRE(okay);
100 
101  xmlNodePtr root = doc.getRoot();
102 
103  xmlNodePtr ein1 = xmlFirstElementChild(root);
104 
105  EinsplineSetBuilder einSet(elec_, ptcl.getPool(), c, ein1);
106  auto spo = einSet.createSPOSetFromXML(ein1);
107  REQUIRE(spo);
108 
109  // for vgl
110  SPOSet::ValueMatrix psiM(elec_.R.size(), spo->getOrbitalSetSize());
111  SPOSet::GradMatrix dpsiM(elec_.R.size(), spo->getOrbitalSetSize());
112  SPOSet::ValueMatrix d2psiM(elec_.R.size(), spo->getOrbitalSetSize());
113  spo->evaluate_notranspose(elec_, 0, elec_.R.size(), psiM, dpsiM, d2psiM);
114 
115 #if !defined(QMC_COMPLEX)
116  // real part
117  // due to the different ordering of bands skip the tests on CUDA+Real builds
118  // checking evaluations, reference values are not independently generated.
119  // value
120  CHECK(std::real(psiM[1][0]) == Approx(-2.6693785191));
121  CHECK(std::real(psiM[1][1]) == Approx(-2.669519186));
122  // grad
123  CHECK(std::real(dpsiM[1][0][0]) == Approx(-0.0002679008));
124  CHECK(std::real(dpsiM[1][0][1]) == Approx(7.622625351));
125  CHECK(std::real(dpsiM[1][0][2]) == Approx(-0.0003001692));
126  CHECK(std::real(dpsiM[1][1][0]) == Approx(-0.0002825252));
127  CHECK(std::real(dpsiM[1][1][1]) == Approx(7.6229834557));
128  CHECK(std::real(dpsiM[1][1][2]) == Approx(-0.0002339484));
129  // lapl
130  CHECK(std::real(d2psiM[1][0]) == Approx(-2.6130394936).epsilon(0.0001));
131  CHECK(std::real(d2psiM[1][1]) == Approx(-2.6067698002).epsilon(0.0001));
132 #endif
133 
134 #if defined(QMC_COMPLEX)
135  // imaginary part
136  // value
137  CHECK(std::imag(psiM[1][0]) == Approx(2.6693942547));
138  CHECK(std::imag(psiM[1][1]) == Approx(-2.6693785191));
139  // grad
140  CHECK(std::imag(dpsiM[1][0][0]) == Approx(-0.000179037));
141  CHECK(std::imag(dpsiM[1][0][1]) == Approx(-7.6221408844));
142  CHECK(std::imag(dpsiM[1][0][2]) == Approx(-0.0000232533));
143  CHECK(std::imag(dpsiM[1][1][0]) == Approx(-0.0002679008));
144  CHECK(std::imag(dpsiM[1][1][1]) == Approx(7.622625351));
145  CHECK(std::imag(dpsiM[1][1][2]) == Approx(-0.0003001692));
146  // lapl
147  CHECK(std::imag(d2psiM[1][0]) == Approx(2.6158618927).epsilon(0.0001));
148  CHECK(std::imag(d2psiM[1][1]) == Approx(-2.6130394936).epsilon(0.0001));
149 #endif
150 
151  // test batched interfaces
152 
153  ParticleSet elec_2(elec_);
154  // interchange positions
155  elec_2.R[0] = elec_.R[1];
156  elec_2.R[1] = elec_.R[0];
157  RefVectorWithLeader<ParticleSet> p_list(elec_);
158  p_list.push_back(elec_);
159  p_list.push_back(elec_2);
160 
161  std::unique_ptr<SPOSet> spo_2(spo->makeClone());
162  RefVectorWithLeader<SPOSet> spo_list(*spo);
163  spo_list.push_back(*spo);
164  spo_list.push_back(*spo_2);
165 
166  ResourceCollection pset_res("test_pset_res");
167  ResourceCollection spo_res("test_spo_res");
168 
169  elec_.createResource(pset_res);
170  spo->createResource(spo_res);
171 
172  ResourceCollectionTeamLock<ParticleSet> mw_pset_lock(pset_res, p_list);
173  ResourceCollectionTeamLock<SPOSet> mw_sposet_lock(spo_res, spo_list);
174 
175  SPOSet::ValueVector psi(spo->getOrbitalSetSize());
176  SPOSet::GradVector dpsi(spo->getOrbitalSetSize());
177  SPOSet::ValueVector d2psi(spo->getOrbitalSetSize());
178  SPOSet::ValueVector psi_2(spo->getOrbitalSetSize());
179  SPOSet::GradVector dpsi_2(spo->getOrbitalSetSize());
180  SPOSet::ValueVector d2psi_2(spo->getOrbitalSetSize());
181 
182  RefVector<SPOSet::ValueVector> psi_v_list;
183  RefVector<SPOSet::GradVector> dpsi_v_list;
184  RefVector<SPOSet::ValueVector> d2psi_v_list;
185 
186  psi_v_list.push_back(psi);
187  psi_v_list.push_back(psi_2);
188  dpsi_v_list.push_back(dpsi);
189  dpsi_v_list.push_back(dpsi_2);
190  d2psi_v_list.push_back(d2psi);
191  d2psi_v_list.push_back(d2psi_2);
192 
193  spo->mw_evaluateVGL(spo_list, p_list, 0, psi_v_list, dpsi_v_list, d2psi_v_list);
194 #if !defined(QMC_COMPLEX)
195  // real part
196  // due to the different ordering of bands skip the tests on CUDA+Real builds
197  // checking evaluations, reference values are not independently generated.
198  // value
199  CHECK(std::real(psi_v_list[1].get()[0]) == Approx(-2.6693785191));
200  CHECK(std::real(psi_v_list[1].get()[1]) == Approx(-2.669519186));
201  // grad
202  CHECK(std::real(dpsi_v_list[1].get()[0][0]) == Approx(-0.0002679008));
203  CHECK(std::real(dpsi_v_list[1].get()[0][1]) == Approx(7.622625351));
204  CHECK(std::real(dpsi_v_list[1].get()[0][2]) == Approx(-0.0003001692));
205  CHECK(std::real(dpsi_v_list[1].get()[1][0]) == Approx(-0.0002825252));
206  CHECK(std::real(dpsi_v_list[1].get()[1][1]) == Approx(7.6229834557));
207  CHECK(std::real(dpsi_v_list[1].get()[1][2]) == Approx(-0.0002339484));
208  // lapl
209  CHECK(std::real(d2psi_v_list[1].get()[0]) == Approx(-2.6130394936).epsilon(0.0001));
210  CHECK(std::real(d2psi_v_list[1].get()[1]) == Approx(-2.6067698002).epsilon(0.0001));
211 #endif
212 
213 #if defined(QMC_COMPLEX)
214  // imaginary part
215  // value
216  CHECK(std::imag(psi_v_list[1].get()[0]) == Approx(2.6693942547));
217  CHECK(std::imag(psi_v_list[1].get()[1]) == Approx(-2.6693785191));
218  // grad
219  CHECK(std::imag(dpsi_v_list[1].get()[0][0]) == Approx(-0.000179037));
220  CHECK(std::imag(dpsi_v_list[1].get()[0][1]) == Approx(-7.6221408844));
221  CHECK(std::imag(dpsi_v_list[1].get()[0][2]) == Approx(-0.0000232533));
222  CHECK(std::imag(dpsi_v_list[1].get()[1][0]) == Approx(-0.0002679008));
223  CHECK(std::imag(dpsi_v_list[1].get()[1][1]) == Approx(7.622625351));
224  CHECK(std::imag(dpsi_v_list[1].get()[1][2]) == Approx(-0.0003001692));
225  // lapl
226  CHECK(std::imag(d2psi_v_list[1].get()[0]) == Approx(2.6158618927).epsilon(0.0001));
227  CHECK(std::imag(d2psi_v_list[1].get()[1]) == Approx(-2.6130394936).epsilon(0.0001));
228 #endif
229 
230  const size_t nw = 2;
231  std::vector<SPOSet::ValueType> ratio_v(nw);
232  std::vector<SPOSet::GradType> grads_v(nw);
233 
234  Vector<SPOSet::ValueType, OffloadPinnedAllocator<SPOSet::ValueType>> inv_row(5);
235  inv_row[0] = 0.1;
236  inv_row[1] = 0.2;
237  inv_row[2] = 0.3;
238  inv_row[3] = 0.4;
239  inv_row[4] = 0.5;
240  inv_row.updateTo();
241 
242  std::vector<const SPOSet::ValueType*> inv_row_ptr(nw, inv_row.device_data());
243 
244  SPOSet::OffloadMWVGLArray phi_vgl_v;
245  phi_vgl_v.resize(QMCTraits::DIM_VGL, nw, 5);
246  spo->mw_evaluateVGLandDetRatioGrads(spo_list, p_list, 0, inv_row_ptr, phi_vgl_v, ratio_v, grads_v);
247  phi_vgl_v.updateFrom();
248 #if !defined(QMC_COMPLEX)
249  CHECK(std::real(ratio_v[0]) == Approx(-0.4838374162));
250  CHECK(std::real(grads_v[0][0]) == Approx(-24.6573209338));
251  CHECK(std::real(grads_v[0][1]) == Approx(24.6573187288));
252  CHECK(std::real(grads_v[0][2]) == Approx(-0.0000002709));
253  CHECK(std::real(phi_vgl_v(0, 0, 0)) == Approx(-1.6125482321));
254  CHECK(std::real(ratio_v[1]) == Approx(-2.4885484445));
255  CHECK(std::real(grads_v[1][0]) == Approx(-0.6732621026));
256  CHECK(std::real(grads_v[1][1]) == Approx(-2.6037152796));
257  CHECK(std::real(grads_v[1][2]) == Approx(-0.0001673779));
258  CHECK(std::real(phi_vgl_v(0, 1, 0)) == Approx(-2.6693785191));
259 #else
260  CHECK(std::imag(ratio_v[0]) == Approx(0.0005530113));
261  CHECK(std::imag(grads_v[0][0]) == Approx(-0.000357729));
262  CHECK(std::imag(grads_v[0][1]) == Approx(-0.0005462062));
263  CHECK(std::imag(grads_v[0][2]) == Approx(0.0004516766));
264  CHECK(std::imag(phi_vgl_v(0, 0, 0)) == Approx(1.6124026775));
265  CHECK(std::imag(ratio_v[1]) == Approx(0.0009635784));
266  CHECK(std::imag(grads_v[1][0]) == Approx(0.0001453869));
267  CHECK(std::imag(grads_v[1][1]) == Approx(-0.0000331457));
268  CHECK(std::imag(grads_v[1][2]) == Approx(0.0000295465));
269  CHECK(std::imag(phi_vgl_v(0, 1, 0)) == Approx(2.6693942547));
270 #endif
271 }
272 } // namespace qmcplusplus
qmcplusplus::CrystalLattice
a class that defines a supercell in D-dimensional Euclean space.
Definition: CrystalLattice.h:55
qmcplusplus::ParticleSet::setName
void setName(const std::string &aname)
Definition: ParticleSet.h:237
qmcplusplus::ParticleSetPool::setSimulationCell
void setSimulationCell(const SimulationCell &simulation_cell)
set simulation cell
Definition: ParticleSetPool.h:98
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::Units::real
QMCTraits::RealType real
Definition: unit_conversion.h:23
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::Vector::device_data
pointer device_data()
Return the device_ptr matching X if this is a vector attached or owning dual space memory...
Definition: OhmmsVector.h:245
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
EinsplineSetBuilder.h
Builder class for einspline-based SPOSet objects.
qmcplusplus::SPOSet::ValueMatrix
OrbitalSetTraits< ValueType >::ValueMatrix ValueMatrix
Definition: SPOSet.h:50
qmcplusplus::ParticleSet::createResource
void createResource(ResourceCollection &collection) const
initialize a shared resource and hand it to a collection
Definition: ParticleSet.cpp:944
EinsplineSpinorSetBuilder.h
Derives EinsplineSetBuilder.
OHMMS::Controller
Communicate * Controller
Global Communicator for a process.
Definition: Communicate.cpp:35
Array::resize
void resize(const std::array< SIZET, D > &dims)
Resize the container.
Definition: OhmmsArray.h:65
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
qmcplusplus::SPOSet::GradMatrix
OrbitalSetTraits< ValueType >::GradMatrix GradMatrix
Definition: SPOSet.h:52
Communicate
Wrapping information on parallelism.
Definition: Communicate.h:68
qmcplusplus::lattice
CrystalLattice< OHMMS_PRECISION, OHMMS_DIM > lattice
Definition: test_SpinDensityNew.cpp:142
qmcplusplus::ParticleSet
Specialized paritlce class for atomistic simulations.
Definition: ParticleSet.h:55
qmcplusplus::imag
float imag(const float &c)
imaginary part of a scalar. Cannot be replaced by std::imag due to AFQMC specific needs...
Definition: complex_help.hpp:91
qmcplusplus::Vector::size
size_type size() const
return the current size
Definition: OhmmsVector.h:162
qmcplusplus::REQUIRE
REQUIRE(std::filesystem::exists(filename))
qmcplusplus::ParticleSetPool
Manage a collection of ParticleSet objects.
Definition: ParticleSetPool.h:34
qmcplusplus::SPOSet::ValueVector
OrbitalSetTraits< ValueType >::ValueVector ValueVector
Definition: SPOSet.h:49
qmcplusplus::ParticleSet::R
ParticlePos R
Position.
Definition: ParticleSet.h:79
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
qmcplusplus::RefVector
std::vector< std::reference_wrapper< T > > RefVector
Definition: template_types.hpp:32
qmcplusplus::SPOSet::GradVector
OrbitalSetTraits< ValueType >::GradVector GradVector
Definition: SPOSet.h:51
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::ResourceCollectionTeamLock
handles acquire/release resource by the consumer (RefVectorWithLeader type).
Definition: ResourceCollection.h:58
WaveFunctionComponent.h
Declaration of WaveFunctionComponent.
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
qmcplusplus::Vector::updateTo
void updateTo(size_type size=0, std::ptrdiff_t offset=0)
Definition: OhmmsVector.h:263
qmcplusplus::EinsplineSetBuilder::createSPOSetFromXML
std::unique_ptr< SPOSet > createSPOSetFromXML(xmlNodePtr cur) override
initialize the Antisymmetric wave function for electrons
Definition: EinsplineSetBuilder_createSPOs.cpp:106
Array
A D-dimensional Array class based on PETE.
Definition: OhmmsArray.h:25
qmcplusplus::ParticleSetPool::getSimulationCell
const auto & getSimulationCell() const
get simulation cell
Definition: ParticleSetPool.h:95
ParticleSetPool.h
Declaration of ParticleSetPool.