Skip to content

[WIP] Add cusp correction based on atomic orbitals #4901

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

Draft
wants to merge 6 commits into
base: develop
Choose a base branch
from
Draft

[WIP] Add cusp correction based on atomic orbitals #4901

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
930ba5d
Add atomic-orbital based cusp correction
markdewing Jan 13, 2024
24556d1
Add comment
markdewing Jan 13, 2024
d759aab
Add script for computing cusp parameters
markdewing Jan 16, 2024
60c644a
Correct all the s-type orbitals
markdewing Jan 22, 2024
e0cda53
Fix bugs
markdewing Jan 29, 2024
6e917d2
Add adaptive rc
markdewing Jan 29, 2024
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
132 changes: 123 additions & 9 deletions src/Numerics/GaussianBasisSet.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -85,6 +85,54 @@ struct GaussianCombo
}
};

bool hasShortRangeCusp = false;
real_type src_rcut = 0.0;
real_type src_alpha;
real_type src_a;
real_type src_c;
real_type src_cp;
real_type src_bp;
real_type src_d0;
real_type src_d1;
real_type src_d2;
real_type src_d3;
real_type src_d4;
real_type src_d5;
real_type src_delta = 0.0;

// The short-ranged cusp is intended to be used with transform="yes" and
// so it only implements the function value and first derivative at the
// smallest grid endpoint (which is assumed to not be in the interpolation region)

real_type evalShortRangeCusp(real_type r)
{
//real_type val = src_a * std::exp(-src_alpha * r) + src_c;
real_type val = src_a * std::exp(-src_alpha * r) + src_bp * r + src_cp;
return val;
}

real_type evalShortRangeInterp(real_type r)
{
//real_type val = src_a * std::exp(-src_alpha * r) + src_c;
real_type d0 = src_d0;
real_type d1 = src_d1;
real_type d2 = src_d2;
real_type d3 = src_d3;
real_type d4 = src_d4;
real_type d5 = src_d5;
real_type rc = src_rcut;
real_type x = r - rc;
real_type val = d5 * std::pow(x, 5) + d4 * std::pow(x, 4) + d3 * x * x * x + d2 * x * x + d1 * x + d0;
return val;
}

real_type evalShortRangeCusp_df(real_type r)
{
//real_type dval = -src_alpha * src_a * std::exp(-src_alpha * r);
real_type dval = -src_alpha * src_a * std::exp(-src_alpha * r) + src_bp;
return dval;
}

///Boolean
bool Normalized;
real_type L;
Expand Down Expand Up @@ -113,10 +161,21 @@ struct GaussianCombo
real_type r2 = r * r;
typename std::vector<BasicGaussian>::const_iterator it(gset.begin());
typename std::vector<BasicGaussian>::const_iterator it_end(gset.end());
while (it != it_end)
if (r < src_rcut)
{
res += (*it).f(r2);
++it;
res = evalShortRangeCusp(r);
}
else if (r < src_rcut + src_delta)
{
res = evalShortRangeInterp(r);
}
else
{
while (it != it_end)
{
res += (*it).f(r2);
++it;
}
}
return res;
}
Expand All @@ -126,10 +185,17 @@ struct GaussianCombo
real_type r2 = r * r;
typename std::vector<BasicGaussian>::const_iterator it(gset.begin());
typename std::vector<BasicGaussian>::const_iterator it_end(gset.end());
while (it != it_end)
if (r < src_rcut)
{
res += (*it).df(r, r2);
++it;
res = evalShortRangeCusp_df(r);
}
else
{
while (it != it_end)
{
res += (*it).df(r, r2);
++it;
}
}
return res;
}
Expand All @@ -139,10 +205,17 @@ struct GaussianCombo
Y = 0.0;
real_type rr = r * r;
typename std::vector<BasicGaussian>::iterator it(gset.begin()), it_end(gset.end());
while (it != it_end)
if (r < src_rcut)
{
Y += (*it).f(rr);
++it;
Y = evalShortRangeCusp(r);
}
else
{
while (it != it_end)
{
Y += (*it).f(rr);
++it;
}
}
return Y;
}
Expand Down Expand Up @@ -274,6 +347,8 @@ bool GaussianCombo<T>::putBasisGroupH5(hdf_archive& hin, Communicate& myComm)
if (myComm.rank() == 0)
{
hin.read(NbRadFunc, "NbRadFunc");
hasShortRangeCusp = hin.is_group("shortrangecusp");
app_log() << "Has short range cusp = " << hasShortRangeCusp << std::endl;
hin.push("radfunctions");
}
myComm.bcast(NbRadFunc);
Expand Down Expand Up @@ -308,6 +383,45 @@ bool GaussianCombo<T>::putBasisGroupH5(hdf_archive& hin, Communicate& myComm)
if (myComm.rank() == 0)
hin.pop();

if (hasShortRangeCusp)
{
if (myComm.rank() == 0)
{
hin.push("shortrangecusp");
hin.read(src_rcut, "rcut");
hin.read(src_alpha, "alpha");
hin.read(src_a, "a");
//hin.read(src_c, "c");
hin.read(src_bp, "bp");
hin.read(src_cp, "cp");
hin.read(src_d0, "d0");
hin.read(src_d1, "d1");
hin.read(src_d2, "d2");
hin.read(src_d3, "d3");
hin.read(src_d4, "d4");
hin.read(src_d5, "d5");
hin.read(src_delta, "delta");
app_log() << " short range cusp rcut = " << src_rcut << " alpha = " << src_alpha << " a = " << src_a
<< " src_cp = " << src_cp << " delta = " << src_delta << std::endl;
app_log() << " d0 = " << src_d0 << " d1 = " << src_d1 << " d2 = " << src_d2 << " d3 = " << src_d3
<< " d4 = " << src_d4 << " d5 = " << src_d5 << std::endl;
hin.pop();
}
myComm.bcast(hasShortRangeCusp);
myComm.bcast(src_rcut);
myComm.bcast(src_alpha);
myComm.bcast(src_a);
myComm.bcast(src_bp);
myComm.bcast(src_cp);
myComm.bcast(src_d0);
myComm.bcast(src_d1);
myComm.bcast(src_d2);
myComm.bcast(src_d3);
myComm.bcast(src_d4);
myComm.bcast(src_d5);
myComm.bcast(src_delta);
}

return true;
}
} // namespace qmcplusplus
Expand Down
Loading