Skip to content

Commit

Permalink
Merge pull request QMCPACK#4935 from markdewing/remove_H2
Browse files Browse the repository at this point in the history 
Remove GEVMethod and H2 from optimizer
  • Loading branch information
@prckent
prckent committed Mar 8, 2024
2 parents ab366c0 + 46f7413 commit 7ca0808
Show file tree
Hide file tree
Showing 6 changed files with 43 additions and 87 deletions.
50 changes: 16 additions & 34 deletions src/QMCDrivers/WFOpt/QMCCostFunction.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -110,7 +110,7 @@ void QMCCostFunction::GradCost(std::vector<Return_rt>& PGradient,
ltz = false;
Return_rt delE = std::pow(std::abs(eloc_new - EtargetEff), PowerE);
Return_rt ddelE = PowerE * std::pow(std::abs(eloc_new - EtargetEff), PowerE - 1);
const Return_t* Dsaved = (*DerivRecords[ip])[iw];
const Return_t* Dsaved = (*DerivRecords[ip])[iw];
const Return_rt* HDsaved = (*HDerivRecords[ip])[iw];
for (int pm = 0; pm < NumOptimizables; pm++)
{
Expand All @@ -137,7 +137,7 @@ void QMCCostFunction::GradCost(std::vector<Return_rt>& PGradient,
Return_rt eloc_new = saved[ENERGY_NEW];
Return_rt delta_l = (eloc_new - curAvg_w);
Return_rt sigma_l = delta_l * delta_l;
const Return_t* Dsaved = (*DerivRecords[ip])[iw];
const Return_t* Dsaved = (*DerivRecords[ip])[iw];
const Return_rt* HDsaved = (*HDerivRecords[ip])[iw];
for (int pm = 0; pm < NumOptimizables; pm++)
{
Expand Down Expand Up @@ -643,28 +643,15 @@ QMCCostFunction::Return_rt QMCCostFunction::fillOverlapHamiltonianMatrices(Matri
{
ScopedTimer tmp_timer(fill_timer_);

RealType b1, b2;
if (GEVType == "H2")
{
b1 = w_beta;
b2 = 0;
}
else
{
b2 = w_beta;
b1 = 0;
}
RealType b2(w_beta);

Right = 0.0;
Left = 0.0;

// resetPsi();
curAvg_w = SumValue[SUM_E_WGT] / SumValue[SUM_WGT];
Return_rt curAvg2_w = SumValue[SUM_ESQ_WGT] / SumValue[SUM_WGT];
// RealType H2_avg = 1.0/curAvg2_w;
RealType H2_avg = 1.0 / (curAvg_w * curAvg_w);
// RealType H2_avg = 1.0/std::sqrt(curAvg_w*curAvg_w*curAvg2_w);
RealType V_avg = curAvg2_w - curAvg_w * curAvg_w;
RealType V_avg = curAvg2_w - curAvg_w * curAvg_w;
std::vector<Return_t> D_avg(getNumParams(), 0.0);
Return_rt wgtinv = 1.0 / SumValue[SUM_WGT];
for (int ip = 0; ip < NumThreads; ip++)
Expand All @@ -674,7 +661,7 @@ QMCCostFunction::Return_rt QMCCostFunction::fillOverlapHamiltonianMatrices(Matri
{
const Return_rt* restrict saved = (*RecordsOnNode[ip])[iw];
Return_rt weight = saved[REWEIGHT] * wgtinv;
const Return_t* Dsaved = (*DerivRecords[ip])[iw];
const Return_t* Dsaved = (*DerivRecords[ip])[iw];
for (int pm = 0; pm < getNumParams(); pm++)
{
D_avg[pm] += Dsaved[pm] * weight;
Expand All @@ -692,19 +679,16 @@ QMCCostFunction::Return_rt QMCCostFunction::fillOverlapHamiltonianMatrices(Matri
const Return_rt* restrict saved = (*RecordsOnNode[ip])[iw];
Return_rt weight = saved[REWEIGHT] * wgtinv;
Return_rt eloc_new = saved[ENERGY_NEW];
const Return_t* Dsaved = (*DerivRecords[ip])[iw];
const Return_t* Dsaved = (*DerivRecords[ip])[iw];
const Return_rt* HDsaved = (*HDerivRecords[ip])[iw];
#pragma omp parallel for
for (int pm = 0; pm < getNumParams(); pm++)
{
Return_t wfe = (HDsaved[pm] + (Dsaved[pm] - D_avg[pm]) * eloc_new) * weight;
Return_t wfd = (Dsaved[pm] - D_avg[pm]) * weight;
Return_t vterm =
HDsaved[pm] * (eloc_new - curAvg_w) + (Dsaved[pm] - D_avg[pm]) * eloc_new * (eloc_new - RealType(2.0) * curAvg_w);
Return_t wfe = (HDsaved[pm] + (Dsaved[pm] - D_avg[pm]) * eloc_new) * weight;
Return_t wfd = (Dsaved[pm] - D_avg[pm]) * weight;
Return_t vterm = HDsaved[pm] * (eloc_new - curAvg_w) +
(Dsaved[pm] - D_avg[pm]) * eloc_new * (eloc_new - RealType(2.0) * curAvg_w);
// Return_t vterm = (HDsaved[pm]+(Dsaved[pm]-D_avg[pm])*eloc_new -curAvg_w)*(eloc_new-curAvg_w);
// H2
Right(0, pm + 1) += b1 * H2_avg * std::real(vterm) * weight;
Right(pm + 1, 0) += b1 * H2_avg * std::real(vterm) * weight;
// Variance
Left(0, pm + 1) += b2 * std::real(vterm) * weight;
Left(pm + 1, 0) += b2 * std::real(vterm) * weight;
Expand All @@ -714,28 +698,26 @@ QMCCostFunction::Return_rt QMCCostFunction::fillOverlapHamiltonianMatrices(Matri
for (int pm2 = 0; pm2 < getNumParams(); pm2++)
{
// Hamiltonian
Left(pm + 1, pm2 + 1) += std::real((1 - b2) * std::conj(wfd) * (HDsaved[pm2] + (Dsaved[pm2] - D_avg[pm2]) * eloc_new));
Left(pm + 1, pm2 + 1) +=
std::real((1 - b2) * std::conj(wfd) * (HDsaved[pm2] + (Dsaved[pm2] - D_avg[pm2]) * eloc_new));
// Overlap
RealType ovlij = std::real(std::conj(wfd) * (Dsaved[pm2] - D_avg[pm2]));
Right(pm + 1, pm2 + 1) += ovlij;
// Variance
RealType varij = weight * std::real((HDsaved[pm] - RealType(2.0) * std::conj(Dsaved[pm] - D_avg[pm]) * eloc_new) *
(HDsaved[pm2] - RealType(2.0) * (Dsaved[pm2] - D_avg[pm2]) * eloc_new));
RealType varij = weight *
std::real((HDsaved[pm] - RealType(2.0) * std::conj(Dsaved[pm] - D_avg[pm]) * eloc_new) *
(HDsaved[pm2] - RealType(2.0) * (Dsaved[pm2] - D_avg[pm2]) * eloc_new));
// RealType varij=weight*(HDsaved[pm] +(Dsaved[pm]-D_avg[pm])*eloc_new-curAvg_w)*
// (HDsaved[pm2] + (Dsaved[pm2]-D_avg[pm2])*eloc_new-curAvg_w);
Left(pm + 1, pm2 + 1) += b2 * (varij + V_avg * ovlij);
// H2
Right(pm + 1, pm2 + 1) += b1 * H2_avg * varij;
}
}
}
}
myComm->allreduce(Right);
myComm->allreduce(Left);
Left(0, 0) = (1 - b2) * curAvg_w + b2 * V_avg;
Right(0, 0) = 1.0 + b1 * H2_avg * V_avg;
if (GEVType == "H2")
return H2_avg;
Right(0, 0) = 1.0;

return 1.0;
}
Expand Down
4 changes: 2 additions & 2 deletions src/QMCDrivers/WFOpt/QMCCostFunctionBase.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,6 @@ QMCCostFunctionBase::QMCCostFunctionBase(ParticleSet& w, TrialWaveFunction& psi,
w_w(0.0),
MaxWeight(1e6),
w_beta(0.0),
GEVType("mixed"),
vmc_or_dmc(2.0),
needGrads(true),
targetExcitedStr("no"),
Expand Down Expand Up @@ -318,6 +317,7 @@ bool QMCCostFunctionBase::put(xmlNodePtr q)
std::string includeNonlocalH;
std::string writeXmlPerStep("no");
std::string computeNLPPderiv;
std::string GEVType;
astring variational_subset_str;
ParameterSet m_param;
m_param.add(writeXmlPerStep, "dumpXML");
Expand All @@ -326,7 +326,7 @@ bool QMCCostFunctionBase::put(xmlNodePtr q)
m_param.add(includeNonlocalH, "nonlocalpp", {}, TagStatus::DEPRECATED);
m_param.add(computeNLPPderiv, "use_nonlocalpp_deriv", {}, TagStatus::DEPRECATED);
m_param.add(w_beta, "beta");
m_param.add(GEVType, "GEVMethod");
m_param.add(GEVType, "GEVMethod", {}, TagStatus::DEPRECATED);
m_param.add(targetExcitedStr, "targetExcited");
m_param.add(omega_shift, "omega");
m_param.add(do_override_output, "output_vp_override", {true});
Expand Down
1 change: 0 additions & 1 deletion src/QMCDrivers/WFOpt/QMCCostFunctionBase.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -221,7 +221,6 @@ class QMCCostFunctionBase : public CostFunctionBase<QMCTraits::RealType>, public
Return_rt curVar_abs;

Return_rt w_beta;
std::string GEVType;
Return_rt vmc_or_dmc;
bool needGrads;
///whether we are targeting an excited state
Expand Down
69 changes: 24 additions & 45 deletions src/QMCDrivers/WFOpt/QMCCostFunctionBatched.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -112,12 +112,12 @@ void QMCCostFunctionBatched::GradCost(std::vector<Return_rt>& PGradient,
ltz = false;
Return_rt delE = std::pow(std::abs(eloc_new - EtargetEff), PowerE);
Return_rt ddelE = PowerE * std::pow(std::abs(eloc_new - EtargetEff), PowerE - 1);
const Return_t* Dsaved = DerivRecords_[iw];
const Return_t* Dsaved = DerivRecords_[iw];
const Return_rt* HDsaved = HDerivRecords_[iw];
for (int pm = 0; pm < NumOptimizables; pm++)
{
//From Toulouse J. Chem. Phys. 126, 084102 (2007), this is H_0j+H_j0, which are independent
//estimates of 1/2 the energy gradient g. So g1+g2 is an estimate of g.
//estimates of 1/2 the energy gradient g. So g1+g2 is an estimate of g.
EDtotals_w[pm] += weight * (HDsaved[pm] + 2.0 * std::real(Dsaved[pm]) * delta_l);
URV[pm] += 2.0 * (eloc_new * HDsaved[pm] - curAvg * HD_avg[pm]);
if (ltz)
Expand All @@ -139,7 +139,7 @@ void QMCCostFunctionBatched::GradCost(std::vector<Return_rt>& PGradient,
Return_rt eloc_new = saved[ENERGY_NEW];
Return_rt delta_l = (eloc_new - curAvg_w);
Return_rt sigma_l = delta_l * delta_l;
const Return_t* Dsaved = DerivRecords_[iw];
const Return_t* Dsaved = DerivRecords_[iw];
const Return_rt* HDsaved = HDerivRecords_[iw];
for (int pm = 0; pm < NumOptimizables; pm++)
{
Expand Down Expand Up @@ -193,7 +193,7 @@ void QMCCostFunctionBatched::getConfigurations(const std::string& aroot)
}
}

/** Compute number of batches and final batch size given the number of samples
/** Compute number of batches and final batch size given the number of samples
* and a batch size.
* \param[in] sample_size number of samples to process.
* \param[in] batch_size process samples in batch_size at a time (typically the number of walkers in a crowd).
Expand Down Expand Up @@ -354,7 +354,7 @@ void QMCCostFunctionBatched::checkConfigurations(EngineHandle& handle)
for (int j = 0; j < nparams; j++)
{
//dlogpsi is in general complex if psi is complex.
DerivRecords[is][j] = dlogpsi_array[ib][j];
DerivRecords[is][j] = dlogpsi_array[ib][j];
//but E_L and d E_L/dc are real if c is real.
HDerivRecords[is][j] = std::real(dhpsioverpsi_array[ib][j]);
}
Expand Down Expand Up @@ -594,7 +594,7 @@ QMCCostFunctionBatched::EffectiveWeight QMCCostFunctionBatched::correlatedSampli
if (optVars.recompute(j))
{
//In general, dlogpsi is complex.
DerivRecords[is][j] = dlogpsi_array[ib][j];
DerivRecords[is][j] = dlogpsi_array[ib][j];
//However, E_L is always real, and so d E_L/dc is real, provided c is real.
HDerivRecords[is][j] = std::real(dhpsioverpsi_array[ib][j]);
}
Expand Down Expand Up @@ -700,36 +700,20 @@ QMCCostFunctionBatched::Return_rt QMCCostFunctionBatched::fillOverlapHamiltonian
{
ScopedTimer tmp_timer(fill_timer_);

RealType b1, b2;
if (GEVType == "H2")
{
b1 = w_beta;
b2 = 0;
}
else
{
b2 = w_beta;
b1 = 0;
}

Right = 0.0;
Left = 0.0;

// resetPsi();
curAvg_w = SumValue[SUM_E_WGT] / SumValue[SUM_WGT];
Return_rt curAvg2_w = SumValue[SUM_ESQ_WGT] / SumValue[SUM_WGT];
// RealType H2_avg = 1.0/curAvg2_w;
RealType H2_avg = 1.0 / (curAvg_w * curAvg_w);
// RealType H2_avg = 1.0/std::sqrt(curAvg_w*curAvg_w*curAvg2_w);
RealType V_avg = curAvg2_w - curAvg_w * curAvg_w;
RealType V_avg = curAvg2_w - curAvg_w * curAvg_w;
std::vector<Return_t> D_avg(getNumParams(), 0.0);
Return_rt wgtinv = 1.0 / SumValue[SUM_WGT];

for (int iw = 0; iw < rank_local_num_samples_; iw++)
{
const Return_rt* restrict saved = RecordsOnNode_[iw];
Return_rt weight = saved[REWEIGHT] * wgtinv;
const Return_t* Dsaved = DerivRecords_[iw];
const Return_t* Dsaved = DerivRecords_[iw];
for (int pm = 0; pm < getNumParams(); pm++)
{
D_avg[pm] += Dsaved[pm] * weight;
Expand All @@ -743,7 +727,7 @@ QMCCostFunctionBatched::Return_rt QMCCostFunctionBatched::fillOverlapHamiltonian
const Return_rt* restrict saved = RecordsOnNode_[iw];
Return_rt weight = saved[REWEIGHT] * wgtinv;
Return_rt eloc_new = saved[ENERGY_NEW];
const Return_t* Dsaved = DerivRecords_[iw];
const Return_t* Dsaved = DerivRecords_[iw];
const Return_rt* HDsaved = HDerivRecords_[iw];

size_t opt_num_crowds = walkers_per_crowd_.size();
Expand All @@ -752,21 +736,18 @@ QMCCostFunctionBatched::Return_rt QMCCostFunctionBatched::fillOverlapHamiltonian


auto constructMatrices = [](int crowd_id, std::vector<int>& crowd_ranges, int numParams, const Return_t* Dsaved,
const Return_rt* HDsaved, Return_rt weight, Return_rt eloc_new, RealType H2_avg,
RealType V_avg, std::vector<Return_t>& D_avg, RealType b1, RealType b2,
RealType curAvg_w, Matrix<Return_rt>& Left, Matrix<Return_rt>& Right) {
const Return_rt* HDsaved, Return_rt weight, Return_rt eloc_new, RealType V_avg,
std::vector<Return_t>& D_avg, RealType b2, RealType curAvg_w, Matrix<Return_rt>& Left,
Matrix<Return_rt>& Right) {
int local_pm_start = crowd_ranges[crowd_id];
int local_pm_end = crowd_ranges[crowd_id + 1];

for (int pm = local_pm_start; pm < local_pm_end; pm++)
{
Return_t wfe = (HDsaved[pm] + (Dsaved[pm] - D_avg[pm]) * eloc_new) * weight;
Return_t wfd = (Dsaved[pm] - D_avg[pm]) * weight;
Return_t vterm =
HDsaved[pm] * (eloc_new - curAvg_w) + (Dsaved[pm] - D_avg[pm]) * eloc_new * (eloc_new - RealType(2.0) * curAvg_w);
// H2
Right(0, pm + 1) += b1 * H2_avg * std::real(vterm) * weight;
Right(pm + 1, 0) += b1 * H2_avg * std::real(vterm) * weight;
Return_t wfe = (HDsaved[pm] + (Dsaved[pm] - D_avg[pm]) * eloc_new) * weight;
Return_t wfd = (Dsaved[pm] - D_avg[pm]) * weight;
Return_t vterm = HDsaved[pm] * (eloc_new - curAvg_w) +
(Dsaved[pm] - D_avg[pm]) * eloc_new * (eloc_new - RealType(2.0) * curAvg_w);
// Variance
Left(0, pm + 1) += b2 * std::real(vterm) * weight;
Left(pm + 1, 0) += b2 * std::real(vterm) * weight;
Expand All @@ -776,30 +757,28 @@ QMCCostFunctionBatched::Return_rt QMCCostFunctionBatched::fillOverlapHamiltonian
for (int pm2 = 0; pm2 < numParams; pm2++)
{
// Hamiltonian
Left(pm + 1, pm2 + 1) += std::real((1 - b2) * std::conj(wfd) * (HDsaved[pm2] + (Dsaved[pm2] - D_avg[pm2]) * eloc_new));
Left(pm + 1, pm2 + 1) +=
std::real((1 - b2) * std::conj(wfd) * (HDsaved[pm2] + (Dsaved[pm2] - D_avg[pm2]) * eloc_new));
// Overlap
RealType ovlij = std::real(std::conj(wfd) * (Dsaved[pm2] - D_avg[pm2]));
Right(pm + 1, pm2 + 1) += ovlij;
// Variance
RealType varij = weight * std::real((HDsaved[pm] - RealType(2.0) * std::conj(Dsaved[pm] - D_avg[pm]) * eloc_new) *
(HDsaved[pm2] - RealType(2.0) * (Dsaved[pm2] - D_avg[pm2]) * eloc_new));
RealType varij = weight *
std::real((HDsaved[pm] - RealType(2.0) * std::conj(Dsaved[pm] - D_avg[pm]) * eloc_new) *
(HDsaved[pm2] - RealType(2.0) * (Dsaved[pm2] - D_avg[pm2]) * eloc_new));
Left(pm + 1, pm2 + 1) += b2 * (varij + V_avg * ovlij);
// H2
Right(pm + 1, pm2 + 1) += b1 * H2_avg * varij;
}
}
};

ParallelExecutor<> crowd_tasks;
crowd_tasks(opt_num_crowds, constructMatrices, params_per_crowd, getNumParams(), Dsaved, HDsaved, weight, eloc_new,
H2_avg, V_avg, D_avg, b1, b2, curAvg_w, Left, Right);
V_avg, D_avg, w_beta, curAvg_w, Left, Right);
}
myComm->allreduce(Right);
myComm->allreduce(Left);
Left(0, 0) = (1 - b2) * curAvg_w + b2 * V_avg;
Right(0, 0) = 1.0 + b1 * H2_avg * V_avg;
if (GEVType == "H2")
return H2_avg;
Left(0, 0) = (1 - w_beta) * curAvg_w + w_beta * V_avg;
Right(0, 0) = 1.0;

return 1.0;
}
Expand Down
4 changes: 0 additions & 4 deletions src/QMCDrivers/WFOpt/QMCFixedSampleLinearOptimize.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,9 +59,7 @@ QMCFixedSampleLinearOptimize::QMCFixedSampleLinearOptimize(const ProjectData& pr
bigChange(50),
exp0(-16),
stepsize(0.25),
GEVtype("mixed"),
StabilizerMethod("best"),
GEVSplit("no"),
bestShift_i(-1.0),
bestShift_s(-1.0),
shift_i_input(0.01),
Expand Down Expand Up @@ -176,8 +174,6 @@ QMCFixedSampleLinearOptimize::QMCFixedSampleLinearOptimize(const ProjectData& pr
// m_param.add(quadstep,"quadstep");
// m_param.add(stepsize,"stepsize");
// m_param.add(exp1,"exp1");
// m_param.add(GEVtype,"GEVMethod");
// m_param.add(GEVSplit,"GEVSplit");
// m_param.add(StabilizerMethod,"StabilizerMethod");
// m_param.add(LambdaMax,"LambdaMax");
//Set parameters for line minimization:
Expand Down
2 changes: 1 addition & 1 deletion src/QMCDrivers/WFOpt/QMCFixedSampleLinearOptimize.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -136,7 +136,7 @@ class QMCFixedSampleLinearOptimize : public QMCDriver, public LinearMethod, priv

int nstabilizers;
RealType stabilizerScale, bigChange, exp0, exp1, stepsize, savedQuadstep;
std::string GEVtype, StabilizerMethod, GEVSplit;
std::string StabilizerMethod;
RealType w_beta;
/// number of previous steps to orthogonalize to.
int eigCG;
Expand Down

0 comments on commit 7ca0808

Please sign in to comment.