diff --git a/doc/src/Commands_compute.rst b/doc/src/Commands_compute.rst index c2753b2afc..f073212524 100644 --- a/doc/src/Commands_compute.rst +++ b/doc/src/Commands_compute.rst @@ -52,6 +52,8 @@ KOKKOS, o = OPENMP, t = OPT. * :doc:`dilatation/atom ` * :doc:`dipole ` * :doc:`dipole/chunk ` + * :doc:`dipole/tip4p ` + * :doc:`dipole/tip4p/chunk ` * :doc:`displace/atom ` * :doc:`dpd ` * :doc:`dpd/atom ` diff --git a/doc/src/Commands_fix.rst b/doc/src/Commands_fix.rst index 89cf3e880a..13bdc34137 100644 --- a/doc/src/Commands_fix.rst +++ b/doc/src/Commands_fix.rst @@ -70,6 +70,7 @@ OPT. * :doc:`dt/reset (k) ` * :doc:`edpd/source ` * :doc:`efield ` + * :doc:`efield/tip4p ` * :doc:`ehex ` * :doc:`electrode/conp (i) ` * :doc:`electrode/conq (i) ` diff --git a/doc/src/compute.rst b/doc/src/compute.rst index 3d74b3884e..880f60a8a6 100644 --- a/doc/src/compute.rst +++ b/doc/src/compute.rst @@ -206,6 +206,8 @@ The individual style names on the :doc:`Commands compute ` pag * :doc:`dilatation/atom ` - Peridynamic dilatation for each atom * :doc:`dipole ` - dipole vector and total dipole * :doc:`dipole/chunk ` - dipole vector and total dipole for each chunk +* :doc:`dipole/tip4p ` - dipole vector and total dipole with TIP4P pair style +* :doc:`dipole/tip4p/chunk ` - dipole vector and total dipole for each chunk with TIP4P pair style * :doc:`displace/atom ` - displacement of each atom * :doc:`dpd ` - total values of internal conductive energy, internal mechanical energy, chemical energy, and harmonic average of internal temperature * :doc:`dpd/atom ` - per-particle values of internal conductive energy, internal mechanical energy, chemical energy, and internal temperature diff --git a/doc/src/compute_dipole.rst b/doc/src/compute_dipole.rst index 95c5e216f0..fa1ca4ce0b 100644 --- a/doc/src/compute_dipole.rst +++ b/doc/src/compute_dipole.rst @@ -1,6 +1,10 @@ .. index:: compute dipole +.. index:: compute dipole/tip4p compute dipole command +====================== + +compute dipole/tip4p command ============================ Syntax @@ -8,10 +12,10 @@ Syntax .. code-block:: LAMMPS - compute ID group-ID dipole arg + compute ID group-ID style arg * ID, group-ID are documented in :doc:`compute ` command -* dipole = style name of this compute command +* style = *dipole* or *dipole/tip4p* * arg = *mass* or *geometry* = use COM or geometric center for charged chunk correction (optional) Examples @@ -21,6 +25,7 @@ Examples compute 1 fluid dipole compute dw water dipole geometry + compute dw water dipole/tip4p Description """"""""""" @@ -28,13 +33,20 @@ Description Define a computation that calculates the dipole vector and total dipole for a group of atoms. -This compute calculates the x,y,z coordinates of the dipole vector -and the total dipole moment for the atoms in the compute group. -This includes all effects due to atoms passing through periodic boundaries. -For a group with a net charge the resulting dipole is made position independent -by subtracting the position vector of the center of mass or geometric center -times the net charge from the computed dipole vector. Both per-atom charges -and per-atom dipole moments, if present, contribute to the computed dipole. +These computes calculate the x,y,z coordinates of the dipole vector and +the total dipole moment for the atoms in the compute group. This +includes all effects due to atoms passing through periodic boundaries. +For a group with a net charge the resulting dipole is made position +independent by subtracting the position vector of the center of mass or +geometric center times the net charge from the computed dipole +vector. Both per-atom charges and per-atom dipole moments, if present, +contribute to the computed dipole. + +.. versionadded:: TBD + +Compute *dipole/tip4p* includes adjustments for the charge carrying +point M in molecules with TIP4P water geometry. The corresponding +parameters are extracted from the pair style. .. note:: @@ -49,10 +61,10 @@ and per-atom dipole moments, if present, contribute to the computed dipole. Output info """"""""""" -This compute calculations a global scalar containing the magnitude of -the computed dipole moment and a global vector of length 3 with the -dipole vector. See the :doc:`Howto output ` page for -an overview of LAMMPS output options. +These computes calculate a global scalar containing the magnitude of the +computed dipole moment and a global vector of length 3 with the dipole +vector. See the :doc:`Howto output ` page for an overview +of LAMMPS output options. The computed values are "intensive". The array values will be in dipole units (i.e., charge :doc:`units ` times distance @@ -60,7 +72,12 @@ dipole units (i.e., charge :doc:`units ` times distance Restrictions """""""""""" - none + +Compute style *dipole/tip4p* is part of the EXTRA-COMPUTE package. It is +only enabled if LAMMPS was built with that package. See the :doc:`Build +package ` page for more info. + +Compute style *dipole/tip4p* can only be used with tip4p pair styles. Related commands """""""""""""""" diff --git a/doc/src/compute_dipole_chunk.rst b/doc/src/compute_dipole_chunk.rst index 504e6f20d0..4e89bb748a 100644 --- a/doc/src/compute_dipole_chunk.rst +++ b/doc/src/compute_dipole_chunk.rst @@ -1,17 +1,21 @@ .. index:: compute dipole/chunk +.. index:: compute dipole/tip4p/chunk compute dipole/chunk command ============================ +compute dipole/tip4p/chunk command +================================== + Syntax """""" .. code-block:: LAMMPS - compute ID group-ID dipole/chunk chunkID arg + compute ID group-ID style chunkID arg * ID, group-ID are documented in :doc:`compute ` command -* dipole/chunk = style name of this compute command +* style = *dipole/chunk* or *dipole/tip4p/chunk* * chunkID = ID of :doc:`compute chunk/atom ` command * arg = *mass* or *geometry* = use COM or geometric center for charged chunk correction (optional) @@ -38,13 +42,20 @@ or atoms in a spatial bin. See the :doc:`compute chunk/atom details of how chunks can be defined and examples of how they can be used to measure properties of a system. -This compute calculates the :math:`(x,y,z)` coordinates of the dipole vector -and the total dipole moment for each chunk, which includes all effects due -to atoms passing through periodic boundaries. For chunks with a net -charge the resulting dipole is made position independent by subtracting -the position vector of the center of mass or geometric center times the -net charge from the computed dipole vector. Both per-atom charges and -per-atom dipole moments, if present, contribute to the computed dipole. +These computes calculate the :math:`(x,y,z)` coordinates of the dipole +vector and the total dipole moment for each chunk, which includes all +effects due to atoms passing through periodic boundaries. For chunks +with a net charge the resulting dipole is made position independent by +subtracting the position vector of the center of mass or geometric +center times the net charge from the computed dipole vector. Both +per-atom charges and per-atom dipole moments, if present, contribute to +the computed dipole. + +.. versionadded:: TBD + +Compute *dipole/tip4p/chunk* includes adjustments for the charge +carrying point M in molecules with TIP4P water geometry. The +corresponding parameters are extracted from the pair style. Note that only atoms in the specified group contribute to the calculation. The :doc:`compute chunk/atom ` command @@ -78,12 +89,12 @@ command, for example: Output info """"""""""" -This compute calculates a global array where the number of rows = the +These computes calculate a global array where the number of rows = the number of chunks *Nchunk* as calculated by the specified :doc:`compute -chunk/atom ` command. The number of columns is 4 for -the :math:`(x,y,z)` dipole vector components and the total dipole of each -chunk. These values can be accessed by any command that uses global -array values from a compute as input. See the :doc:`Howto output +chunk/atom ` command. The number of columns is 4 +for the :math:`(x,y,z)` dipole vector components and the total dipole of +each chunk. These values can be accessed by any command that uses +global array values from a compute as input. See the :doc:`Howto output ` page for an overview of LAMMPS output options. The array values are "intensive". The array values will be in @@ -92,7 +103,13 @@ dipole units (i.e., charge :doc:`units ` times distance Restrictions """""""""""" - none + +Compute style *dipole/tip4p/chunk* is part of the EXTRA-COMPUTE +package. It is only enabled if LAMMPS was built with that package. See +the :doc:`Build package ` page for more info. + +Compute style *dipole/tip4p/chunk* can only be used with tip4p pair +styles. Related commands """""""""""""""" diff --git a/doc/src/fix.rst b/doc/src/fix.rst index 7d505fec63..9b7ada9369 100644 --- a/doc/src/fix.rst +++ b/doc/src/fix.rst @@ -222,6 +222,7 @@ accelerated styles exist. * :doc:`dt/reset ` - reset the timestep based on velocity, forces * :doc:`edpd/source ` - add heat source to eDPD simulations * :doc:`efield ` - impose electric field on system +* :doc:`efield/tip4p ` - impose electric field on system with TIP4P molecules * :doc:`ehex ` - enhanced heat exchange algorithm * :doc:`electrode/conp ` - impose electric potential * :doc:`electrode/conq ` - impose total electric charge diff --git a/doc/src/fix_efield.rst b/doc/src/fix_efield.rst index a44f8ddb9d..eebfaba8c6 100644 --- a/doc/src/fix_efield.rst +++ b/doc/src/fix_efield.rst @@ -1,17 +1,21 @@ .. index:: fix efield +.. index:: fix efield/tip4p fix efield command ================== +fix efield/tip4p command +======================== + Syntax """""" .. parsed-literal:: - fix ID group-ID efield ex ey ez keyword value ... + fix ID group-ID style ex ey ez keyword value ... * ID, group-ID are documented in :doc:`fix ` command -* efield = style name of this fix command +* style = *efield* or *efield/tip4p* * ex,ey,ez = E-field component values (electric field units) * any of ex,ey,ez can be a variable (see below) * zero or more keyword/value pairs may be appended to args @@ -31,27 +35,36 @@ Examples fix kick external-field efield 1.0 0.0 0.0 fix kick external-field efield 0.0 0.0 v_oscillate + fix kick external-field efield/tip4p 1.0 0.0 0.0 Description """"""""""" -Add a force F = qE to each charged atom in the group due to an +Add a force :math:`\vec{F} = q\vec{E}` to each charged atom in the group due to an external electric field being applied to the system. If the system contains point-dipoles, also add a torque on the dipoles due to the external electric field. -For charges, any of the 3 quantities defining the E-field components -can be specified as an equal-style or atom-style -:doc:`variable `, namely *ex*, *ey*, *ez*\ . If the value is a -variable, it should be specified as v_name, where name is the variable -name. In this case, the variable will be evaluated each timestep, and -its value used to determine the E-field component. +.. versionadded:: TBD + +When the *efield/tip4p* style is used, the E-field will be applied to +the position of the virtual charge site M of a TIP4P molecule instead of +the oxygen position as it is defined by a corresponding :doc:`TIP4P pair +style `. The forces on the M site due to the +external field are projected on the oxygen and hydrogen atoms of the +TIP4P molecules. + +For charges, any of the 3 quantities defining the E-field components can +be specified as an equal-style or atom-style :doc:`variable `, +namely *ex*, *ey*, *ez*\ . If the value is a variable, it should be +specified as v_name, where name is the variable name. In this case, the +variable will be evaluated each timestep, and its value used to +determine the E-field component. For point-dipoles, equal-style variables can be used, but atom-style variables are not currently supported, since they imply a spatial gradient in the electric field which means additional terms with -gradients of the field are required for the force and torque on -dipoles. +gradients of the field are required for the force and torque on dipoles. Equal-style variables can specify formulas with various mathematical functions, and include :doc:`thermo_style ` command @@ -81,10 +94,18 @@ self-consistent minimization problem (see below). The *energy* keyword is not allowed if the added field is a constant vector (ex,ey,ez), with all components defined as numeric constants and not as variables. This is because LAMMPS can compute the energy -for each charged particle directly as E = -x dot qE = -q (x\*ex + y\*ey -+ z\*ez), so that -Grad(E) = F. Similarly for point-dipole particles -the energy can be computed as E = -mu dot E = -(mux\*ex + muy\*ey + -muz\*ez). +for each charged particle directly as + +.. math:: + + U_{efield} = -\vec{x} \cdot q\vec{E} = -q (x\cdot E_x + y\cdot E_y + z\cdot Ez), + +so that :math:`-\nabla U_{efield} = \vec{F}`. Similarly for point-dipole particles +the energy can be computed as + +.. math:: + + U_{efield} = -\vec{\mu} \cdot \vec{E} = -\mu_x\cdot E_x + \mu_y\cdot E_y + \mu_z\cdot E_z The *energy* keyword is optional if the added force is defined with one or more variables, and if you are performing dynamics via the @@ -120,29 +141,28 @@ Restart, fix_modify, output, run start/stop, minimize info No information about this fix is written to :doc:`binary restart files `. -The :doc:`fix_modify ` *energy* option is supported by -this fix to add the potential energy inferred by the added force due -to the electric field to the global potential energy of the system as -part of :doc:`thermodynamic output `. The default -setting for this fix is :doc:`fix_modify energy no `. -Note that this energy is a fictitious quantity but is needed so that -the :doc:`minimize ` command can include the forces added by -this fix in a consistent manner. I.e. there is a decrease in -potential energy when atoms move in the direction of the added force -due to the electric field. +The :doc:`fix_modify ` *energy* option is supported by this +fix to add the potential energy inferred by the added force due to the +electric field to the global potential energy of the system as part of +:doc:`thermodynamic output `. The default setting for +this fix is :doc:`fix_modify energy no `. Note that this +energy is a fictitious quantity but is needed so that the :doc:`minimize +` command can include the forces added by this fix in a +consistent manner. I.e. there is a decrease in potential energy when +atoms move in the direction of the added force due to the electric +field. -The :doc:`fix_modify ` *virial* option is supported by -this fix to add the contribution due to the added forces on atoms to -both the global pressure and per-atom stress of the system via the -:doc:`compute pressure ` and :doc:`compute -stress/atom ` commands. The former can be -accessed by :doc:`thermodynamic output `. The default -setting for this fix is :doc:`fix_modify virial no `. +The :doc:`fix_modify ` *virial* option is supported by this +fix to add the contribution due to the added forces on atoms to both the +global pressure and per-atom stress of the system via the :doc:`compute +pressure ` and :doc:`compute stress/atom +` commands. The former can be accessed by +:doc:`thermodynamic output `. The default setting for +this fix is :doc:`fix_modify virial no `. The :doc:`fix_modify ` *respa* option is supported by this -fix. This allows to set at which level of the :doc:`r-RESPA -` integrator the fix adding its forces. Default is the -outermost level. +fix. This allows to set at which level of the :doc:`r-RESPA ` +integrator the fix adding its forces. Default is the outermost level. This fix computes a global scalar and a global 3-vector of forces, which can be accessed by various :doc:`output commands @@ -169,7 +189,11 @@ the iteration count during the minimization. Restrictions """""""""""" -None +Fix style *efield/tip4p* is part of the EXTRA-FIX package. It is only +enabled if LAMMPS was built with that package. See the :doc:`Build +package ` page for more info. + +Fix style *efield/tip4p* can only be used with tip4p pair styles. Related commands """""""""""""""" diff --git a/src/.gitignore b/src/.gitignore index aa57927669..6eefcb5b34 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -536,6 +536,10 @@ /compute_damage_atom.h /compute_dilatation_atom.cpp /compute_dilatation_atom.h +/compute_dipole_tip4p.cpp +/compute_dipole_tip4p.h +/compute_dipole_tip4p_chunk.cpp +/compute_dipole_tip4p_chunk.h /compute_dpd.cpp /compute_dpd.h /compute_dpd_atom.cpp @@ -754,6 +758,8 @@ /fix_damping_cundall.h /fix_dpd_energy.cpp /fix_dpd_energy.h +/fix_efield_tip4p.cpp +/fix_efield_tip4p.h /fix_electron_stopping.cpp /fix_electron_stopping.h /fix_electron_stopping_fit.cpp diff --git a/src/EXTRA-COMPUTE/compute_dipole_tip4p.cpp b/src/EXTRA-COMPUTE/compute_dipole_tip4p.cpp new file mode 100644 index 0000000000..f57fe6a9fb --- /dev/null +++ b/src/EXTRA-COMPUTE/compute_dipole_tip4p.cpp @@ -0,0 +1,220 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "compute_dipole_tip4p.h" + +#include "angle.h" +#include "atom.h" +#include "bond.h" +#include "domain.h" +#include "error.h" +#include "force.h" +#include "pair.h" +#include "update.h" + +#include +#include + +using namespace LAMMPS_NS; + +enum { MASSCENTER, GEOMCENTER }; + +/* ---------------------------------------------------------------------- */ + +ComputeDipoleTIP4P::ComputeDipoleTIP4P(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) +{ + if ((narg < 3) || (narg > 4)) error->all(FLERR, "Illegal compute dipole/tip4p command"); + + scalar_flag = 1; + vector_flag = 1; + size_vector = 3; + extscalar = 0; + extvector = 0; + + vector = new double[size_vector]; + vector[0] = vector[1] = vector[2] = 0.0; + usecenter = MASSCENTER; + + if (narg == 4) { + if (utils::strmatch(arg[3], "^geom")) + usecenter = GEOMCENTER; + else if (strcmp(arg[3], "mass") == 0) + usecenter = MASSCENTER; + else + error->all(FLERR, "Illegal compute dipole/tip4p command"); + } +} + +/* ---------------------------------------------------------------------- */ + +ComputeDipoleTIP4P::~ComputeDipoleTIP4P() +{ + delete[] vector; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeDipoleTIP4P::init() +{ + if (!force->pair) error->all(FLERR, "Pair style must be defined for compute dipole/ti4p"); + + int itmp; + double *p_qdist = (double *) force->pair->extract("qdist", itmp); + int *p_typeO = (int *) force->pair->extract("typeO", itmp); + int *p_typeH = (int *) force->pair->extract("typeH", itmp); + int *p_typeA = (int *) force->pair->extract("typeA", itmp); + int *p_typeB = (int *) force->pair->extract("typeB", itmp); + if (!p_qdist || !p_typeO || !p_typeH || !p_typeA || !p_typeB) + error->all(FLERR, "Pair style is incompatible with compute dipole/tip4p"); + typeO = *p_typeO; + typeH = *p_typeH; + int typeA = *p_typeA; + int typeB = *p_typeB; + + if (!force->angle || !force->bond || !force->angle->setflag || !force->bond->setflag) + error->all(FLERR, "Bond and angle potentials must be defined for compute dipole/tip4p"); + if ((typeA < 1) || (typeA > atom->nangletypes) || (force->angle->setflag[typeA] == 0)) + error->all(FLERR, "Bad TIP4P angle type for compute dipole/tip4p"); + if ((typeB < 1) || (typeB > atom->nbondtypes) || (force->bond->setflag[typeB] == 0)) + error->all(FLERR, "Bad TIP4P bond type for compute dipole/tip4p"); + double theta = force->angle->equilibrium_angle(typeA); + double blen = force->bond->equilibrium_distance(typeB); + alpha = *p_qdist / (cos(0.5 * theta) * blen); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeDipoleTIP4P::compute_vector() +{ + invoked_vector = update->ntimestep; + + const auto x = atom->x; + const auto mask = atom->mask; + const auto type = atom->type; + const auto image = atom->image; + const auto mass = atom->mass; + const auto rmass = atom->rmass; + const auto q = atom->q; + const auto mu = atom->mu; + const auto nlocal = atom->nlocal; + + double dipole[3] = {0.0, 0.0, 0.0}; + double comproc[3] = {0.0, 0.0, 0.0}; + double com[3] = {0.0, 0.0, 0.0}; + double masstotal = 0.0; + double chrgtotal = 0.0; + double massproc = 0.0; + double chrgproc = 0.0; + double unwrap[3], xM[3]; + double *xi; + + for (int i = 0; i < nlocal; ++i) { + if (mask[i] & groupbit) { + double massone = 1.0; // for usecenter == GEOMCENTER + if (usecenter == MASSCENTER) { + if (rmass) + massone = rmass[i]; + else + massone = mass[type[i]]; + } + massproc += massone; + if (atom->q_flag) chrgproc += q[i]; + domain->unmap(x[i], image[i], unwrap); + comproc[0] += unwrap[0] * massone; + comproc[1] += unwrap[1] * massone; + comproc[2] += unwrap[2] * massone; + } + } + MPI_Allreduce(&massproc, &masstotal, 1, MPI_DOUBLE, MPI_SUM, world); + MPI_Allreduce(&chrgproc, &chrgtotal, 1, MPI_DOUBLE, MPI_SUM, world); + MPI_Allreduce(comproc, com, 3, MPI_DOUBLE, MPI_SUM, world); + + if (masstotal > 0.0) { + com[0] /= masstotal; + com[1] /= masstotal; + com[2] /= masstotal; + } + + // compute dipole moment + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + if (type[i] == typeO) { + find_M(i, xM); + xi = xM; + } else { + xi = x[i]; + } + domain->unmap(xi, image[i], unwrap); + if (atom->q_flag) { + dipole[0] += q[i] * unwrap[0]; + dipole[1] += q[i] * unwrap[1]; + dipole[2] += q[i] * unwrap[2]; + } + if (atom->mu_flag) { + dipole[0] += mu[i][0]; + dipole[1] += mu[i][1]; + dipole[2] += mu[i][2]; + } + } + } + + MPI_Allreduce(dipole, vector, 3, MPI_DOUBLE, MPI_SUM, world); + + // correct for position dependence with a net charged group + vector[0] -= chrgtotal * com[0]; + vector[1] -= chrgtotal * com[1]; + vector[2] -= chrgtotal * com[2]; +} + +/* ---------------------------------------------------------------------- */ + +double ComputeDipoleTIP4P::compute_scalar() +{ + if (invoked_vector != update->ntimestep) compute_vector(); + + invoked_scalar = update->ntimestep; + scalar = sqrt(vector[0] * vector[0] + vector[1] * vector[1] + vector[2] * vector[2]); + return scalar; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeDipoleTIP4P::find_M(int i, double *xM) +{ + double **x = atom->x; + + int iH1 = atom->map(atom->tag[i] + 1); + int iH2 = atom->map(atom->tag[i] + 2); + + if ((iH1 == -1) || (iH2 == -1)) error->one(FLERR, "TIP4P hydrogen is missing"); + if ((atom->type[iH1] != typeH) || (atom->type[iH2] != typeH)) + error->one(FLERR, "TIP4P hydrogen has incorrect atom type"); + + // set iH1,iH2 to index of closest image to O + + iH1 = domain->closest_image(i, iH1); + iH2 = domain->closest_image(i, iH2); + + double delx1 = x[iH1][0] - x[i][0]; + double dely1 = x[iH1][1] - x[i][1]; + double delz1 = x[iH1][2] - x[i][2]; + + double delx2 = x[iH2][0] - x[i][0]; + double dely2 = x[iH2][1] - x[i][1]; + double delz2 = x[iH2][2] - x[i][2]; + + xM[0] = x[i][0] + alpha * 0.5 * (delx1 + delx2); + xM[1] = x[i][1] + alpha * 0.5 * (dely1 + dely2); + xM[2] = x[i][2] + alpha * 0.5 * (delz1 + delz2); +} diff --git a/src/EXTRA-COMPUTE/compute_dipole_tip4p.h b/src/EXTRA-COMPUTE/compute_dipole_tip4p.h new file mode 100644 index 0000000000..75784b83ec --- /dev/null +++ b/src/EXTRA-COMPUTE/compute_dipole_tip4p.h @@ -0,0 +1,46 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef COMPUTE_CLASS +// clang-format off +ComputeStyle(dipole/tip4p,ComputeDipoleTIP4P); +// clang-format on +#else + +#ifndef LMP_COMPUTE_DIPOLE_TIP4P_H +#define LMP_COMPUTE_DIPOLE_TIP4P_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeDipoleTIP4P : public Compute { + public: + ComputeDipoleTIP4P(class LAMMPS *, int, char **); + ~ComputeDipoleTIP4P() override; + void init() override; + void compute_vector() override; + double compute_scalar() override; + + private: + int usecenter; + + int typeO, typeH; + double alpha; + void find_M(int i, double *xM); +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/EXTRA-COMPUTE/compute_dipole_tip4p_chunk.cpp b/src/EXTRA-COMPUTE/compute_dipole_tip4p_chunk.cpp new file mode 100644 index 0000000000..651811bba5 --- /dev/null +++ b/src/EXTRA-COMPUTE/compute_dipole_tip4p_chunk.cpp @@ -0,0 +1,365 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "compute_dipole_tip4p_chunk.h" + +#include "angle.h" +#include "atom.h" +#include "bond.h" +#include "comm.h" +#include "compute_chunk_atom.h" +#include "domain.h" +#include "error.h" +#include "force.h" +#include "pair.h" +#include "math_special.h" +#include "memory.h" +#include "modify.h" +#include "update.h" + +#include +#include + +using namespace LAMMPS_NS; +using namespace MathSpecial; + +enum { MASSCENTER, GEOMCENTER }; + +/* ---------------------------------------------------------------------- */ + +ComputeDipoleTIP4PChunk::ComputeDipoleTIP4PChunk(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg), + idchunk(nullptr), massproc(nullptr), masstotal(nullptr), chrgproc(nullptr), + chrgtotal(nullptr), com(nullptr), + comall(nullptr), dipole(nullptr), dipoleall(nullptr) +{ + if ((narg != 4) && (narg != 5)) + error->all(FLERR,"Illegal compute dipole/tip4p/chunk command"); + + array_flag = 1; + size_array_cols = 4; + size_array_rows = 0; + size_array_rows_variable = 1; + extarray = 0; + + // ID of compute chunk/atom + + idchunk = utils::strdup(arg[3]); + + usecenter = MASSCENTER; + + if (narg == 5) { + if (strncmp(arg[4],"geom",4) == 0) usecenter = GEOMCENTER; + else if (strcmp(arg[4],"mass") == 0) usecenter = MASSCENTER; + else error->all(FLERR,"Illegal compute dipole/tip4p/chunk command"); + } + + ComputeDipoleTIP4PChunk::init(); + + // chunk-based data + + nchunk = 1; + maxchunk = 0; + allocate(); +} + +/* ---------------------------------------------------------------------- */ + +ComputeDipoleTIP4PChunk::~ComputeDipoleTIP4PChunk() +{ + delete[] idchunk; + memory->destroy(massproc); + memory->destroy(masstotal); + memory->destroy(chrgproc); + memory->destroy(chrgtotal); + memory->destroy(com); + memory->destroy(comall); + memory->destroy(dipole); + memory->destroy(dipoleall); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeDipoleTIP4PChunk::init() +{ + int icompute = modify->find_compute(idchunk); + if (icompute < 0) + error->all(FLERR,"Chunk/atom compute does not exist for compute dipole/tip4p/chunk"); + cchunk = dynamic_cast(modify->compute[icompute]); + if (strcmp(cchunk->style,"chunk/atom") != 0) + error->all(FLERR,"Compute dipole/tip4p/chunk does not use chunk/atom compute"); + + if (!force->pair) error->all(FLERR, "Pair style must be defined for compute dipole/tip4p/chunk"); + + int itmp; + double *p_qdist = (double *) force->pair->extract("qdist", itmp); + int *p_typeO = (int *) force->pair->extract("typeO", itmp); + int *p_typeH = (int *) force->pair->extract("typeH", itmp); + int *p_typeA = (int *) force->pair->extract("typeA", itmp); + int *p_typeB = (int *) force->pair->extract("typeB", itmp); + if (!p_qdist || !p_typeO || !p_typeH || !p_typeA || !p_typeB) + error->all(FLERR, "Pair style is incompatible with compute dipole/tip4p/chunk"); + typeO = *p_typeO; + typeH = *p_typeH; + int typeA = *p_typeA; + int typeB = *p_typeB; + + if (!force->angle || !force->bond || !force->angle->setflag || !force->bond->setflag) + error->all(FLERR, "Bond and angle potentials must be defined for compute dipole/tip4p/chunk"); + if ((typeA < 1) || (typeA > atom->nangletypes) || (force->angle->setflag[typeA] == 0)) + error->all(FLERR, "Bad TIP4P angle type for compute dipole/tip4p/chunk"); + if ((typeB < 1) || (typeB > atom->nbondtypes) || (force->bond->setflag[typeB] == 0)) + error->all(FLERR, "Bad TIP4P bond type for compute dipole/tip4p/chunk"); + double theta = force->angle->equilibrium_angle(typeA); + double blen = force->bond->equilibrium_distance(typeB); + alpha = *p_qdist / (cos(0.5 * theta) * blen); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeDipoleTIP4PChunk::compute_array() +{ + int i,index; + double massone; + double unwrap[3]; + + invoked_array = update->ntimestep; + + // compute chunk/atom assigns atoms to chunk IDs + // extract ichunk index vector from compute + // ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms + + nchunk = cchunk->setup_chunks(); + cchunk->compute_ichunk(); + int *ichunk = cchunk->ichunk; + + if (nchunk > maxchunk) allocate(); + size_array_rows = nchunk; + + // zero local per-chunk values + + for (i = 0; i < nchunk; i++) { + massproc[i] = chrgproc[i] = 0.0; + com[i][0] = com[i][1] = com[i][2] = 0.0; + dipole[i][0] = dipole[i][1] = dipole[i][2] = dipole[i][3] = 0.0; + } + + // compute COM for each chunk + + double **x = atom->x; + int *mask = atom->mask; + int *type = atom->type; + imageint *image = atom->image; + double *mass = atom->mass; + double *rmass = atom->rmass; + double *q = atom->q; + double **mu = atom->mu; + double xM[3]; + double *xi; + + int nlocal = atom->nlocal; + + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + index = ichunk[i]-1; + if (index < 0) continue; + if (usecenter == MASSCENTER) { + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; + } else massone = 1.0; // usecenter == GEOMCENTER + domain->unmap(x[i],image[i],unwrap); + massproc[index] += massone; + if (atom->q_flag) chrgproc[index] += q[i]; + com[index][0] += unwrap[0] * massone; + com[index][1] += unwrap[1] * massone; + com[index][2] += unwrap[2] * massone; + } + + MPI_Allreduce(massproc,masstotal,nchunk,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(chrgproc,chrgtotal,nchunk,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&com[0][0],&comall[0][0],3*nchunk,MPI_DOUBLE,MPI_SUM,world); + + for (i = 0; i < nchunk; i++) { + if (masstotal[i] > 0.0) { + comall[i][0] /= masstotal[i]; + comall[i][1] /= masstotal[i]; + comall[i][2] /= masstotal[i]; + } + } + + // compute dipole for each chunk + + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + index = ichunk[i]-1; + if (index < 0) continue; + + if (type[i] == typeO) { + find_M(i,xM); + xi = xM; + } else { + xi = x[i]; + } + domain->unmap(xi,image[i],unwrap); + if (atom->q_flag) { + dipole[index][0] += q[i]*unwrap[0]; + dipole[index][1] += q[i]*unwrap[1]; + dipole[index][2] += q[i]*unwrap[2]; + } + if (atom->mu_flag) { + dipole[index][0] += mu[i][0]; + dipole[index][1] += mu[i][1]; + dipole[index][2] += mu[i][2]; + } + } + } + + MPI_Allreduce(&dipole[0][0],&dipoleall[0][0],4*nchunk,MPI_DOUBLE,MPI_SUM,world); + + for (i = 0; i < nchunk; i++) { + // correct for position dependence with charged chunks + dipoleall[i][0] -= chrgtotal[i]*comall[i][0]; + dipoleall[i][1] -= chrgtotal[i]*comall[i][1]; + dipoleall[i][2] -= chrgtotal[i]*comall[i][2]; + // compute total dipole moment + dipoleall[i][3] = sqrt(square(dipoleall[i][0]) + + square(dipoleall[i][1]) + + square(dipoleall[i][2])); + } +} + +/* ---------------------------------------------------------------------- + lock methods: called by fix ave/time + these methods ensure vector/array size is locked for Nfreq epoch + by passing lock info along to compute chunk/atom +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + increment lock counter +------------------------------------------------------------------------- */ + +void ComputeDipoleTIP4PChunk::lock_enable() +{ + cchunk->lockcount++; +} + +/* ---------------------------------------------------------------------- + decrement lock counter in compute chunk/atom, it if still exists +------------------------------------------------------------------------- */ + +void ComputeDipoleTIP4PChunk::lock_disable() +{ + int icompute = modify->find_compute(idchunk); + if (icompute >= 0) { + cchunk = dynamic_cast(modify->compute[icompute]); + cchunk->lockcount--; + } +} + +/* ---------------------------------------------------------------------- + calculate and return # of chunks = length of vector/array +------------------------------------------------------------------------- */ + +int ComputeDipoleTIP4PChunk::lock_length() +{ + nchunk = cchunk->setup_chunks(); + return nchunk; +} + +/* ---------------------------------------------------------------------- + set the lock from startstep to stopstep +------------------------------------------------------------------------- */ + +void ComputeDipoleTIP4PChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep) +{ + cchunk->lock(fixptr,startstep,stopstep); +} + +/* ---------------------------------------------------------------------- + unset the lock +------------------------------------------------------------------------- */ + +void ComputeDipoleTIP4PChunk::unlock(Fix *fixptr) +{ + cchunk->unlock(fixptr); +} + +/* ---------------------------------------------------------------------- + free and reallocate per-chunk arrays +------------------------------------------------------------------------- */ + +void ComputeDipoleTIP4PChunk::allocate() +{ + memory->destroy(massproc); + memory->destroy(masstotal); + memory->destroy(chrgproc); + memory->destroy(chrgtotal); + memory->destroy(com); + memory->destroy(comall); + memory->destroy(dipole); + memory->destroy(dipoleall); + maxchunk = nchunk; + memory->create(massproc,maxchunk,"dipole/tip4p/chunk:massproc"); + memory->create(masstotal,maxchunk,"dipole/tip4p/chunk:masstotal"); + memory->create(chrgproc,maxchunk,"dipole/tip4p/chunk:chrgproc"); + memory->create(chrgtotal,maxchunk,"dipole/tip4p/chunk:chrgtotal"); + memory->create(com,maxchunk,3,"dipole/tip4p/chunk:com"); + memory->create(comall,maxchunk,3,"dipole/tip4p/chunk:comall"); + memory->create(dipole,maxchunk,4,"dipole/tip4p/chunk:dipole"); + memory->create(dipoleall,maxchunk,4,"dipole/tip4p/chunk:dipoleall"); + array = dipoleall; +} + +/* ---------------------------------------------------------------------- + memory usage of local data +------------------------------------------------------------------------- */ + +double ComputeDipoleTIP4PChunk::memory_usage() +{ + double bytes = (bigint) maxchunk * 2 * sizeof(double); + bytes += (double)maxchunk * 2*3 * sizeof(double); + bytes += (double)maxchunk * 2*4 * sizeof(double); + return bytes; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeDipoleTIP4PChunk::find_M(int i, double *xM) +{ + double **x = atom->x; + + int iH1 = atom->map(atom->tag[i] + 1); + int iH2 = atom->map(atom->tag[i] + 2); + + if ((iH1 == -1) || (iH2 == -1)) error->one(FLERR,"TIP4P hydrogen is missing"); + if ((atom->type[iH1] != typeH) || (atom->type[iH2] != typeH)) + error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); + + // set iH1,iH2 to index of closest image to O + + iH1 = domain->closest_image(i,iH1); + iH2 = domain->closest_image(i,iH2); + + double delx1 = x[iH1][0] - x[i][0]; + double dely1 = x[iH1][1] - x[i][1]; + double delz1 = x[iH1][2] - x[i][2]; + + double delx2 = x[iH2][0] - x[i][0]; + double dely2 = x[iH2][1] - x[i][1]; + double delz2 = x[iH2][2] - x[i][2]; + + xM[0] = x[i][0] + alpha * 0.5 * (delx1 + delx2); + xM[1] = x[i][1] + alpha * 0.5 * (dely1 + dely2); + xM[2] = x[i][2] + alpha * 0.5 * (delz1 + delz2); +} diff --git a/src/EXTRA-COMPUTE/compute_dipole_tip4p_chunk.h b/src/EXTRA-COMPUTE/compute_dipole_tip4p_chunk.h new file mode 100644 index 0000000000..63d6e8401e --- /dev/null +++ b/src/EXTRA-COMPUTE/compute_dipole_tip4p_chunk.h @@ -0,0 +1,64 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef COMPUTE_CLASS +// clang-format off +ComputeStyle(dipole/tip4p/chunk,ComputeDipoleTIP4PChunk); +// clang-format on +#else + +#ifndef LMP_COMPUTE_DIPOLE_TIP4P_CHUNK_H +#define LMP_COMPUTE_DIPOLE_TIP4P_CHUNK_H + +#include "compute.h" + +namespace LAMMPS_NS { +class Fix; + +class ComputeDipoleTIP4PChunk : public Compute { + public: + ComputeDipoleTIP4PChunk(class LAMMPS *, int, char **); + ~ComputeDipoleTIP4PChunk() override; + void init() override; + void compute_array() override; + + void lock_enable() override; + void lock_disable() override; + int lock_length() override; + void lock(Fix *, bigint, bigint) override; + void unlock(Fix *) override; + + double memory_usage() override; + + private: + int nchunk, maxchunk; + char *idchunk; + class ComputeChunkAtom *cchunk; + + double *massproc, *masstotal; + double *chrgproc, *chrgtotal; + double **com, **comall; + double **dipole, **dipoleall; + int usecenter; + + int typeO, typeH; + double alpha; + void find_M(int i, double *xM); + + void allocate(); +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/EXTRA-FIX/fix_efield_tip4p.cpp b/src/EXTRA-FIX/fix_efield_tip4p.cpp new file mode 100644 index 0000000000..a204fe6f7d --- /dev/null +++ b/src/EXTRA-FIX/fix_efield_tip4p.cpp @@ -0,0 +1,644 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "fix_efield_tip4p.h" + +#include "angle.h" +#include "atom.h" +#include "bond.h" +#include "comm.h" +#include "domain.h" +#include "error.h" +#include "force.h" +#include "input.h" +#include "memory.h" +#include "modify.h" +#include "pair.h" +#include "region.h" +#include "respa.h" +#include "update.h" +#include "variable.h" + +#include + +using namespace LAMMPS_NS; +using namespace FixConst; + +/* ---------------------------------------------------------------------- */ + +FixEfieldTIP4P::FixEfieldTIP4P(LAMMPS *_lmp, int narg, char **arg) : FixEfield(_lmp, narg, arg) {} + +/* ---------------------------------------------------------------------- */ + +void FixEfieldTIP4P::init() +{ + FixEfield::init(); + + if (atom->tag_enable == 0) error->all(FLERR, "Fix efield/tip4p requires atom IDs"); + if (!atom->q_flag) error->all(FLERR, "Fix efield/tip4p requires atom attribute q"); + + int itmp; + double *p_qdist = (double *) force->pair->extract("qdist", itmp); + int *p_typeO = (int *) force->pair->extract("typeO", itmp); + int *p_typeH = (int *) force->pair->extract("typeH", itmp); + int *p_typeA = (int *) force->pair->extract("typeA", itmp); + int *p_typeB = (int *) force->pair->extract("typeB", itmp); + if (!p_qdist || !p_typeO || !p_typeH || !p_typeA || !p_typeB) + error->all(FLERR, "Pair style is incompatible with compute {}", style); + typeO = *p_typeO; + typeH = *p_typeH; + int typeA = *p_typeA; + int typeB = *p_typeB; + + if (force->angle == nullptr || force->bond == nullptr || force->angle->setflag == nullptr || + force->bond->setflag == nullptr) + error->all(FLERR, "Bond and angle potentials must be defined for compute {}", style); + if ((typeA < 1) || (typeA > atom->nangletypes) || (force->angle->setflag[typeA] == 0)) + error->all(FLERR, "Bad TIP4P angle type for compute {}", style); + if ((typeB < 1) || (typeB > atom->nbondtypes) || (force->bond->setflag[typeB] == 0)) + error->all(FLERR, "Bad TIP4P bond type for PPPM/TIP4P"); + + // determine alpha parameter + + const double theta = force->angle->equilibrium_angle(typeA); + const double blen = force->bond->equilibrium_distance(typeB); + alpha = *p_qdist / (cos(0.5 * theta) * blen); +} + +/* ---------------------------------------------------------------------- + apply F = qE +------------------------------------------------------------------------- */ + +void FixEfieldTIP4P::post_force(int vflag) +{ + double **f = atom->f; + double *q = atom->q; + tagint *tag = atom->tag; + int *type = atom->type; + int *mask = atom->mask; + imageint *image = atom->image; + int nlocal = atom->nlocal; + + // virial setup + + v_init(vflag); + + // reallocate efield array if necessary + + if ((varflag == ATOM) && (atom->nmax > maxatom)) { + maxatom = atom->nmax; + memory->destroy(efield); + memory->create(efield, maxatom, 4, "efield:efield"); + } + + // update region if necessary + + if (region) region->prematch(); + + // fsum[0] = "potential energy" for added force + // fsum[123] = extra force added to atoms + + fsum[0] = fsum[1] = fsum[2] = fsum[3] = 0.0; + force_flag = 0; + + double **x = atom->x; + double fx, fy, fz, xM[3]; + double v[6], unwrap[3]; + int iO, iH1, iH2; + + // constant efield + + if (varflag == CONSTANT) { + + // charge interactions + // force = qE, potential energy = F dot x in unwrapped coords + + if (qflag) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + + // process *all* atoms belonging to a TIP4P molecule since we need + // to consider molecules where parts of the molecule are ghost atoms + + if ((type[i] == typeO) || (type[i] == typeH)) { + + if (type[i] == typeO) { + iO = i; + iH1 = atom->map(tag[i] + 1); + iH2 = atom->map(tag[i] + 2); + } else { + + // set indices for first or second hydrogen + + iO = atom->map(tag[i] - 1); + if ((iO != -1) && (type[iO] == typeO)) { + iH1 = i; + iH2 = atom->map(tag[i] + 1); + } else { + iO = atom->map(tag[i] - 2); + iH1 = atom->map(tag[i] - 1); + iH2 = i; + } + } + if ((iH1 == -1) || (iH2 == -1)) error->one(FLERR, "TIP4P hydrogen is missing"); + if (iO == -1) error->one(FLERR, "TIP4P oxygen is missing"); + if ((type[iH1] != typeH) || (type[iH2] != typeH)) + error->one(FLERR, "TIP4P hydrogen has incorrect atom type"); + if (type[iO] != typeO) error->one(FLERR, "TIP4P oxygen has incorrect atom type"); + iH1 = domain->closest_image(i, iH1); + iH2 = domain->closest_image(i, iH2); + + find_M(x[iO], x[iH1], x[iH2], xM); + + // M site contributions + + if (!region || region->match(xM[0], xM[1], xM[2])) { + + fx = q[iO] * ex; + fy = q[iO] * ey; + fz = q[iO] * ez; + + // distribute and apply forces + + if (i == iO) { + f[iO][0] += fx * (1.0 - alpha); + f[iO][1] += fy * (1.0 - alpha); + f[iO][2] += fz * (1.0 - alpha); + if (iH1 < nlocal) { + f[iH1][0] += 0.5 * alpha * fx; + f[iH1][1] += 0.5 * alpha * fy; + f[iH1][2] += 0.5 * alpha * fz; + } + if (iH2 < nlocal) { + f[iH2][0] += 0.5 * alpha * fx; + f[iH2][1] += 0.5 * alpha * fy; + f[iH2][2] += 0.5 * alpha * fz; + } + } else { + if (iO >= nlocal) { + if (i == iH1) { + f[iH1][0] += 0.5 * alpha * fx; + f[iH1][1] += 0.5 * alpha * fy; + f[iH1][2] += 0.5 * alpha * fz; + } + if (i == iH2) { + f[iH2][0] += 0.5 * alpha * fx; + f[iH2][1] += 0.5 * alpha * fy; + f[iH2][2] += 0.5 * alpha * fz; + } + } + } + + if (i == iO) { + domain->unmap(xM, image[iO], unwrap); + fsum[0] -= fx * unwrap[0] + fy * unwrap[1] + fz * unwrap[2]; + fsum[1] += fx; + fsum[2] += fy; + fsum[3] += fz; + + // tally virial contribution with Oxygen + + if (evflag) { + v[0] = fx * unwrap[0]; + v[1] = fy * unwrap[1]; + v[2] = fz * unwrap[2]; + v[3] = fx * unwrap[1]; + v[4] = fx * unwrap[2]; + v[5] = fy * unwrap[2]; + v_tally(iO, v); + } + } + } + + // H1 site contributions + + if (!region || region->match(x[iH1][0], x[iH1][1], x[iH1][2])) { + + fx = q[iH1] * ex; + fy = q[iH1] * ey; + fz = q[iH1] * ez; + + if (i == iH1) { + f[iH1][0] += fx; + f[iH1][1] += fy; + f[iH1][2] += fz; + + // tally global force + + domain->unmap(x[iH1], image[iH1], unwrap); + fsum[0] -= fx * unwrap[0] + fy * unwrap[1] + fz * unwrap[2]; + fsum[1] += fx; + fsum[2] += fy; + fsum[3] += fz; + + // tally virial contributions + + if (evflag) { + v[0] = fx * unwrap[0]; + v[1] = fy * unwrap[1]; + v[2] = fz * unwrap[2]; + v[3] = fx * unwrap[1]; + v[4] = fx * unwrap[2]; + v[5] = fy * unwrap[2]; + v_tally(iH1, v); + } + } + } + + // H2 site contributions + + if (!region || region->match(x[iH2][0], x[iH2][1], x[iH2][2])) { + + fx = q[iH2] * ex; + fy = q[iH2] * ey; + fz = q[iH2] * ez; + + if (i == iH2) { + f[iH2][0] += fx; + f[iH2][1] += fy; + f[iH2][2] += fz; + + // tally global force + + domain->unmap(x[iH2], image[iH2], unwrap); + fsum[0] -= fx * unwrap[0] + fy * unwrap[1] + fz * unwrap[2]; + fsum[1] += fx; + fsum[2] += fy; + fsum[3] += fz; + + // tally virial contributions + + if (evflag) { + v[0] = fx * unwrap[0]; + v[1] = fy * unwrap[1]; + v[2] = fz * unwrap[2]; + v[3] = fx * unwrap[1]; + v[4] = fx * unwrap[2]; + v[5] = fy * unwrap[2]; + v_tally(iH2, v); + } + } + } + + // non-TIP4P atoms + + } else { + + if (region && !region->match(x[i][0], x[i][1], x[i][2])) continue; + fx = q[i] * ex; + fy = q[i] * ey; + fz = q[i] * ez; + f[i][0] += fx; + f[i][1] += fy; + f[i][2] += fz; + + domain->unmap(x[i], image[i], unwrap); + fsum[0] -= fx * unwrap[0] + fy * unwrap[1] + fz * unwrap[2]; + fsum[1] += fx; + fsum[2] += fy; + fsum[3] += fz; + if (evflag) { + v[0] = fx * unwrap[0]; + v[1] = fy * unwrap[1]; + v[2] = fz * unwrap[2]; + v[3] = fx * unwrap[1]; + v[4] = fx * unwrap[2]; + v[5] = fy * unwrap[2]; + v_tally(i, v); + } + } + } + } + } + + // dipole interactions, no special TIP4P treatment needed + // no force, torque = mu cross E, potential energy = -mu dot E + + if (muflag) { + double **mu = atom->mu; + double **t = atom->torque; + double tx, ty, tz; + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + if (region && !region->match(x[i][0], x[i][1], x[i][2])) continue; + tx = ez * mu[i][1] - ey * mu[i][2]; + ty = ex * mu[i][2] - ez * mu[i][0]; + tz = ey * mu[i][0] - ex * mu[i][1]; + t[i][0] += tx; + t[i][1] += ty; + t[i][2] += tz; + fsum[0] -= mu[i][0] * ex + mu[i][1] * ey + mu[i][2] * ez; + } + } + + } else { + + // variable efield, wrap with clear/add + // potential energy = evar if defined, else 0.0 + + modify->clearstep_compute(); + + if (xstyle == EQUAL) { + ex = qe2f * input->variable->compute_equal(xvar); + } else if (xstyle == ATOM) { + input->variable->compute_atom(xvar, igroup, &efield[0][0], 4, 0); + } + if (ystyle == EQUAL) { + ey = qe2f * input->variable->compute_equal(yvar); + } else if (ystyle == ATOM) { + input->variable->compute_atom(yvar, igroup, &efield[0][1], 4, 0); + } + if (zstyle == EQUAL) { + ez = qe2f * input->variable->compute_equal(zvar); + } else if (zstyle == ATOM) { + input->variable->compute_atom(zvar, igroup, &efield[0][2], 4, 0); + } + if (estyle == ATOM) input->variable->compute_atom(evar, igroup, &efield[0][3], 4, 0); + + modify->addstep_compute(update->ntimestep + 1); + + // charge interactions + // force = qE + + if (qflag) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + + // process *all* atoms belonging to a TIP4P molecule since we need + // to consider molecules where parts of the molecule are ghost atoms + + if ((type[i] == typeO) || (type[i] == typeH)) { + + if (type[i] == typeO) { + iO = i; + iH1 = atom->map(tag[i] + 1); + iH2 = atom->map(tag[i] + 2); + } else { + + // set indices for first or second hydrogen + + iO = atom->map(tag[i] - 1); + if ((iO != -1) && (type[iO] == typeO)) { + iH1 = i; + iH2 = atom->map(tag[i] + 1); + } else { + iO = atom->map(tag[i] - 2); + iH1 = atom->map(tag[i] - 1); + iH2 = i; + } + } + if ((iH1 == -1) || (iH2 == -1)) error->one(FLERR, "TIP4P hydrogen is missing"); + if (iO == -1) error->one(FLERR, "TIP4P oxygen is missing"); + if ((type[iH1] != typeH) || (type[iH2] != typeH)) + error->one(FLERR, "TIP4P hydrogen has incorrect atom type"); + if (type[iO] != typeO) error->one(FLERR, "TIP4P oxygen has incorrect atom type"); + iH1 = domain->closest_image(i, iH1); + iH2 = domain->closest_image(i, iH2); + + find_M(x[iO], x[iH1], x[iH2], xM); + + // M site contributions + + if (!region || region->match(xM[0], xM[1], xM[2])) { + + if (xstyle == ATOM) { + fx = qe2f * q[iO] * efield[iO][0]; + } else { + fx = q[iO] * ex; + } + if (ystyle == ATOM) { + fy = qe2f * q[iO] * efield[iO][1]; + } else { + fy = q[iO] * ey; + } + if (zstyle == ATOM) { + fz = qe2f * q[iO] * efield[iO][2]; + } else { + fz = q[iO] * ez; + } + + // distribute and apply forces, but only to local atoms + + if (i == iO) { + f[iO][0] += fx * (1.0 - alpha); + f[iO][1] += fy * (1.0 - alpha); + f[iO][2] += fz * (1.0 - alpha); + if (iH1 < nlocal) { + f[iH1][0] += 0.5 * alpha * fx; + f[iH1][1] += 0.5 * alpha * fy; + f[iH1][2] += 0.5 * alpha * fz; + } + if (iH2 < nlocal) { + f[iH2][0] += 0.5 * alpha * fx; + f[iH2][1] += 0.5 * alpha * fy; + f[iH2][2] += 0.5 * alpha * fz; + } + } else { + if (iO >= nlocal) { + if (i == iH1) { + f[iH1][0] += 0.5 * alpha * fx; + f[iH1][1] += 0.5 * alpha * fy; + f[iH1][2] += 0.5 * alpha * fz; + } + if (i == iH2) { + f[iH2][0] += 0.5 * alpha * fx; + f[iH2][1] += 0.5 * alpha * fy; + f[iH2][2] += 0.5 * alpha * fz; + } + } + } + + if (i == iO) { + if (estyle == ATOM) fsum[0] += efield[0][3]; + fsum[1] += fx; + fsum[2] += fy; + fsum[3] += fz; + + // tally virial contribution for point M with Oxygen + + if (evflag) { + v[0] = fx * unwrap[0]; + v[1] = fy * unwrap[1]; + v[2] = fz * unwrap[2]; + v[3] = fx * unwrap[1]; + v[4] = fx * unwrap[2]; + v[5] = fy * unwrap[2]; + v_tally(iO, v); + } + } + } + + // H1 site contributions + + if (!region || region->match(x[iH1][0], x[iH1][1], x[iH1][2])) { + + if (xstyle == ATOM) { + fx = qe2f * q[iH1] * efield[iH1][0]; + } else { + fx = q[iH1] * ex; + } + if (ystyle == ATOM) { + fy = qe2f * q[iH1] * efield[iH1][1]; + } else { + fy = q[iH1] * ey; + } + if (zstyle == ATOM) { + fz = qe2f * q[iH1] * efield[iH1][2]; + } else { + fz = q[iH1] * ez; + } + + if (i == iH1) { + f[iH1][0] += fx; + f[iH1][1] += fy; + f[iH1][2] += fz; + + // tally global force + + if (estyle == ATOM) fsum[0] += efield[0][3]; + fsum[1] += fx; + fsum[2] += fy; + fsum[3] += fz; + + // tally virial contribution + + if (evflag) { + v[0] = fx * unwrap[0]; + v[1] = fy * unwrap[1]; + v[2] = fz * unwrap[2]; + v[3] = fx * unwrap[1]; + v[4] = fx * unwrap[2]; + v[5] = fy * unwrap[2]; + v_tally(iH1, v); + } + } + } + + // H2 site contributions + + if (!region || region->match(x[iH2][0], x[iH2][1], x[iH2][2])) { + + if (xstyle == ATOM) { + fx = qe2f * q[iH2] * efield[iH2][0]; + } else { + fx = q[iH2] * ex; + } + if (ystyle == ATOM) { + fy = qe2f * q[iH2] * efield[iH2][1]; + } else { + fy = q[iH2] * ey; + } + if (zstyle == ATOM) { + fz = qe2f * q[iH2] * efield[iH2][2]; + } else { + fz = q[iH2] * ez; + } + + if (i == iH2) { + f[iH2][0] += fx; + f[iH2][1] += fy; + f[iH2][2] += fz; + + // tally global force + + if (estyle == ATOM) fsum[0] += efield[0][3]; + fsum[1] += fx; + fsum[2] += fy; + fsum[3] += fz; + + // tally virial contribution + + if (evflag) { + v[0] = fx * unwrap[0]; + v[1] = fy * unwrap[1]; + v[2] = fz * unwrap[2]; + v[3] = fx * unwrap[1]; + v[4] = fx * unwrap[2]; + v[5] = fy * unwrap[2]; + v_tally(iH2, v); + } + } + } + + } else { + + // non-TIP4P charge interactions + // force = qE + + if (region && !region->match(x[i][0], x[i][1], x[i][2])) continue; + + if (xstyle == ATOM) { + fx = qe2f * q[i] * efield[i][0]; + } else { + fx = q[i] * ex; + } + f[i][0] += fx; + fsum[1] += fx; + if (ystyle == ATOM) { + fy = qe2f * q[i] * efield[i][1]; + } else { + fy = q[i] * ey; + } + f[i][1] += fy; + fsum[2] += fy; + if (zstyle == ATOM) { + fz = qe2f * q[i] * efield[i][2]; + } else { + fz = q[i] * ez; + } + f[i][2] += fz; + fsum[3] += fz; + + if (estyle == ATOM) fsum[0] += efield[0][3]; + } + } + } + } + + // dipole interactions + // no force, torque = mu cross E + + if (muflag) { + double **mu = atom->mu; + double **t = atom->torque; + double tx, ty, tz; + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + if (region && !region->match(x[i][0], x[i][1], x[i][2])) continue; + tx = ez * mu[i][1] - ey * mu[i][2]; + ty = ex * mu[i][2] - ez * mu[i][0]; + tz = ey * mu[i][0] - ex * mu[i][1]; + t[i][0] += tx; + t[i][1] += ty; + t[i][2] += tz; + } + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixEfieldTIP4P::find_M(double *xO, double *xH1, double *xH2, double *xM) +{ + double delx1 = xH1[0] - xO[0]; + double dely1 = xH1[1] - xO[1]; + double delz1 = xH1[2] - xO[2]; + + double delx2 = xH2[0] - xO[0]; + double dely2 = xH2[1] - xO[1]; + double delz2 = xH2[2] - xO[2]; + + xM[0] = xO[0] + alpha * 0.5 * (delx1 + delx2); + xM[1] = xO[1] + alpha * 0.5 * (dely1 + dely2); + xM[2] = xO[2] + alpha * 0.5 * (delz1 + delz2); +} diff --git a/src/EXTRA-FIX/fix_efield_tip4p.h b/src/EXTRA-FIX/fix_efield_tip4p.h new file mode 100644 index 0000000000..6c3cfa58ab --- /dev/null +++ b/src/EXTRA-FIX/fix_efield_tip4p.h @@ -0,0 +1,43 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS +// clang-format off +FixStyle(efield/tip4p,FixEfieldTIP4P); +// clang-format on +#else + +#ifndef LMP_FIX_EFIELD_TIP4P_H +#define LMP_FIX_EFIELD_TIP4P_H + +#include "fix_efield.h" + +namespace LAMMPS_NS { + +class FixEfieldTIP4P : public FixEfield { + + public: + FixEfieldTIP4P(class LAMMPS *, int, char **); + + void init() override; + void post_force(int) override; + + protected: + double alpha; + int typeO, typeH; // atom types for TIP4P molecule + void find_M(double *, double *, double *, double *); +}; +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/MOLECULE/pair_lj_cut_tip4p_cut.cpp b/src/MOLECULE/pair_lj_cut_tip4p_cut.cpp index 19b4d6585f..5be1de82a9 100644 --- a/src/MOLECULE/pair_lj_cut_tip4p_cut.cpp +++ b/src/MOLECULE/pair_lj_cut_tip4p_cut.cpp @@ -18,20 +18,20 @@ #include "pair_lj_cut_tip4p_cut.h" -#include -#include -#include "atom.h" -#include "force.h" -#include "neighbor.h" -#include "neigh_list.h" -#include "domain.h" #include "angle.h" +#include "atom.h" #include "bond.h" #include "comm.h" +#include "domain.h" +#include "error.h" +#include "force.h" #include "math_const.h" #include "memory.h" -#include "error.h" +#include "neigh_list.h" +#include "neighbor.h" +#include +#include using namespace LAMMPS_NS; using namespace MathConst; @@ -743,12 +743,18 @@ void PairLJCutTIP4PCut::compute_newsite(double *xO, double *xH1, void *PairLJCutTIP4PCut::extract(const char *str, int &dim) { dim = 0; + if (strcmp(str,"qdist") == 0) return (void *) &qdist; + if (strcmp(str,"typeO") == 0) return (void *) &typeO; + if (strcmp(str,"typeH") == 0) return (void *) &typeH; + if (strcmp(str,"typeA") == 0) return (void *) &typeA; + if (strcmp(str,"typeB") == 0) return (void *) &typeB; if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul; dim = 2; if (strcmp(str,"epsilon") == 0) return (void *) epsilon; if (strcmp(str,"sigma") == 0) return (void *) sigma; return nullptr; } + /* ---------------------------------------------------------------------- memory usage of hneigh ------------------------------------------------------------------------- */ diff --git a/src/MOLECULE/pair_tip4p_cut.cpp b/src/MOLECULE/pair_tip4p_cut.cpp index 25c9b82468..6d27c1a164 100644 --- a/src/MOLECULE/pair_tip4p_cut.cpp +++ b/src/MOLECULE/pair_tip4p_cut.cpp @@ -562,3 +562,17 @@ double PairTIP4PCut::memory_usage() bytes += (double)2 * nmax * sizeof(double); return bytes; } + +/* ---------------------------------------------------------------------- */ + +void *PairTIP4PCut::extract(const char *str, int &dim) +{ + dim = 0; + if (strcmp(str,"qdist") == 0) return (void *) &qdist; + if (strcmp(str,"typeO") == 0) return (void *) &typeO; + if (strcmp(str,"typeH") == 0) return (void *) &typeH; + if (strcmp(str,"typeA") == 0) return (void *) &typeA; + if (strcmp(str,"typeB") == 0) return (void *) &typeB; + if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul; + return nullptr; +} diff --git a/src/MOLECULE/pair_tip4p_cut.h b/src/MOLECULE/pair_tip4p_cut.h index 85f9e4a72b..73dced7432 100644 --- a/src/MOLECULE/pair_tip4p_cut.h +++ b/src/MOLECULE/pair_tip4p_cut.h @@ -38,6 +38,7 @@ class PairTIP4PCut : public Pair { void write_restart(FILE *) override; void read_restart(FILE *) override; double memory_usage() override; + void *extract(const char *, int &) override; protected: double cut_coul_global; diff --git a/src/fix_efield.cpp b/src/fix_efield.cpp index 3514ad1c34..3046f21ef3 100644 --- a/src/fix_efield.cpp +++ b/src/fix_efield.cpp @@ -42,7 +42,7 @@ FixEfield::FixEfield(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), xstr(nullptr), ystr(nullptr), zstr(nullptr), estr(nullptr), idregion(nullptr), region(nullptr), efield(nullptr) { - if (narg < 6) utils::missing_cmd_args(FLERR, "fix efield", error); + if (narg < 6) utils::missing_cmd_args(FLERR, std::string("fix ") + style, error); dynamic_group_allow = 1; vector_flag = 1; @@ -85,20 +85,23 @@ FixEfield::FixEfield(LAMMPS *lmp, int narg, char **arg) : int iarg = 6; while (iarg < narg) { if (strcmp(arg[iarg], "region") == 0) { - if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix efield region", error); + if (iarg + 2 > narg) + utils::missing_cmd_args(FLERR, std::string("fix ") + style + " region", error); region = domain->get_region_by_id(arg[iarg + 1]); if (!region) error->all(FLERR, "Region {} for fix efield does not exist", arg[iarg + 1]); idregion = utils::strdup(arg[iarg + 1]); iarg += 2; } else if (strcmp(arg[iarg], "energy") == 0) { - if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix efield energy", error); + if (iarg + 2 > narg) + utils::missing_cmd_args(FLERR, std::string("fix ") + style + "energy", error); if (utils::strmatch(arg[iarg + 1], "^v_")) { estr = utils::strdup(arg[iarg + 1] + 2); } else - error->all(FLERR, "Illegal fix efield energy value argument"); + error->all(FLERR, "Unsupported argument for fix {} energy command: {}", style, arg[iarg]); iarg += 2; - } else - error->all(FLERR, "Unknown fix efield keyword: {}", arg[iarg]); + } else { + error->all(FLERR, "Unknown keyword for fix {} command: {}", style, arg[iarg]); + } } force_flag = 0; @@ -140,47 +143,47 @@ void FixEfield::init() qflag = muflag = 0; if (atom->q_flag) qflag = 1; if (atom->mu_flag && atom->torque_flag) muflag = 1; - if (!qflag && !muflag) error->all(FLERR, "Fix efield requires atom attribute q or mu"); + if (!qflag && !muflag) error->all(FLERR, "Fix {} requires atom attribute q or mu", style); // check variables if (xstr) { xvar = input->variable->find(xstr); - if (xvar < 0) error->all(FLERR, "Variable {} for fix efield does not exist", xstr); + if (xvar < 0) error->all(FLERR, "Variable {} for fix {} does not exist", xstr, style); if (input->variable->equalstyle(xvar)) xstyle = EQUAL; else if (input->variable->atomstyle(xvar)) xstyle = ATOM; else - error->all(FLERR, "Variable {} for fix efield is invalid style", xstr); + error->all(FLERR, "Variable {} for fix {} is invalid style", xstr, style); } if (ystr) { yvar = input->variable->find(ystr); - if (yvar < 0) error->all(FLERR, "Variable {} for fix efield does not exist", ystr); + if (yvar < 0) error->all(FLERR, "Variable {} for fix {} does not exist", ystr, style); if (input->variable->equalstyle(yvar)) ystyle = EQUAL; else if (input->variable->atomstyle(yvar)) ystyle = ATOM; else - error->all(FLERR, "Variable {} for fix efield is invalid style", ystr); + error->all(FLERR, "Variable {} for fix {} is invalid style", ystr, style); } if (zstr) { zvar = input->variable->find(zstr); - if (zvar < 0) error->all(FLERR, "Variable {} for fix efield does not exist", zstr); + if (zvar < 0) error->all(FLERR, "Variable {} for fix {} does not exist", zstr, style); if (input->variable->equalstyle(zvar)) zstyle = EQUAL; else if (input->variable->atomstyle(zvar)) zstyle = ATOM; else - error->all(FLERR, "Variable {} for fix efield is invalid style", zstr); + error->all(FLERR, "Variable {} for fix {} is invalid style", zstr, style); } if (estr) { evar = input->variable->find(estr); - if (evar < 0) error->all(FLERR, "Variable {} for fix efield does not exist", estr); + if (evar < 0) error->all(FLERR, "Variable {} for fix {} does not exist", estr, style); if (input->variable->atomstyle(evar)) estyle = ATOM; else - error->all(FLERR, "Variable {} for fix efield is invalid style", estr); + error->all(FLERR, "Variable {} for fix {} is invalid style", estr, style); } else estyle = NONE; @@ -188,7 +191,7 @@ void FixEfield::init() if (idregion) { region = domain->get_region_by_id(idregion); - if (!region) error->all(FLERR, "Region {} for fix efield does not exist", idregion); + if (!region) error->all(FLERR, "Region {} for fix {} does not exist", idregion, style); } if (xstyle == ATOM || ystyle == ATOM || zstyle == ATOM) @@ -199,15 +202,15 @@ void FixEfield::init() varflag = CONSTANT; if (muflag && varflag == ATOM) - error->all(FLERR, "Fix efield with dipoles cannot use atom-style variables"); + error->all(FLERR, "Fix {} with dipoles cannot use atom-style variables", style); if (muflag && update->whichflag == 2 && comm->me == 0) - error->warning(FLERR, "The minimizer does not re-orient dipoles when using fix efield"); + error->warning(FLERR, "The minimizer does not re-orient dipoles when using fix {}", style); if (varflag == CONSTANT && estyle != NONE) - error->all(FLERR, "Cannot use variable energy with constant efield in fix efield"); + error->all(FLERR, "Cannot use variable energy with constant efield in fix {}", style); if ((varflag == EQUAL || varflag == ATOM) && update->whichflag == 2 && estyle == NONE) - error->all(FLERR, "Must use variable energy with fix efield"); + error->all(FLERR, "Must use variable energy with fix {}", style); if (utils::strmatch(update->integrate_style, "^respa")) { ilevel_respa = (dynamic_cast(update->integrate))->nlevels - 1; @@ -254,7 +257,7 @@ void FixEfield::post_force(int vflag) // reallocate efield array if necessary - if (varflag == ATOM && atom->nmax > maxatom) { + if ((varflag == ATOM) && (atom->nmax > maxatom)) { maxatom = atom->nmax; memory->destroy(efield); memory->create(efield, maxatom, 4, "efield:efield"); @@ -272,12 +275,12 @@ void FixEfield::post_force(int vflag) double **x = atom->x; double fx, fy, fz; - double v[6]; + double v[6], unwrap[3]; + ; // constant efield if (varflag == CONSTANT) { - double unwrap[3]; // charge interactions // force = qE, potential energy = F dot x in unwrapped coords diff --git a/src/fix_efield.h b/src/fix_efield.h index 2ebdcbfaed..52c827bb50 100644 --- a/src/fix_efield.h +++ b/src/fix_efield.h @@ -60,8 +60,6 @@ class FixEfield : public Fix { int force_flag; double fsum[4], fsum_all[4]; }; - } // namespace LAMMPS_NS - #endif #endif diff --git a/unittest/force-styles/tests/fix-timestep-efield_tip4p_const.yaml b/unittest/force-styles/tests/fix-timestep-efield_tip4p_const.yaml new file mode 100644 index 0000000000..06d9ac6788 --- /dev/null +++ b/unittest/force-styles/tests/fix-timestep-efield_tip4p_const.yaml @@ -0,0 +1,85 @@ +--- +lammps_version: 8 Feb 2023 +date_generated: Sat Mar 11 18:36:57 2023 +epsilon: 2e-13 +skip_tests: +prerequisites: ! | + atom full + fix efield/tip4p + pair tip4p/cut +pre_commands: ! | + variable newton_pair delete + variable newton_pair index on +post_commands: ! | + pair_style tip4p/cut 5 2 5 1 0.15 8 + pair_coeff * * + set mol 1*3 type 1 + pair_modify compute off + fix move all nve + fix test all efield/tip4p 0.0 -0.1 0.1 +input_file: in.fourmol +natoms: 29 +global_scalar: 4.491995770754199 +global_vector: ! |- + 3 0 3.3306690738754696e-16 -3.3306690738754696e-16 +run_pos: ! |2 + 1 -2.7552415093234162e-01 2.4719487659017276e+00 -1.7659260155592149e-01 + 2 3.0667295195551686e-01 2.9595207006091298e+00 -8.4948873046248308e-01 + 3 -6.9501841828317501e-01 1.2527375852468967e+00 -6.1999248956326514e-01 + 4 -1.5797583064222933e+00 1.4872506539147059e+00 -1.2557393847172438e+00 + 5 -9.0059081341694203e-01 9.3371332191055079e-01 4.0123828894262992e-01 + 6 2.9466004201322965e-01 2.3098154203185278e-01 -1.2847818233212565e+00 + 7 3.3883194064643740e-01 -9.4939443955616154e-03 -2.4648068462973325e+00 + 8 1.1652661823392052e+00 -4.8788892293481145e-01 -6.7051084726427845e-01 + 9 1.3792181599117415e+00 -2.5785144806547511e-01 2.7127116166812570e-01 + 10 2.0205755944752375e+00 -1.4248083043774298e+00 -9.7168038781705335e-01 + 11 1.7908052895357973e+00 -1.9906451925361470e+00 -1.8892937481668459e+00 + 12 3.0048165232344948e+00 -4.9103749217400644e-01 -1.6223692885574461e+00 + 13 4.0522477613596308e+00 -8.9351481123719190e-01 -1.6399777605636197e+00 + 14 2.6043344097441490e+00 -4.1669249548232629e-01 -2.6631499467387134e+00 + 15 2.9678538772616037e+00 5.5587250140887356e-01 -1.2363630117859896e+00 + 16 2.6503290482597039e+00 -2.3948066095776737e+00 3.7940961923343369e-02 + 17 2.2332493442069845e+00 -2.1021976986931690e+00 1.1486155813700776e+00 + 18 2.1368656301495856e+00 3.0159850430539912e+00 -3.5184857372030947e+00 + 19 1.5342191259863922e+00 2.6247080014171758e+00 -4.2375259746262710e+00 + 20 2.7733620959877663e+00 3.6930851299036145e+00 -3.9339459598190301e+00 + 21 4.9045879801671548e+00 -4.0745123882892225e+00 -3.6229125886346081e+00 + 22 4.3617321654820183e+00 -4.2118993857238047e+00 -4.4557088010084982e+00 + 23 5.7381475002461277e+00 -3.5857943422556979e+00 -3.8747045101649773e+00 + 24 2.0686006913404285e+00 3.1531330464503946e+00 3.1540340249574279e+00 + 25 1.3094870938232261e+00 3.2656011662469178e+00 2.5154941909251844e+00 + 26 2.5772659745607247e+00 4.0054947408934565e+00 3.2209131657742485e+00 + 27 -1.9615465740567124e+00 -4.3544512503622474e+00 2.1090363824362317e+00 + 28 -2.7415808583638364e+00 -4.0226893662371292e+00 1.5870458306784345e+00 + 29 -1.3167312777782281e+00 -3.6011902209042037e+00 2.2737597954628499e+00 +run_vel: ! |2 + 1 3.5866086122019378e-03 -1.2154347473708673e-03 -4.2392200362722891e-03 + 2 2.0082524022385501e-03 3.3511787841678979e-03 3.7688429839581980e-03 + 3 6.2192542069958688e-04 6.9301276762936742e-03 2.5708729584915588e-03 + 4 -1.6431765497648366e-03 -2.9960364454768916e-03 -2.8741662319631094e-03 + 5 -4.3853806891633691e-03 -2.6901128825407325e-03 -1.1121197568886760e-03 + 6 8.7148031038297540e-04 2.4954313330053800e-03 -1.2309066226600864e-03 + 7 -1.2568415210119106e-03 4.0807509738615430e-04 -6.9712338524110520e-04 + 8 6.8689620421805110e-04 -3.2780884533380726e-03 5.8272396977482131e-03 + 9 2.3640248488616234e-04 -1.5224658663571485e-03 1.4175297300497880e-03 + 10 1.6082301926220160e-03 2.9963750033118689e-03 -4.0623208158007441e-03 + 11 -1.7960053899953208e-03 -1.9848960612000211e-03 -2.2845602155362049e-03 + 12 7.1704667350825442e-04 -1.1633952760482067e-03 -2.2137502152034124e-03 + 13 3.4262167231009551e-03 4.2909145176285319e-03 -7.4448972505535764e-04 + 14 1.1671343597943564e-03 -4.6984783616368505e-03 -3.7471832948834287e-03 + 15 -3.4509552609576417e-03 -4.2688014665352289e-03 4.1647084114789104e-03 + 16 -1.2920232746597238e-03 1.3898109477904978e-03 3.4080654373162835e-03 + 17 1.5549230284425913e-03 -2.4052541544723738e-04 -1.3250375236016142e-03 + 18 -8.8068025374161475e-04 -7.0649672396011162e-04 -1.9975231159692470e-03 + 19 -2.1408231218136301e-05 -3.2882182579271491e-03 5.3923401875054571e-03 + 20 4.3713031453335733e-03 4.5512662640716187e-03 2.3532301247409359e-03 + 21 -9.2462567997913224e-04 4.1213962901802100e-04 -1.0839649018049458e-03 + 22 -3.7082492847362717e-03 -3.5677336805362207e-03 5.5992552065310119e-03 + 23 5.1261812506730889e-04 -4.8292196668287962e-03 4.4061275828068932e-03 + 24 -8.5477751371944852e-06 7.2834656439669328e-04 -8.2723276308350297e-04 + 25 4.9711072680197884e-03 -5.2985686094745897e-03 3.3843456016457627e-03 + 26 -1.4738757253221410e-03 -3.2290060494108487e-03 4.0272837484112388e-03 + 27 7.3001515368261834e-05 7.6655436532142407e-04 -5.3634821214295335e-04 + 28 -7.6962825083727625e-04 -1.0470804470619696e-03 1.9904866134623045e-03 + 29 -3.3162166593954292e-03 -1.5880017630017177e-03 2.8796190980974416e-03 +... diff --git a/unittest/force-styles/tests/fix-timestep-efield_tip4p_variable.yaml b/unittest/force-styles/tests/fix-timestep-efield_tip4p_variable.yaml new file mode 100644 index 0000000000..e364de7d07 --- /dev/null +++ b/unittest/force-styles/tests/fix-timestep-efield_tip4p_variable.yaml @@ -0,0 +1,90 @@ +--- +lammps_version: 8 Feb 2023 +date_generated: Sat Mar 11 18:52:03 2023 +epsilon: 2e-13 +skip_tests: +prerequisites: ! | + atom full + fix efield/tip4p + pair tip4p/cut +pre_commands: ! | + variable newton_pair delete + variable newton_pair index on +post_commands: ! | + pair_style tip4p/cut 5 2 5 1 0.15 8.0 + pair_coeff * * + set mol 1*3 type 1 + pair_modify compute off + fix init all store/state 0 x y z + variable xforce delete + variable xforce atom f_init[1] + variable zforce delete + variable zforce equal 0.1*step/10.0 + fix move all nve + fix test all efield/tip4p v_xforce 0.0 v_zforce +input_file: in.fourmol +natoms: 29 +global_scalar: 0 +global_vector: ! |- + 3 9.596230365887404 0 -4.440892098500626e-16 +run_pos: ! |2 + 1 -2.7531286450623932e-01 2.4718738394938558e+00 -1.7653750622874620e-01 + 2 3.0682277741052061e-01 2.9595699761707506e+00 -8.4952487409977362e-01 + 3 -6.9499689283442667e-01 1.2527341715135745e+00 -6.1999076705548650e-01 + 4 -1.5799847405798819e+00 1.4872652200436742e+00 -1.2557496756887883e+00 + 5 -9.0071967709813461e-01 9.3372768633609737e-01 4.0122756859740766e-01 + 6 2.9489677035093576e-01 2.3106805370014455e-01 -1.2848415584173214e+00 + 7 3.3855383234893327e-01 -9.5765527451944828e-03 -2.4647480539199424e+00 + 8 1.1643962409868904e+00 -4.8797342487030382e-01 -6.7045108599319103e-01 + 9 1.3799023150910010e+00 -2.5780046066835277e-01 2.7123202628826903e-01 + 10 2.0207975889886196e+00 -1.4247994958278709e+00 -9.7168706713719422e-01 + 11 1.7910648474147319e+00 -1.9906312942531232e+00 -1.8893042032132690e+00 + 12 3.0035425420205333e+00 -4.9107995377108410e-01 -1.6223383760887220e+00 + 13 4.0528201523622283e+00 -8.9349969755442493e-01 -1.6399873319065361e+00 + 14 2.6047057144491697e+00 -4.1667644840131191e-01 -2.6631622713789995e+00 + 15 2.9682779280011724e+00 5.5588692229168801e-01 -1.2363744303219071e+00 + 16 2.6524868604870657e+00 -2.3947127635091343e+00 3.7898084963605735e-02 + 17 2.2314319778419343e+00 -2.1022876263990393e+00 1.1486586486933228e+00 + 18 2.1339901615202725e+00 3.0158791355914092e+00 -3.5183828730903643e+00 + 19 1.5352491765678440e+00 2.6247671297589257e+00 -4.2375921810526602e+00 + 20 2.7752225351538069e+00 3.6931319090244483e+00 -3.9339826175053725e+00 + 21 4.9009803165421442e+00 -4.0745404736656203e+00 -3.6228510538925023e+00 + 22 4.3679206782268585e+00 -4.2117777340324549e+00 -4.4559421729556803e+00 + 23 5.7469587988417672e+00 -3.5858038792020333e+00 -3.8747167803563598e+00 + 24 2.0670666352241178e+00 3.1530720129788077e+00 3.1541041395003924e+00 + 25 1.3109853135785989e+00 3.2657703855053795e+00 2.5153388842970363e+00 + 26 2.5813689868123637e+00 4.0055691627387287e+00 3.2207885803334761e+00 + 27 -1.9601054593669249e+00 -4.3545249904253041e+00 2.1090702834664450e+00 + 28 -2.7460065388962032e+00 -4.0225877485635015e+00 1.5870129241176147e+00 + 29 -1.3183210257363591e+00 -3.6009974737008963e+00 2.2736573716188233e+00 +run_vel: ! |2 + 1 3.7977419185848097e-03 -1.2897669714584856e-03 -4.1950151251349604e-03 + 2 2.1575172581734408e-03 3.3999203430656075e-03 3.7400839475233913e-03 + 3 6.4262848510624609e-04 6.9265031484497313e-03 2.5714895877023557e-03 + 4 -1.8679239811759064e-03 -2.9813578087670117e-03 -2.8820882510653990e-03 + 5 -4.5136571372885112e-03 -2.6758289003017102e-03 -1.1209155388337711e-03 + 6 1.1037671134296492e-03 2.5866105050619246e-03 -1.2789549161066559e-03 + 7 -1.5342861422661686e-03 3.2481824192138846e-04 -6.5106555914167202e-04 + 8 -1.7354144798367277e-04 -3.3718731447888424e-03 5.8810120566168342e-03 + 9 9.1842435073134841e-04 -1.4702773441316094e-03 1.3826250682568268e-03 + 10 1.8250546263153188e-03 3.0026348558799894e-03 -4.0657750868727510e-03 + 11 -1.5361372713327986e-03 -1.9715983670545788e-03 -2.2928062204248456e-03 + 12 -5.2715705670034873e-04 -1.2049239009499940e-03 -2.1898118814179450e-03 + 13 3.9856626950166771e-03 4.3066930552004591e-03 -7.5097419634071252e-04 + 14 1.5331088183544960e-03 -4.6807759413229257e-03 -3.7592753709259140e-03 + 15 -3.0333470556333213e-03 -4.2544006359050701e-03 4.1544659083419620e-03 + 16 8.5032527220227317e-04 1.4960209341250981e-03 3.3943991963756425e-03 + 17 -2.5150987291715812e-04 -3.3871968539122063e-04 -1.3110122226154578e-03 + 18 -3.7194157533833472e-03 -7.8097728820252952e-04 -1.9105012016605539e-03 + 19 9.9357069481510198e-04 -3.2384951918997540e-03 5.3187556083604745e-03 + 20 6.2100808371219670e-03 4.5760237622866403e-03 2.3397927895772242e-03 + 21 -4.4361083278321254e-03 4.3281730793372350e-04 -1.0265740062531993e-03 + 22 2.3759782136059185e-03 -3.4766457991458125e-03 5.2602144818325288e-03 + 23 9.0442546211972809e-03 -5.0028513128573702e-03 4.5160680978488929e-03 + 24 -1.5158245979151304e-03 6.8210911564513992e-04 -7.5228991658315980e-04 + 25 6.4198421777612721e-03 -5.1099391053400440e-03 3.2035258023780651e-03 + 26 2.5717204131090240e-03 -3.2330590784595391e-03 3.9089372482089858e-03 + 27 1.4748018001575585e-03 6.9455594382931529e-04 -5.3514527521761964e-04 + 28 -5.0752470104177567e-03 -9.9571139633829468e-04 2.0584364913419587e-03 + 29 -4.8690861937706323e-03 -1.3519584410798830e-03 2.8068671849969034e-03 +... diff --git a/unittest/force-styles/tests/mol-pair-lj_cut_tip4p_cut.yaml b/unittest/force-styles/tests/mol-pair-lj_cut_tip4p_cut.yaml index 36607e3a22..95080fcdc6 100644 --- a/unittest/force-styles/tests/mol-pair-lj_cut_tip4p_cut.yaml +++ b/unittest/force-styles/tests/mol-pair-lj_cut_tip4p_cut.yaml @@ -1,5 +1,6 @@ --- lammps_version: 17 Feb 2022 +tags: unstable date_generated: Fri Mar 18 22:17:32 2022 epsilon: 5e-10 skip_tests: @@ -21,6 +22,11 @@ pair_coeff: ! | 4 4 0.015 3.1 5 5 0.015 3.1 extract: ! | + qdist 0 + typeO 0 + typeH 0 + typeA 0 + typeB 0 epsilon 2 sigma 2 cut_coul 0 diff --git a/unittest/force-styles/tests/mol-pair-lj_cut_tip4p_long.yaml b/unittest/force-styles/tests/mol-pair-lj_cut_tip4p_long.yaml index 28961cd04e..00ff7f3842 100644 --- a/unittest/force-styles/tests/mol-pair-lj_cut_tip4p_long.yaml +++ b/unittest/force-styles/tests/mol-pair-lj_cut_tip4p_long.yaml @@ -1,5 +1,6 @@ --- lammps_version: 17 Feb 2022 +tags: unstable date_generated: Fri Mar 18 22:17:32 2022 epsilon: 2e-11 skip_tests: diff --git a/unittest/force-styles/tests/mol-pair-lj_cut_tip4p_long_soft.yaml b/unittest/force-styles/tests/mol-pair-lj_cut_tip4p_long_soft.yaml index 29cd06c4a9..55331e09e9 100644 --- a/unittest/force-styles/tests/mol-pair-lj_cut_tip4p_long_soft.yaml +++ b/unittest/force-styles/tests/mol-pair-lj_cut_tip4p_long_soft.yaml @@ -1,5 +1,6 @@ --- lammps_version: 17 Feb 2022 +tags: unstable date_generated: Fri Mar 18 22:17:32 2022 epsilon: 1.0e-12 skip_tests: diff --git a/unittest/force-styles/tests/mol-pair-lj_cut_tip4p_table.yaml b/unittest/force-styles/tests/mol-pair-lj_cut_tip4p_table.yaml index 34a48e8f5f..5e6d460b65 100644 --- a/unittest/force-styles/tests/mol-pair-lj_cut_tip4p_table.yaml +++ b/unittest/force-styles/tests/mol-pair-lj_cut_tip4p_table.yaml @@ -1,5 +1,6 @@ --- lammps_version: 17 Feb 2022 +tags: unstable date_generated: Fri Mar 18 22:17:32 2022 epsilon: 5e-13 skip_tests: gpu diff --git a/unittest/force-styles/tests/mol-pair-lj_long_cut_tip4p_long.yaml b/unittest/force-styles/tests/mol-pair-lj_long_cut_tip4p_long.yaml index ef1fb4f4cd..80b844e7c5 100644 --- a/unittest/force-styles/tests/mol-pair-lj_long_cut_tip4p_long.yaml +++ b/unittest/force-styles/tests/mol-pair-lj_long_cut_tip4p_long.yaml @@ -1,5 +1,6 @@ --- lammps_version: 17 Feb 2022 +tags: unstable date_generated: Fri Mar 18 22:17:33 2022 epsilon: 2.5e-09 skip_tests: diff --git a/unittest/force-styles/tests/mol-pair-lj_long_tip4p_long.yaml b/unittest/force-styles/tests/mol-pair-lj_long_tip4p_long.yaml index 0a4639b1c9..72fc77efa4 100644 --- a/unittest/force-styles/tests/mol-pair-lj_long_tip4p_long.yaml +++ b/unittest/force-styles/tests/mol-pair-lj_long_tip4p_long.yaml @@ -1,6 +1,6 @@ --- lammps_version: 17 Feb 2022 -tags: slow +tags: slow unstable date_generated: Fri Mar 18 22:17:33 2022 epsilon: 2.5e-09 skip_tests: diff --git a/unittest/force-styles/tests/mol-pair-lj_table_tip4p_long.yaml b/unittest/force-styles/tests/mol-pair-lj_table_tip4p_long.yaml index 817a1fcae8..6a20bf7b73 100644 --- a/unittest/force-styles/tests/mol-pair-lj_table_tip4p_long.yaml +++ b/unittest/force-styles/tests/mol-pair-lj_table_tip4p_long.yaml @@ -1,6 +1,6 @@ --- lammps_version: 17 Feb 2022 -tags: slow +tags: slow unstable date_generated: Fri Mar 18 22:17:34 2022 epsilon: 2.5e-09 skip_tests: diff --git a/unittest/force-styles/tests/mol-pair-lj_table_tip4p_table.yaml b/unittest/force-styles/tests/mol-pair-lj_table_tip4p_table.yaml index c64b01d34a..02a58d0e95 100644 --- a/unittest/force-styles/tests/mol-pair-lj_table_tip4p_table.yaml +++ b/unittest/force-styles/tests/mol-pair-lj_table_tip4p_table.yaml @@ -1,6 +1,6 @@ --- lammps_version: 17 Feb 2022 -tags: slow +tags: slow unstable date_generated: Fri Mar 18 22:17:35 2022 epsilon: 2.5e-09 skip_tests: diff --git a/unittest/force-styles/tests/mol-pair-tip4p_cut.yaml b/unittest/force-styles/tests/mol-pair-tip4p_cut.yaml index bcd16233b8..09e96c1a7f 100644 --- a/unittest/force-styles/tests/mol-pair-tip4p_cut.yaml +++ b/unittest/force-styles/tests/mol-pair-tip4p_cut.yaml @@ -1,5 +1,6 @@ --- lammps_version: 17 Feb 2022 +tags: unstable date_generated: Fri Mar 18 22:17:36 2022 epsilon: 5.0e-10 skip_tests: @@ -15,7 +16,13 @@ input_file: in.fourmol pair_style: tip4p/cut 5 2 5 1 0.15 10.0 pair_coeff: ! | * * -extract: ! "" +extract: ! | + qdist 0 + typeO 0 + typeH 0 + typeA 0 + typeB 0 + cut_coul 0 natoms: 29 init_vdwl: 0 init_coul: -128.97793216803515 diff --git a/unittest/force-styles/tests/mol-pair-tip4p_long.yaml b/unittest/force-styles/tests/mol-pair-tip4p_long.yaml index e79a216951..072642c471 100644 --- a/unittest/force-styles/tests/mol-pair-tip4p_long.yaml +++ b/unittest/force-styles/tests/mol-pair-tip4p_long.yaml @@ -1,5 +1,6 @@ --- lammps_version: 17 Feb 2022 +tags: unstable date_generated: Fri Mar 18 22:17:36 2022 epsilon: 5e-13 skip_tests: diff --git a/unittest/force-styles/tests/mol-pair-tip4p_long_soft.yaml b/unittest/force-styles/tests/mol-pair-tip4p_long_soft.yaml index cf96b906c2..ec4590c0f8 100644 --- a/unittest/force-styles/tests/mol-pair-tip4p_long_soft.yaml +++ b/unittest/force-styles/tests/mol-pair-tip4p_long_soft.yaml @@ -1,5 +1,6 @@ --- lammps_version: 17 Feb 2022 +tags: unstable date_generated: Fri Mar 18 22:17:36 2022 epsilon: 1e-13 skip_tests: diff --git a/unittest/force-styles/tests/mol-pair-tip4p_table.yaml b/unittest/force-styles/tests/mol-pair-tip4p_table.yaml index eeef0d30ad..e76ecdf0f0 100644 --- a/unittest/force-styles/tests/mol-pair-tip4p_table.yaml +++ b/unittest/force-styles/tests/mol-pair-tip4p_table.yaml @@ -1,5 +1,6 @@ --- lammps_version: 17 Feb 2022 +tags: unstable date_generated: Fri Mar 18 22:17:36 2022 epsilon: 2.5e-13 skip_tests: