Skip to content

Commit

Permalink
Replace 1 nel det with 2 nel/2 dets.
Browse files Browse the repository at this point in the history
  • Loading branch information
@ye-luo
ye-luo committed Jan 9, 2018
1 parent 66e8ff5 commit 026e8f7
Show file tree
Hide file tree
Showing 5 changed files with 84 additions and 34 deletions.
13 changes: 8 additions & 5 deletions src/QMCWaveFunctions/Determinant.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -178,7 +178,8 @@ void checkDiff(const MT1 &a, const MT2 &b, const std::string &tag)

struct DiracDeterminant : public WaveFunctionComponentBase
{
DiracDeterminant(int nels, RandomGenerator<RealType> RNG)
DiracDeterminant(int nels, RandomGenerator<RealType> RNG, int First=0)
: FirstIndex(First)
{
psiMinv.resize(nels, nels);
psiV.resize(nels);
Expand Down Expand Up @@ -222,7 +223,7 @@ struct DiracDeterminant : public WaveFunctionComponentBase

GradType evalGrad(ParticleSet &P, int iat) {}

ValueType ratioGrad(ParticleSet &P, int iat, GradType &grad) {}
ValueType ratioGrad(ParticleSet &P, int iat, GradType &grad) { return ratio(P, iat); }

void evaluateGL(ParticleSet &P, ParticleSet::ParticleGradient_t &G,
ParticleSet::ParticleLaplacian_t &L, bool fromscratch = false) {}
Expand All @@ -246,16 +247,16 @@ struct DiracDeterminant : public WaveFunctionComponentBase
constexpr double czero(0);
for (int j = 0; j < nels; ++j)
psiV[j] = myRandom() - shift;
curRatio = inner_product_n(psiV.data(), psiMinv[iel], nels, czero);
curRatio = inner_product_n(psiV.data(), psiMinv[iel-FirstIndex], nels, czero);
return curRatio;
}

/** accept the row and update the inverse */
inline void acceptMove(ParticleSet &P, int iel)
{
const int nels = psiV.size();
updateRow(psiMinv.data(), psiV.data(), nels, nels, iel, curRatio);
std::copy_n(psiV.data(), nels, psiMsave[iel]);
updateRow(psiMinv.data(), psiV.data(), nels, nels, iel-FirstIndex, curRatio);
std::copy_n(psiV.data(), nels, psiMsave[iel-FirstIndex]);
}

/** accessor functions for checking */
Expand All @@ -269,6 +270,8 @@ struct DiracDeterminant : public WaveFunctionComponentBase
double curRatio;
/// workspace size
int LWork;
/// initial particle index
const int FirstIndex;
/// inverse matrix to be update
Matrix<RealType> psiMinv;
/// a SPO set for the row update
Expand Down
13 changes: 8 additions & 5 deletions src/QMCWaveFunctions/DeterminantRef.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -182,7 +182,8 @@ struct DiracDeterminantRef : public qmcplusplus::WaveFunctionComponentBase
{
using ParticleSet = qmcplusplus::ParticleSet;

DiracDeterminantRef(int nels, qmcplusplus::RandomGenerator<RealType> RNG)
DiracDeterminantRef(int nels, qmcplusplus::RandomGenerator<RealType> RNG, int First=0)
: FirstIndex(First)
{
psiMinv.resize(nels, nels);
psiV.resize(nels);
Expand Down Expand Up @@ -225,7 +226,7 @@ struct DiracDeterminantRef : public qmcplusplus::WaveFunctionComponentBase
// FIXME do we want remainder of evaluateLog?
}
GradType evalGrad(ParticleSet &P, int iat) {}
ValueType ratioGrad(ParticleSet &P, int iat, GradType &grad) {}
ValueType ratioGrad(ParticleSet &P, int iat, GradType &grad) { return ratio(P, iat); }
void evaluateGL(ParticleSet &P, ParticleSet::ParticleGradient_t &G,
ParticleSet::ParticleLaplacian_t &L, bool fromscratch = false) {}

Expand All @@ -248,16 +249,16 @@ struct DiracDeterminantRef : public qmcplusplus::WaveFunctionComponentBase
constexpr double czero(0);
for (int j = 0; j < nels; ++j)
psiV[j] = myRandom() - shift;
curRatio = inner_product_n(psiV.data(), psiMinv[iel], nels, czero);
curRatio = inner_product_n(psiV.data(), psiMinv[iel-FirstIndex], nels, czero);
return curRatio;
}

/** accept the row and update the inverse */
inline void acceptMove(ParticleSet &P, int iel)
{
const int nels = psiV.size();
updateRow(psiMinv.data(), psiV.data(), nels, nels, iel, curRatio);
std::copy_n(psiV.data(), nels, psiMsave[iel]);
updateRow(psiMinv.data(), psiV.data(), nels, nels, iel-FirstIndex, curRatio);
std::copy_n(psiV.data(), nels, psiMsave[iel-FirstIndex]);
}

/** accessor functions for checking */
Expand All @@ -271,6 +272,8 @@ struct DiracDeterminantRef : public qmcplusplus::WaveFunctionComponentBase
double curRatio;
/// workspace size
int LWork;
/// initial particle index
const int FirstIndex;
/// inverse matrix to be update
qmcplusplus::Matrix<RealType> psiMinv;
/// a SPO set for the row update
Expand Down
42 changes: 31 additions & 11 deletions src/QMCWaveFunctions/WaveFunction.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,9 @@ WaveFunction::WaveFunction(ParticleSet &ions, ParticleSet &els,
d_ie = DistanceTable::add(ions, els, DT_SOA);

// determinant component
Det = new DetType(els.getTotalNum(), RNG);
nelup = els.getTotalNum()/2;
Det_up = new DetType(nelup, RNG, 0);
Det_dn = new DetType(els.getTotalNum()-nelup, RNG, nelup);

// J1 component
J1 = new J1OrbType(ions, els);
Expand All @@ -48,7 +50,14 @@ WaveFunction::WaveFunction(ParticleSet &ions, ParticleSet &els,
buildJeeI(*J3, els.Lattice.WignerSeitzRadius);
}

WaveFunction::~WaveFunction() { delete J2; }
WaveFunction::~WaveFunction()
{
delete Det_up;
delete Det_dn;
delete J1;
delete J2;
delete J3;
}

void WaveFunction::evaluateLog(ParticleSet &P)
{
Expand All @@ -57,7 +66,8 @@ void WaveFunction::evaluateLog(ParticleSet &P)
{
P.G = czero;
P.L = czero;
LogValue = Det->evaluateLog(P, P.G, P.L);
LogValue = Det_up->evaluateLog(P, P.G, P.L);
LogValue += Det_dn->evaluateLog(P, P.G, P.L);
LogValue += J1->evaluateLog(P, P.G, P.L);
LogValue += J2->evaluateLog(P, P.G, P.L);
LogValue += J3->evaluateLog(P, P.G, P.L);
Expand All @@ -67,26 +77,35 @@ void WaveFunction::evaluateLog(ParticleSet &P)

WaveFunctionBase::posT WaveFunction::evalGrad(ParticleSet &P, int iat)
{
return Det->evalGrad(P, iat) + J1->evalGrad(P, iat) + J2->evalGrad(P, iat) +
J3->evalGrad(P, iat);
return ( iat<nelup ? Det_up->evalGrad(P, iat) : Det_dn->evalGrad(P, iat) )
+ J1->evalGrad(P, iat)
+ J2->evalGrad(P, iat)
+ J3->evalGrad(P, iat);
}

WaveFunctionBase::valT WaveFunction::ratioGrad(ParticleSet &P, int iat,
posT &grad)
{
return Det->ratioGrad(P, iat, grad) + J1->ratioGrad(P, iat, grad) +
J2->ratioGrad(P, iat, grad) + J3->ratioGrad(P, iat, grad);
return ( iat<nelup ? Det_up->ratioGrad(P, iat, grad) : Det_dn->ratioGrad(P, iat, grad) )
* J1->ratioGrad(P, iat, grad)
* J2->ratioGrad(P, iat, grad)
* J3->ratioGrad(P, iat, grad);
}

WaveFunctionBase::valT WaveFunction::ratio(ParticleSet &P, int iat)
{
return Det->ratio(P, iat) * J1->ratio(P, iat) * J2->ratio(P, iat) *
J3->ratio(P, iat);
return ( iat<nelup ? Det_up->ratio(P, iat) : Det_dn->ratio(P, iat) )
* J1->ratio(P, iat)
* J2->ratio(P, iat)
* J3->ratio(P, iat);
}

void WaveFunction::acceptMove(ParticleSet &P, int iat)
{
Det->acceptMove(P, iat);
if(iat<nelup)
Det_up->acceptMove(P, iat);
else
Det_dn->acceptMove(P, iat);
J1->acceptMove(P, iat);
J2->acceptMove(P, iat);
J3->acceptMove(P, iat);
Expand All @@ -99,7 +118,8 @@ void WaveFunction::evaluateGL(ParticleSet &P)
constexpr valT czero(0);
P.G = czero;
P.L = czero;
Det->evaluateGL(P, P.G, P.L);
Det_up->evaluateGL(P, P.G, P.L);
Det_dn->evaluateGL(P, P.G, P.L);
J1->evaluateGL(P, P.G, P.L);
J2->evaluateGL(P, P.G, P.L);
J3->evaluateGL(P, P.G, P.L);
Expand Down
8 changes: 6 additions & 2 deletions src/QMCWaveFunctions/WaveFunction.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,7 +73,9 @@ struct WaveFunction : public WaveFunctionBase
J1OrbType *J1;
J2OrbType *J2;
J3OrbType *J3;
DetType *Det;
int nelup;
DetType *Det_up;
DetType *Det_dn;

WaveFunction(ParticleSet &ions, ParticleSet &els,
RandomGenerator<RealType> RNG);
Expand Down Expand Up @@ -106,7 +108,9 @@ struct WaveFunctionRef : public qmcplusplus::WaveFunctionBase
J1OrbType *J1;
J2OrbType *J2;
J3OrbType *J3;
DetType *Det;
int nelup;
DetType *Det_up;
DetType *Det_dn;
PooledData<valT> Buffer;

WaveFunctionRef(ParticleSet &ions, ParticleSet &els,
Expand Down
42 changes: 31 additions & 11 deletions src/QMCWaveFunctions/WaveFunctionRef.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,9 @@ WaveFunctionRef::WaveFunctionRef(ParticleSet &ions, ParticleSet &els,
d_ie = DistanceTable::add(ions, els, DT_SOA);

// determinant component
Det = new DetType(els.getTotalNum(), RNG);
nelup = els.getTotalNum()/2;
Det_up = new DetType(nelup, RNG, 0);
Det_dn = new DetType(els.getTotalNum()-nelup, RNG, nelup);

// J1 component
J1 = new J1OrbType(ions, els);
Expand All @@ -54,7 +56,14 @@ WaveFunctionRef::WaveFunctionRef(ParticleSet &ions, ParticleSet &els,
buildJeeI(*J3, els.Lattice.WignerSeitzRadius);
}

WaveFunctionRef::~WaveFunctionRef() { delete J2; }
WaveFunctionRef::~WaveFunctionRef()
{
delete Det_up;
delete Det_dn;
delete J1;
delete J2;
delete J3;
}

void WaveFunctionRef::evaluateLog(ParticleSet &P)
{
Expand All @@ -63,7 +72,8 @@ void WaveFunctionRef::evaluateLog(ParticleSet &P)
{
P.G = czero;
P.L = czero;
LogValue = Det->evaluateLog(P, P.G, P.L);
LogValue = Det_up->evaluateLog(P, P.G, P.L);
LogValue += Det_dn->evaluateLog(P, P.G, P.L);
LogValue += J1->evaluateLog(P, P.G, P.L);
LogValue += J2->evaluateLog(P, P.G, P.L);
LogValue += J3->evaluateLog(P, P.G, P.L);
Expand All @@ -73,26 +83,35 @@ void WaveFunctionRef::evaluateLog(ParticleSet &P)

WaveFunctionBase::posT WaveFunctionRef::evalGrad(ParticleSet &P, int iat)
{
return Det->evalGrad(P, iat) + J1->evalGrad(P, iat) + J2->evalGrad(P, iat) +
J3->evalGrad(P, iat);
return ( iat<nelup ? Det_up->evalGrad(P, iat) : Det_dn->evalGrad(P, iat) )
+ J1->evalGrad(P, iat)
+ J2->evalGrad(P, iat)
+ J3->evalGrad(P, iat);
}

WaveFunctionBase::valT WaveFunctionRef::ratioGrad(ParticleSet &P, int iat,
posT &grad)
{
return Det->ratioGrad(P, iat, grad) + J1->ratioGrad(P, iat, grad) +
J2->ratioGrad(P, iat, grad) + J3->ratioGrad(P, iat, grad);
return ( iat<nelup ? Det_up->ratioGrad(P, iat, grad) : Det_dn->ratioGrad(P, iat, grad) )
* J1->ratioGrad(P, iat, grad)
* J2->ratioGrad(P, iat, grad)
* J3->ratioGrad(P, iat, grad);
}

WaveFunctionBase::valT WaveFunctionRef::ratio(ParticleSet &P, int iat)
{
return Det->ratio(P, iat) * J1->ratio(P, iat) * J2->ratio(P, iat) *
J3->ratio(P, iat);
return ( iat<nelup ? Det_up->ratio(P, iat) : Det_dn->ratio(P, iat) )
* J1->ratio(P, iat)
* J2->ratio(P, iat)
* J3->ratio(P, iat);
}

void WaveFunctionRef::acceptMove(ParticleSet &P, int iat)
{
Det->acceptMove(P, iat);
if(iat<nelup)
Det_up->acceptMove(P, iat);
else
Det_dn->acceptMove(P, iat);
J1->acceptMove(P, iat);
J2->acceptMove(P, iat);
J3->acceptMove(P, iat);
Expand All @@ -105,7 +124,8 @@ void WaveFunctionRef::evaluateGL(ParticleSet &P)
constexpr valT czero(0);
P.G = czero;
P.L = czero;
Det->evaluateGL(P, P.G, P.L);
Det_up->evaluateGL(P, P.G, P.L);
Det_dn->evaluateGL(P, P.G, P.L);
J1->evaluateGL(P, P.G, P.L);
J2->evaluateGL(P, P.G, P.L);
J3->evaluateGL(P, P.G, P.L);
Expand Down

0 comments on commit 026e8f7

Please sign in to comment.