jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

avoid division by zero in ewald for empty and uncharged systems. require kspace_modify gewald

This commit is contained in:
Axel Kohlmeyer 2018-01-05 20:14:25 -05:00
parent 91993b236d
commit 6c058fb56c
2 changed files with 8 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
src/KSPACE/ewald.cpp
View File

@ -154,6 +154,8 @@ void Ewald::init()
if (!gewaldflag) { if (!gewaldflag) {
if (accuracy <= 0.0) if (accuracy <= 0.0)
error->all(FLERR,"KSpace accuracy must be > 0"); error->all(FLERR,"KSpace accuracy must be > 0");
if (q2 == 0.0)
error->all(FLERR,"Must use 'kspace_modify gewald' for uncharged system");
g_ewald = accuracy*sqrt(natoms*cutoff*xprd*yprd*zprd) / (2.0*q2); g_ewald = accuracy*sqrt(natoms*cutoff*xprd*yprd*zprd) / (2.0*q2);
if (g_ewald >= 1.0) g_ewald = (1.35 - 0.15*log(accuracy))/cutoff; if (g_ewald >= 1.0) g_ewald = (1.35 - 0.15*log(accuracy))/cutoff;
else g_ewald = sqrt(-log(g_ewald)) / cutoff; else g_ewald = sqrt(-log(g_ewald)) / cutoff;
@ -340,6 +342,7 @@ void Ewald::setup()
double Ewald::rms(int km, double prd, bigint natoms, double q2) double Ewald::rms(int km, double prd, bigint natoms, double q2)
{ {
if (natoms == 0) natoms = 1; // avoid division by zero
double value = 2.0*q2*g_ewald/prd * double value = 2.0*q2*g_ewald/prd *
sqrt(1.0/(MY_PI*km*natoms)) * sqrt(1.0/(MY_PI*km*natoms)) *
exp(-MY_PI*MY_PI*km*km/(g_ewald*g_ewald*prd*prd)); exp(-MY_PI*MY_PI*km*km/(g_ewald*g_ewald*prd*prd));

5
src/KSPACE/ewald_disp.cpp
View File

@ -193,6 +193,10 @@ void EwaldDisp::init()
if (!gewaldflag) { if (!gewaldflag) {
if (function[0]) { if (function[0]) {
if (accuracy <= 0.0)
error->all(FLERR,"KSpace accuracy must be > 0");
if (q2 == 0.0)
error->all(FLERR,"Must use 'kspace_modify gewald' for uncharged system");
g_ewald = accuracy*sqrt(natoms*(*cutoff)*shape_det(domain->h)) / (2.0*q2); g_ewald = accuracy*sqrt(natoms*(*cutoff)*shape_det(domain->h)) / (2.0*q2);
if (g_ewald >= 1.0) g_ewald = (1.35 - 0.15*log(accuracy))/(*cutoff); if (g_ewald >= 1.0) g_ewald = (1.35 - 0.15*log(accuracy))/(*cutoff);
else g_ewald = sqrt(-log(g_ewald)) / (*cutoff); else g_ewald = sqrt(-log(g_ewald)) / (*cutoff);
@ -298,6 +302,7 @@ double EwaldDisp::rms(int km, double prd, bigint natoms,
double q2, double b2, double M2) double q2, double b2, double M2)
{ {
double value = 0.0; double value = 0.0;
if (natoms == 0) natoms = 1; // avoid division by zero
// Coulombic // Coulombic