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

Changes from 1 commit
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
Next Next commit
Add atomic-orbital based cusp correction 
This is from Manten and Luchow JCP 115 5362 (2001).
The actual parameters need to be computed and added to the HDF file by a
separate script.
  • Loading branch information
@markdewing
markdewing committed Jan 13, 2024
commit 930ba5d6b736b070d67d37b90057cd117e526d47
126 changes: 117 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,50 @@ struct GaussianCombo
}
};

bool hasShortRangeCusp = false;
real_type src_rcut;
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;

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 +157,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 +181,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 +201,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 +343,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 +379,43 @@ 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_d5, "delta");
app_log() << " short range cusp rcut = " << src_rcut << " alpha = " << src_alpha << " a = " << src_a
<< " src_c = " << src_c << 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