Skip to content

Commit

Permalink
Merge pull request QMCPACK#4686 from ye-luo/reduce-app_error
Browse files Browse the repository at this point in the history 
Correct app_error usage
  • Loading branch information
@PDoakORNL
PDoakORNL authored Jul 25, 2023
2 parents aea5dda + 84208e9 commit 19c5834
Show file tree
Hide file tree
Showing 8 changed files with 63 additions and 138 deletions.
4 changes: 2 additions & 2 deletions src/Concurrency/ParallelExecutorOPENMP.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,9 +19,9 @@

#include <stdexcept>
#include <string>
#include <iostream>

#include "Concurrency/ParallelExecutor.hpp"
#include "Platforms/Host/OutputManager.h"
#include "Concurrency/OpenMP.h"

namespace qmcplusplus
Expand Down Expand Up @@ -54,7 +54,7 @@ void ParallelExecutor<Executor::OPENMP>::operator()(int num_tasks, F&& f, Args&&
++nested_throw_count;
else
{
app_error() << re.what() << std::flush;
std::cerr << re.what() << std::flush;
++throw_count;
}
}
Expand Down
2 changes: 1 addition & 1 deletion src/QMCApp/QMCMain.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -243,7 +243,7 @@ bool QMCMain::execute()
if (qmc_common.dryrun)
{
app_log() << " dryrun == 1 Ignore qmc/loop elements " << std::endl;
APP_ABORT("QMCMain::execute");
myComm->barrier_and_abort("QMCMain::execute");
}
t3.stop();
Timer t1;
Expand Down
14 changes: 11 additions & 3 deletions src/QMCHamiltonians/CoulombPotentialFactory.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,15 @@ void HamiltonianFactory::addMPCPotential(xmlNodePtr cur, bool isphysical)
app_summary() << " Name: " << title << " Physical : " << physical << std::endl;
app_summary() << std::endl;

std::unique_ptr<MPC> mpc = std::make_unique<MPC>(targetPtcl, cutoff);
if (targetPtcl.Density_G.size() == 0)
myComm->barrier_and_abort("HamiltonianFactory::addMPCPotential\n"
"************************\n"
"** Error in MPC setup **\n"
"************************\n"
" The electron density was not setup by the "
"wave function builder.\n");

auto mpc = std::make_unique<MPC>(targetPtcl, cutoff);
targetH->addOperator(std::move(mpc), "MPC", isphysical);
#else
APP_ABORT(
Expand Down Expand Up @@ -251,8 +259,8 @@ void HamiltonianFactory::addPseudoPotential(xmlNodePtr cur)
{
if (psiPool.empty())
return;
app_error() << " Cannot find " << wfname << " in the Wavefunction pool. Using the first wavefunction."
<< std::endl;
app_warning() << " Cannot find " << wfname << " in the Wavefunction pool. Using the first wavefunction."
<< std::endl;
psi = psiPool.begin()->second.get();
}
else
Expand Down
68 changes: 1 addition & 67 deletions src/QMCHamiltonians/ECPComponentBuilder.1.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,10 +45,6 @@ void ECPComponentBuilder::addSemiLocal(xmlNodePtr cur)
{
pp_nonloc->add(l, createVrWithBasisGroup(cur1, grid_semilocal.get()));
}
//else if(cname1 == "data")
//{
// pp_nonloc->add(l,createVrWithData(cur1,grid_semilocal));
//}
cur1 = cur1->next;
}
NumNonLocal++;
Expand Down Expand Up @@ -143,12 +139,6 @@ void ECPComponentBuilder::buildLocal(xmlNodePtr cur)
vr.putBasisGroup(cur, vPowerCorrection);
bareCoulomb = false;
}
else if (cname == "data")
{
pp_loc = std::unique_ptr<RadialPotentialType>(createVrWithData(cur, grid_local_inp.get(), vPowerCorrection));
app_log() << " Local pseduopotential in a <data/>" << std::endl;
return;
}
cur = cur->next;
}
if (grid_local_inp == nullptr)
Expand Down Expand Up @@ -193,9 +183,7 @@ void ECPComponentBuilder::buildLocal(xmlNodePtr cur)
--last;
}
if (last == 0)
{
app_error() << " Illegal Local Pseudopotential " << std::endl;
}
myComm->barrier_and_abort("ECPComponentBuilder::buildLocal. Illegal Local Pseudopotential");
//Add the reset values here
int ng = static_cast<int>(r / 1e-3) + 1;
app_log() << " Use a Linear Grid: [0," << r << "] Number of points = " << ng << std::endl;
Expand Down Expand Up @@ -298,58 +286,4 @@ std::unique_ptr<ECPComponentBuilder::mGridType> ECPComponentBuilder::createGrid(
return agrid;
}

