Skip to content

[WIP] Anisotropic Ewald code #2182

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 13 commits into
base: develop
Choose a base branch
from
Open

[WIP] Anisotropic Ewald code #2182

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
5a0ff83
add anisotropic ewald code
jtkrogel Jan 3, 2020
5415831
remove commented code
jtkrogel Jan 3, 2020
1c927f0
fix typos in comments
jtkrogel Jan 3, 2020
fd926a2
add namespace, const, and light function docs
jtkrogel Jan 6, 2020
f9db5f4
split off cpp file for EwaldRef
jtkrogel Jan 6, 2020
8b06e0a
comparing aniso and qmcpack ewald
jtkrogel Jan 6, 2020
00ef7f7
merge with develop
jtkrogel Jan 7, 2020
51a77f7
start optimized ewald
jtkrogel Jan 8, 2020
634e8c7
cleanup state machine
jtkrogel Jan 9, 2020
1aa8058
rename functions
jtkrogel Jan 9, 2020
9162af3
add SR and LR error bounds
jtkrogel Jan 13, 2020
3ef40f0
add search to enforce bounds
jtkrogel Jan 13, 2020
d55eeba
enable accuracy vs time comparisons
jtkrogel Jan 13, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
start optimized ewald
  • Loading branch information
@jtkrogel
jtkrogel committed Jan 8, 2020
commit 51a77f73553582f24df767776076e5704228cb5d
51 changes: 40 additions & 11 deletions src/QMCHamiltonians/CoulombPBCAA.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -62,15 +62,18 @@ CoulombPBCAA::CoulombPBCAA(ParticleSet& ref, bool active, bool computeForces)
Q[i] = Zat[i];


ewald.initialize(A,Q);

//RealType qmcpack_kappa = AA->LR_kc;
RealType qmcpack_sigma = std::sqrt(AA->LR_kc / (2.0 * AA->LR_rc));
RealType qmcpack_kappa = 1./(std::sqrt(2.0)*qmcpack_sigma);
ewaldtools::AnisotropicEwald ewald_qp(A,Q,1e-10,qmcpack_kappa);

//ewaldref::IntVec nmax = 0;
//ewald.setNmax(nmax);
//ewald.initialize(A,Q);
ewald.initialize(A,Q,1e-10,qmcpack_kappa);


ewaldref::IntVec nmax = 20;
ewald.setNmax(nmax);

if (!is_active)
{
Expand Down Expand Up @@ -102,9 +105,11 @@ CoulombPBCAA::CoulombPBCAA(ParticleSet& ref, bool active, bool computeForces)
//RealType lr_qp = ewald_qp.ewaldEnergyLRDT(dt);
//RealType sr0_qp = ewald_qp.ewaldEnergySR0DT(dt);

//RealType c_qmc = myConst;
//RealType sr_qmc = evalSR(Ps);
//RealType lr_qmc = evalLR(Ps);
RealType c_qmc = myConst;
RealType sr_qmc = evalSR(Ps);
RealType lr_qmc = evalLR(Ps);
RealType E_qmc = c_qmc + sr_qmc + lr_qmc;


app_log()<<std::setprecision(14);
app_log()<<std::endl;
Expand Down Expand Up @@ -135,12 +140,29 @@ CoulombPBCAA::CoulombPBCAA(ParticleSet& ref, bool active, bool computeForces)
app_log()<<" aniso ewald nmax : "<<ewald_qp.getNmax()<<std::endl;
app_log()<<std::setprecision(14);


//ewaldtools::IntVec nmax = 20;
//ewald.setupOpt(nmax);
ewald.setupOpt();
RealType lr_opt = ewald.ewaldEnergyLROpt(Ps.R);
app_log()<<" ewald LR : "<<lr<<std::endl;

//auto& ewald_opt = ewald_qp;
auto& ewald_opt = ewald;

RealType c_ref = ewald_opt.ewaldEnergyConst();
RealType sr_ref = ewald_opt.ewaldEnergySR(Ps.R);
RealType lr_ref = ewald_opt.ewaldEnergyLR(Ps.R);
//ewald_opt.setupOpt();
ewald_opt.setupOpt(nmax);
RealType c_opt = ewald_opt.ewaldEnergyConst();
RealType sr_opt = ewald_opt.ewaldEnergySROpt(dt);
RealType lr_opt = ewald_opt.ewaldEnergyLROpt(Ps.R);
app_log()<<" ewald opt energy : "<<c_opt+sr_opt+lr_opt<<std::endl;
app_log()<<" ewald opt energy2 : "<<ewald.ewaldEnergyOpt(Ps.R,dt)<<std::endl;
app_log()<<" qmcpack energy : "<<E_qmc<<std::endl;
app_log()<<" ewald LR : "<<lr_ref<<std::endl;
app_log()<<" ewald opt LR : "<<lr_opt<<std::endl;
app_log()<<" ewald SR : "<<sr_ref<<std::endl;
app_log()<<" ewald opt SR : "<<sr_opt<<std::endl;
app_log()<<" qmcpack SR : "<<sr_qmc<<std::endl;

app_log()<<std::endl;
app_log()<<std::endl;
Expand Down Expand Up @@ -175,7 +197,7 @@ CoulombPBCAA::CoulombPBCAA(ParticleSet& ref, bool active, bool computeForces)
app_log() << " Check passed." << std::endl;
}

