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 #include "OhmmsData/Libxml2Doc.h" 15 #include "Lattice/CrystalLattice.h" 16 #include "LongRange/StructFact.h" 17 #include "Particle/ParticleSet.h" 18 #include "Particle/DistanceTableData.h" 19 #include "QMCHamiltonians/SkAllEstimator.h" 20 #include "Particle/ParticleSetPool.h" 21 #include <stdio.h> 22 #include <string> 23 using std::string; 24 25 26 /* 27 -- jpt 25/03/2020 -- 28 Deterministic unit test for the SkAllEstimator class. 29 Despite the plumbing, it seems that _for the moment_ (March 2020) 30 that only the e-e structure factor is computed and stored. 31 32 This is the important piece, as we use it for finite size corrections. 33 However, it looks like one might be able to generalize it for e-i and i-i 34 at a later date. If that happens, this test should be updated! 35 */ 36 37 38 namespace qmcplusplus 39 { 40 TEST_CASE("SkAll", "[hamiltonian]") 41 { 42 // Boiler plate setup 43 std::cout << std::fixed; 44 std::cout << std::setprecision(8); 45 typedef QMCTraits::RealType RealType; 46 47 Communicate* c; 48 c = OHMMS::Controller; 49 50 // XML parser 51 Libxml2Document doc; 52 53 54 /* 55 Minimal xml input blocks 56 */ 57 // Lattice block 58 const char* lat_xml = "<simulationcell> \ 59 <parameter name=\"bconds\"> \ 60 p p p \ 61 </parameter> \ 62 <parameter name=\"LR_dim_cutoff\" > 6 </parameter>\ 63 <parameter name=\"rs\" > 1.0 </parameter>\ 64 <parameter name=\"nparticles\" > 8 </parameter>\ 65 </simulationcell>"; 66 67 // Particleset block for the electrons 68 const char* elec_pset_xml = "<particleset name=\"e\" random=\"yes\"> \ 69 <group name=\"u\" size=\"4\"> \ 70 <parameter name=\"charge\" > -1 </parameter>\ 71 <parameter name=\"mass\" > 1.0 </parameter>\ 72 </group> \ 73 <group name=\"d\" size=\"4\"> \ 74 <parameter name=\"charge\" > -1 </parameter>\ 75 <parameter name=\"mass\" > 1.0 </parameter>\ 76 </group> \ 77 </particleset>"; 78 79 // Particleset block for the ions 80 const char* ion_pset_xml = "<particleset name=\"i\" size=\"4\"> \ 81 <group name=\"He\"> \ 82 <parameter name=\"charge\" > 2 </parameter>\ 83 <parameter name=\"valence\" > 2 </parameter>\ 84 <parameter name=\"atomicnumber\"> 2 </parameter>\ 85 </group> \ 86 <attrib name=\"position\" datatype=\"posArray\" \ 87 condition=\"1\"> \ 88 0.00 0.00 0.00 \ 89 0.25 0.25 0.25 \ 90 0.50 0.50 0.50 \ 91 0.75 0.75 0.75 \ 92 </attrib> \ 93 <attrib name=\"ionid\" datatype=\"stringArray\"> \ 94 He He He He \ 95 </attrib> \ 96 </particleset>"; 97 98 // SKAllEstimator block, seems that writeionion does nothing? 99 const char* skall_xml = "<estimator \ 100 name=\"Skall\" type=\"skall\" \ 101 source=\"i\" target=\"e\" hdf5=\"yes\" \ 102 />"; 103 104 // Build a Lattice 105 CrystalLattice<OHMMS_PRECISION, OHMMS_DIM> lattice; 106 lattice.BoxBConds = true; // periodic 107 lattice.R.diagonal(2.0); 108 lattice.reset(); 109 110 // Read in xml, add to pset_builder below 111 bool lat_okay = doc.parseFromString(lat_xml); 112 REQUIRE(lat_okay); 113 xmlNodePtr lat_xml_root = doc.getRoot(); 114 115 std::cout << "\n\n\ntest_SkAllEstimator: START\n"; 116 117 // Build a ParticleSetPool - makes ParticleSets 118 ParticleSetPool pset_builder(c, "pset_builder"); 119 120 // First attach the Lattice defined above 121 pset_builder.putLattice(lat_xml_root); 122 123 // Now build the elec ParticleSet 124 bool elec_pset_okay = doc.parseFromString(elec_pset_xml); 125 REQUIRE(elec_pset_okay); 126 xmlNodePtr elec_pset_xml_root = doc.getRoot(); 127 pset_builder.put(elec_pset_xml_root); 128 129 // Get the (now assembled) elec ParticleSet, sanity check, report 130 ParticleSet* elec = pset_builder.getParticleSet("e"); 131 elec->Lattice = lattice; // copy in the new Lattice 132 REQUIRE(elec->SameMass); 133 REQUIRE(elec->getName() == "e"); 134 135 // Move the particles manually onto B1 lattice 136 // NB: Spins are grouped contiguously 137 // Up spins 138 elec->R[0][0] = 0.0; 139 elec->R[0][1] = 0.0; 140 elec->R[0][2] = 0.0; 141 elec->R[1][0] = 1.0; 142 elec->R[1][1] = 1.0; 143 elec->R[1][2] = 0.0; 144 elec->R[2][0] = 1.0; 145 elec->R[2][1] = 0.0; 146 elec->R[2][2] = 1.0; 147 elec->R[3][0] = 0.0; 148 elec->R[3][1] = 1.0; 149 elec->R[3][2] = 1.0; 150 151 // Down spins 152 elec->R[4][0] = 1.0; 153 elec->R[4][1] = 0.0; 154 elec->R[4][2] = 0.0; 155 elec->R[5][0] = 0.0; 156 elec->R[5][1] = 1.0; 157 elec->R[5][2] = 0.0; 158 elec->R[6][0] = 0.0; 159 elec->R[6][1] = 0.0; 160 elec->R[6][2] = 1.0; 161 elec->R[7][0] = 1.0; 162 elec->R[7][1] = 1.0; 163 elec->R[7][2] = 1.0; 164 165 elec->get(std::cout); // print particleset info to stdout 166 167 168 // Get the (now assembled) ion ParticleSet, sanity check, report 169 bool ion_pset_okay = doc.parseFromString(ion_pset_xml); 170 REQUIRE(ion_pset_okay); 171 xmlNodePtr ion_pset_xml_root = doc.getRoot(); 172 pset_builder.put(ion_pset_xml_root); 173 174 // Finally, make the ion ParticleSet, sanity check, report 175 // It seems that we need this only to construct skall, but it 176 // is never used to evaluate the estimator in this test. 177 ParticleSet* ion = pset_builder.getParticleSet("i"); 178 ion->Lattice = lattice; // copy in the new Lattice 179 REQUIRE(ion->SameMass); 180 REQUIRE(ion->getName() == "i"); 181 ion->get(std::cout); // print particleset info to stdout 182 183 184 // Set up the distance table, match expected layout 185 const int ee_table_id = elec->addTable(*elec); 186 187 const auto& dii(elec->getDistTable(ee_table_id)); 188 elec->update(); // distance table evaluation here 189 190 // Check that the skall xml block is valid 191 bool skall_okay = doc.parseFromString(skall_xml); 192 REQUIRE(skall_okay); 193 xmlNodePtr skall_xml_root = doc.getRoot(); 194 195 // Make a SkAllEstimator, call put() to set up internals 196 SkAllEstimator skall(*ion, *elec); 197 skall.put(skall_xml_root); 198 skall.addObservables(elec->PropertyList, elec->Collectables); 199 skall.get(app_log()); // pretty print settings 200 201 // Hack to make a walker so that tWalker points to something 202 // Only used to set tWalker->Weight = 1 so that skall->evaluate() 203 // doesn't segfault. 204 // NB: setHistories(dummy) attaches dummy to tWalker 205 ParticleSet::Walker_t dummy = ParticleSet::Walker_t(1); 206 skall.setHistories(dummy); 207 skall.evaluate(*elec); 208 209 // s(k) is computed by & held by the ParticleSet. SkAll just 210 // takes that pre-computed s(k) and pulls out rho(k).. 211 // In order to compare to analytic result, need the list 212 // of k-vectors in cartesian coordinates. 213 // Luckily, ParticleSet stores that in SK->KLists.kpts_cart 214 int nkpts = elec->SK->KLists.numk; 215 std::cout << "\n"; 216 std::cout << "SkAll results:\n"; 217 std::cout << std::fixed; 218 std::cout << std::setprecision(6); 219 std::cout << std::setw(4) << "i" 220 << " " << std::setw(8) << "kx" 221 << " " << std::setw(8) << "ky" 222 << " " << std::setw(8) << "kz" 223 << " " << std::setw(8) << "rhok_r" 224 << " " << std::setw(8) << "rhok_i" 225 << " " << std::setw(8) << "c.c." 226 << "\n"; 227 std::cout << "================================================================\n"; 228 229 // Extract rhok out of Collectables, print values 230 auto rhok = elec->Collectables; 231 std::cout << std::fixed; 232 std::cout << std::setprecision(5); 233 for (int k = 0; k < nkpts; k++) 234 { 235 auto kvec = elec->SK->KLists.kpts_cart[k]; 236 RealType kx = kvec[0]; 237 RealType ky = kvec[1]; 238 RealType kz = kvec[2]; 239 RealType rk_r = rhok[nkpts + k]; 240 RealType rk_i = rhok[2 * nkpts + k]; 241 RealType rk_cc = rhok[k]; 242 std::cout << std::setw(4) << k << " " << std::setw(8) << kx << " " << std::setw(8) << ky << " " << std::setw(8) 243 << kz << " " << std::setw(8) << rk_r << " " << std::setw(8) << rk_i << " " << std::setw(8) << rk_cc 244 << "\n"; 245 } 246 247 /* 248 Verify against analytic result: 249 NB: MUST match xml input!!! 250 rho(-1.94889, 0.00000, 0.00000)= 2.52341 + -3.71748i 251 rho( 0.00000, -1.94889, 0.00000)= 2.52341 + -3.71748i 252 rho( 0.00000, 0.00000, -1.94889)= 2.52341 + -3.71748i 253 rho( 0.00000, 0.00000, 1.94889)= 2.52341 + 3.71748i 254 rho( 0.00000, 1.94889, 0.00000)= 2.52341 + 3.71748i 255 rho( 1.94889, 0.00000, 0.00000)= 2.52341 + 3.71748i 256 rho(-1.94889, -1.94889, 0.00000)= -0.93151 + -2.34518i 257 rho(-1.94889, 0.00000, -1.94889)= -0.93151 + -2.34518i 258 rho(-1.94889, 0.00000, 1.94889)= 2.52341 + 0.00000i 259 rho(-1.94889, 1.94889, 0.00000)= 2.52341 + 0.00000i 260 rho( 0.00000, -1.94889, -1.94889)= -0.93151 + -2.34518i 261 rho( 0.00000, -1.94889, 1.94889)= 2.52341 + 0.00000i 262 rho( 0.00000, 1.94889, -1.94889)= 2.52341 + 0.00000i 263 rho( 0.00000, 1.94889, 1.94889)= -0.93151 + 2.34518i 264 rho( 1.94889, -1.94889, 0.00000)= 2.52341 + 0.00000i 265 rho( 1.94889, 0.00000, -1.94889)= 2.52341 + 0.00000i 266 rho( 1.94889, 0.00000, 1.94889)= -0.93151 + 2.34518i 267 rho( 1.94889, 1.94889, 0.00000)= -0.93151 + 2.34518i 268 rho(-1.94889, -1.94889, -1.94889)= -1.38359 + -0.30687i 269 rho(-1.94889, -1.94889, 1.94889)= 0.79595 + -1.17259i 270 rho(-1.94889, 1.94889, -1.94889)= 0.79595 + -1.17259i 271 rho(-1.94889, 1.94889, 1.94889)= 0.79595 + 1.17259i 272 rho( 1.94889, -1.94889, -1.94889)= 0.79595 + -1.17259i 273 rho( 1.94889, -1.94889, 1.94889)= 0.79595 + 1.17259i 274 rho( 1.94889, 1.94889, -1.94889)= 0.79595 + 1.17259i 275 rho( 1.94889, 1.94889, 1.94889)= -1.38359 + 0.30687i 276 */ 277 278 const RealType eps = 1E-04; // tolerance 279 REQUIRE(std::fabs(rhok[26] - 2.52341) < eps); 280 REQUIRE(std::fabs(rhok[52] + 3.71748) < eps); 281 REQUIRE(std::fabs(rhok[26 + 1] - 2.52341) < eps); 282 REQUIRE(std::fabs(rhok[52 + 1] + 3.71748) < eps); 283 REQUIRE(std::fabs(rhok[26 + 2] - 2.52341) < eps); 284 REQUIRE(std::fabs(rhok[52 + 2] + 3.71748) < eps); 285 286 REQUIRE(std::fabs(rhok[26 + 3] - 2.52341) < eps); 287 REQUIRE(std::fabs(rhok[52 + 3] - 3.71748) < eps); 288 REQUIRE(std::fabs(rhok[26 + 4] - 2.52341) < eps); 289 REQUIRE(std::fabs(rhok[52 + 4] - 3.71748) < eps); 290 REQUIRE(std::fabs(rhok[26 + 5] - 2.52341) < eps); 291 REQUIRE(std::fabs(rhok[52 + 5] - 3.71748) < eps); 292 293 REQUIRE(std::fabs(rhok[26 + 6] + 0.93151) < eps); 294 REQUIRE(std::fabs(rhok[52 + 6] + 2.34518) < eps); 295 REQUIRE(std::fabs(rhok[26 + 7] + 0.93151) < eps); 296 REQUIRE(std::fabs(rhok[52 + 7] + 2.34518) < eps); 297 REQUIRE(std::fabs(rhok[26 + 8] - 2.52341) < eps); 298 REQUIRE(std::fabs(rhok[52 + 8] - 0.00000) < eps); 299 300 REQUIRE(std::fabs(rhok[26 + 9] - 2.52341) < eps); 301 REQUIRE(std::fabs(rhok[52 + 9] - 0.00000) < eps); 302 REQUIRE(std::fabs(rhok[26 + 10] + 0.93151) < eps); 303 REQUIRE(std::fabs(rhok[52 + 10] + 2.34518) < eps); 304 REQUIRE(std::fabs(rhok[26 + 11] - 2.52341) < eps); 305 REQUIRE(std::fabs(rhok[52 + 11] - 0.00000) < eps); 306 307 REQUIRE(std::fabs(rhok[26 + 12] - 2.52341) < eps); 308 REQUIRE(std::fabs(rhok[52 + 12] - 0.00000) < eps); 309 REQUIRE(std::fabs(rhok[26 + 13] + 0.93151) < eps); 310 REQUIRE(std::fabs(rhok[52 + 13] - 2.34518) < eps); 311 REQUIRE(std::fabs(rhok[26 + 14] - 2.52341) < eps); 312 REQUIRE(std::fabs(rhok[52 + 14] - 0.00000) < eps); 313 314 REQUIRE(std::fabs(rhok[26 + 15] - 2.52341) < eps); 315 REQUIRE(std::fabs(rhok[52 + 15] - 0.00000) < eps); 316 REQUIRE(std::fabs(rhok[26 + 16] + 0.93151) < eps); 317 REQUIRE(std::fabs(rhok[52 + 16] - 2.34518) < eps); 318 REQUIRE(std::fabs(rhok[26 + 17] + 0.93151) < eps); 319 REQUIRE(std::fabs(rhok[52 + 17] - 2.34518) < eps); 320 321 REQUIRE(std::fabs(rhok[26 + 18] + 1.38359) < eps); 322 REQUIRE(std::fabs(rhok[52 + 18] + 0.30687) < eps); 323 REQUIRE(std::fabs(rhok[26 + 19] - 0.79595) < eps); 324 REQUIRE(std::fabs(rhok[52 + 19] + 1.17259) < eps); 325 REQUIRE(std::fabs(rhok[26 + 20] - 0.79595) < eps); 326 REQUIRE(std::fabs(rhok[52 + 20] + 1.17259) < eps); 327 328 REQUIRE(std::fabs(rhok[26 + 21] - 0.79595) < eps); 329 REQUIRE(std::fabs(rhok[52 + 21] - 1.17259) < eps); 330 REQUIRE(std::fabs(rhok[26 + 22] - 0.79595) < eps); 331 REQUIRE(std::fabs(rhok[52 + 22] + 1.17259) < eps); 332 REQUIRE(std::fabs(rhok[26 + 23] - 0.79595) < eps); 333 REQUIRE(std::fabs(rhok[52 + 23] - 1.17259) < eps); 334 335 REQUIRE(std::fabs(rhok[26 + 24] - 0.79595) < eps); 336 REQUIRE(std::fabs(rhok[52 + 24] - 1.17259) < eps); 337 REQUIRE(std::fabs(rhok[26 + 25] + 1.38359) < eps); 338 REQUIRE(std::fabs(rhok[52 + 25] - 0.30687) < eps); 339 340 std::cout << "test_SkAllEstimator:: STOP\n"; 341 } 342 } // namespace qmcplusplus 343