Skip to content

Commit

Permalink
Merge pull request #4504 from camelto2/mw_APIs_for_SOECP
Browse files Browse the repository at this point in the history 
mw_ APIs for SOECPotential and SOECPComponent
  • Loading branch information
@ye-luo
ye-luo authored Mar 10, 2023
2 parents f0571b3 + f5480e3 commit 98c8f2d
Show file tree
Hide file tree
Showing 9 changed files with 790 additions and 81 deletions.
53 changes: 52 additions & 1 deletion src/Particle/VirtualParticleSet.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,8 @@ struct VPMultiWalkerMem : public Resource
Resource* makeClone() const override { return new VPMultiWalkerMem(*this); }
};

VirtualParticleSet::VirtualParticleSet(const ParticleSet& p, int nptcl, size_t dt_count_limit) : ParticleSet(p.getSimulationCell())
VirtualParticleSet::VirtualParticleSet(const ParticleSet& p, int nptcl, size_t dt_count_limit)
: ParticleSet(p.getSimulationCell())
{
setName("virtual");

Expand Down Expand Up @@ -194,4 +195,54 @@ void VirtualParticleSet::mw_makeMoves(const RefVectorWithLeader<VirtualParticleS
ParticleSet::mw_update(p_list);
}

void VirtualParticleSet::mw_makeMovesWithSpin(const RefVectorWithLeader<VirtualParticleSet>& vp_list,
const RefVectorWithLeader<ParticleSet>& refp_list,
const RefVector<const std::vector<PosType>>& deltaV_list,
const RefVector<const std::vector<RealType>>& deltaS_list,
const RefVector<const NLPPJob<RealType>>& joblist,
bool sphere)
{
auto& vp_leader = vp_list.getLeader();
if (!vp_leader.isSpinor())
throw std::runtime_error(
"VirtualParticleSet::mw_makeMovesWithSpin should not be called if particle sets aren't spionor types");
vp_leader.onSphere = sphere;
vp_leader.refPS = refp_list.getLeader();

const size_t nVPs = countVPs(vp_list);
auto& mw_refPctls = vp_leader.getMultiWalkerRefPctls();
mw_refPctls.resize(nVPs);

RefVectorWithLeader<ParticleSet> p_list(vp_leader);
p_list.reserve(vp_list.size());

size_t ivp = 0;
for (int iw = 0; iw < vp_list.size(); iw++)
{
VirtualParticleSet& vp(vp_list[iw]);
const std::vector<PosType>& deltaV(deltaV_list[iw]);
const std::vector<RealType>& deltaS(deltaS_list[iw]);
const NLPPJob<RealType>& job(joblist[iw]);

vp.onSphere = sphere;
vp.refPS = refp_list[iw];
vp.refPtcl = job.electron_id;
vp.refSourcePtcl = job.ion_id;
assert(vp.R.size() == deltaV.size());
assert(vp.spins.size() == deltaS.size());
assert(vp.R.size() == vp.spins.size());
for (size_t k = 0; k < vp.R.size(); k++, ivp++)
{
vp.R[k] = refp_list[iw].R[vp.refPtcl] + deltaV[k];
vp.spins[k] = refp_list[iw].spins[vp.refPtcl] + deltaS[k];
mw_refPctls[ivp] = vp.refPtcl;
}
p_list.push_back(vp);
}
assert(ivp == nVPs);

mw_refPctls.updateTo();
ParticleSet::mw_update(p_list);
}

} // namespace qmcplusplus
13 changes: 12 additions & 1 deletion src/Particle/VirtualParticleSet.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -88,7 +88,11 @@ class VirtualParticleSet : public ParticleSet
* @param sphere set true if VP are on a sphere around the reference source particle
* @param iat reference source particle
*/
void makeMoves(const ParticleSet& refp, int jel, const std::vector<PosType>& deltaV, bool sphere = false, int iat = -1);
void makeMoves(const ParticleSet& refp,
int jel,
const std::vector<PosType>& deltaV,
bool sphere = false,
int iat = -1);

