diff --git a/src/Estimators/OneBodyDensityMatrices.cpp b/src/Estimators/OneBodyDensityMatrices.cpp index 82583750f1..0790bf707c 100644 --- a/src/Estimators/OneBodyDensityMatrices.cpp +++ b/src/Estimators/OneBodyDensityMatrices.cpp @@ -36,6 +36,7 @@ OneBodyDensityMatrices::OneBodyDensityMatrices(OneBodyDensityMatricesInput&& obd lattice_(lattice), species_(species), basis_functions_("OneBodyDensityMatrices::basis"), + is_spinor_(pset_target.isSpinor()), timers_("OneBodyDensityMatrix") { my_name_ = "OneBodyDensityMatrices"; @@ -81,6 +82,8 @@ OneBodyDensityMatrices::OneBodyDensityMatrices(OneBodyDensityMatricesInput&& obd break; } rsamples_.resize(samples_); + if (is_spinor_) + ssamples_.resize(samples_); // get the sposets that form the basis auto& sposets = input_.get_basis_sets(); @@ -122,7 +125,6 @@ OneBodyDensityMatrices::OneBodyDensityMatrices(OneBodyDensityMatricesInput&& obd for (int i = 0; i < basis_size_; ++i) basis_norms_[i] = bn_standard; - rsamples_.resize(samples_); samples_weights_.resize(samples_); psi_ratios_.resize(nparticles); @@ -147,6 +149,8 @@ OneBodyDensityMatrices::OneBodyDensityMatrices(OneBodyDensityMatricesInput&& obd { basis_gradients_.resize(basis_size_); basis_laplacians_.resize(basis_size_); + if (is_spinor_) + basis_spin_gradients_.resize(basis_size_); } // so if the input is not normalized, normalize it. @@ -220,9 +224,12 @@ void OneBodyDensityMatrices::generateSamples(const Real weight, case Integrator::UNIFORM: generateUniformSamples(rng); break; - case Integrator::DENSITY: { - generateDensitySamples(save, steps, rng, pset_target); - } + case Integrator::DENSITY: + if (is_spinor_) + generateDensitySamplesWithSpin(save, steps, rng, pset_target); + else + generateDensitySamples(save, steps, rng, pset_target); + break; } if (save) @@ -354,6 +361,78 @@ inline void OneBodyDensityMatrices::generateDensitySamples(bool save, rhocur_ = rho; } +inline void OneBodyDensityMatrices::generateDensitySamplesWithSpin(bool save, + int steps, + RandomBase& rng, + ParticleSet& pset_target) +{ + const auto timestep = input_.get_timestep(); + Real sqt = std::sqrt(timestep); + Real ot = 1.0 / timestep; + Position r = rpcur_; //current position + Position d = dpcur_; //current drift + Real rho = rhocur_; //current density + for (int s = 0; s < steps; ++s) + { + nmoves_++; + Position rp; // trial pos + Position dp; // trial drift + Position ds; // drift sum + Real rhop; // trial density + Real ratio; // dens ratio + Real Pacc; // acc prob + Position diff = diffuse(sqt, rng); // get diffusion + + //now do spin variables + Real spinp; // trial spin + Real dspinp; // trial spindrifty + Real spin_ds; // spin drift sum + Real sdiff = diffuseSpin(sqt, rng); //spin diffusion + if (input_.get_use_drift()) + { + rp = r + diff + d; //update trial position + spinp = spcur_ + sdiff + dspcur_; //update trial spin + calcDensityDriftWithSpin(rp, spinp, rhop, dp, dspinp, pset_target); //get trial drift and density + ratio = rhop / rho; //density ratio + ds = dp + d; //drift sum + spin_ds = dspinp + dspcur_; //spin drift sum + Pacc = ratio * std::exp(-ot * (dot(diff, ds) + .5 * dot(ds, ds))) * + std::exp(-ot * (sdiff * spin_ds) + 0.5 * spin_ds * spin_ds); //acceptance probability + } + else + { + rp = r + diff; //update trial position + spinp = spcur_ + sdiff; //update trial spin + calcDensityWithSpin(rp, spinp, rhop, pset_target); //get trial density + ratio = rhop / rho; //density ratio + Pacc = ratio; //acceptance probability + } + if (rng() < Pacc) + { //accept move + r = rp; + d = dp; + spcur_ = spinp; + dspcur_ = dspinp; + rho = rhop; + naccepted_++; + } + if (save) + { + rsamples_[s] = r; + samples_weights_[s] = 1.0 / rho; + ssamples_[s] = spcur_; + } + } + acceptance_ratio_ = Real(naccepted_) / nmoves_; + + if (input_.get_write_acceptance_ratio() && omp_get_thread_num() == 0) + app_log() << "dm1b acceptance_ratio = " << acceptance_ratio_ << std::endl; + + rpcur_ = r; + dpcur_ = d; + rhocur_ = rho; +} + OneBodyDensityMatrices::Position OneBodyDensityMatrices::diffuse(const Real sqt, RandomBase& rng) { Position diff; @@ -362,6 +441,13 @@ OneBodyDensityMatrices::Position OneBodyDensityMatrices::diffuse(const Real sqt, return diff; } +OneBodyDensityMatrices::Real OneBodyDensityMatrices::diffuseSpin(const Real sqt, RandomBase& rng) +{ + Real diff; + assignGaussRand(&diff, 1, rng); + diff *= sqt; + return diff; +} inline void OneBodyDensityMatrices::calcDensity(const Position& r, Real& dens, ParticleSet& pset_target) { @@ -375,6 +461,20 @@ inline void OneBodyDensityMatrices::calcDensity(const Position& r, Real& dens, P dens /= basis_size_; } +inline void OneBodyDensityMatrices::calcDensityWithSpin(const Position& r, + const Real& s, + Real& dens, + ParticleSet& pset_target) +{ + updateBasisWithSpin(r, s, pset_target); + dens = 0.0; + for (int i = 0; i < basis_size_; ++i) + { + Value b = basis_values_[i]; + dens += std::abs(qmcplusplus::conj(b) * b); + } + dens /= basis_size_; +} void OneBodyDensityMatrices::calcDensityDrift(const Position& r, Real& dens, Position& drift, ParticleSet& pset_target) { @@ -394,6 +494,33 @@ void OneBodyDensityMatrices::calcDensityDrift(const Position& r, Real& dens, Pos dens /= basis_size_; } +void OneBodyDensityMatrices::calcDensityDriftWithSpin(const Position& r, + const Real& s, + Real& dens, + Position& drift, + Real& sdrift, + ParticleSet& pset_target) +{ + updateBasisD012WithSpin(r, s, pset_target); + dens = 0.0; + drift = 0.0; + sdrift = 0.0; + for (int i = 0; i < basis_size_; ++i) + { + const Grad& bg = basis_gradients_[i]; + const Value& bsg = basis_spin_gradients_[i]; + Value b = basis_values_[i]; + Value bc = qmcplusplus::conj(b); + dens += std::abs(bc * b); + for (int d = 0; d < OHMMS_DIM; ++d) + drift[d] += std::real(bc * bg[d]); + sdrift += std::real(bc * bsg); + } + drift *= input_.get_timestep() / dens; + sdrift *= input_.get_timestep() / dens; + dens /= basis_size_; +} + void OneBodyDensityMatrices::accumulate(const RefVector& walkers, const RefVector& psets, const RefVector& wfns, @@ -479,7 +606,10 @@ void OneBodyDensityMatrices::generateParticleBasis(ParticleSet& pset_target, std Matrix& P_nb = phi_nb[s]; for (int n = 0; n < species_sizes_[s]; ++n, ++p) { - updateBasis(pset_target.R[p], pset_target); + if (is_spinor_) + updateBasisWithSpin(pset_target.R[p], pset_target.spins[p], pset_target); + else + updateBasis(pset_target.R[p], pset_target); for (int b = 0; b < basis_size_; ++b, ++nb) P_nb(nb) = qmcplusplus::conj(basis_values_[b]); } @@ -494,7 +624,10 @@ void OneBodyDensityMatrices::generateSampleBasis(Matrix& Phi_mb, int mb = 0; for (int m = 0; m < samples_; ++m) { - updateBasis(rsamples_[m], pset_target); + if (is_spinor_) + updateBasisWithSpin(rsamples_[m], ssamples_[m], pset_target); + else + updateBasis(rsamples_[m], pset_target); for (int b = 0; b < basis_size_; ++b, ++mb) Phi_mb(mb) = basis_values_[b]; } @@ -508,8 +641,25 @@ void OneBodyDensityMatrices::generateSampleRatios(ParticleSet& pset_target, for (int m = 0; m < samples_; ++m) { // get N ratios for the current sample point - pset_target.makeVirtualMoves(rsamples_[m]); - psi_target.evaluateRatiosAlltoOne(pset_target, psi_ratios_); + if (!is_spinor_) + { + pset_target.makeVirtualMoves(rsamples_[m]); + psi_target.evaluateRatiosAlltoOne(pset_target, psi_ratios_); + } + else + { + //note: makeVirtualMoves updates distance tables + //There is not a corresponding "distance table" for the spins, so we need to brute force this by moving each electron and using makeMoveWithSpin, calcRatio, and rejectMove, and resetPhaseDiff, similar to how NonLocalECP works for evaluating quadrature points without virtual particles + int p = 0; + for (int s = 0; s < species_.size(); ++s) + for (int n = 0; n < species_sizes_[s]; ++n, ++p) + { + pset_target.makeMoveWithSpin(p, rsamples_[m] - pset_target.R[p], ssamples_[m] - pset_target.spins[p]); + psi_ratios_[p] = psi_target.calcRatio(pset_target, p); + pset_target.rejectMove(p); + psi_target.resetPhaseDiff(); + } + } // collect ratios into per-species matrices int p = 0; @@ -534,6 +684,15 @@ inline void OneBodyDensityMatrices::updateBasis(const Position& r, ParticleSet& basis_values_[i] *= basis_norms_[i]; } +inline void OneBodyDensityMatrices::updateBasisWithSpin(const Position& r, const Real& s, ParticleSet& pset_target) +{ + // This is ridiculous in the case of splines, still necessary for hybrid/LCAO + pset_target.makeMoveWithSpin(0, r - pset_target.R[0], s - pset_target.spins[0]); + basis_functions_.evaluateValue(pset_target, 0, basis_values_); + pset_target.rejectMove(0); + for (int i = 0; i < basis_size_; ++i) + basis_values_[i] *= basis_norms_[i]; +} inline void OneBodyDensityMatrices::updateBasisD012(const Position& r, ParticleSet& pset_target) { @@ -548,6 +707,21 @@ inline void OneBodyDensityMatrices::updateBasisD012(const Position& r, ParticleS basis_laplacians_[i] *= basis_norms_[i]; } +inline void OneBodyDensityMatrices::updateBasisD012WithSpin(const Position& r, const Real& s, ParticleSet& pset_target) +{ + pset_target.makeMoveWithSpin(0, r - pset_target.R[0], s - pset_target.spins[0]); + basis_functions_.evaluateVGL_spin(pset_target, 0, basis_values_, basis_gradients_, basis_laplacians_, + basis_spin_gradients_); + pset_target.rejectMove(0); + for (int i = 0; i < basis_size_; ++i) + { + basis_values_[i] *= basis_norms_[i]; + basis_gradients_[i] *= basis_norms_[i]; + basis_laplacians_[i] *= basis_norms_[i]; + basis_spin_gradients_[i] *= basis_norms_[i]; + } +} + void OneBodyDensityMatrices::warmupSampling(ParticleSet& pset_target, RandomBase& rng) { if (sampling_ == Sampling::METROPOLIS) @@ -556,7 +730,13 @@ void OneBodyDensityMatrices::warmupSampling(ParticleSet& pset_target, RandomBase { rpcur_ = diffuse(std::sqrt(input_.get_timestep()), rng); rpcur_ += center_; - calcDensityDrift(rpcur_, rhocur_, dpcur_, pset_target); + if (!is_spinor_) + calcDensityDrift(rpcur_, rhocur_, dpcur_, pset_target); + else + { + spcur_ = diffuseSpin(std::sqrt(input_.get_timestep()), rng); + calcDensityDriftWithSpin(rpcur_, spcur_, rhocur_, dpcur_, dspcur_, pset_target); + } } generateSamples(1.0, pset_target, rng, input_.get_warmup_samples()); warmed_up_ = true; diff --git a/src/Estimators/OneBodyDensityMatrices.h b/src/Estimators/OneBodyDensityMatrices.h index 422226c72b..9d8923d899 100644 --- a/src/Estimators/OneBodyDensityMatrices.h +++ b/src/Estimators/OneBodyDensityMatrices.h @@ -90,7 +90,9 @@ class OneBodyDensityMatrices : public OperatorEstBase Vector basis_norms_; Vector basis_gradients_; Vector basis_laplacians_; + Vector basis_spin_gradients_; std::vector rsamples_; + std::vector ssamples_; Vector samples_weights_; int basis_size_; std::vector species_sizes_; @@ -145,6 +147,11 @@ class OneBodyDensityMatrices : public OperatorEstBase /// current density Real rhocur_ = 0.0; + ///spin related variables + Real spcur_; + Real dspcur_; + const bool is_spinor_; + public: /** Standard Constructor * Call this to make a new OBDM this is what you should be calling @@ -234,11 +241,19 @@ class OneBodyDensityMatrices : public OperatorEstBase * * */ void generateDensitySamples(bool save, int steps, RandomBase& rng, ParticleSet& pset_target); + + /// same as above, but with spin variables included + void generateDensitySamplesWithSpin(bool save, int steps, RandomBase& rng, ParticleSet& pset_target); + void generateSampleRatios(ParticleSet& pset_target, TrialWaveFunction& psi_target, std::vector>& Psi_nm); /// produce a position difference vector from timestep Position diffuse(const Real sqt, RandomBase& rng); + + /// spin diffusion + Real diffuseSpin(const Real sqt, RandomBase& rng); + /** calculate density based on r * \param[in] r position * \param[out] dens density @@ -249,6 +264,8 @@ class OneBodyDensityMatrices : public OperatorEstBase * * updateBasis is called */ void calcDensity(const Position& r, Real& dens, ParticleSet& pset_target); + /// same as above, but with spin move + void calcDensityWithSpin(const Position& r, const Real& s, Real& dens, ParticleSet& pset_target); /** calculate density and drift bashed on r * \param[in] r position * \param[out] dens density @@ -261,6 +278,13 @@ class OneBodyDensityMatrices : public OperatorEstBase * * updateBasisD012 is called */ void calcDensityDrift(const Position& r, Real& dens, Position& drift, ParticleSet& pset_target); + /// same as above, but with spin move + void calcDensityDriftWithSpin(const Position& r, + const Real& s, + Real& dens, + Position& drift, + Real& sdrift, + ParticleSet& pset_target); // basis & wavefunction ratio matrix construction /** set Phi_mp to basis vaules per sample @@ -274,14 +298,18 @@ class OneBodyDensityMatrices : public OperatorEstBase */ void generateParticleBasis(ParticleSet& pset_target, std::vector>& phi_nb); - // basis set updates + /// basis set updates void updateBasis(const Position& r, ParticleSet& pset_target); + /// basis set updates with spin + void updateBasisWithSpin(const Position& r, const Real& s, ParticleSet& pset_target); /** evaluates vgl on basis_functions_ for r * sideeffects: * * sets basis_values_, basis_gradients_, basis_laplacians_ * all are normalized by basis norms_ */ void updateBasisD012(const Position& r, ParticleSet& pset_target); + /// same as above, but includes spin gradients + void updateBasisD012WithSpin(const Position& r, const Real& s, ParticleSet& pset_target); /** does some warmup sampling i.e. samples but throws away the results * Only when integrator_ = Integrator::DENSITY * sets rpcur_ initial rpcur + one diffusion step diff --git a/src/Estimators/tests/CMakeLists.txt b/src/Estimators/tests/CMakeLists.txt index 0d97d93a27..b294cb29eb 100644 --- a/src/Estimators/tests/CMakeLists.txt +++ b/src/Estimators/tests/CMakeLists.txt @@ -14,8 +14,10 @@ set(SRC_DIR estimators) set(UTEST_EXE test_${SRC_DIR}) set(UTEST_NAME deterministic-unit_test_${SRC_DIR}) set(UTEST_HDF_INPUT ${qmcpack_SOURCE_DIR}/tests/solids/diamondC_1x1x1_pp/pwscf.pwscf.h5) +set(UTEST_HDF_INPUT_SPINOR ${qmcpack_SOURCE_DIR}/tests/solids/monoO_noncollinear_1x1x1_pp/o2_45deg_spins.pwscf.h5) set(UTEST_DIR ${CMAKE_CURRENT_BINARY_DIR}) maybe_symlink(${UTEST_HDF_INPUT} ${UTEST_DIR}/diamondC_1x1x1.pwscf.h5) +maybe_symlink(${UTEST_HDF_INPUT_SPINOR} ${UTEST_DIR}/o2_45deg_spins.pwscf.h5) set(SRCS test_accumulator.cpp diff --git a/src/Estimators/tests/test_OneBodyDensityMatrices.cpp b/src/Estimators/tests/test_OneBodyDensityMatrices.cpp index 28b57182d5..407fb2e3ed 100644 --- a/src/Estimators/tests/test_OneBodyDensityMatrices.cpp +++ b/src/Estimators/tests/test_OneBodyDensityMatrices.cpp @@ -125,6 +125,17 @@ class OneBodyDensityMatricesTests } } + void testGenerateSamplesForSpinor(onebodydensitymatrices::Inputs input, + OneBodyDensityMatrices& obdm, + ParticleSet& pset_target, + StdRandom& rng) + { + testGenerateSamples(input, obdm, pset_target, rng); + //confirm that the spins were sampled, will be zero if not + CHECK(obdm.is_spinor_ == true); + CHECK(std::abs(obdm.spcur_) > 1e-8); + } + /** Checking approximate equality for complex valued data as * two reals is not consistent with testing practices elsewhere in the code. * Values that are slightly off may now fall in approximation limit properly. @@ -329,6 +340,48 @@ TEST_CASE("OneBodyDensityMatrices::generateSamples", "[estimators]") samplingCaseRunner(valid_obdm_input_grid); } +#ifdef QMC_COMPLEX +TEST_CASE("OneBodyDensityMatrices::generateSamplesForSpinor", "[estimators]") +{ + using namespace testing; + using namespace onebodydensitymatrices; + + using MCPWalker = OperatorEstBase::MCPWalker; + + ProjectData test_project("test", ProjectData::DriverVersion::BATCH); + Communicate* comm; + comm = OHMMS::Controller; + + auto particle_pool = MinimalParticlePool::make_O2_spinor(comm); + auto wavefunction_pool = + MinimalWaveFunctionPool::make_O2_spinor(test_project.getRuntimeOptions(), comm, particle_pool); + auto& pset_target = *(particle_pool.getParticleSet("e")); + auto& species_set = pset_target.getSpeciesSet(); + auto& spo_map = wavefunction_pool.getWaveFunction("wavefunction")->getSPOMap(); + + auto samplingCaseRunner = [&pset_target, &species_set, &spo_map](Inputs test_case) { + Libxml2Document doc; + + bool okay = doc.parseFromString(valid_one_body_density_matrices_input_sections[test_case]); + if (!okay) + throw std::runtime_error("cannot parse OneBodyDensitMatricesInput section"); + xmlNodePtr node = doc.getRoot(); + OneBodyDensityMatricesInput obdmi(node); + + OneBodyDensityMatrices obDenMat(std::move(obdmi), pset_target.getLattice(), species_set, spo_map, pset_target); + + OneBodyDensityMatricesTests obdmt; + //Get control over which rng is used. + //we don't want FakeRandom. + StdRandom rng; + obdmt.testGenerateSamplesForSpinor(test_case, obDenMat, pset_target, rng); + }; + + //Spin sampling only added for density sampling + samplingCaseRunner(valid_obdm_input); +} +#endif + TEST_CASE("OneBodyDensityMatrices::spawnCrowdClone()", "[estimators]") { using namespace testing; diff --git a/src/Particle/tests/MinimalParticlePool.h b/src/Particle/tests/MinimalParticlePool.h index 2903076d8a..2d3471b30d 100644 --- a/src/Particle/tests/MinimalParticlePool.h +++ b/src/Particle/tests/MinimalParticlePool.h @@ -55,6 +55,40 @@ class MinimalParticlePool +)"; + + static constexpr const char* const particles_xml_spinor = R"( + + + + 5.10509515 -3.23993545 0.00000000 + 5.10509515 3.23993545 -0.00000000 + -6.49690625 0.00000000 7.08268015 + + + p p p + + 15 + + + + 6 + 6 + 8 + 29164.3928678 + + -0.00000000 -0.00000000 1.08659253 + 0.00000000 0.00000000 -1.08659253 + + + + + + -1 + 1.0 + + + )"; public: @@ -80,6 +114,29 @@ class MinimalParticlePool return pp; } + + static ParticleSetPool make_O2_spinor(Communicate* c) + { + Libxml2Document doc; + + doc.parseFromString(particles_xml_spinor); + + xmlNodePtr root = doc.getRoot(); + xmlNodePtr sim_cell = xmlFirstElementChild(root); + + ParticleSetPool pp(c); + + // Need to set up simulation cell lattice before reading particle sets + pp.readSimulationCellXML(sim_cell); + + xmlNodePtr part_ion = xmlNextElementSibling(sim_cell); + pp.put(part_ion); + xmlNodePtr part_elec = xmlNextElementSibling(part_ion); + pp.put(part_elec); + pp.randomize(); + + return pp; + } }; } // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/CompositeSPOSet.cpp b/src/QMCWaveFunctions/CompositeSPOSet.cpp index 7110a831b1..6fd4c98303 100644 --- a/src/QMCWaveFunctions/CompositeSPOSet.cpp +++ b/src/QMCWaveFunctions/CompositeSPOSet.cpp @@ -65,6 +65,7 @@ void CompositeSPOSet::add(std::unique_ptr component) component_values.emplace_back(norbs); component_gradients.emplace_back(norbs); component_laplacians.emplace_back(norbs); + component_spin_gradients.emplace_back(norbs); OrbitalSetSize += norbs; component_offsets.push_back(OrbitalSetSize); @@ -114,6 +115,30 @@ void CompositeSPOSet::evaluateVGL(const ParticleSet& P, int iat, ValueVector& ps } } +void CompositeSPOSet::evaluateVGL_spin(const ParticleSet& P, + int iat, + ValueVector& psi, + GradVector& dpsi, + ValueVector& d2psi, + ValueVector& dspin_psi) +{ + int n = 0; + for (int c = 0; c < components.size(); ++c) + { + SPOSet& component = *components[c]; + ValueVector& values = component_values[c]; + GradVector& gradients = component_gradients[c]; + ValueVector& laplacians = component_laplacians[c]; + ValueVector& spin_gradients = component_spin_gradients[c]; + component.evaluateVGL_spin(P, iat, values, gradients, laplacians, spin_gradients); + std::copy(values.begin(), values.end(), psi.begin() + n); + std::copy(gradients.begin(), gradients.end(), dpsi.begin() + n); + std::copy(laplacians.begin(), laplacians.end(), d2psi.begin() + n); + std::copy(spin_gradients.begin(), spin_gradients.end(), dspin_psi.begin() + n); + n += component.size(); + } +} + void CompositeSPOSet::evaluate_notranspose(const ParticleSet& P, int first, int last, diff --git a/src/QMCWaveFunctions/CompositeSPOSet.h b/src/QMCWaveFunctions/CompositeSPOSet.h index 1c03eb356f..ff6761eb4f 100644 --- a/src/QMCWaveFunctions/CompositeSPOSet.h +++ b/src/QMCWaveFunctions/CompositeSPOSet.h @@ -33,6 +33,8 @@ class CompositeSPOSet : public SPOSet std::vector component_gradients; ///temporary storage for laplacians std::vector component_laplacians; + ///temporary storage for spin gradients + std::vector component_spin_gradients; ///store the precomputed offsets std::vector component_offsets; @@ -58,6 +60,13 @@ class CompositeSPOSet : public SPOSet void evaluateVGL(const ParticleSet& P, int iat, ValueVector& psi, GradVector& dpsi, ValueVector& d2psi) override; + void evaluateVGL_spin(const ParticleSet& P, + int iat, + ValueVector& psi, + GradVector& dpsi, + ValueVector& d2psi, + ValueVector& dspin_psi) override; + ///unimplemented functions call this to abort inline void not_implemented(const std::string& method) { diff --git a/src/QMCWaveFunctions/tests/MinimalWaveFunctionPool.h b/src/QMCWaveFunctions/tests/MinimalWaveFunctionPool.h index bb340ab751..6d895dcc84 100644 --- a/src/QMCWaveFunctions/tests/MinimalWaveFunctionPool.h +++ b/src/QMCWaveFunctions/tests/MinimalWaveFunctionPool.h @@ -39,6 +39,20 @@ class MinimalWaveFunctionPool )"; + static constexpr const char* const wf_input_spinor = R"( + + + + + + + + + + + + )"; + public: static WaveFunctionPool make_diamondC_1x1x1(const RuntimeOptions& runtime_options, Communicate* comm, @@ -56,6 +70,23 @@ class MinimalWaveFunctionPool return wp; } + + static WaveFunctionPool make_O2_spinor(const RuntimeOptions& runtime_options, + Communicate* comm, + ParticleSetPool& particle_pool) + { + WaveFunctionPool wp(runtime_options, particle_pool, comm); + + Libxml2Document doc; + bool okay = doc.parseFromString(wf_input_spinor); + REQUIRE(okay); + + xmlNodePtr root = doc.getRoot(); + + wp.put(root); + + return wp; + } }; } // namespace qmcplusplus