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

Merge pull request #3685 from akohlmey/more-tip4p-support

Browse Source 
Add TIP4P versions for compute dipole, compute dipole/chunk and fix efield
This commit is contained in:
Axel Kohlmeyer
2023-03-16 17:59:58 -04:00
committed by GitHub
parent 6069e8cb14 32c695a7d4
commit 224b1369a1
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:`dilatation/atom <compute_dilatation_atom>`
* :doc:`dipole <compute_dipole>` * :doc:`dipole <compute_dipole>`
* :doc:`dipole/chunk <compute_dipole_chunk>` * :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:`displace/atom <compute_displace_atom>`
* :doc:`dpd <compute_dpd>` * :doc:`dpd <compute_dpd>`
* :doc:`dpd/atom <compute_dpd_atom>` * :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:`dt/reset (k) <fix_dt_reset>`
* :doc:`edpd/source <fix_dpd_source>` * :doc:`edpd/source <fix_dpd_source>`
* :doc:`efield <fix_efield>` * :doc:`efield <fix_efield>`
* :doc:`efield/tip4p <fix_efield>`
* :doc:`ehex <fix_ehex>` * :doc:`ehex <fix_ehex>`
* :doc:`electrode/conp (i) <fix_electrode>` * :doc:`electrode/conp (i) <fix_electrode>`
* :doc:`electrode/conq (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:`dilatation/atom <compute_dilatation_atom>` - Peridynamic dilatation for each atom
* :doc:`dipole <compute_dipole>` - dipole vector and total dipole * :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/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:`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 <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 * :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
.. index:: compute dipole/tip4p
compute dipole command compute dipole command
======================
compute dipole/tip4p command
============================ ============================
Syntax Syntax
@ -8,10 +12,10 @@ Syntax
.. code-block:: LAMMPS .. 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 * 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) * arg = *mass* or *geometry* = use COM or geometric center for charged chunk correction (optional)
Examples Examples
@ -21,6 +25,7 @@ Examples
compute 1 fluid dipole compute 1 fluid dipole
compute dw water dipole geometry compute dw water dipole geometry
compute dw water dipole/tip4p
Description Description
""""""""""" """""""""""
@ -28,13 +33,20 @@ Description
Define a computation that calculates the dipole vector and total dipole Define a computation that calculates the dipole vector and total dipole
for a group of atoms. for a group of atoms.
This compute calculates the x,y,z coordinates of the dipole vector These computes calculate the x,y,z coordinates of the dipole vector and
and the total dipole moment for the atoms in the compute group. the total dipole moment for the atoms in the compute group. This
This includes all effects due to atoms passing through periodic boundaries. includes all effects due to atoms passing through periodic boundaries.
For a group with a net charge the resulting dipole is made position independent For a group with a net charge the resulting dipole is made position
by subtracting the position vector of the center of mass or geometric center independent by subtracting the position vector of the center of mass or
times the net charge from the computed dipole vector. Both per-atom charges geometric center times the net charge from the computed dipole
and per-atom dipole moments, if present, contribute to 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:: .. note::
@ -49,10 +61,10 @@ and per-atom dipole moments, if present, contribute to the computed dipole.
Output info Output info
""""""""""" """""""""""
This compute calculations a global scalar containing the magnitude of These computes calculate a global scalar containing the magnitude of the
the computed dipole moment and a global vector of length 3 with the computed dipole moment and a global vector of length 3 with the dipole
dipole vector. See the :doc:`Howto output <Howto_output>` page for vector. See the :doc:`Howto output <Howto_output>` page for an overview
an overview of LAMMPS output options. of LAMMPS output options.
The computed values are "intensive". The array values will be in The computed values are "intensive". The array values will be in
dipole units (i.e., charge :doc:`units <units>` times distance 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 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 Related commands
"""""""""""""""" """"""""""""""""

47
doc/src/compute_dipole_chunk.rst
View File

@ -1,17 +1,21 @@
.. index:: compute dipole/chunk .. index:: compute dipole/chunk
.. index:: compute dipole/tip4p/chunk
compute dipole/chunk command compute dipole/chunk command
============================ ============================
compute dipole/tip4p/chunk command
==================================
Syntax Syntax
"""""" """"""
.. code-block:: LAMMPS .. 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 * 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 * 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) * 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 details of how chunks can be defined and examples of how they can be
used to measure properties of a system. used to measure properties of a system.
This compute calculates the :math:`(x,y,z)` coordinates of the dipole vector These computes calculate the :math:`(x,y,z)` coordinates of the dipole
and the total dipole moment for each chunk, which includes all effects due vector and the total dipole moment for each chunk, which includes all
to atoms passing through periodic boundaries. For chunks with a net effects due to atoms passing through periodic boundaries. For chunks
charge the resulting dipole is made position independent by subtracting with a net charge the resulting dipole is made position independent by
the position vector of the center of mass or geometric center times the subtracting the position vector of the center of mass or geometric
net charge from the computed dipole vector. Both per-atom charges and center times the net charge from the computed dipole vector. Both
per-atom dipole moments, if present, contribute to the computed dipole. 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 Note that only atoms in the specified group contribute to the
calculation. The :doc:`compute chunk/atom <compute_chunk_atom>` command calculation. The :doc:`compute chunk/atom <compute_chunk_atom>` command
@ -78,12 +89,12 @@ command, for example:
Output info 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 number of chunks *Nchunk* as calculated by the specified :doc:`compute
chunk/atom <compute_chunk_atom>` command. The number of columns is 4 for chunk/atom <compute_chunk_atom>` command. The number of columns is 4
the :math:`(x,y,z)` dipole vector components and the total dipole of each for the :math:`(x,y,z)` dipole vector components and the total dipole of
chunk. These values can be accessed by any command that uses global each chunk. These values can be accessed by any command that uses
array values from a compute as input. See the :doc:`Howto output global array values from a compute as input. See the :doc:`Howto output
<Howto_output>` page for an overview of LAMMPS output options. <Howto_output>` page for an overview of LAMMPS output options.
The array values are "intensive". The array values will be in 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 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 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:`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:`edpd/source <fix_dpd_source>` - add heat source to eDPD simulations
* :doc:`efield <fix_efield>` - impose electric field on system * :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:`ehex <fix_ehex>` - enhanced heat exchange algorithm
* :doc:`electrode/conp <fix_electrode>` - impose electric potential * :doc:`electrode/conp <fix_electrode>` - impose electric potential
* :doc:`electrode/conq <fix_electrode>` - impose total electric charge * :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
.. index:: fix efield/tip4p
fix efield command fix efield command
================== ==================
fix efield/tip4p command
========================
Syntax Syntax
"""""" """"""
.. parsed-literal:: .. 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 * 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) * ex,ey,ez = E-field component values (electric field units)
* any of ex,ey,ez can be a variable (see below) * any of ex,ey,ez can be a variable (see below)
* zero or more keyword/value pairs may be appended to args * 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 1.0 0.0 0.0
fix kick external-field efield 0.0 0.0 v_oscillate fix kick external-field efield 0.0 0.0 v_oscillate
fix kick external-field efield/tip4p 1.0 0.0 0.0
Description 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 external electric field being applied to the system. If the system
contains point-dipoles, also add a torque on the dipoles due to the contains point-dipoles, also add a torque on the dipoles due to the
external electric field. external electric field.
For charges, any of the 3 quantities defining the E-field components .. versionadded:: TBD
can be specified as an equal-style or atom-style
:doc:`variable <variable>`, namely *ex*, *ey*, *ez*\ . If the value is a When the *efield/tip4p* style is used, the E-field will be applied to
variable, it should be specified as v_name, where name is the variable the position of the virtual charge site M of a TIP4P molecule instead of
name. In this case, the variable will be evaluated each timestep, and the oxygen position as it is defined by a corresponding :doc:`TIP4P pair
its value used to determine the E-field component. 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 For point-dipoles, equal-style variables can be used, but atom-style
variables are not currently supported, since they imply a spatial variables are not currently supported, since they imply a spatial
gradient in the electric field which means additional terms with gradient in the electric field which means additional terms with
gradients of the field are required for the force and torque on gradients of the field are required for the force and torque on dipoles.
dipoles.
Equal-style variables can specify formulas with various mathematical Equal-style variables can specify formulas with various mathematical
functions, and include :doc:`thermo_style <thermo_style>` command 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 The *energy* keyword is not allowed if the added field is a constant
vector (ex,ey,ez), with all components defined as numeric constants vector (ex,ey,ez), with all components defined as numeric constants
and not as variables. This is because LAMMPS can compute the energy 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 for each charged particle directly as
+ 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 + .. math::
muz\*ez).
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 The *energy* keyword is optional if the added force is defined with
one or more variables, and if you are performing dynamics via the 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 No information about this fix is written to :doc:`binary restart files
<restart>`. <restart>`.
The :doc:`fix_modify <fix_modify>` *energy* option is supported by The :doc:`fix_modify <fix_modify>` *energy* option is supported by this
this fix to add the potential energy inferred by the added force due fix to add the potential energy inferred by the added force due to the
to the electric field to the global potential energy of the system as electric field to the global potential energy of the system as part of
part of :doc:`thermodynamic output <thermo_style>`. The default :doc:`thermodynamic output <thermo_style>`. The default setting for
setting for this fix is :doc:`fix_modify energy no <fix_modify>`. this fix is :doc:`fix_modify energy no <fix_modify>`. Note that this
Note that this energy is a fictitious quantity but is needed so that energy is a fictitious quantity but is needed so that the :doc:`minimize
the :doc:`minimize <minimize>` command can include the forces added by <minimize>` command can include the forces added by this fix in a
this fix in a consistent manner. I.e. there is a decrease in consistent manner. I.e. there is a decrease in potential energy when
potential energy when atoms move in the direction of the added force atoms move in the direction of the added force due to the electric
due to the electric field. field.
The :doc:`fix_modify <fix_modify>` *virial* option is supported by The :doc:`fix_modify <fix_modify>` *virial* option is supported by this
this fix to add the contribution due to the added forces on atoms to fix to add the contribution due to the added forces on atoms to both the
both the global pressure and per-atom stress of the system via the global pressure and per-atom stress of the system via the :doc:`compute
:doc:`compute pressure <compute_pressure>` and :doc:`compute pressure <compute_pressure>` and :doc:`compute stress/atom
stress/atom <compute_stress_atom>` commands. The former can be <compute_stress_atom>` commands. The former can be accessed by
accessed by :doc:`thermodynamic output <thermo_style>`. The default :doc:`thermodynamic output <thermo_style>`. The default setting for
setting for this fix is :doc:`fix_modify virial no <fix_modify>`. this fix is :doc:`fix_modify virial no <fix_modify>`.
The :doc:`fix_modify <fix_modify>` *respa* option is supported by this 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 fix. This allows to set at which level of the :doc:`r-RESPA <run_style>`
<run_style>` integrator the fix adding its forces. Default is the integrator the fix adding its forces. Default is the outermost level.
outermost level.
This fix computes a global scalar and a global 3-vector of forces, This fix computes a global scalar and a global 3-vector of forces,
which can be accessed by various :doc:`output commands which can be accessed by various :doc:`output commands
@ -169,7 +189,11 @@ the iteration count during the minimization.
Restrictions 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 Related commands
"""""""""""""""" """"""""""""""""

6
src/.gitignore vendored
View File

@ -536,6 +536,10 @@
/compute_damage_atom.h /compute_damage_atom.h
/compute_dilatation_atom.cpp /compute_dilatation_atom.cpp
/compute_dilatation_atom.h /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.cpp
/compute_dpd.h /compute_dpd.h
/compute_dpd_atom.cpp /compute_dpd_atom.cpp
@ -754,6 +758,8 @@
/fix_damping_cundall.h /fix_damping_cundall.h
/fix_dpd_energy.cpp /fix_dpd_energy.cpp
/fix_dpd_energy.h /fix_dpd_energy.h
/fix_efield_tip4p.cpp
/fix_efield_tip4p.h
/fix_electron_stopping.cpp /fix_electron_stopping.cpp
/fix_electron_stopping.h /fix_electron_stopping.h
/fix_electron_stopping_fit.cpp /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 "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 "angle.h"
#include "atom.h"
#include "bond.h" #include "bond.h"
#include "comm.h" #include "comm.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "math_const.h" #include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "neigh_list.h"
#include "neighbor.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst; using namespace MathConst;
@ -743,12 +743,18 @@ void PairLJCutTIP4PCut::compute_newsite(double *xO, double *xH1,
void *PairLJCutTIP4PCut::extract(const char *str, int &dim) void *PairLJCutTIP4PCut::extract(const char *str, int &dim)
{ {
dim = 0; 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; if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul;
dim = 2; dim = 2;
if (strcmp(str,"epsilon") == 0) return (void *) epsilon; if (strcmp(str,"epsilon") == 0) return (void *) epsilon;
if (strcmp(str,"sigma") == 0) return (void *) sigma; if (strcmp(str,"sigma") == 0) return (void *) sigma;
return nullptr; return nullptr;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
memory usage of hneigh 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); bytes += (double)2 * nmax * sizeof(double);
return bytes; 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 write_restart(FILE *) override;
void read_restart(FILE *) override; void read_restart(FILE *) override;
double memory_usage() override; double memory_usage() override;
void *extract(const char *, int &) override;
protected: protected:
double cut_coul_global; 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), Fix(lmp, narg, arg), xstr(nullptr), ystr(nullptr), zstr(nullptr), estr(nullptr),
idregion(nullptr), region(nullptr), efield(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; dynamic_group_allow = 1;
vector_flag = 1; vector_flag = 1;
@ -85,20 +85,23 @@ FixEfield::FixEfield(LAMMPS *lmp, int narg, char **arg) :
int iarg = 6; int iarg = 6;
while (iarg < narg) { while (iarg < narg) {
if (strcmp(arg[iarg], "region") == 0) { 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]); region = domain->get_region_by_id(arg[iarg + 1]);
if (!region) error->all(FLERR, "Region {} for fix efield does not exist", arg[iarg + 1]); if (!region) error->all(FLERR, "Region {} for fix efield does not exist", arg[iarg + 1]);
idregion = utils::strdup(arg[iarg + 1]); idregion = utils::strdup(arg[iarg + 1]);
iarg += 2; iarg += 2;
} else if (strcmp(arg[iarg], "energy") == 0) { } 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_")) { if (utils::strmatch(arg[iarg + 1], "^v_")) {
estr = utils::strdup(arg[iarg + 1] + 2); estr = utils::strdup(arg[iarg + 1] + 2);
} else } else
error->all(FLERR, "Illegal fix efield energy value argument"); error->all(FLERR, "Unsupported argument for fix {} energy command: {}", style, arg[iarg]);
iarg += 2; iarg += 2;
} else } else {
error->all(FLERR, "Unknown fix efield keyword: {}", arg[iarg]); error->all(FLERR, "Unknown keyword for fix {} command: {}", style, arg[iarg]);
}
} }
force_flag = 0; force_flag = 0;
@ -140,47 +143,47 @@ void FixEfield::init()
qflag = muflag = 0; qflag = muflag = 0;
if (atom->q_flag) qflag = 1; if (atom->q_flag) qflag = 1;
if (atom->mu_flag && atom->torque_flag) muflag = 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 // check variables
if (xstr) { if (xstr) {
xvar = input->variable->find(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)) if (input->variable->equalstyle(xvar))
xstyle = EQUAL; xstyle = EQUAL;
else if (input->variable->atomstyle(xvar)) else if (input->variable->atomstyle(xvar))
xstyle = ATOM; xstyle = ATOM;
else 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) { if (ystr) {
yvar = input->variable->find(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)) if (input->variable->equalstyle(yvar))
ystyle = EQUAL; ystyle = EQUAL;
else if (input->variable->atomstyle(yvar)) else if (input->variable->atomstyle(yvar))
ystyle = ATOM; ystyle = ATOM;
else 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) { if (zstr) {
zvar = input->variable->find(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)) if (input->variable->equalstyle(zvar))
zstyle = EQUAL; zstyle = EQUAL;
else if (input->variable->atomstyle(zvar)) else if (input->variable->atomstyle(zvar))
zstyle = ATOM; zstyle = ATOM;
else 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) { if (estr) {
evar = input->variable->find(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)) if (input->variable->atomstyle(evar))
estyle = ATOM; estyle = ATOM;
else else
error->all(FLERR, "Variable {} for fix efield is invalid style", estr); error->all(FLERR, "Variable {} for fix {} is invalid style", estr, style);
} else } else
estyle = NONE; estyle = NONE;
@ -188,7 +191,7 @@ void FixEfield::init()
if (idregion) { if (idregion) {
region = domain->get_region_by_id(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) if (xstyle == ATOM || ystyle == ATOM || zstyle == ATOM)
@ -199,15 +202,15 @@ void FixEfield::init()
varflag = CONSTANT; varflag = CONSTANT;
if (muflag && varflag == ATOM) 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) 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) 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) 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")) { if (utils::strmatch(update->integrate_style, "^respa")) {
ilevel_respa = (dynamic_cast<Respa *>(update->integrate))->nlevels - 1; ilevel_respa = (dynamic_cast<Respa *>(update->integrate))->nlevels - 1;
@ -254,7 +257,7 @@ void FixEfield::post_force(int vflag)
// reallocate efield array if necessary // reallocate efield array if necessary
if (varflag == ATOM && atom->nmax > maxatom) { if ((varflag == ATOM) && (atom->nmax > maxatom)) {
maxatom = atom->nmax; maxatom = atom->nmax;
memory->destroy(efield); memory->destroy(efield);
memory->create(efield, maxatom, 4, "efield:efield"); memory->create(efield, maxatom, 4, "efield:efield");
@ -272,12 +275,12 @@ void FixEfield::post_force(int vflag)
double **x = atom->x; double **x = atom->x;
double fx, fy, fz; double fx, fy, fz;
double v[6]; double v[6], unwrap[3];
;
// constant efield // constant efield
if (varflag == CONSTANT) { if (varflag == CONSTANT) {
double unwrap[3];
// charge interactions // charge interactions
// force = qE, potential energy = F dot x in unwrapped coords // 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; int force_flag;
double fsum[4], fsum_all[4]; double fsum[4], fsum_all[4];
}; };
} // namespace LAMMPS_NS } // namespace LAMMPS_NS
#endif #endif
#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 lammps_version: 17 Feb 2022
tags: unstable
date_generated: Fri Mar 18 22:17:32 2022 date_generated: Fri Mar 18 22:17:32 2022
epsilon: 5e-10 epsilon: 5e-10
skip_tests: skip_tests:
@ -21,6 +22,11 @@ pair_coeff: ! |
4 4 0.015 3.1 4 4 0.015 3.1
5 5 0.015 3.1 5 5 0.015 3.1
extract: ! | extract: ! |
qdist 0
typeO 0
typeH 0
typeA 0
typeB 0
epsilon 2 epsilon 2
sigma 2 sigma 2
cut_coul 0 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 lammps_version: 17 Feb 2022
tags: unstable
date_generated: Fri Mar 18 22:17:32 2022 date_generated: Fri Mar 18 22:17:32 2022
epsilon: 2e-11 epsilon: 2e-11
skip_tests: 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 lammps_version: 17 Feb 2022
tags: unstable
date_generated: Fri Mar 18 22:17:32 2022 date_generated: Fri Mar 18 22:17:32 2022
epsilon: 1.0e-12 epsilon: 1.0e-12
skip_tests: 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 lammps_version: 17 Feb 2022
tags: unstable
date_generated: Fri Mar 18 22:17:32 2022 date_generated: Fri Mar 18 22:17:32 2022
epsilon: 5e-13 epsilon: 5e-13
skip_tests: gpu 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 lammps_version: 17 Feb 2022
tags: unstable
date_generated: Fri Mar 18 22:17:33 2022 date_generated: Fri Mar 18 22:17:33 2022
epsilon: 2.5e-09 epsilon: 2.5e-09
skip_tests: 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 lammps_version: 17 Feb 2022
tags: slow tags: slow unstable
date_generated: Fri Mar 18 22:17:33 2022 date_generated: Fri Mar 18 22:17:33 2022
epsilon: 2.5e-09 epsilon: 2.5e-09
skip_tests: 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 lammps_version: 17 Feb 2022
tags: slow tags: slow unstable
date_generated: Fri Mar 18 22:17:34 2022 date_generated: Fri Mar 18 22:17:34 2022
epsilon: 2.5e-09 epsilon: 2.5e-09
skip_tests: 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 lammps_version: 17 Feb 2022
tags: slow tags: slow unstable
date_generated: Fri Mar 18 22:17:35 2022 date_generated: Fri Mar 18 22:17:35 2022
epsilon: 2.5e-09 epsilon: 2.5e-09
skip_tests: skip_tests:

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

@ -1,5 +1,6 @@
--- ---
lammps_version: 17 Feb 2022 lammps_version: 17 Feb 2022
tags: unstable
date_generated: Fri Mar 18 22:17:36 2022 date_generated: Fri Mar 18 22:17:36 2022
epsilon: 5.0e-10 epsilon: 5.0e-10
skip_tests: skip_tests:
@ -15,7 +16,13 @@ input_file: in.fourmol
pair_style: tip4p/cut 5 2 5 1 0.15 10.0 pair_style: tip4p/cut 5 2 5 1 0.15 10.0
pair_coeff: ! | pair_coeff: ! |
* * * *
extract: ! "" extract: ! |
qdist 0
typeO 0
typeH 0
typeA 0
typeB 0
cut_coul 0
natoms: 29 natoms: 29
init_vdwl: 0 init_vdwl: 0
init_coul: -128.97793216803515 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 lammps_version: 17 Feb 2022
tags: unstable
date_generated: Fri Mar 18 22:17:36 2022 date_generated: Fri Mar 18 22:17:36 2022
epsilon: 5e-13 epsilon: 5e-13
skip_tests: skip_tests:

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

@ -1,5 +1,6 @@
--- ---
lammps_version: 17 Feb 2022 lammps_version: 17 Feb 2022
tags: unstable
date_generated: Fri Mar 18 22:17:36 2022 date_generated: Fri Mar 18 22:17:36 2022
epsilon: 1e-13 epsilon: 1e-13
skip_tests: skip_tests:

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

@ -1,5 +1,6 @@
--- ---
lammps_version: 17 Feb 2022 lammps_version: 17 Feb 2022
tags: unstable
date_generated: Fri Mar 18 22:17:36 2022 date_generated: Fri Mar 18 22:17:36 2022
epsilon: 2.5e-13 epsilon: 2.5e-13
skip_tests: skip_tests: