ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

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

Browse Source
This commit is contained in:
Axel Kohlmeyer
2018-01-05 19:52:51 -05:00
parent 5ecc3ce366
commit 91993b236d
2 changed files with 6 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
src/KOKKOS/pppm_kokkos.cpp
View File

@ -969,6 +969,8 @@ void PPPMKokkos<DeviceType>::set_grid_global()
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;
@ -1179,6 +1181,7 @@ double PPPMKokkos<DeviceType>::final_accuracy()
double yprd = domain->yprd; double yprd = domain->yprd;
double zprd = domain->zprd; double zprd = domain->zprd;
bigint natoms = atomKK->natoms; bigint natoms = atomKK->natoms;
if (natoms == 0) natoms = 1; // avoid division by zero
double df_kspace = compute_df_kspace(); double df_kspace = compute_df_kspace();
double q2_over_sqrt = q2 / sqrt(natoms*cutoff*xprd*yprd*zprd); double q2_over_sqrt = q2 / sqrt(natoms*cutoff*xprd*yprd*zprd);

3
src/KSPACE/pppm.cpp
View File

@ -1004,6 +1004,8 @@ void PPPM::set_grid_global()
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;
@ -1346,6 +1348,7 @@ double PPPM::final_accuracy()
double yprd = domain->yprd; double yprd = domain->yprd;
double zprd = domain->zprd; double zprd = domain->zprd;
bigint natoms = atom->natoms; bigint natoms = atom->natoms;
if (natoms == 0) natoms = 1; // avoid division by zero
double df_kspace = compute_df_kspace(); double df_kspace = compute_df_kspace();
double q2_over_sqrt = q2 / sqrt(natoms*cutoff*xprd*yprd*zprd); double q2_over_sqrt = q2 / sqrt(natoms*cutoff*xprd*yprd*zprd);