/** move virtual particles to new postions and update distance tables
* @param refp reference particle set
Expand All @@ -111,6 +115,13 @@ class VirtualParticleSet : public ParticleSet
const RefVector<const NLPPJob<RealType>>& joblist,
bool sphere);

static void mw_makeMovesWithSpin(const RefVectorWithLeader<VirtualParticleSet>& vp_list,
const RefVectorWithLeader<ParticleSet>& p_list,
const RefVector<const std::vector<PosType>>& deltaV_list,
const RefVector<const std::vector<RealType>>& deltaS_list,
const RefVector<const NLPPJob<RealType>>& joblist,
bool sphere);

static RefVectorWithLeader<ParticleSet> RefVectorWithLeaderParticleSet(
const RefVectorWithLeader<VirtualParticleSet>& vp_list)
{
Expand Down
6 changes: 3 additions & 3 deletions src/QMCHamiltonians/ECPComponentBuilder.2.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -135,7 +135,7 @@ void ECPComponentBuilder::buildSemiLocalAndLocal(std::vector<xmlNodePtr>& semiPt
// we may not know which one is local yet.

std::vector<int> angList;
std::vector<int> angListSO; //For spin-orbit, if it exists
std::vector<int> angListSO; //For spin-orbit, if it exists
std::vector<xmlNodePtr> vpsPtr;
std::vector<xmlNodePtr> vpsoPtr; //For spin-orbit, if it exists.
Lmax = -1;
Expand Down Expand Up @@ -364,8 +364,8 @@ void ECPComponentBuilder::buildSO(const std::vector<int>& angList,
app->spline();
pp_so->add(angList[l], app);
}
NumSO = angList.size();
pp_so->Rmax_ = rmax;
NumSO = angList.size();
pp_so->setRmax(rmax);
}

bool ECPComponentBuilder::parseCasino(const std::string& fname, xmlNodePtr cur)
Expand Down
121 changes: 101 additions & 20 deletions src/QMCHamiltonians/SOECPComponent.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,36 +15,37 @@
#include "SOECPComponent.h"
#include "Numerics/Ylm.h"
#include "CPU/BLAS.hpp"
#include "NLPPJob.h"

namespace qmcplusplus
{
SOECPComponent::SOECPComponent()
: lmax_(0), nchannel_(0), nknot_(0), sknot_(0), total_knots_(0), Rmax_(-1), VP_(nullptr)
: lmax_(0), nchannel_(0), nknot_(0), sknot_(0), total_knots_(0), rmax_(-1), vp_(nullptr)
{}

SOECPComponent::~SOECPComponent()
{
for (int i = 0; i < sopp_m_.size(); i++)
delete sopp_m_[i];
if (VP_)
delete VP_;
if (vp_)
delete vp_;
}

void SOECPComponent::print(std::ostream& os) {}

void SOECPComponent::initVirtualParticle(const ParticleSet& qp)
{
assert(VP_ == nullptr);
assert(vp_ == nullptr);
outputManager.pause();
VP_ = new VirtualParticleSet(qp, total_knots_);
vp_ = new VirtualParticleSet(qp, total_knots_);
outputManager.resume();
}

void SOECPComponent::deleteVirtualParticle()
{
if (VP_)
delete VP_;
VP_ = nullptr;
if (vp_)
delete vp_;
vp_ = nullptr;
}

void SOECPComponent::add(int l, RadialPotentialType* pp)
Expand All @@ -58,8 +59,8 @@ SOECPComponent* SOECPComponent::makeClone(const ParticleSet& qp)
SOECPComponent* myclone = new SOECPComponent(*this);
for (int i = 0; i < sopp_m_.size(); i++)
myclone->sopp_m_[i] = sopp_m_[i]->makeClone();
if (VP_)
myclone->VP_ = new VirtualParticleSet(qp, total_knots_);
if (vp_)
myclone->vp_ = new VirtualParticleSet(qp, total_knots_);
return myclone;
}

Expand Down Expand Up @@ -141,17 +142,13 @@ SOECPComponent::RealType SOECPComponent::evaluateOne(ParticleSet& W,
RealType r,
const PosType& dr)
{

for (int ip = 0; ip < nchannel_; ip++)
vrad_[ip] = sopp_m_[ip]->splint(r);

RealType sold = W.spins[iel];
buildTotalQuadrature(r, dr, sold);

if (VP_)
if (vp_)
{
VP_->makeMovesWithSpin(W, iel, deltaV_, deltaS_, true, iat);
Psi.evaluateRatios(*VP_, psiratio_);
vp_->makeMovesWithSpin(W, iel, deltaV_, deltaS_, true, iat);
Psi.evaluateRatios(*vp_, psiratio_);
}
else
for (int iq = 0; iq < total_knots_; iq++)
Expand All @@ -162,6 +159,11 @@ SOECPComponent::RealType SOECPComponent::evaluateOne(ParticleSet& W,
Psi.resetPhaseDiff();
}

return calculateProjector(r, dr, sold);
}

SOECPComponent::RealType SOECPComponent::calculateProjector(RealType r, const PosType& dr, RealType sold)
{
ComplexType pairpot;
for (int iq = 0; iq < total_knots_; iq++)
{
Expand Down Expand Up @@ -190,6 +192,81 @@ SOECPComponent::RealType SOECPComponent::evaluateOne(ParticleSet& W,
return std::real(pairpot);
}

void SOECPComponent::mw_evaluateOne(const RefVectorWithLeader<SOECPComponent>& soecp_component_list,
const RefVectorWithLeader<ParticleSet>& p_list,
const RefVectorWithLeader<TrialWaveFunction>& psi_list,
const RefVector<const NLPPJob<RealType>>& joblist,
std::vector<RealType>& pairpots,
ResourceCollection& collection)
{
auto& soecp_component_leader = soecp_component_list.getLeader();
if (soecp_component_leader.vp_)
{
// Compute ratios with VP
RefVectorWithLeader<VirtualParticleSet> vp_list(*soecp_component_leader.vp_);
RefVectorWithLeader<const VirtualParticleSet> const_vp_list(*soecp_component_leader.vp_);
RefVector<const std::vector<PosType>> deltaV_list;
RefVector<const std::vector<RealType>> deltaS_list;
RefVector<std::vector<ValueType>> psiratios_list;
vp_list.reserve(soecp_component_list.size());
const_vp_list.reserve(soecp_component_list.size());
deltaV_list.reserve(soecp_component_list.size());
deltaS_list.reserve(soecp_component_list.size());
psiratios_list.reserve(soecp_component_list.size());

for (size_t i = 0; i < soecp_component_list.size(); i++)
{
SOECPComponent& component(soecp_component_list[i]);
const NLPPJob<RealType>& job = joblist[i];
const RealType sold = p_list[i].spins[job.electron_id];

component.buildTotalQuadrature(job.ion_elec_dist, job.ion_elec_displ, sold);

vp_list.push_back(*component.vp_);
const_vp_list.push_back(*component.vp_);
deltaV_list.push_back(component.deltaV_);
deltaS_list.push_back(component.deltaS_);
psiratios_list.push_back(component.psiratio_);
}

ResourceCollectionTeamLock<VirtualParticleSet> vp_res_lock(collection, vp_list);

VirtualParticleSet::mw_makeMovesWithSpin(vp_list, p_list, deltaV_list, deltaS_list, joblist, true);

TrialWaveFunction::mw_evaluateRatios(psi_list, const_vp_list, psiratios_list);
}
else
{
// Compute ratios without VP. Slow
for (size_t i = 0; i < p_list.size(); i++)
{
SOECPComponent& component(soecp_component_list[i]);
ParticleSet& W(p_list[i]);
TrialWaveFunction& psi(psi_list[i]);
const NLPPJob<RealType>& job = joblist[i];

const RealType sold = W.spins[job.electron_id];
component.buildTotalQuadrature(job.ion_elec_dist, job.ion_elec_displ, sold);

for (int j = 0; j < component.total_knots_; j++)
{
W.makeMoveWithSpin(job.electron_id, component.deltaV_[j], component.deltaS_[j]);
component.psiratio_[j] = psi.calcRatio(W, job.electron_id);
W.rejectMove(job.electron_id);
psi.resetPhaseDiff();
}
}
}

for (size_t i = 0; i < p_list.size(); i++)
{
SOECPComponent& component(soecp_component_list[i]);
const NLPPJob<RealType>& job = joblist[i];
const RealType sold = p_list[i].spins[job.electron_id];
pairpots[i] = component.calculateProjector(job.ion_elec_dist, job.ion_elec_displ, sold);
}
}

SOECPComponent::RealType SOECPComponent::evaluateValueAndDerivatives(ParticleSet& W,
int iat,
TrialWaveFunction& Psi,
Expand All @@ -214,11 +291,11 @@ SOECPComponent::RealType SOECPComponent::evaluateValueAndDerivatives(ParticleSet

//Now we have all the spin and spatial quadrature points acculated to use in evaluation
//Now we need to obtain dlogpsi and dlogpsi_vp
if (VP_)
if (vp_)
{
// Compute ratios with VP
VP_->makeMovesWithSpin(W, iel, deltaV_, deltaS_, true, iat);
Psi.evaluateDerivRatios(*VP_, optvars, psiratio_, dratio_);
vp_->makeMovesWithSpin(W, iel, deltaV_, deltaS_, true, iat);
Psi.evaluateDerivRatios(*vp_, optvars, psiratio_, dratio_);
}
else
for (int iq = 0; iq < total_knots_; iq++)
Expand Down Expand Up @@ -317,6 +394,10 @@ void SOECPComponent::buildTotalQuadrature(const RealType r, const PosType& dr, c
addSpatialQuadrature(sknot_, r, dr, smax - sold, RealType(1. / 3.) * dS);

assert(count == total_knots_);

//also set the radial function
for (int ip = 0; ip < nchannel_; ip++)
vrad_[ip] = sopp_m_[ip]->splint(r);
}

void SOECPComponent::rotateQuadratureGrid(const TensorType& rmat)
Expand Down
32 changes: 26 additions & 6 deletions src/QMCHamiltonians/SOECPComponent.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,13 +16,19 @@
#include "QMCHamiltonians/OperatorBase.h"
#include "QMCHamiltonians/RandomRotationMatrix.h"
#include "QMCWaveFunctions/TrialWaveFunction.h"
#include <ResourceCollection.h>
#include "Numerics/OneDimGridBase.h"
#include "Numerics/OneDimGridFunctor.h"
#include "Numerics/OneDimLinearSpline.h"
#include "Numerics/OneDimCubicSpline.h"

namespace qmcplusplus
{

namespace testing
{
class TestSOECPotential;
}
/** class SOECPComponent
** brief Computes the nonlocal spin-orbit interaction \f$\Delta V_SO(r) |ljm_j><ljm_j|\f$.
** details This computes the nonlocal spin-orbit interaction between a single ion species and
Expand All @@ -46,7 +52,7 @@ class SOECPComponent : public QMCTraits
int sknot_;
int total_knots_; //spin + spatial knots
///Maximum cutoff the non-local pseudopotential
RealType Rmax_;
RealType rmax_;
///Angular momentum map
aligned_vector<int> angpp_m_;
///Non-Local part of the pseudo-potential
Expand All @@ -70,9 +76,9 @@ class SOECPComponent : public QMCTraits
//scratch spaces used by evaluateValueAndDerivative
Matrix<ValueType> dratio_;
Vector<ValueType> dlogpsi_vp_;
VirtualParticleSet* VP_;
VirtualParticleSet* vp_;

//This builds the full quadrature grid for the Simpsons rule used for spin integrals as well as
//This builds the full quadrature grid for the Simpsons rule used for spin integrals as well as
//the spatial quadrature. In this function, it specifies the deltaS_ and deltaV_ for all the quadrature points and sets the interal weights
//in spin_quad_weights
//If there are s0,s1,...sN spin integral points and q0,q1,...qM spatial quadrature points, the order is
Expand Down Expand Up @@ -108,6 +114,16 @@ class SOECPComponent : public QMCTraits
*/
RealType evaluateOne(ParticleSet& W, int iat, TrialWaveFunction& Psi, int iel, RealType r, const PosType& dr);

