diff --git a/doc/src/pair_lj_cut_tip4p.rst b/doc/src/pair_lj_cut_tip4p.rst index 4f55f3b14e..5f3b4e48e7 100644 --- a/doc/src/pair_lj_cut_tip4p.rst +++ b/doc/src/pair_lj_cut_tip4p.rst @@ -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 `. For the same + reason, when using one of these pair styles with + :doc:`pair_style 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 diff --git a/src/compute_dipole_chunk.cpp b/src/compute_dipole_chunk.cpp index cb0bb8d914..4ca991f90b 100644 --- a/src/compute_dipole_chunk.cpp +++ b/src/compute_dipole_chunk.cpp @@ -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"); } /* ---------------------------------------------------------------------- */ diff --git a/src/pair_hybrid.cpp b/src/pair_hybrid.cpp index ecbb8ffa0c..7c5976e65e 100644 --- a/src/pair_hybrid.cpp +++ b/src/pair_hybrid.cpp @@ -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();