/** Disable pseudo/semilocal/vps/data */
ECPComponentBuilder::RadialPotentialType* ECPComponentBuilder::createVrWithData(xmlNodePtr cur,
mGridType* agrid,
int rCorrection)
{
return nullptr;
// RealType rcIn = agrid->rmax();
// //use the maximum value of the grid
// if(RcutMax<0) RcutMax=rcIn;
// //create a new linear grid if the input grid is not good enough
// GridType *newgrid=0;
// if(agrid->GridTag != LINEAR_1DGRID || RcutMax < rcIn)
// {
// const RealType delta=1000.; // use 1/000
// newgrid = new LinearGrid<RealType>;
// newgrid->set(0.0,RcutMax,static_cast<int>(RcutMax*delta)+1);
// }
// //read the numerical data
// std::vector<RealType> pdata;
// putContent(pdata,cur);
// if(pdata.size() != agrid->size())
// {
// app_error() << " ECPComponentBuilder::createVrWithData vsp/data size does not match." << std::endl;
// abort(); //FIXABORT
// }
// if(rCorrection == 1)
// {
// for(int i=0; i<agrid->size(); i++) pdata[i] *= (*agrid)[i];
// }
// if(newgrid)
// {
// OneDimCubicSpline<RealType> inFunc(grid_global,pdata);
// inFunc.spline();
// int ng=newgrid->size();
// pdata.resize(ng);
// for(int i=0; i<ng; i++)
// {
// RealType r((*agrid)[i]);
// pdata[i]=inFunc.splint(r);
// }
// if(agrid->rmin()>0.0) pdata[0]=pdata[1];
// RadialPotentialType *app = new RadialPotentialType(newgrid,pdata);
// app->spline();
// return app;
// }
// else
// {//use Radial potential with the input grid
// RadialPotentialType *app = new RadialPotentialType(agrid,pdata);
// app->spline();
// return app;
// }
}


} // namespace qmcplusplus
1 change: 0 additions & 1 deletion src/QMCHamiltonians/ECPComponentBuilder.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,7 +74,6 @@ struct ECPComponentBuilder : public MPIObjectBase, public QMCTraits

std::unique_ptr<mGridType> createGrid(xmlNodePtr cur, bool useLinear = false);
RadialPotentialType* createVrWithBasisGroup(xmlNodePtr cur, mGridType* agrid);
RadialPotentialType* createVrWithData(xmlNodePtr cur, mGridType* agrid, int rCorrection = 0);

void doBreakUp(const std::vector<int>& angList,
const Matrix<mRealType>& vnn,
Expand Down
6 changes: 3 additions & 3 deletions src/QMCHamiltonians/ECPComponentBuilder_L2.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ void ECPComponentBuilder::buildL2(xmlNodePtr cur)
if (grid_global == 0)
{
app_error() << " Global grid needs to be defined." << std::endl;
APP_ABORT("ECPComponentBuilder::buildL2");
myComm->barrier_and_abort("ECPComponentBuilder::buildL2");
}

// read <L2> attributes
Expand Down Expand Up @@ -52,12 +52,12 @@ void ECPComponentBuilder::buildL2(xmlNodePtr cur)
else
{
app_error() << " Unrecognized format \"" << format << "\" in PP file." << std::endl;
APP_ABORT("ECPComponentBuilder::buildL2");
myComm->barrier_and_abort("ECPComponentBuilder::buildL2");
}
if (rcut < 0.0)
{
app_error() << " L2 attribute cutoff is missing or negative. Cutoff must be a positive real number." << std::endl;
APP_ABORT("ECPComponentBuilder::buildL2");
myComm->barrier_and_abort("ECPComponentBuilder::buildL2");
}
int npts = grid_global->size();

Expand Down
87 changes: 36 additions & 51 deletions src/QMCHamiltonians/MPC.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,30 +33,27 @@ namespace qmcplusplus
void MPC::resetTargetParticleSet(ParticleSet& ptcl) {}

MPC::MPC(ParticleSet& ptcl, double cutoff)
: Ecut(cutoff),
d_aa_ID(ptcl.addTable(ptcl, DTModes::NEED_FULL_TABLE_ON_HOST_AFTER_DONEPBYP)),
PtclRef(&ptcl),
FirstTime(true)
: Ecut(cutoff), d_aa_ID(ptcl.addTable(ptcl, DTModes::NEED_FULL_TABLE_ON_HOST_AFTER_DONEPBYP)), FirstTime(true)
{
initBreakup();
initBreakup(ptcl);
}

MPC::~MPC() = default;

void MPC::init_gvecs()
void MPC::init_gvecs(const ParticleSet& ptcl)
{
TinyVector<int, OHMMS_DIM> maxIndex(0);
PosType b[OHMMS_DIM];
for (int j = 0; j < OHMMS_DIM; j++)
b[j] = static_cast<RealType>(2.0 * M_PI) * PtclRef->getLattice().b(j);
int numG1 = PtclRef->Density_G.size();
int numG2 = PtclRef->DensityReducedGvecs.size();
assert(PtclRef->Density_G.size() == PtclRef->DensityReducedGvecs.size());
b[j] = static_cast<RealType>(2.0 * M_PI) * ptcl.getLattice().b(j);
int numG1 = ptcl.Density_G.size();
int numG2 = ptcl.DensityReducedGvecs.size();
assert(ptcl.Density_G.size() == ptcl.DensityReducedGvecs.size());
// Loop through all the G-vectors, and find the largest
// indices in each direction with energy less than the cutoff
for (int iG = 0; iG < PtclRef->DensityReducedGvecs.size(); iG++)
for (int iG = 0; iG < ptcl.DensityReducedGvecs.size(); iG++)
{
TinyVector<int, OHMMS_DIM> gint = PtclRef->DensityReducedGvecs[iG];
TinyVector<int, OHMMS_DIM> gint = ptcl.DensityReducedGvecs[iG];
PosType G = (double)gint[0] * b[0];
for (int j = 1; j < OHMMS_DIM; j++)
G += (double)gint[j] * b[j];
Expand All @@ -66,7 +63,7 @@ void MPC::init_gvecs()
maxIndex[j] = std::max(maxIndex[j], std::abs(gint[j]));
Gvecs.push_back(G);
Gints.push_back(gint);
Rho_G.push_back(PtclRef->Density_G[iG]);
Rho_G.push_back(ptcl.Density_G[iG]);
}
}
for (int idim = 0; idim < OHMMS_DIM; idim++)
Expand All @@ -78,17 +75,17 @@ void MPC::init_gvecs()
}


