diff --git a/src/KSPACE/ewald.cpp b/src/KSPACE/ewald.cpp index f5fb680f51..7a383d8472 100644 --- a/src/KSPACE/ewald.cpp +++ b/src/KSPACE/ewald.cpp @@ -176,6 +176,8 @@ void Ewald::init() double lpr = sqrt(lprx*lprx + lpry*lpry + lprz*lprz) / sqrt(3.0); double spr = 2.0*q2 * exp(-g_ewald*g_ewald*cutoff*cutoff) / sqrt(natoms*cutoff*xprd*yprd*zprd_slab); + double tpr = estimate_table_accuracy(spr); + double accuracy = sqrt(lpr*lpr + spr*spr + tpr*tpr); // stats @@ -183,18 +185,18 @@ void Ewald::init() if (screen) { fprintf(screen," G vector (1/distance) = %g\n",g_ewald); fprintf(screen," estimated absolute RMS force accuracy = %g\n", - MAX(lpr,spr)); + accuracy); fprintf(screen," estimated relative force accuracy = %g\n", - MAX(lpr,spr)/two_charge_force); + accuracy/two_charge_force); fprintf(screen," KSpace vectors: actual max1d max3d = %d %d %d\n", kcount,kmax,kmax3d); } if (logfile) { fprintf(logfile," G vector (1/distnace) = %g\n",g_ewald); fprintf(logfile," estimated absolute RMS force accuracy = %g\n", - MAX(lpr,spr)); + accuracy); fprintf(logfile," estimated relative force accuracy = %g\n", - MAX(lpr,spr)/two_charge_force); + accuracy/two_charge_force); fprintf(logfile," KSpace vectors: actual max1d max3d = %d %d %d\n", kcount,kmax,kmax3d); } diff --git a/src/KSPACE/pppm.cpp b/src/KSPACE/pppm.cpp index 74479551d4..966d3a98a8 100644 --- a/src/KSPACE/pppm.cpp +++ b/src/KSPACE/pppm.cpp @@ -1091,6 +1091,8 @@ void PPPM::set_grid() double lpr = sqrt(lprx*lprx + lpry*lpry + lprz*lprz) / sqrt(3.0); double spr = 2.0*q2 * exp(-g_ewald*g_ewald*cutoff*cutoff) / sqrt(natoms*cutoff*xprd*yprd*zprd_slab); + double tpr = estimate_table_accuracy(spr); + double accuracy = sqrt(lpr*lpr + spr*spr + tpr*tpr); // free local memory @@ -1109,9 +1111,9 @@ void PPPM::set_grid() fprintf(screen," grid = %d %d %d\n",nx_pppm,ny_pppm,nz_pppm); fprintf(screen," stencil order = %d\n",order); fprintf(screen," estimated absolute RMS force accuracy = %g\n", - MAX(lpr,spr)); + accuracy); fprintf(screen," estimated relative force accuracy = %g\n", - MAX(lpr,spr)/two_charge_force); + accuracy/two_charge_force); fprintf(screen," using %s precision FFTs\n",fft_prec); } if (logfile) { @@ -1119,9 +1121,9 @@ void PPPM::set_grid() fprintf(logfile," grid = %d %d %d\n",nx_pppm,ny_pppm,nz_pppm); fprintf(logfile," stencil order = %d\n",order); fprintf(logfile," estimated absolute RMS force accuracy = %g\n", - MAX(lpr,spr)); + accuracy); fprintf(logfile," estimated relative force accuracy = %g\n", - MAX(lpr,spr)/two_charge_force); + accuracy/two_charge_force); fprintf(logfile," using %s precision FFTs\n",fft_prec); } } diff --git a/src/kspace.cpp b/src/kspace.cpp index b0b345bdfc..397c2067c7 100644 --- a/src/kspace.cpp +++ b/src/kspace.cpp @@ -17,6 +17,7 @@ #include "atom.h" #include "comm.h" #include "force.h" +#include "pair.h" #include "memory.h" #include "error.h" #include "suffix.h" @@ -124,6 +125,37 @@ void KSpace::ev_setup(int eflag, int vflag) } } +/* ---------------------------------------------------------------------- + estimate the accuracy of the short-range coulomb tables +------------------------------------------------------------------------- */ + +double KSpace::estimate_table_accuracy(double spr) +{ + double table_accuracy = 0.0; + int nctb = force->pair->ncoultablebits; + if (nctb) { + double empirical_precision[17]; + empirical_precision[6] = 2.12E-04; + empirical_precision[7] = 4.97E-05; + empirical_precision[8] = 1.24E-05; + empirical_precision[9] = 3.04E-06; + empirical_precision[10] = 8.51E-07; + empirical_precision[11] = 1.85E-07; + empirical_precision[12] = 5.87E-08; + empirical_precision[13] = 2.81E-08; + empirical_precision[14] = 2.20E-08; + empirical_precision[15] = 2.13E-08; + empirical_precision[16] = 2.11E-08; + if (nctb <= 6) table_accuracy = empirical_precision[6]; + else if (nctb <= 16) table_accuracy = empirical_precision[nctb]; + else table_accuracy = empirical_precision[16]; + table_accuracy *= two_charge_force; + if (table_accuracy > spr) + error->warning(FLERR,"For better accuracy use 'pair_modify table 0'"); + } + return table_accuracy; +} + /* ---------------------------------------------------------------------- modify parameters of the KSpace style ------------------------------------------------------------------------- */ diff --git a/src/kspace.h b/src/kspace.h index c5ea971d4d..eeb586a287 100644 --- a/src/kspace.h +++ b/src/kspace.h @@ -70,6 +70,7 @@ class KSpace : protected Pointers { int maxeatom,maxvatom; void ev_setup(int, int); + double estimate_table_accuracy(double); }; } @@ -93,4 +94,9 @@ W: Kspace_modify slab param < 2.0 may cause unphysical behavior The kspace_modify slab parameter should be larger to insure periodic grids padded with empty space do not overlap. +W: For better accuracy use 'pair_modify table 0' + +The user-specified force accuracy cannot be achieved unless the table +feature is disabled by using 'pair_modify table 0'. + */ diff --git a/src/pair.h b/src/pair.h index 2c67e9bdb0..e854041d22 100644 --- a/src/pair.h +++ b/src/pair.h @@ -57,6 +57,8 @@ class Pair : protected Pointers { int evflag; // energy,virial settings int eflag_either,eflag_global,eflag_atom; int vflag_either,vflag_global,vflag_atom; + + int ncoultablebits; // size of Coulomb table int nextra; // # of extra quantities pair style calculates double *pvector; // vector of extra pair quantities @@ -146,7 +148,6 @@ class Pair : protected Pointers { // pair_modify settings int offset_flag,mix_flag; // flags for offset and mixing - int ncoultablebits; // size of Coulomb table double tabinner; // inner cutoff for Coulomb table // custom data type for accessing Coulomb tables