Merge pull request #2775 from akohlmey/tip4p-coulomb-warn

Print warning when a tip4p pair style may cause incorrect results
This commit is contained in:
Axel Kohlmeyer 2021-05-19 20:06:05 -04:00 committed by GitHub
commit b2641a4836
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
3 changed files with 38 additions and 0 deletions

View File

@ -145,6 +145,22 @@ specified since a Coulombic cutoff cannot be specified for an individual I,J
type pair. All type pairs use the same global Coulombic cutoff specified in
the pair_style command.
.. warning::
Because of how these pair styles implement the coulomb interactions
by implicitly defining a fourth site for the negative charge
of the TIP4P and similar water models, special care must be taken
when using these pair styles with other computations that also use
charges. Unless they are specially set up to also handle the implicit
definition of the 4th site, results are likely incorrect. Example:
:doc:`compute dipole/chunk <compute_dipole_chunk>`. For the same
reason, when using one of these pair styles with
:doc:`pair_style hybrid <pair_hybrid>`, **all** coulomb interactions
should be handled by a single sub-style with TIP4P support. All other
instances and styles will "see" the M point charges at the position
of the Oxygen atom and thus compute incorrect forces and energies.
LAMMPS will print a warning when it detects one of these issues.
----------
A version of these styles with a soft core, *lj/cut/tip4p/long/soft*\ , suitable

View File

@ -15,9 +15,11 @@
#include "compute_dipole_chunk.h"
#include "atom.h"
#include "comm.h"
#include "compute_chunk_atom.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "math_special.h"
#include "memory.h"
#include "modify.h"
@ -95,6 +97,10 @@ void ComputeDipoleChunk::init()
cchunk = (ComputeChunkAtom *) modify->compute[icompute];
if (strcmp(cchunk->style,"chunk/atom") != 0)
error->all(FLERR,"Compute dipole/chunk does not use chunk/atom compute");
if ((force->pair_match("/tip4p/",0) != nullptr) && (comm->me == 0))
error->warning(FLERR,"Computed dipole moments may be incorrect when "
"using a tip4p pair style");
}
/* ---------------------------------------------------------------------- */

View File

@ -341,6 +341,8 @@ void PairHybrid::settings(int narg, char **arg)
// multiple[i] = 1 to M if sub-style used multiple times, else 0
int num_tip4p = 0, num_coul = 0; // count sub-styles with tip4p and coulomb
for (int i = 0; i < nstyles; i++) {
int count = 0;
for (int j = 0; j < nstyles; j++) {
@ -348,8 +350,22 @@ void PairHybrid::settings(int narg, char **arg)
if (j == i) multiple[i] = count;
}
if (count == 1) multiple[i] = 0;
if (utils::strmatch(keywords[i],"/tip4p/")) ++num_tip4p;
if (utils::strmatch(keywords[i],"/coul/")
|| utils::strmatch(keywords[i],"^comb")
|| utils::strmatch(keywords[i],"^reax/c")) ++num_coul;
}
if ((num_tip4p > 1) && (comm->me == 0))
error->warning(FLERR,"Using multiple tip4p sub-styles can result in "
"inconsistent calculation of coulomb interactions");
if ((num_tip4p > 0) && (num_coul > 0) && (comm->me == 0))
error->warning(FLERR,"Using a tip4p sub-style with other sub-styles "
"that include coulomb interactions can result in "
"inconsistent calculation of the coulomb interactions");
// set pair flags from sub-style flags
flags();