void MPC::compute_g_G(double& g_0, std::vector<double>& g_G, int N)
void MPC::compute_g_G(const ParticleSet& ptcl, double& g_0, std::vector<double>& g_G, int N)
{
double L = PtclRef->getLattice().WignerSeitzRadius;
double L = ptcl.getLattice().WignerSeitzRadius;
double Linv = 1.0 / L;
double Linv3 = Linv * Linv * Linv;
// create an FFTW plan
Array<std::complex<double>, 3> rBox(N, N, N);
Array<std::complex<double>, 3> GBox(N, N, N);
// app_log() << "Doing " << N << " x " << N << " x " << N << " FFT.\n";
//create BC handler
DTD_BConds<RealType, 3, SUPERCELL_BULK> mybc(PtclRef->getLattice());
DTD_BConds<RealType, 3, SUPERCELL_BULK> mybc(ptcl.getLattice());
// Fill the real-space array with f(r)
double Ninv = 1.0 / (double)N;
TinyVector<RealType, 3> u, r;
Expand All @@ -101,8 +98,8 @@ void MPC::compute_g_G(double& g_0, std::vector<double>& g_G, int N)
for (int iz = 0; iz < N; iz++)
{
u[2] = Ninv * iz;
r = PtclRef->getLattice().toCart(u);
//DTD_BConds<double,3,SUPERCELL_BULK>::apply (PtclRef->getLattice(), r);
r = ptcl.getLattice().toCart(u);
//DTD_BConds<double,3,SUPERCELL_BULK>::apply (ptcl.getLattice(), r);
//double rmag = std::sqrt(dot(r,r));
double rmag = std::sqrt(mybc.apply_bc(r));
if (rmag < L)
Expand Down Expand Up @@ -165,21 +162,21 @@ inline double extrap(int N, TinyVector<double, 2> g_12)
}


void MPC::init_f_G()
void MPC::init_f_G(const ParticleSet& ptcl)
{
int numG = Gints.size();
f_G.resize(numG);
int N = std::max(64, 2 * MaxDim + 1);
std::vector<double> g_G_N(numG), g_G_2N(numG), g_G_4N(numG);
double g_0_N, g_0_2N, g_0_4N;
compute_g_G(g_0_N, g_G_N, 1 * N);
compute_g_G(g_0_2N, g_G_2N, 2 * N);
compute_g_G(g_0_4N, g_G_4N, 4 * N);
compute_g_G(ptcl, g_0_N, g_G_N, 1 * N);
compute_g_G(ptcl, g_0_2N, g_G_2N, 2 * N);
compute_g_G(ptcl, g_0_4N, g_G_4N, 4 * N);
// fprintf (stderr, "g_G_1N[0] = %18.14e\n", g_G_N[0]);
// fprintf (stderr, "g_G_2N[0] = %18.14e\n", g_G_2N[0]);
// fprintf (stderr, "g_G_4N[0] = %18.14e\n", g_G_4N[0]);
double volInv = 1.0 / PtclRef->getLattice().Volume;
double L = PtclRef->getLattice().WignerSeitzRadius;
double volInv = 1.0 / ptcl.getLattice().Volume;
double L = ptcl.getLattice().WignerSeitzRadius;
TinyVector<double, 2> g0_12(g_0_2N, g_0_4N);
TinyVector<double, 3> g0_124(g_0_N, g_0_2N, g_0_4N);
f_0 = extrap(N, g0_124);
Expand Down Expand Up @@ -221,15 +218,15 @@ void MPC::init_f_G()
}


void MPC::init_spline()
void MPC::init_spline(const ParticleSet& ptcl)
{
Array<std::complex<double>, 3> rBox(SplineDim), GBox(SplineDim);
Array<double, 3> splineData(SplineDim);
GBox = std::complex<double>();
Vconst = 0.0;
// Now fill in elements of GBox
const RealType vol = PtclRef->getLattice().Volume;
const RealType volInv = 1.0 / PtclRef->getLattice().Volume;
const RealType vol = ptcl.getLattice().Volume;
const RealType volInv = 1.0 / ptcl.getLattice().Volume;
const RealType halfvol = vol / 2.0;
for (int iG = 0; iG < Gvecs.size(); iG++)
{
Expand Down Expand Up @@ -284,38 +281,29 @@ void MPC::init_spline()
VlongSpline =
std::shared_ptr<UBspline_3d_d>(create_UBspline_3d_d(grid0, grid1, grid2, bc0, bc1, bc2, splineData.data()),
destroy_Bspline);
// grid0.num = PtclRef->Density_r.size(0);
// grid1.num = PtclRef->Density_r.size(1);
// grid2.num = PtclRef->Density_r.size(2);
// grid0.num = ptcl.Density_r.size(0);
// grid1.num = ptcl.Density_r.size(1);
// grid2.num = ptcl.Density_r.size(2);
// DensitySpline = create_UBspline_3d_d (grid0, grid1, grid2, bc0, bc1, bc2,
// PtclRef->Density_r.data());
// ptcl.Density_r.data());
}

void MPC::initBreakup()
void MPC::initBreakup(const ParticleSet& ptcl)
{
NParticles = PtclRef->getTotalNum();
NParticles = ptcl.getTotalNum();
app_log() << "\n === Initializing MPC interaction === " << std::endl;
if (PtclRef->Density_G.size() == 0)
{
app_error() << "************************\n"
<< "** Error in MPC setup **\n"
<< "************************\n"
<< " The electron density was not setup by the "
<< "wave function builder.\n";
abort();
}
init_gvecs();
init_f_G();
init_spline();
init_gvecs(ptcl);
init_f_G(ptcl);
init_spline(ptcl);
// FILE *fout = fopen ("MPC.dat", "w");
// double vol = PtclRef->getLattice().Volume;
// double vol = ptcl.getLattice().Volume;
// PosType r0 (0.0, 0.0, 0.0);
// PosType r1 (10.26499236, 10.26499236, 10.26499236);
// int nPoints=1001;
// for (int i=0; i<nPoints; i++) {
// double s = (double)i/(double)(nPoints-1);
// PosType r = (1.0-s)*r0 + s*r1;
// PosType u = PtclRef->getLattice().toUnit(r);
// PosType u = ptcl.getLattice().toUnit(r);
// double V, rho(0.0);
// eval_UBspline_3d_d (VlongSpline, u[0], u[1], u[2], &V);
// // eval_UBspline_3d_d (DensitySpline, u[0], u[1], u[2], &rho);
Expand All @@ -327,9 +315,7 @@ void MPC::initBreakup()

std::unique_ptr<OperatorBase> MPC::makeClone(ParticleSet& qp, TrialWaveFunction& psi)
{
// return new MPC(qp, Ecut);
std::unique_ptr<MPC> newMPC = std::make_unique<MPC>(*this);
newMPC->resetTargetParticleSet(qp);
auto newMPC = std::make_unique<MPC>(*this);
return newMPC;
}

Expand Down Expand Up @@ -368,7 +354,6 @@ MPC::Return_t MPC::evalLR(ParticleSet& P) const

MPC::Return_t MPC::evaluate(ParticleSet& P)
{
//if (FirstTime || P.tag() == PtclRef->tag())
value_ = evalSR(P) + evalLR(P) + Vconst;
return value_;
}
Expand Down
Loading

0 comments on commit 19c5834

Please sign in to comment.