diff --git a/src/Estimators/SelfHealingOverlap.cpp b/src/Estimators/SelfHealingOverlap.cpp index 2ac8328380..5c76c8a769 100644 --- a/src/Estimators/SelfHealingOverlap.cpp +++ b/src/Estimators/SelfHealingOverlap.cpp @@ -11,6 +11,7 @@ #include "SelfHealingOverlap.h" #include "TrialWaveFunction.h" #include "QMCWaveFunctions/Fermion/MultiSlaterDetTableMethod.h" +#include "QMCWaveFunctions/Fermion/SlaterDet.h" #include #include @@ -19,19 +20,47 @@ namespace qmcplusplus { SelfHealingOverlap::SelfHealingOverlap(SelfHealingOverlapInput&& inp_, const TrialWaveFunction& wfn, DataLocality dl) - : OperatorEstBase(dl), input_(std::move(inp_)) + : OperatorEstBase(dl), + input_(std::move(inp_)), + wf_type(no_wf), + use_param_deriv(input_.input_section_.get("param_deriv")) { - //my_name_ = input_.get_name(); - - auto& inp = this->input_.input_section_; - auto msd_refvec = wfn.findMSD(); - if (msd_refvec.size() != 1) + auto sd_refvec = wfn.findSD(); + + auto nsd = sd_refvec.size(); + auto nmsd = msd_refvec.size(); + + size_t nparams; + if (nmsd == 1 && nsd == 0) + { // multi-slater-det wavefunction + wf_type = msd_wf; + const MultiSlaterDetTableMethod& msd = msd_refvec[0]; + if (!use_param_deriv) + nparams = msd.getLinearExpansionCoefs().size(); + else + { + throw std::runtime_error("SelfHealingOverlap: use_param_deriv implementation incomplete, needs access to param " + "count from wavefunction component myVars"); + } + if (nparams == 0) + throw std::runtime_error("SelfHealingOverlap: multidet wavefunction has no parameters."); + } + else if (nmsd == 0 && nsd == 1) + { // slater-det wavefunction + throw std::runtime_error("SelfHealingOverlap: slaterdet wavefunction implementation incomplete"); + } + else + { throw std::runtime_error( - "SelfHealingOverlap requires one and only one multi slater determinant component in the trial wavefunction."); + "SelfHealingOverlap requires a single slater or multi-slater determinant component in the trial wavefunction."); + } - const MultiSlaterDetTableMethod& msd = msd_refvec[0]; - const size_t data_size = msd.getLinearExpansionCoefs().size(); +#ifndef QMC_COMPLEX + const size_t data_size = nparams; +#else + const size_t data_size = 2 * nparams; +#endif data_.resize(data_size, 0.0); } @@ -78,31 +107,38 @@ void SelfHealingOverlap::accumulate(const RefVector& walkers, RealType weight = walker.Weight; auto& wcs = psi.getOrbitals(); - // separate jastrow and fermi wavefunction components + // find jastrow wavefunction components std::vector wcs_jastrow; - std::vector wcs_fermi; for (auto& wc : wcs) - if (wc->isFermionic()) - wcs_fermi.push_back(wc.get()); - else + if (!wc->isFermionic()) wcs_jastrow.push_back(wc.get()); - // fermionic must have only one component, and must be multideterminant - assert(wcs_fermi.size() == 1); - WaveFunctionComponent& wf = *wcs_fermi[0]; - if (!wf.isMultiDet()) - throw std::runtime_error("SelfHealingOverlap estimator requires use of multideterminant wavefunction"); - auto msd_refvec = psi.findMSD(); - MultiSlaterDetTableMethod& msd = msd_refvec[0]; - - // collect parameter derivatives: (dpsi/dc_i)/psi - msd.calcIndividualDetRatios(det_ratios); + if (wf_type == msd_wf) + { + auto msd_refvec = psi.findMSD(); + MultiSlaterDetTableMethod& msd = msd_refvec[0]; + // collect parameter derivatives: (dpsi/dc_i)/psi + if (!use_param_deriv) + msd.calcIndividualDetRatios(det_ratios); + else + { + throw std::runtime_error("SelfHealingOverlap: use_param_deriv implementation incomplete, needs call to " + "msd.evaluateDerivatives with correct myVars"); + } + } + else if (wf_type == sd_rot_wf) + { + throw std::runtime_error("SelfHealingOverlap: slaterdet wavefunction implementation incomplete"); + auto sd_refvec = psi.findSD(); + } + else + throw std::runtime_error("SelfHealingOverlap: impossible branch reached, contact the developers"); // collect jastrow prefactor WaveFunctionComponent::LogValue Jval = 0.0; for (auto& wc : wcs_jastrow) Jval += wc->get_log_value(); - auto Jprefactor = std::real(std::exp(-2. * Jval)); + RealType Jprefactor = std::exp(std::real(-2. * Jval)); // accumulate weight (required by all estimators, otherwise inf results) walkers_weight_ += weight; @@ -110,7 +146,15 @@ void SelfHealingOverlap::accumulate(const RefVector& walkers, // accumulate data assert(det_ratios.size() == data_.size()); for (int ic = 0; ic < det_ratios.size(); ++ic) - data_[ic] += weight * Jprefactor * std::real(det_ratios[ic]); // only real supported for now + { +#ifndef QMC_COMPLEX + data_[ic] += weight * Jprefactor * det_ratios[ic]; +#else + auto value = weight * Jprefactor * std::conj(det_ratios[ic]); + data_[2 * ic] += std::real(value); + data_[2 * ic + 1] += std::imag(value); +#endif + } } } diff --git a/src/Estimators/SelfHealingOverlap.h b/src/Estimators/SelfHealingOverlap.h index f3422b2395..a761244db8 100644 --- a/src/Estimators/SelfHealingOverlap.h +++ b/src/Estimators/SelfHealingOverlap.h @@ -33,6 +33,14 @@ class SelfHealingOverlap : public OperatorEstBase using ValueType = QMCTraits::ValueType; using PosType = QMCTraits::PosType; + + enum wf_types + { + msd_wf = 0, + sd_rot_wf, + no_wf + }; + //data members set only during construction const SelfHealingOverlapInput input_; @@ -40,6 +48,13 @@ class SelfHealingOverlap : public OperatorEstBase */ Vector det_ratios; + /// wavefunction type + wf_types wf_type; + + /// use direct parameter derivative for MSD or not + const bool use_param_deriv; + + public: /** Constructor for SelfHealingOverlapInput */ diff --git a/src/Estimators/SelfHealingOverlapInput.h b/src/Estimators/SelfHealingOverlapInput.h index 4b0106f3e0..6987571644 100644 --- a/src/Estimators/SelfHealingOverlapInput.h +++ b/src/Estimators/SelfHealingOverlapInput.h @@ -23,7 +23,7 @@ class SelfHealingOverlapInput { public: using Consumer = SelfHealingOverlap; - using Real = QMCTraits::RealType; + using Real = QMCTraits::RealType; class SelfHealingOverlapInputSection : public InputSection { @@ -32,9 +32,12 @@ class SelfHealingOverlapInput SelfHealingOverlapInputSection() { section_name = "SelfHealingOverlap"; - attributes = {"type", "name"}; + attributes = {"type", "name", "param_deriv"}; strings = {"type", "name"}; - default_values = {{"type", std::string("sh_overlap")},{"name", std::string("sh_overlap")}}; + bools = {"param_deriv"}; + default_values = {{"type", std::string("sh_overlap")}, + {"name", std::string("sh_overlap")}, + {"param_deriv", false}}; } // clang-format: on }; diff --git a/src/QMCWaveFunctions/TrialWaveFunction.cpp b/src/QMCWaveFunctions/TrialWaveFunction.cpp index faf97db21f..0f8e928d9e 100644 --- a/src/QMCWaveFunctions/TrialWaveFunction.cpp +++ b/src/QMCWaveFunctions/TrialWaveFunction.cpp @@ -24,6 +24,7 @@ #include "Concurrency/Info.hpp" #include "type_traits/ConvertToReal.h" #include "NaNguard.h" +#include "Fermion/SlaterDet.h" #include "Fermion/MultiSlaterDetTableMethod.h" namespace qmcplusplus @@ -107,6 +108,15 @@ const SPOSet& TrialWaveFunction::getSPOSet(const std::string& name) const return *spoit->second; } +RefVector TrialWaveFunction::findSD() const +{ + RefVector refs; + for (auto& component : Z) + if (auto* comp_ptr = dynamic_cast(component.get()); comp_ptr) + refs.push_back(*comp_ptr); + return refs; +} + RefVector TrialWaveFunction::findMSD() const { RefVector refs; diff --git a/src/QMCWaveFunctions/TrialWaveFunction.h b/src/QMCWaveFunctions/TrialWaveFunction.h index f3eaa1824a..3d46052709 100644 --- a/src/QMCWaveFunctions/TrialWaveFunction.h +++ b/src/QMCWaveFunctions/TrialWaveFunction.h @@ -38,6 +38,7 @@ namespace qmcplusplus { +class SlaterDet; class MultiSlaterDetTableMethod; /** @ingroup MBWfs @@ -539,6 +540,9 @@ class TrialWaveFunction /// spomap_ reference accessor const SPOMap& getSPOMap() const { return *spomap_; } + /// find SD WFCs if exist + RefVector findSD() const; + /// find MSD WFCs if exist RefVector findMSD() const;