RealType calculateProjector(RealType r, const PosType& dr, RealType sold);


static void mw_evaluateOne(const RefVectorWithLeader<SOECPComponent>& soecp_component_list,
const RefVectorWithLeader<ParticleSet>& p_list,
const RefVectorWithLeader<TrialWaveFunction>& psi_list,
const RefVector<const NLPPJob<RealType>>& joblist,
std::vector<RealType>& pairpots,
ResourceCollection& collection);

RealType evaluateValueAndDerivatives(ParticleSet& P,
int iat,
TrialWaveFunction& psi,
Expand All @@ -123,15 +139,19 @@ class SOECPComponent : public QMCTraits
void initVirtualParticle(const ParticleSet& qp);
void deleteVirtualParticle();

inline void setRmax(RealType rmax) { Rmax_ = rmax; }
inline RealType getRmax() const { return Rmax_; }
inline void setLmax(int Lmax) { lmax_ = Lmax; }
inline void setRmax(RealType rmax) { rmax_ = rmax; }
inline RealType getRmax() const { return rmax_; }
inline void setLmax(int lmax) { lmax_ = lmax; }
inline int getLmax() const { return lmax_; }
inline int getNknot() const { return nknot_; }
inline int getSknot() const { return sknot_; }

const VirtualParticleSet* getVP() const { return vp_; };

friend struct ECPComponentBuilder;
friend void copyGridUnrotatedForTest(SOECPComponent& nlpp);

friend class testing::TestSOECPotential;
};

} // namespace qmcplusplus
Expand Down
Loading

0 comments on commit 98c8f2d

Please sign in to comment.