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

Merge branch 'develop' into collected-small-changes

Browse Source
This commit is contained in:
Axel Kohlmeyer
2023-03-16 18:00:36 -04:00
parent d3a23817f1 224b1369a1
commit 5e6004961d
33 changed files with 1771 additions and 102 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
doc/src/Commands_compute.rst
View File

@ -52,6 +52,8 @@ KOKKOS, o = OPENMP, t = OPT.
* :doc:`dilatation/atom <compute_dilatation_atom>`
* :doc:`dipole <compute_dipole>`
* :doc:`dipole/chunk <compute_dipole_chunk>`
* :doc:`dipole/tip4p <compute_dipole>`
* :doc:`dipole/tip4p/chunk <compute_dipole_chunk>`
* :doc:`displace/atom <compute_displace_atom>`
* :doc:`dpd <compute_dpd>`
* :doc:`dpd/atom <compute_dpd_atom>`

1
doc/src/Commands_fix.rst
View File

@ -70,6 +70,7 @@ OPT.
* :doc:`dt/reset (k) <fix_dt_reset>`
* :doc:`edpd/source <fix_dpd_source>`
* :doc:`efield <fix_efield>`
* :doc:`efield/tip4p <fix_efield>`
* :doc:`ehex <fix_ehex>`
* :doc:`electrode/conp (i) <fix_electrode>`
* :doc:`electrode/conq (i) <fix_electrode>`

2
doc/src/compute.rst
View File

