Skip to content

Commit

Permalink
Merge pull request QMCPACK#4807 from camelto2/1rdm_spinor
Browse files Browse the repository at this point in the history 
1rdm with spinor support
  • Loading branch information
@PDoakORNL
PDoakORNL authored Dec 6, 2023
2 parents 150dd85 + 55483da commit b0faa12
Show file tree
Hide file tree
Showing 8 changed files with 395 additions and 10 deletions.
198 changes: 189 additions & 9 deletions src/Estimators/OneBodyDensityMatrices.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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";
Expand Down Expand Up @@ -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();
Expand Down Expand Up @@ -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);

Expand All @@ -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.
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -354,6 +361,78 @@ inline void OneBodyDensityMatrices::generateDensitySamples(bool save,
rhocur_ = rho;
}

inline void OneBodyDensityMatrices::generateDensitySamplesWithSpin(bool save,
int steps,
RandomBase<FullPrecReal>& 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<FullPrecReal>& rng)
{
Position diff;
Expand All @@ -362,6 +441,13 @@ OneBodyDensityMatrices::Position OneBodyDensityMatrices::diffuse(const Real sqt,
return diff;
}

OneBodyDensityMatrices::Real OneBodyDensityMatrices::diffuseSpin(const Real sqt, RandomBase<FullPrecReal>& rng)
{
Real diff;
assignGaussRand(&diff, 1, rng);
diff *= sqt;
return diff;
}

inline void OneBodyDensityMatrices::calcDensity(const Position& r, Real& dens, ParticleSet& pset_target)
{
Expand All @@ -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)
{
Expand All @@ -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<MCPWalker>& walkers,
const RefVector<ParticleSet>& psets,
const RefVector<TrialWaveFunction>& wfns,
Expand Down Expand Up @@ -479,7 +606,10 @@ void OneBodyDensityMatrices::generateParticleBasis(ParticleSet& pset_target, std
Matrix<Value>& 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]);
}
Expand All @@ -494,7 +624,10 @@ void OneBodyDensityMatrices::generateSampleBasis(Matrix<Value>& 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];
}
Expand All @@ -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;
Expand All @@ -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)
{
Expand All @@ -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<FullPrecReal>& rng)
{
if (sampling_ == Sampling::METROPOLIS)
Expand All @@ -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;
Expand Down
30 changes: 29 additions & 1 deletion src/Estimators/OneBodyDensityMatrices.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -90,7 +90,9 @@ class OneBodyDensityMatrices : public OperatorEstBase
Vector<Value> basis_norms_;
Vector<Grad> basis_gradients_;
Vector<Value> basis_laplacians_;
Vector<Value> basis_spin_gradients_;
std::vector<Position> rsamples_;
std::vector<Real> ssamples_;
Vector<Real> samples_weights_;
int basis_size_;
std::vector<int> species_sizes_;
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -234,11 +241,19 @@ class OneBodyDensityMatrices : public OperatorEstBase
* *
*/
void generateDensitySamples(bool save, int steps, RandomBase<FullPrecReal>& rng, ParticleSet& pset_target);

/// same as above, but with spin variables included
void generateDensitySamplesWithSpin(bool save, int steps, RandomBase<FullPrecReal>& rng, ParticleSet& pset_target);

void generateSampleRatios(ParticleSet& pset_target,
TrialWaveFunction& psi_target,
std::vector<Matrix<Value>>& Psi_nm);
/// produce a position difference vector from timestep
Position diffuse(const Real sqt, RandomBase<FullPrecReal>& rng);

/// spin diffusion
Real diffuseSpin(const Real sqt, RandomBase<FullPrecReal>& rng);

/** calculate density based on r
* \param[in] r position
* \param[out] dens density
Expand All @@ -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
Expand All @@ -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
Expand All @@ -274,14 +298,18 @@ class OneBodyDensityMatrices : public OperatorEstBase
*/
void generateParticleBasis(ParticleSet& pset_target, std::vector<Matrix<Value>>& 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
Expand Down
2 changes: 2 additions & 0 deletions src/Estimators/tests/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Loading

0 comments on commit b0faa12

Please sign in to comment.