APP_ABORT("explore")
//APP_ABORT("explore")

}
prefix = "F_AA";
Expand Down Expand Up @@ -247,6 +269,8 @@ CoulombPBCAA::Return_t CoulombPBCAA::evaluate(ParticleSet& P)
{
ScopedTimer evaltimer(TimerManager.createTimer("CoulombPBCAA_eval_total"));

const DistanceTableData& dt(P.getDistTable(d_aa_ID));

if (is_active)
{
{
Expand All @@ -263,12 +287,17 @@ CoulombPBCAA::Return_t CoulombPBCAA::evaluate(ParticleSet& P)
RealType eval;
{
ScopedTimer evaltimer(TimerManager.createTimer("Aniso eval"));
eval = ewald.fixedGridEwaldEnergy(P.R);
//eval = ewald.fixedGridEwaldEnergy(P.R);
eval = ewald.ewaldEnergyOpt(P.R,dt);
}

app_log()<<std::endl;
app_log()<<" QMCPACK value: "<<Value<<std::endl;
app_log()<<" Aniso value: "<<eval<<std::endl;
app_log()<<" QMCPACK SR: "<<evalSR(P)<<std::endl;
app_log()<<" Aniso SR: "<<ewald.ewaldEnergySROpt(dt)<<std::endl;
app_log()<<" QMCPACK LR: "<<evalLR(P)+myConst<<std::endl;
app_log()<<" Aniso LR: "<<ewald.ewaldEnergyLROpt(P.R)<<std::endl;
app_log()<<std::endl;
}
return Value;
Expand Down
86 changes: 51 additions & 35 deletions src/QMCHamiltonians/EwaldTools.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -203,6 +203,24 @@ real_t fixedGridSum(T& function, IntVec nmax, bool zero=true)



/// Functor for short range part of e-e pair potential
class RspaceEwaldPairPotential
{
private:
/// The constant 1/(\sqrt{2}\kappa) in Drummond 2008 formula 6
const real_t rexponent;

public:
RspaceEwaldPairPotential(real_t rexponent_in) : rexponent(rexponent_in) {}

real_t operator()(const RealVec& r) const
{
real_t R = std::sqrt(dot(r, r));
return std::erfc(rexponent * R) / R;
}
};


/// Functor for Fourier transform of long range part of e-e pair potential
class KspaceEwaldPairPotential
{
Expand Down Expand Up @@ -553,8 +571,31 @@ class AnisotropicEwald
setupOpt(nmax_anisotropic);
}


template<typename PA, typename DT>
real_t ewaldEnergyOpt(const PA& R, const DT& dt)
{
real_t Vc = ewaldEnergyConst();
real_t Vsr = ewaldEnergySROpt(dt);
real_t Vlr = ewaldEnergyLROpt(R);
return Vc + Vsr + Vlr;
}


template<typename DT>
real_t ewaldEnergySROpt(const DT& dt)
{
const auto& dr = dt.getDisplacements();
real_t ee = 0.0;
RspaceEwaldPairPotential vsr(rexponent_ewald);
for (size_t i = 0; i < N; ++i)
for (size_t j = 0; j < i; ++j)
ee += Q[i]*Q[j]*vsr(dr[i][j]);
return ee;
}


/// LR (k-space) part of total energy computed optimally
/// LR (k-space) part of total energy computed optimally
template<typename PA>
real_t ewaldEnergyLROpt(const PA& R)
{
Expand Down Expand Up @@ -606,27 +647,23 @@ class AnisotropicEwald





template<typename PA>
real_t ewaldEnergySR0(const PA& R)
template<typename DT>
real_t ewaldEnergySR0(const DT& dt)
{
const auto& dr = dt.getDisplacements();
real_t ee = 0.0;
size_t Npairs = (N*(N-1))/2;
IntVec nmax;
RspaceEwaldPairPotential vsr(rexponent_ewald);
for (size_t i = 0; i < N; ++i)
for (size_t j = 0; j < i; ++j)
{
real_t tol_ewald = tol/(Q[i]*Q[j])/Npairs;
RspaceEwaldTerm rfunc(R[i]-R[j], A, rexponent_ewald);
IntVec iv = 0;
real_t rval = rfunc(iv);
ee += Q[i]*Q[j]*rval;
}
ee += Q[i]*Q[j]*vsr(dr[i][j]);
return ee;
}






template<typename DT>
real_t ewaldEnergySRDT(const DT& dt)
{
Expand Down Expand Up @@ -667,27 +704,6 @@ class AnisotropicEwald



template<typename DT>
real_t ewaldEnergySR0DT(const DT& dt)
{
auto& dr = dt.getDisplacements();
real_t ee = 0.0;
size_t Npairs = (N*(N-1))/2;
IntVec nmax;
for (size_t i = 0; i < N; ++i)
for (size_t j = 0; j < i; ++j)
{
real_t tol_ewald = tol/(Q[i]*Q[j])/Npairs;
RspaceEwaldTerm rfunc(dr[i][j], A, rexponent_ewald);
IntVec iv = 0;
real_t rval = rfunc(iv);
ee += Q[i]*Q[j]*rval;
}
return ee;
}





};
Expand Down