@ -206,6 +206,8 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` pag
* :doc:`dilatation/atom <compute_dilatation_atom>` - Peridynamic dilatation for each atom
* :doc:`dipole <compute_dipole>` - dipole vector and total dipole
* :doc:`dipole/chunk <compute_dipole_chunk>` - dipole vector and total dipole for each chunk
* :doc:`dipole/tip4p <compute_dipole>` - dipole vector and total dipole with TIP4P pair style
* :doc:`dipole/tip4p/chunk <compute_dipole_chunk>` - dipole vector and total dipole for each chunk with TIP4P pair style
* :doc:`displace/atom <compute_displace_atom>` - displacement of each atom
* :doc:`dpd <compute_dpd>` - total values of internal conductive energy, internal mechanical energy, chemical energy, and harmonic average of internal temperature
* :doc:`dpd/atom <compute_dpd_atom>` - per-particle values of internal conductive energy, internal mechanical energy, chemical energy, and internal temperature

45
doc/src/compute_dipole.rst
View File

@ -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 <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 <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 <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 <units>` times distance
@ -60,7 +72,12 @@ dipole units (i.e., charge :doc:`units <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 <Build_package>` page for more info.
Compute style *dipole/tip4p* can only be used with tip4p pair styles.
Related commands
""""""""""""""""

47
doc/src/compute_dipole_chunk.rst
View File

@ -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 <compute>` command
* dipole/chunk = style name of this compute command
* style = *dipole/chunk* or *dipole/tip4p/chunk*
* chunkID = ID of :doc:`compute chunk/atom <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 <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 <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 <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
<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 <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 <Build_package>` page for more info.
Compute style *dipole/tip4p/chunk* can only be used with tip4p pair
styles.
Related commands
""""""""""""""""

1
doc/src/fix.rst
View File

@ -222,6 +222,7 @@ accelerated styles exist.
* :doc:`dt/reset <fix_dt_reset>` - reset the timestep based on velocity, forces
* :doc:`edpd/source <fix_dpd_source>` - add heat source to eDPD simulations
* :doc:`efield <fix_efield>` - impose electric field on system
* :doc:`efield/tip4p <fix_efield>` - impose electric field on system with TIP4P molecules
* :doc:`ehex <fix_ehex>` - enhanced heat exchange algorithm
* :doc:`electrode/conp <fix_electrode>` - impose electric potential
* :doc:`electrode/conq <fix_electrode>` - impose total electric charge

96
doc/src/fix_efield.rst
View File

@ -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 <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 <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 <pair_lj_cut_tip4p>`. 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 <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 <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
<restart>`.
The :doc:`fix_modify <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 <thermo_style>`. The default
setting for this fix is :doc:`fix_modify energy no <fix_modify>`.
Note that this energy is a fictitious quantity but is needed so that
the :doc:`minimize <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 <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 <thermo_style>`. The default setting for
this fix is :doc:`fix_modify energy no <fix_modify>`. Note that this
energy is a fictitious quantity but is needed so that the :doc:`minimize
<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 <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 <compute_pressure>` and :doc:`compute
stress/atom <compute_stress_atom>` commands. The former can be
accessed by :doc:`thermodynamic output <thermo_style>`. The default
setting for this fix is :doc:`fix_modify virial no <fix_modify>`.
The :doc:`fix_modify <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 <compute_pressure>` and :doc:`compute stress/atom
<compute_stress_atom>` commands. The former can be accessed by
:doc:`thermodynamic output <thermo_style>`. The default setting for
this fix is :doc:`fix_modify virial no <fix_modify>`.
The :doc:`fix_modify <fix_modify>` *respa* option is supported by this
fix. This allows to set at which level of the :doc:`r-RESPA
<run_style>` integrator the fix adding its forces. Default is the
outermost level.
fix. This allows to set at which level of the :doc:`r-RESPA <run_style>`
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 <Build_package>` page for more info.
Fix style *efield/tip4p* can only be used with tip4p pair styles.
Related commands
""""""""""""""""

6
src/.gitignore vendored
View File

@ -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

220
src/EXTRA-COMPUTE/compute_dipole_tip4p.cpp Normal file
View File

@ -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 <cmath>
#include <cstring>
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);
}

46
src/EXTRA-COMPUTE/compute_dipole_tip4p.h Normal file
View File

@ -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

365
src/EXTRA-COMPUTE/compute_dipole_tip4p_chunk.cpp Normal file
View File

@ -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 <cmath>
#include <cstring>
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<ComputeChunkAtom *>(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<ComputeChunkAtom *>(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);
}

64
src/EXTRA-COMPUTE/compute_dipole_tip4p_chunk.h Normal file
View File

@ -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

644
src/EXTRA-FIX/fix_efield_tip4p.cpp Normal file
View File

@ -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 <cstring>
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);
}

43
src/EXTRA-FIX/fix_efield_tip4p.h Normal file
View File

@ -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

22
src/MOLECULE/pair_lj_cut_tip4p_cut.cpp
View File

@ -18,20 +18,20 @@
#include "pair_lj_cut_tip4p_cut.h"
#include <cmath>
#include <cstring>
#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 <cmath>
#include <cstring>
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
------------------------------------------------------------------------- */

14
src/MOLECULE/pair_tip4p_cut.cpp
View File

@ -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;
}

1
src/MOLECULE/pair_tip4p_cut.h
View File

@ -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;

49
src/fix_efield.cpp
View File

@ -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<Respa *>(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

2
src/fix_efield.h
View File

@ -60,8 +60,6 @@ class FixEfield : public Fix {
int force_flag;
double fsum[4], fsum_all[4];
};
} // namespace LAMMPS_NS
#endif
#endif

85
unittest/force-styles/tests/fix-timestep-efield_tip4p_const.yaml Normal file
View File

@ -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
...

90
unittest/force-styles/tests/fix-timestep-efield_tip4p_variable.yaml Normal file
View File

@ -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
...

6
unittest/force-styles/tests/mol-pair-lj_cut_tip4p_cut.yaml
View File

@ -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

1
unittest/force-styles/tests/mol-pair-lj_cut_tip4p_long.yaml
View File

@ -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:

1
unittest/force-styles/tests/mol-pair-lj_cut_tip4p_long_soft.yaml
View File

@ -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:

1
unittest/force-styles/tests/mol-pair-lj_cut_tip4p_table.yaml
View File

@ -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

1
unittest/force-styles/tests/mol-pair-lj_long_cut_tip4p_long.yaml
View File

@ -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:

2
unittest/force-styles/tests/mol-pair-lj_long_tip4p_long.yaml
View File

@ -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:

2
unittest/force-styles/tests/mol-pair-lj_table_tip4p_long.yaml
View File

@ -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:

2
unittest/force-styles/tests/mol-pair-lj_table_tip4p_table.yaml
View File

@ -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:

9
unittest/force-styles/tests/mol-pair-tip4p_cut.yaml
View File

@ -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

1
unittest/force-styles/tests/mol-pair-tip4p_long.yaml
View File

@ -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:

1
unittest/force-styles/tests/mol-pair-tip4p_long_soft.yaml
View File

@ -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:

1
unittest/force-styles/tests/mol-pair-tip4p_table.yaml
View File

@ -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: