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

Merge branch 'alphataubio-kokkos-fixes' of https://github.com/alphataubio/lammps-alphataubio into alphataubio-kokkos-fixes

Browse Source
This commit is contained in:
alphataubio
2024-10-24 12:20:14 -04:00
parent ce830fcdbd 56ac4aeedf
commit e1a129fbf2
28 changed files with 822 additions and 709 deletions
Download Patch File Download Diff File Expand all files Collapse all files

41
doc/src/Howto_bioFF.rst
View File

@ -1,5 +1,5 @@
CHARMM, AMBER, COMPASS, and DREIDING force fields
=================================================
CHARMM, AMBER, COMPASS, DREIDING, and OPLS force fields
=======================================================
A compact summary of the concepts, definitions, and properties of
force fields with explicit bonded interactions (like the ones discussed
@ -236,6 +236,40 @@ documentation for the formula it computes.
* :doc:`special_bonds <special_bonds>` dreiding
OPLS
----
OPLS (Optimized Potentials for Liquid Simulations) is a general force
field for atomistic simulation of organic molecules in solvent. It was
developed by the `Jorgensen group
<https://traken.chem.yale.edu/oplsaam.html>`_ at Purdue University and
later at Yale University. Multiple versions of the OPLS parameters
exist for united atom representations (OPLS-UA) and for all-atom
representations (OPLS-AA).
This force field is based on atom types mapped to specific functional
groups in organic and biological molecules. Each atom includes a
static, partial atomic charge reflecting the oxidation state of the
element derived from its bonded neighbors :ref:`(Jorgensen)
<howto-jorgensen>` and computed based on increments determined by the
atom type of the atoms bond to it.
The interaction styles listed below compute force field formulas that
are fully or in part consistent with the OPLS style force fields. See
each command's documentation for the formula it computes. Some are only
compatible with a subset of OPLS interactions.
* :doc:`bond_style <bond_harmonic>` harmonic
* :doc:`angle_style <angle_harmonic>` harmonic
* :doc:`dihedral_style <dihedral_opls>` opls
* :doc:`improper_style <improper_cvff>` cvff
* :doc:`improper_style <improper_fourier>` fourier
* :doc:`improper_style <improper_harmonic>` harmonic
* :doc:`pair_style <pair_lj_cut_coul>` lj/cut/coul/cut
* :doc:`pair_style <pair_lj_cut_coul>` lj/cut/coul/long
* :doc:`pair_modify <pair_modify>` geometric
* :doc:`special_bonds <special_bonds>` lj/coul 0.0 0.0 0.5
----------
.. _Typelabel2:
@ -266,3 +300,6 @@ documentation for the formula it computes.
**(Mayo)** Mayo, Olfason, Goddard III (1990). J Phys Chem, 94, 8897-8909. https://doi.org/10.1021/j100389a010
.. _howto-Jorgensen:
**(Jorgensen)** Jorgensen, Tirado-Rives (1988). J Am Chem Soc, 110, 1657-1666. https://doi.org/10.1021/ja00214a001

67
doc/src/Howto_chunk.rst
View File

@ -58,28 +58,30 @@ chunk ID for an individual atom can also be static (e.g. a molecule
ID), or dynamic (e.g. what spatial bin an atom is in as it moves).
Note that this compute allows the per-atom output of other
:doc:`computes <compute>`, :doc:`fixes <fix>`, and
:doc:`variables <variable>` to be used to define chunk IDs for each
atom. This means you can write your own compute or fix to output a
per-atom quantity to use as chunk ID. See the :doc:`Modify <Modify>`
doc pages for info on how to do this. You can also define a :doc:`per-atom variable <variable>` in the input script that uses a formula to
generate a chunk ID for each atom.
:doc:`computes <compute>`, :doc:`fixes <fix>`, and :doc:`variables
<variable>` to be used to define chunk IDs for each atom. This means
you can write your own compute or fix to output a per-atom quantity to
use as chunk ID. See the :doc:`Modify <Modify>` doc pages for info on
how to do this. You can also define a :doc:`per-atom variable
<variable>` in the input script that uses a formula to generate a chunk
ID for each atom.
Fix ave/chunk command:
----------------------
This fix takes the ID of a :doc:`compute chunk/atom <compute_chunk_atom>` command as input. For each chunk,
it then sums one or more specified per-atom values over the atoms in
each chunk. The per-atom values can be any atom property, such as
velocity, force, charge, potential energy, kinetic energy, stress,
etc. Additional keywords are defined for per-chunk properties like
density and temperature. More generally any per-atom value generated
by other :doc:`computes <compute>`, :doc:`fixes <fix>`, and :doc:`per-atom variables <variable>`, can be summed over atoms in each chunk.
This fix takes the ID of a :doc:`compute chunk/atom
<compute_chunk_atom>` command as input. For each chunk, it then sums
one or more specified per-atom values over the atoms in each chunk. The
per-atom values can be any atom property, such as velocity, force,
charge, potential energy, kinetic energy, stress, etc. Additional
keywords are defined for per-chunk properties like density and
temperature. More generally any per-atom value generated by other
:doc:`computes <compute>`, :doc:`fixes <fix>`, and :doc:`per-atom
variables <variable>`, can be summed over atoms in each chunk.
Similar to other averaging fixes, this fix allows the summed per-chunk
values to be time-averaged in various ways, and output to a file. The
fix produces a global array as output with one row of values per
chunk.
fix produces a global array as output with one row of values per chunk.
Compute \*/chunk commands:
--------------------------
@ -97,17 +99,20 @@ category:
* :doc:`compute torque/chunk <compute_vcm_chunk>`
* :doc:`compute vcm/chunk <compute_vcm_chunk>`
They each take the ID of a :doc:`compute chunk/atom <compute_chunk_atom>` command as input. As their names
indicate, they calculate the center-of-mass, radius of gyration,
moments of inertia, mean-squared displacement, temperature, torque,
and velocity of center-of-mass for each chunk of atoms. The :doc:`compute property/chunk <compute_property_chunk>` command can tally the
count of atoms in each chunk and extract other per-chunk properties.
They each take the ID of a :doc:`compute chunk/atom
<compute_chunk_atom>` command as input. As their names indicate, they
calculate the center-of-mass, radius of gyration, moments of inertia,
mean-squared displacement, temperature, torque, and velocity of
center-of-mass for each chunk of atoms. The :doc:`compute
property/chunk <compute_property_chunk>` command can tally the count of
atoms in each chunk and extract other per-chunk properties.
The reason these various calculations are not part of the :doc:`fix ave/chunk command <fix_ave_chunk>`, is that each requires a more
The reason these various calculations are not part of the :doc:`fix
ave/chunk command <fix_ave_chunk>`, is that each requires a more
complicated operation than simply summing and averaging over per-atom
values in each chunk. For example, many of them require calculation
of a center of mass, which requires summing mass\*position over the
atoms and then dividing by summed mass.
values in each chunk. For example, many of them require calculation of
a center of mass, which requires summing mass\*position over the atoms
and then dividing by summed mass.
All of these computes produce a global vector or global array as
output, with one or more values per chunk. The output can be used in
@ -118,9 +123,10 @@ various ways:
* As input to the :doc:`fix ave/histo <fix_ave_histo>` command to
histogram values across chunks. E.g. a histogram of cluster sizes or
molecule diffusion rates.
* As input to special functions of :doc:`equal-style variables <variable>`, like sum() and max() and ave(). E.g. to
find the largest cluster or fastest diffusing molecule or average
radius-of-gyration of a set of molecules (chunks).
* As input to special functions of :doc:`equal-style variables
<variable>`, like sum() and max() and ave(). E.g. to find the largest
cluster or fastest diffusing molecule or average radius-of-gyration of
a set of molecules (chunks).
Other chunk commands:
---------------------
@ -138,9 +144,10 @@ spatially average per-chunk values calculated by a per-chunk compute.
The :doc:`compute reduce/chunk <compute_reduce_chunk>` command reduces a
peratom value across the atoms in each chunk to produce a value per
chunk. When used with the :doc:`compute chunk/spread/atom <compute_chunk_spread_atom>` command it can
create peratom values that induce a new set of chunks with a second
:doc:`compute chunk/atom <compute_chunk_atom>` command.
chunk. When used with the :doc:`compute chunk/spread/atom
<compute_chunk_spread_atom>` command it can create peratom values that
induce a new set of chunks with a second :doc:`compute chunk/atom
<compute_chunk_atom>` command.
Example calculations with chunks
--------------------------------

2
doc/src/compute_angle_local.rst
View File

@ -78,7 +78,7 @@ system and output the statistics in various ways:
compute 2 all angle/local eng theta v_cos v_cossq set theta t
dump 1 all local 100 tmp.dump c_1[*] c_2[*]
compute 3 all reduce ave c_2[*]
compute 3 all reduce ave c_2[*] inputs local
thermo_style custom step temp press c_3[*]
fix 10 all ave/histo 10 10 100 -1 1 20 c_2[3] mode vector file tmp.histo

2
doc/src/compute_bond_local.rst
View File

@ -139,7 +139,7 @@ output the statistics in various ways:
compute 2 all bond/local engpot dist v_dsq set dist d
dump 1 all local 100 tmp.dump c_1[*] c_2[*]
compute 3 all reduce ave c_2[*]
compute 3 all reduce ave c_2[*] inputs local
thermo_style custom step temp press c_3[*]
fix 10 all ave/histo 10 10 100 0 6 20 c_2[3] mode vector file tmp.histo

2
doc/src/compute_dihedral_local.rst
View File

@ -76,7 +76,7 @@ angle in the system and output the statistics in various ways:
compute 2 all dihedral/local phi v_cos v_cossq set phi p
dump 1 all local 100 tmp.dump c_1[*] c_2[*]
compute 3 all reduce ave c_2[*]
compute 3 all reduce ave c_2[*] inputs local
thermo_style custom step temp press c_3[*]
fix 10 all ave/histo 10 10 100 -1 1 20 c_2[2] mode vector file tmp.histo

10
doc/src/compute_reduce.rst
View File

@ -206,11 +206,13 @@ IDs and the bond stretch will be printed with thermodynamic output.
The *inputs* keyword allows selection of whether all the inputs are
per-atom or local quantities. As noted above, all the inputs must be
the same kind (per-atom or local). Per-atom is the default setting.
If a compute or fix is specified as an input, it must produce per-atom
or local data to match this setting. If it produces both, e.g. for
the same kind (per-atom or local). Per-atom is the default setting. If
a compute or fix is specified as an input, it must produce per-atom or
local data to match this setting. If it produces both, like for example
the :doc:`compute voronoi/atom <compute_voronoi_atom>` command, then
this keyword selects between them.
this keyword selects between them. If a compute *only* produces local
data, like for example the :doc:`compute bond/local command
<compute_bond_local>`, the setting "inputs local" is *required*.
----------

93
doc/src/compute_reduce_chunk.rst
View File

@ -37,55 +37,57 @@ Description
Define a calculation that reduces one or more per-atom vectors into
per-chunk values. This can be useful for diagnostic output. Or when
used in conjunction with the :doc:`compute chunk/spread/atom <compute_chunk_spread_atom>` command it can be
used to create per-atom values that induce a new set of chunks with a
second :doc:`compute chunk/atom <compute_chunk_atom>` command. An
example is given below.
used in conjunction with the :doc:`compute chunk/spread/atom
<compute_chunk_spread_atom>` command it can be used to create per-atom
values that induce a new set of chunks with a second :doc:`compute
chunk/atom <compute_chunk_atom>` command. An example is given below.
In LAMMPS, chunks are collections of atoms defined by a :doc:`compute chunk/atom <compute_chunk_atom>` command, which assigns each atom
to a single chunk (or no chunk). The ID for this command is specified
as chunkID. For example, a single chunk could be the atoms in a
molecule or atoms in a spatial bin. See the :doc:`compute chunk/atom <compute_chunk_atom>` and :doc:`Howto chunk <Howto_chunk>`
doc pages for details of how chunks can be defined and examples of how
they can be used to measure properties of a system.
In LAMMPS, chunks are collections of atoms defined by a :doc:`compute
chunk/atom <compute_chunk_atom>` command, which assigns each atom to a
single chunk (or no chunk). The ID for this command is specified as
chunkID. For example, a single chunk could be the atoms in a molecule
or atoms in a spatial bin. See the :doc:`compute chunk/atom
<compute_chunk_atom>` and :doc:`Howto chunk <Howto_chunk>` doc pages for
details of how chunks can be defined and examples of how they can be
used to measure properties of a system.
For each atom, this compute accesses its chunk ID from the specified
*chunkID* compute. The per-atom value from an input contributes
to a per-chunk value corresponding the the chunk ID.
*chunkID* compute. The per-atom value from an input contributes to a
per-chunk value corresponding the chunk ID.
The reduction operation is specified by the *mode* setting and is
performed over all the per-atom values from the atoms in each chunk.
The *sum* option adds the pre-atom values to a per-chunk total. The
*min* or *max* options find the minimum or maximum value of the
per-atom values for each chunk.
The *sum* option adds the per-atom values to a per-chunk total. The
*min* or *max* options find the minimum or maximum value of the per-atom
values for each chunk.
Note that only atoms in the specified group contribute to the
reduction operation. If the *chunkID* compute returns a 0 for the
chunk ID of an atom (i.e., the atom is not in a chunk defined by the
:doc:`compute chunk/atom <compute_chunk_atom>` command), that atom will
also not contribute to the reduction operation. An input that is a
compute or fix may define its own group which affects the quantities
it returns. For example, a compute with return a zero value for atoms
that are not in the group specified for that compute.
Note that only atoms in the specified group contribute to the reduction
operation. If the *chunkID* compute returns a 0 for the chunk ID of an
atom (i.e., the atom is not in a chunk defined by the :doc:`compute
chunk/atom <compute_chunk_atom>` command), that atom will also not
contribute to the reduction operation. An input that is a compute or
fix may define its own group which affects the quantities it returns.
For example, a compute will return a zero value for atoms that are not
in the group specified for that compute.
Each listed input is operated on independently. Each input can be the
result of a :doc:`compute <compute>` or :doc:`fix <fix>` or the evaluation
of an atom-style :doc:`variable <variable>`.
result of a :doc:`compute <compute>` or :doc:`fix <fix>` or the
evaluation of an atom-style :doc:`variable <variable>`.
Note that for values from a compute or fix, the bracketed index I can
be specified using a wildcard asterisk with the index to effectively
Note that for values from a compute or fix, the bracketed index I can be
specified using a wildcard asterisk with the index to effectively
specify multiple values. This takes the form "\*" or "\*n" or "m\*" or
"m\*n". If :math:`N` is the size of the vector (for *mode* = scalar) or the
number of columns in the array (for *mode* = vector), then an asterisk
with no numeric values means all indices from 1 to :math:`N`. A leading
asterisk means all indices from 1 to n (inclusive). A trailing
asterisk means all indices from n to :math:`N` (inclusive). A middle asterisk
means all indices from m to n (inclusive).
"m\*n". If :math:`N` is the size of the vector (for *mode* = scalar) or
the number of columns in the array (for *mode* = vector), then an
asterisk with no numeric values means all indices from 1 to :math:`N`.
A leading asterisk means all indices from 1 to n (inclusive). A
trailing asterisk means all indices from n to :math:`N` (inclusive). A
middle asterisk means all indices from m to n (inclusive).
Using a wildcard is the same as if the individual columns of the array
had been listed one by one. For example, the following two compute reduce/chunk
commands are equivalent, since the
:doc:`compute property/chunk <compute_property_chunk>` command creates a per-atom
had been listed one by one. For example, the following two compute
reduce/chunk commands are equivalent, since the :doc:`compute
property/chunk <compute_property_chunk>` command creates a per-atom
array with 3 columns:
.. code-block:: LAMMPS
@ -164,13 +166,14 @@ Output info
"""""""""""
This compute calculates a global vector if a single input value is
specified, otherwise a global array is output. The number of columns
in the array is the number of inputs provided. The length of the
vector or the number of vector elements or array rows = the number of
chunks *Nchunk* as calculated by the specified :doc:`compute chunk/atom <compute_chunk_atom>` command. The vector or array can
be accessed by any command that uses global values from a compute as
input. See the :doc:`Howto output <Howto_output>` page for an
overview of LAMMPS output options.
specified, otherwise a global array is output. The number of columns in
the array is the number of inputs provided. The length of the vector or
the number of vector elements or array rows = the number of chunks
*Nchunk* as calculated by the specified :doc:`compute chunk/atom
<compute_chunk_atom>` command. The vector or array can be accessed by
any command that uses global values from a compute as input. See the
:doc:`Howto output <Howto_output>` page for an overview of LAMMPS output
options.
The per-atom values for the vector or each column of the array will be
in whatever :doc:`units <units>` the corresponding input value is in.
@ -183,7 +186,9 @@ Restrictions
Related commands
""""""""""""""""
:doc:`compute chunk/atom <compute_chunk_atom>`, :doc:`compute reduce <compute_reduce>`, :doc:`compute chunk/spread/atom <compute_chunk_spread_atom>`
:doc:`compute chunk/atom <compute_chunk_atom>`,
:doc:`compute reduce <compute_reduce>`,
:doc:`compute chunk/spread/atom <compute_chunk_spread_atom>`
Default
"""""""

3
doc/src/compute_stress_mop.rst
View File

@ -129,6 +129,9 @@ package <Build_package>` doc page on for more info.
The method is implemented for orthogonal simulation boxes whose
size does not change in time, and axis-aligned planes.
Contributions from bonds, angles, and dihedrals are not compatible
with MPI parallel runs.
The method only works with two-body pair interactions, because it
requires the class method ``Pair::single()`` to be implemented, which is
not possible for manybody potentials. In particular, compute

13
doc/src/fix_atom_swap.rst
View File

@ -119,15 +119,14 @@ groups of atoms that have different charges, these charges will not be
changed when the atom types change.
Since this fix computes total potential energies before and after
proposed swaps, so even complicated potential energy calculations are
OK, including the following:
proposed swaps, even complicated potential energy calculations are
acceptable, including the following:
* long-range electrostatics (:math:`k`-space)
* many body pair styles
* hybrid pair styles
* eam pair styles
* hybrid pair styles (with restrictions)
* EAM pair styles
* triclinic systems
* need to include potential energy contributions from other fixes
Some fixes have an associated potential energy. Examples of such fixes
include: :doc:`efield <fix_efield>`, :doc:`gravity <fix_gravity>`,
@ -181,6 +180,10 @@ This fix is part of the MC package. It is only enabled if LAMMPS was
built with that package. See the :doc:`Build package <Build_package>`
doc page for more info.
When this fix is used with a :doc:`hybrid pair style <pair_hybrid>`
system, only swaps between atom types of the same sub-style (or
combination of sub-styles) are permitted.
This fix cannot be used with systems that do not have per-type masses
(e.g. atom style sphere) since the implemented algorithm pre-computes
velocity rescaling factors from per-type masses and ignores any per-atom

2
doc/src/fix_halt.rst
View File

@ -101,7 +101,7 @@ hstyle = bondmax option.
.. code-block:: LAMMPS
compute bdist all bond/local dist
compute bmax all reduce max c_bdist
compute bmax all reduce max c_bdist inputs local
variable bondmax equal c_bmax
Thus these two versions of a fix halt command will do the same thing:

6
doc/src/fix_langevin.rst
View File

@ -231,12 +231,6 @@ the particles. As described below, this energy can then be printed
out or added to the potential energy of the system to monitor energy
conservation.
.. note::
This accumulated energy does NOT include kinetic energy removed
by the *zero* flag. LAMMPS will print a warning when both options are
active.
The keyword *zero* can be used to eliminate drift due to the
thermostat. Because the random forces on different atoms are
independent, they do not sum exactly to zero. As a result, this fix

4
doc/src/group.rst
View File

@ -162,7 +162,7 @@ potential energy is above the threshold value :math:`-3.0`.
compute 1 all pe/atom
compute 2 all reduce sum c_1
thermo_style custom step temp pe c_2
run 0
run 0 post no
variable eatom atom "c_1 > -3.0"
group hienergy variable eatom
@ -173,7 +173,7 @@ Note that these lines
compute 2 all reduce sum c_1
thermo_style custom step temp pe c_2
run 0
run 0 post no
are necessary to ensure that the "eatom" variable is current when the
group command invokes it. Because the eatom variable computes the

2
doc/utils/requirements.txt
View File

@ -1,4 +1,4 @@
Sphinx >= 5.3.0, <9.0
Sphinx >= 5.3.0, <8.2.0
sphinxcontrib-spelling
sphinxcontrib-jquery
sphinx-design

2
doc/utils/sphinx-config/false_positives.txt
View File

@ -1286,6 +1286,7 @@ gcc
gcmc
gdot
GeC
Geesthacht
Geier
gencode
Geocomputing
@ -4219,6 +4220,7 @@ zeeman
Zeeman
Zemer
zenodo
Zentrum
Zepeda
zflag
Zhang

2
src/DIELECTRIC/pppm_disp_dielectric.cpp
View File

@ -715,7 +715,7 @@ void PPPMDispDielectric::fieldforce_c_ad()
{
int i,l,m,n,nx,ny,nz,mx,my,mz;
FFT_SCALAR dx,dy,dz;
FFT_SCALAR ekx,eky,ekz,u;
FFT_SCALAR ekx,eky,ekz;
double s1,s2,s3;
double sf = 0.0;

156
src/EXTRA-COMPUTE/compute_stress_mop.cpp
View File

@ -39,7 +39,7 @@
using namespace LAMMPS_NS;
static constexpr double SMALL = 0.001;
static constexpr double SMALL = 0.001;
enum { X, Y, Z };
enum { TOTAL, CONF, KIN, PAIR, BOND, ANGLE, DIHEDRAL };
@ -63,7 +63,7 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : Compute(
} else if (strcmp(arg[3], "z") == 0) {
dir = Z;
} else
error->all(FLERR, "Illegal compute stress/mop command");
error->all(FLERR, "Illegal compute stress/mop plane direction {}", arg[3]);
// Position of the plane
@ -89,7 +89,7 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : Compute(
error->all(FLERR, "Plane for compute stress/mop is out of bounds");
}
if (pos < (domain->boxlo[dir] + domain->prd_half[dir])) {
if (pos < (domain->boxlo[dir] + domain->prd_half[dir])) {
pos1 = pos + domain->prd[dir];
} else {
pos1 = pos - domain->prd[dir];
@ -133,27 +133,17 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : Compute(
which[nvalues] = ANGLE;
nvalues++;
}
} else if (strcmp(arg[iarg],"dihedral") == 0) {
for (i=0; i<3; i++) {
} else if (strcmp(arg[iarg], "dihedral") == 0) {
for (i = 0; i < 3; i++) {
which[nvalues] = DIHEDRAL;
nvalues++;
}
} else
error->all(FLERR, "Illegal compute stress/mop command"); //break;
error->all(FLERR, "Unknown compute stress/mop keyword {}", arg[iarg]);
iarg++;
}
// Error checks:
// orthogonal simulation box
if (domain->triclinic != 0)
error->all(FLERR, "Compute stress/mop is incompatible with triclinic simulation box");
// 2D and pressure calculation in the Z coordinate
if (domain->dimension == 2 && dir == Z)
error->all(FLERR, "Compute stress/mop is incompatible with Z in 2d system");
// Initialize some variables
values_local = values_global = vector = nullptr;
@ -173,8 +163,8 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : Compute(
memory->create(bond_global, nvalues, "stress/mop:bond_global");
memory->create(angle_local, nvalues, "stress/mop:angle_local");
memory->create(angle_global, nvalues, "stress/mop:angle_global");
memory->create(dihedral_local,nvalues,"stress/mop:dihedral_local");
memory->create(dihedral_global,nvalues,"stress/mop:dihedral_global");
memory->create(dihedral_local, nvalues, "stress/mop:dihedral_local");
memory->create(dihedral_global, nvalues, "stress/mop:dihedral_global");
size_vector = nvalues;
vector_flag = 1;
@ -202,6 +192,16 @@ ComputeStressMop::~ComputeStressMop()
void ComputeStressMop::init()
{
// Error checks:
// 2D and pressure calculation in the Z coordinate
if (domain->dimension == 2 && dir == Z)
error->all(FLERR, "Compute stress/mop is incompatible with Z in 2d system");
// orthogonal simulation box
if (domain->triclinic != 0)
error->all(FLERR, "Compute stress/mop is incompatible with triclinic simulation box");
// Conversion constants
@ -224,11 +224,12 @@ void ComputeStressMop::init()
dt = update->dt;
// Error check
// Error checks:
// Compute stress/mop requires fixed simulation box
if (domain->box_change_size || domain->box_change_shape || domain->deform_flag)
error->all(FLERR, "Compute stress/mop requires a fixed size simulation box");
error->all(FLERR, "Compute stress/mop requires a fixed simulation box geometry");
// This compute requires a pair style with pair_single method implemented
@ -242,34 +243,46 @@ void ComputeStressMop::init()
// issue an error for unimplemented intramolecular potentials or Kspace.
if (force->bond) bondflag = 1;
if (force->bond) {
bondflag = 1;
if (comm->nprocs > 1)
error->one(FLERR, "compute stress/mop with bonds does not (yet) support MPI parallel runs");
}
if (force->angle) {
if (force->angle->born_matrix_enable == 0) {
if ((strcmp(force->angle_style, "zero") != 0) && (strcmp(force->angle_style, "none") != 0))
error->all(FLERR, "compute stress/mop does not account for angle potentials");
error->one(FLERR, "compute stress/mop does not account for angle potentials");
} else {
angleflag = 1;
if (comm->nprocs > 1)
error->one(FLERR,
"compute stress/mop with angles does not (yet) support MPI parallel runs");
}
}
if (force->dihedral) {
if (force->dihedral->born_matrix_enable == 0) {
if ((strcmp(force->dihedral_style, "zero") != 0) &&
(strcmp(force->dihedral_style, "none") != 0))
error->all(FLERR, "compute stress/mop does not account for dihedral potentials");
error->one(FLERR, "compute stress/mop does not account for dihedral potentials");
} else {
dihedralflag = 1;
if (comm->nprocs > 1)
error->one(FLERR,
"compute stress/mop with dihedrals does not (yet) support MPI parallel runs");
}
}
if (force->improper) {
if ((strcmp(force->improper_style, "zero") != 0) &&
(strcmp(force->improper_style, "none") != 0))
error->all(FLERR, "compute stress/mop does not account for improper potentials");
error->one(FLERR, "compute stress/mop does not account for improper potentials");
}
if (force->kspace)
error->warning(FLERR, "compute stress/mop does not account for kspace contributions");
}
// need an occasional half neighbor list
neighbor->add_request(this, NeighConst::REQ_OCCASIONAL);
}
@ -324,11 +337,11 @@ void ComputeStressMop::compute_vector()
//Compute dihedral contribution on separate procs
compute_dihedrals();
} else {
for (int i=0; i<nvalues; i++) dihedral_local[i] = 0.0;
for (int i = 0; i < nvalues; i++) dihedral_local[i] = 0.0;
}
// sum dihedral contribution over all procs
MPI_Allreduce(dihedral_local,dihedral_global,nvalues,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(dihedral_local, dihedral_global, nvalues, MPI_DOUBLE, MPI_SUM, world);
for (int m = 0; m < nvalues; m++) {
vector[m] = values_global[m] + bond_global[m] + angle_global[m] + dihedral_global[m];
@ -425,7 +438,7 @@ void ComputeStressMop::compute_pairs()
values_local[m + 1] += fpair * (xi[1] - xj[1]) / area * nktv2p;
values_local[m + 2] += fpair * (xi[2] - xj[2]) / area * nktv2p;
} else if (((xi[dir] < pos) && (xj[dir] > pos)) ||
((xi[dir] < pos1) && (xj[dir] > pos1))) {
((xi[dir] < pos1) && (xj[dir] > pos1))) {
pair->single(i, j, itype, jtype, rsq, factor_coul, factor_lj, fpair);
values_local[m] -= fpair * (xi[0] - xj[0]) / area * nktv2p;
values_local[m + 1] -= fpair * (xi[1] - xj[1]) / area * nktv2p;
@ -869,9 +882,7 @@ void ComputeStressMop::compute_dihedrals()
double x_atom_4[3] = {0.0, 0.0, 0.0};
// initialization
for (int i = 0; i < nvalues; i++) {
dihedral_local[i] = 0.0;
}
for (int i = 0; i < nvalues; i++) { dihedral_local[i] = 0.0; }
double local_contribution[3] = {0.0, 0.0, 0.0};
for (atom2 = 0; atom2 < nlocal; atom2++) {
@ -889,9 +900,9 @@ void ComputeStressMop::compute_dihedrals()
for (i = 0; i < nd; i++) {
if (molecular == 1) {
if (tag[atom2] != dihedral_atom2[atom2][i]) continue;
atom1 = atom->map(dihedral_atom1[atom2][i]);
atom3 = atom->map(dihedral_atom3[atom2][i]);
atom4 = atom->map(dihedral_atom4[atom2][i]);
atom1 = atom->map(dihedral_atom1[atom2][i]);
atom3 = atom->map(dihedral_atom3[atom2][i]);
atom4 = atom->map(dihedral_atom4[atom2][i]);
} else {
if (tag[atom2] != onemols[imol]->dihedral_atom2[atom2][i]) continue;
tagprev = tag[atom2] - iatom - 1;
@ -943,9 +954,9 @@ void ComputeStressMop::compute_dihedrals()
double tau_right = (x_atom_2[dir] - pos) / (x_atom_2[dir] - x_atom_1[dir]);
double tau_middle = (x_atom_3[dir] - pos) / (x_atom_3[dir] - x_atom_2[dir]);
double tau_left = (x_atom_4[dir] - pos) / (x_atom_4[dir] - x_atom_3[dir]);
bool right_cross = ((tau_right >=0) && (tau_right <= 1));
bool right_cross = ((tau_right >= 0) && (tau_right <= 1));
bool middle_cross = ((tau_middle >= 0) && (tau_middle <= 1));
bool left_cross = ((tau_left >=0) && (tau_left <= 1));
bool left_cross = ((tau_left >= 0) && (tau_left <= 1));
// no bonds crossing the plane
if (!right_cross && !middle_cross && !left_cross) continue;
@ -972,45 +983,45 @@ void ComputeStressMop::compute_dihedrals()
vb3z = x_atom_4[2] - x_atom_3[2];
// c0 calculation
sb1 = 1.0 / (vb1x*vb1x + vb1y*vb1y + vb1z*vb1z);
sb2 = 1.0 / (vb2x*vb2x + vb2y*vb2y + vb2z*vb2z);
sb3 = 1.0 / (vb3x*vb3x + vb3y*vb3y + vb3z*vb3z);
sb1 = 1.0 / (vb1x * vb1x + vb1y * vb1y + vb1z * vb1z);
sb2 = 1.0 / (vb2x * vb2x + vb2y * vb2y + vb2z * vb2z);
sb3 = 1.0 / (vb3x * vb3x + vb3y * vb3y + vb3z * vb3z);
rb1 = sqrt(sb1);
rb3 = sqrt(sb3);
c0 = (vb1x*vb3x + vb1y*vb3y + vb1z*vb3z) * rb1*rb3;
c0 = (vb1x * vb3x + vb1y * vb3y + vb1z * vb3z) * rb1 * rb3;
// 1st and 2nd angle
b1mag2 = vb1x*vb1x + vb1y*vb1y + vb1z*vb1z;
b1mag2 = vb1x * vb1x + vb1y * vb1y + vb1z * vb1z;
b1mag = sqrt(b1mag2);
b2mag2 = vb2x*vb2x + vb2y*vb2y + vb2z*vb2z;
b2mag2 = vb2x * vb2x + vb2y * vb2y + vb2z * vb2z;
b2mag = sqrt(b2mag2);
b3mag2 = vb3x*vb3x + vb3y*vb3y + vb3z*vb3z;
b3mag2 = vb3x * vb3x + vb3y * vb3y + vb3z * vb3z;
b3mag = sqrt(b3mag2);
ctmp = vb1x*vb2x + vb1y*vb2y + vb1z*vb2z;
r12c1 = 1.0 / (b1mag*b2mag);
ctmp = vb1x * vb2x + vb1y * vb2y + vb1z * vb2z;
r12c1 = 1.0 / (b1mag * b2mag);
c1mag = ctmp * r12c1;
ctmp = vb2xm*vb3x + vb2ym*vb3y + vb2zm*vb3z;
r12c2 = 1.0 / (b2mag*b3mag);
ctmp = vb2xm * vb3x + vb2ym * vb3y + vb2zm * vb3z;
r12c2 = 1.0 / (b2mag * b3mag);
c2mag = ctmp * r12c2;
// cos and sin of 2 angles and final c
sin2 = MAX(1.0 - c1mag*c1mag,0.0);
sin2 = MAX(1.0 - c1mag * c1mag, 0.0);
sc1 = sqrt(sin2);
if (sc1 < SMALL) sc1 = SMALL;
sc1 = 1.0/sc1;
sc1 = 1.0 / sc1;
sin2 = MAX(1.0 - c2mag*c2mag,0.0);
sin2 = MAX(1.0 - c2mag * c2mag, 0.0);
sc2 = sqrt(sin2);
if (sc2 < SMALL) sc2 = SMALL;
sc2 = 1.0/sc2;
sc2 = 1.0 / sc2;
s1 = sc1 * sc1;
s2 = sc2 * sc2;
s12 = sc1 * sc2;
c = (c0 + c1mag*c2mag) * s12;
c = (c0 + c1mag * c2mag) * s12;
// error check
if (c > 1.0) c = 1.0;
@ -1020,28 +1031,28 @@ void ComputeStressMop::compute_dihedrals()
double a = dudih;
c = c * a;
s12 = s12 * a;
a11 = c*sb1*s1;
a22 = -sb2 * (2.0*c0*s12 - c*(s1+s2));
a33 = c*sb3*s2;
a12 = -r12c1 * (c1mag*c*s1 + c2mag*s12);
a13 = -rb1*rb3*s12;
a23 = r12c2 * (c2mag*c*s2 + c1mag*s12);
a11 = c * sb1 * s1;
a22 = -sb2 * (2.0 * c0 * s12 - c * (s1 + s2));
a33 = c * sb3 * s2;
a12 = -r12c1 * (c1mag * c * s1 + c2mag * s12);
a13 = -rb1 * rb3 * s12;
a23 = r12c2 * (c2mag * c * s2 + c1mag * s12);
sx2 = a12*vb1x + a22*vb2x + a23*vb3x;
sy2 = a12*vb1y + a22*vb2y + a23*vb3y;
sz2 = a12*vb1z + a22*vb2z + a23*vb3z;
sx2 = a12 * vb1x + a22 * vb2x + a23 * vb3x;
sy2 = a12 * vb1y + a22 * vb2y + a23 * vb3y;
sz2 = a12 * vb1z + a22 * vb2z + a23 * vb3z;
f1[0] = a11*vb1x + a12*vb2x + a13*vb3x;
f1[1] = a11*vb1y + a12*vb2y + a13*vb3y;
f1[2] = a11*vb1z + a12*vb2z + a13*vb3z;
f1[0] = a11 * vb1x + a12 * vb2x + a13 * vb3x;
f1[1] = a11 * vb1y + a12 * vb2y + a13 * vb3y;
f1[2] = a11 * vb1z + a12 * vb2z + a13 * vb3z;
f2[0] = -sx2 - f1[0];
f2[1] = -sy2 - f1[1];
f2[2] = -sz2 - f1[2];
f4[0] = a13*vb1x + a23*vb2x + a33*vb3x;
f4[1] = a13*vb1y + a23*vb2y + a33*vb3y;
f4[2] = a13*vb1z + a23*vb2z + a33*vb3z;
f4[0] = a13 * vb1x + a23 * vb2x + a33 * vb3x;
f4[1] = a13 * vb1y + a23 * vb2y + a33 * vb3y;
f4[2] = a13 * vb1z + a23 * vb2z + a33 * vb3z;
f3[0] = sx2 - f4[0];
f3[1] = sy2 - f4[1];
@ -1106,9 +1117,9 @@ void ComputeStressMop::compute_dihedrals()
df[1] = sgn * (f1[1] + f3[1]);
df[2] = sgn * (f1[2] + f3[2]);
}
local_contribution[0] += df[0]/area*nktv2p;
local_contribution[1] += df[1]/area*nktv2p;
local_contribution[2] += df[2]/area*nktv2p;
local_contribution[0] += df[0] / area * nktv2p;
local_contribution[1] += df[1] / area * nktv2p;
local_contribution[2] += df[2] / area * nktv2p;
}
}
@ -1116,11 +1127,10 @@ void ComputeStressMop::compute_dihedrals()
int m = 0;
while (m < nvalues) {
if ((which[m] == CONF) || (which[m] == TOTAL) || (which[m] == DIHEDRAL)) {
dihedral_local[m] = local_contribution[0];
dihedral_local[m+1] = local_contribution[1];
dihedral_local[m+2] = local_contribution[2];
dihedral_local[m] = local_contribution[0];
dihedral_local[m + 1] = local_contribution[1];
dihedral_local[m + 2] = local_contribution[2];
}
m += 3;
}
}

282
src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp
View File

@ -39,7 +39,7 @@
using namespace LAMMPS_NS;
static constexpr double SMALL = 0.001;
static constexpr double SMALL = 0.001;
enum { X, Y, Z };
enum { TOTAL, CONF, KIN, PAIR, BOND, ANGLE, DIHEDRAL };
@ -64,7 +64,7 @@ ComputeStressMopProfile::ComputeStressMopProfile(LAMMPS *lmp, int narg, char **a
} else if (strcmp(arg[3], "z") == 0) {
dir = Z;
} else
error->all(FLERR, "Illegal compute stress/mop/profile command");
error->all(FLERR, "Illegal compute stress/mop/profile plane direction {}", arg[3]);
// bin parameters
@ -118,30 +118,18 @@ ComputeStressMopProfile::ComputeStressMopProfile(LAMMPS *lmp, int narg, char **a
which[nvalues] = ANGLE;
nvalues++;
}
} else if (strcmp(arg[iarg],"dihedral") == 0) {
for (i=0; i<3; i++) {
} else if (strcmp(arg[iarg], "dihedral") == 0) {
for (i = 0; i < 3; i++) {
which[nvalues] = DIHEDRAL;
nvalues++;
}
} else
error->all(FLERR, "Illegal compute stress/mop/profile command"); //break;
error->all(FLERR, "Unknown compute stress/mop/profile keyword {}", arg[iarg]);
iarg++;
}
// check domain related errors
// 3D only
if (domain->dimension == 2 && dir == Z)
error->all(FLERR, "Compute stress/mop/profile incompatible with Z in 2d system");
// orthogonal simulation box
if (domain->triclinic != 0)
error->all(FLERR, "Compute stress/mop/profile incompatible with triclinic simulation box");
// initialize some variables
// Initialize some variables
nbins = 0;
coord = coordp = nullptr;
@ -192,12 +180,25 @@ ComputeStressMopProfile::~ComputeStressMopProfile()
void ComputeStressMopProfile::init()
{
// conversion constants
// check domain related errors
// 3D only
if (domain->dimension == 2 && dir == Z)
error->all(FLERR, "Compute stress/mop/profile is incompatible with Z in 2d system");
// check for orthogonal simulation box
if (domain->triclinic != 0)
error->all(FLERR, "Compute stress/mop/profile is incompatible with triclinic simulation box");
// Conversion constants
nktv2p = force->nktv2p;
ftm2v = force->ftm2v;
// plane area
// Plane area
if (domain->dimension == 3) {
area = 1;
int i;
@ -208,7 +209,7 @@ void ComputeStressMopProfile::init()
area = (dir == X) ? domain->prd[1] : domain->prd[0];
}
// timestep Value
// Timestep Value
dt = update->dt;
@ -217,46 +218,59 @@ void ComputeStressMopProfile::init()
// Compute stress/mop/profile requires fixed simulation box
if (domain->box_change_size || domain->box_change_shape || domain->deform_flag)
error->all(FLERR, "Compute stress/mop/profile requires a fixed simulation box");
error->all(FLERR, "Compute stress/mop/profile requires a fixed simulation box geometry");
// This compute requires a pair style with pair_single method implemented
if (!force->pair) error->all(FLERR, "No pair style is defined for compute stress/mop/profile");
if (force->pair->single_enable == 0)
error->all(FLERR, "Pair style does not support compute stress/mop/profile");
error->all(FLERR, "Pair style {} does not support compute stress/mop/profile",
force->pair_style);
// Errors
if (comm->me == 0) {
// Compute stress/mop/profile only accounts for pair interactions.
// issue an error if any intramolecular potential or Kspace is defined.
// issue an error for unimplemented intramolecular potentials or Kspace.
if (force->bond) bondflag = 1;
if (force->bond) {
bondflag = 1;
if (comm->nprocs > 1)
error->one(
FLERR,
"compute stress/mop/profile with bonds does not (yet) support MPI parallel runs");
}
if (force->angle) {
if (force->angle->born_matrix_enable == 0) {
if ((strcmp(force->angle_style, "zero") != 0) && (strcmp(force->angle_style, "none") != 0))
error->all(FLERR,"compute stress/mop/profile does not account for angle potentials");
error->one(FLERR, "compute stress/mop/profile does not account for angle potentials");
} else {
angleflag = 1;
if (comm->nprocs > 1)
error->one(
FLERR,
"compute stress/mop/profile with angles does not (yet) support MPI parallel runs");
}
}
if (force->dihedral) {
if (force->dihedral->born_matrix_enable == 0) {
if ((strcmp(force->dihedral_style, "zero") != 0) &&
(strcmp(force->dihedral_style, "none") != 0))
error->all(FLERR, "compute stress/mop/profile does not account for dihedral potentials");
error->one(FLERR, "compute stress/mop/profile does not account for dihedral potentials");
} else {
dihedralflag = 1;
if (comm->nprocs > 1)
error->one(
FLERR,
"compute stress/mop/profile with dihedrals does not (yet) support MPI parallel runs");
}
}
if (force->improper)
if (force->improper) {
if ((strcmp(force->improper_style, "zero") != 0) &&
(strcmp(force->improper_style, "none") != 0))
error->all(FLERR, "compute stress/mop/profile does not account for improper potentials");
error->one(FLERR, "compute stress/mop/profile does not account for improper potentials");
}
if (force->kspace)
error->warning(FLERR, "compute stress/mop/profile does not account for kspace contributions");
}
@ -308,28 +322,26 @@ void ComputeStressMopProfile::compute_array()
compute_angles();
} else {
for (int m = 0; m < nbins; m++) {
for (int i = 0; i < nvalues; i++) {
angle_local[m][i] = 0.0;
}
for (int i = 0; i < nvalues; i++) { angle_local[m][i] = 0.0; }
}
}
// sum angle contribution over all procs
MPI_Allreduce(&angle_local[0][0],&angle_global[0][0],nbins*nvalues,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&angle_local[0][0], &angle_global[0][0], nbins * nvalues, MPI_DOUBLE, MPI_SUM,
world);
if (dihedralflag) {
//Compute dihedral contribution on separate procs
compute_dihedrals();
} else {
for (int m = 0; m < nbins; m++) {
for (int i = 0; i < nvalues; i++) {
dihedral_local[m][i] = 0.0;
}
for (int i = 0; i < nvalues; i++) { dihedral_local[m][i] = 0.0; }
}
}
// sum dihedral contribution over all procs
MPI_Allreduce(&dihedral_local[0][0],&dihedral_global[0][0],nbins*nvalues,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&dihedral_local[0][0], &dihedral_global[0][0], nbins * nvalues, MPI_DOUBLE, MPI_SUM,
world);
for (int ibin = 0; ibin < nbins; ibin++) {
array[ibin][0] = coord[ibin];
@ -337,7 +349,8 @@ void ComputeStressMopProfile::compute_array()
int mo = 1;
int m = 0;
while (m < nvalues) {
array[ibin][m + mo] = values_global[ibin][m] + bond_global[ibin][m] + angle_global[ibin][m] + dihedral_global[ibin][m];
array[ibin][m + mo] = values_global[ibin][m] + bond_global[ibin][m] + angle_global[ibin][m] +
dihedral_global[ibin][m];
m++;
}
}
@ -726,15 +739,12 @@ void ComputeStressMopProfile::compute_angles()
// initialization
for (int m = 0; m < nbins; m++) {
for (int i = 0; i < nvalues; i++) {
angle_local[m][i] = 0.0;
}
for (int i = 0; i < nvalues; i++) { angle_local[m][i] = 0.0; }
local_contribution[m][0] = 0.0;
local_contribution[m][1] = 0.0;
local_contribution[m][2] = 0.0;
}
for (atom2 = 0; atom2 < nlocal; atom2++) {
if (!(mask[atom2] & groupbit)) continue;
@ -765,7 +775,7 @@ void ComputeStressMopProfile::compute_angles()
if (atom3 < 0 || !(mask[atom3] & groupbit)) continue;
if (atype <= 0) continue;
for (int ibin = 0; ibin<nbins; ibin++) {
for (int ibin = 0; ibin < nbins; ibin++) {
double pos = coord[ibin];
// minimum image of atom1 with respect to the plane of interest
@ -800,59 +810,59 @@ void ComputeStressMopProfile::compute_angles()
// check if any bond vector crosses the plane of interest
double tau_right = (x_angle_right[dir] - pos) / (x_angle_right[dir] - x_angle_middle[dir]);
double tau_left = (x_angle_middle[dir] - pos) / (x_angle_middle[dir] - x_angle_left[dir]);
bool right_cross = ((tau_right >=0) && (tau_right <= 1));
bool left_cross = ((tau_left >=0) && (tau_left <= 1));
bool right_cross = ((tau_right >= 0) && (tau_right <= 1));
bool left_cross = ((tau_left >= 0) && (tau_left <= 1));
// no bonds crossing the plane
if (!right_cross && !left_cross) continue;
// compute the cos(theta) of the angle
r1 = sqrt(dx_left[0]*dx_left[0] + dx_left[1]*dx_left[1] + dx_left[2]*dx_left[2]);
r2 = sqrt(dx_right[0]*dx_right[0] + dx_right[1]*dx_right[1] + dx_right[2]*dx_right[2]);
cos_theta = -(dx_right[0]*dx_left[0] + dx_right[1]*dx_left[1] + dx_right[2]*dx_left[2])/(r1*r2);
r1 = sqrt(dx_left[0] * dx_left[0] + dx_left[1] * dx_left[1] + dx_left[2] * dx_left[2]);
r2 =
sqrt(dx_right[0] * dx_right[0] + dx_right[1] * dx_right[1] + dx_right[2] * dx_right[2]);
cos_theta =
-(dx_right[0] * dx_left[0] + dx_right[1] * dx_left[1] + dx_right[2] * dx_left[2]) /
(r1 * r2);
if (cos_theta > 1.0) cos_theta = 1.0;
if (cos_theta > 1.0) cos_theta = 1.0;
if (cos_theta < -1.0) cos_theta = -1.0;
// The method returns derivative with regards to cos(theta)
angle->born_matrix(atype, atom1, atom2, atom3, duang, du2ang);
// only right bond crossing the plane
if (right_cross && !left_cross)
{
if (right_cross && !left_cross) {
double sgn = copysign(1.0, x_angle_right[dir] - pos);
dcos_theta[0] = sgn*(dx_right[0]*cos_theta/r2 + dx_left[0]/r1)/r2;
dcos_theta[1] = sgn*(dx_right[1]*cos_theta/r2 + dx_left[1]/r1)/r2;
dcos_theta[2] = sgn*(dx_right[2]*cos_theta/r2 + dx_left[2]/r1)/r2;
dcos_theta[0] = sgn * (dx_right[0] * cos_theta / r2 + dx_left[0] / r1) / r2;
dcos_theta[1] = sgn * (dx_right[1] * cos_theta / r2 + dx_left[1] / r1) / r2;
dcos_theta[2] = sgn * (dx_right[2] * cos_theta / r2 + dx_left[2] / r1) / r2;
}
// only left bond crossing the plane
if (!right_cross && left_cross)
{
if (!right_cross && left_cross) {
double sgn = copysign(1.0, x_angle_left[dir] - pos);
dcos_theta[0] = -sgn*(dx_left[0]*cos_theta/r1 + dx_right[0]/r2)/r1;
dcos_theta[1] = -sgn*(dx_left[1]*cos_theta/r1 + dx_right[1]/r2)/r1;
dcos_theta[2] = -sgn*(dx_left[2]*cos_theta/r1 + dx_right[2]/r2)/r1;
dcos_theta[0] = -sgn * (dx_left[0] * cos_theta / r1 + dx_right[0] / r2) / r1;
dcos_theta[1] = -sgn * (dx_left[1] * cos_theta / r1 + dx_right[1] / r2) / r1;
dcos_theta[2] = -sgn * (dx_left[2] * cos_theta / r1 + dx_right[2] / r2) / r1;
}
// both bonds crossing the plane
if (right_cross && left_cross)
{
if (right_cross && left_cross) {
// due to right bond
double sgn = copysign(1.0, x_angle_middle[dir] - pos);
dcos_theta[0] = -sgn*(dx_right[0]*cos_theta/r2 + dx_left[0]/r1)/r2;
dcos_theta[1] = -sgn*(dx_right[1]*cos_theta/r2 + dx_left[1]/r1)/r2;
dcos_theta[2] = -sgn*(dx_right[2]*cos_theta/r2 + dx_left[2]/r1)/r2;
dcos_theta[0] = -sgn * (dx_right[0] * cos_theta / r2 + dx_left[0] / r1) / r2;
dcos_theta[1] = -sgn * (dx_right[1] * cos_theta / r2 + dx_left[1] / r1) / r2;
dcos_theta[2] = -sgn * (dx_right[2] * cos_theta / r2 + dx_left[2] / r1) / r2;
// due to left bond
dcos_theta[0] += sgn*(dx_left[0]*cos_theta/r1 + dx_right[0]/r2)/r1;
dcos_theta[1] += sgn*(dx_left[1]*cos_theta/r1 + dx_right[1]/r2)/r1;
dcos_theta[2] += sgn*(dx_left[2]*cos_theta/r1 + dx_right[2]/r2)/r1;
dcos_theta[0] += sgn * (dx_left[0] * cos_theta / r1 + dx_right[0] / r2) / r1;
dcos_theta[1] += sgn * (dx_left[1] * cos_theta / r1 + dx_right[1] / r2) / r1;
dcos_theta[2] += sgn * (dx_left[2] * cos_theta / r1 + dx_right[2] / r2) / r1;
}
// final contribution of the given angle term
local_contribution[ibin][0] += duang*dcos_theta[0]/area*nktv2p;
local_contribution[ibin][1] += duang*dcos_theta[1]/area*nktv2p;
local_contribution[ibin][2] += duang*dcos_theta[2]/area*nktv2p;
local_contribution[ibin][0] += duang * dcos_theta[0] / area * nktv2p;
local_contribution[ibin][1] += duang * dcos_theta[1] / area * nktv2p;
local_contribution[ibin][2] += duang * dcos_theta[2] / area * nktv2p;
}
}
}
@ -863,8 +873,8 @@ void ComputeStressMopProfile::compute_angles()
if (which[m] == CONF || which[m] == TOTAL || which[m] == ANGLE) {
for (int ibin = 0; ibin < nbins; ibin++) {
angle_local[ibin][m] = local_contribution[ibin][0];
angle_local[ibin][m+1] = local_contribution[ibin][1];
angle_local[ibin][m+2] = local_contribution[ibin][2];
angle_local[ibin][m + 1] = local_contribution[ibin][1];
angle_local[ibin][m + 2] = local_contribution[ibin][2];
}
}
m += 3;
@ -917,9 +927,7 @@ void ComputeStressMopProfile::compute_dihedrals()
// initialization
for (int m = 0; m < nbins; m++) {
for (int i = 0; i < nvalues; i++) {
dihedral_local[m][i] = 0.0;
}
for (int i = 0; i < nvalues; i++) { dihedral_local[m][i] = 0.0; }
local_contribution[m][0] = 0.0;
local_contribution[m][1] = 0.0;
local_contribution[m][2] = 0.0;
@ -940,9 +948,9 @@ void ComputeStressMopProfile::compute_dihedrals()
for (i = 0; i < nd; i++) {
if (molecular == 1) {
if (tag[atom2] != dihedral_atom2[atom2][i]) continue;
atom1 = atom->map(dihedral_atom1[atom2][i]);
atom3 = atom->map(dihedral_atom3[atom2][i]);
atom4 = atom->map(dihedral_atom4[atom2][i]);
atom1 = atom->map(dihedral_atom1[atom2][i]);
atom3 = atom->map(dihedral_atom3[atom2][i]);
atom4 = atom->map(dihedral_atom4[atom2][i]);
} else {
if (tag[atom2] != onemols[imol]->dihedral_atom2[atom2][i]) continue;
tagprev = tag[atom2] - iatom - 1;
@ -955,7 +963,7 @@ void ComputeStressMopProfile::compute_dihedrals()
if (atom3 < 0 || !(mask[atom3] & groupbit)) continue;
if (atom4 < 0 || !(mask[atom4] & groupbit)) continue;
for (int ibin = 0; ibin<nbins; ibin++) {
for (int ibin = 0; ibin < nbins; ibin++) {
double pos = coord[ibin];
// minimum image of atom1 with respect to the plane of interest
@ -997,9 +1005,9 @@ void ComputeStressMopProfile::compute_dihedrals()
double tau_right = (x_atom_2[dir] - pos) / (x_atom_2[dir] - x_atom_1[dir]);
double tau_middle = (x_atom_3[dir] - pos) / (x_atom_3[dir] - x_atom_2[dir]);
double tau_left = (x_atom_4[dir] - pos) / (x_atom_4[dir] - x_atom_3[dir]);
bool right_cross = ((tau_right >=0) && (tau_right <= 1));
bool right_cross = ((tau_right >= 0) && (tau_right <= 1));
bool middle_cross = ((tau_middle >= 0) && (tau_middle <= 1));
bool left_cross = ((tau_left >=0) && (tau_left <= 1));
bool left_cross = ((tau_left >= 0) && (tau_left <= 1));
// no bonds crossing the plane
if (!right_cross && !middle_cross && !left_cross) continue;
@ -1026,45 +1034,45 @@ void ComputeStressMopProfile::compute_dihedrals()
vb3z = x_atom_4[2] - x_atom_3[2];
// c0 calculation
sb1 = 1.0 / (vb1x*vb1x + vb1y*vb1y + vb1z*vb1z);
sb2 = 1.0 / (vb2x*vb2x + vb2y*vb2y + vb2z*vb2z);
sb3 = 1.0 / (vb3x*vb3x + vb3y*vb3y + vb3z*vb3z);
sb1 = 1.0 / (vb1x * vb1x + vb1y * vb1y + vb1z * vb1z);
sb2 = 1.0 / (vb2x * vb2x + vb2y * vb2y + vb2z * vb2z);
sb3 = 1.0 / (vb3x * vb3x + vb3y * vb3y + vb3z * vb3z);
rb1 = sqrt(sb1);
rb3 = sqrt(sb3);
c0 = (vb1x*vb3x + vb1y*vb3y + vb1z*vb3z) * rb1*rb3;
c0 = (vb1x * vb3x + vb1y * vb3y + vb1z * vb3z) * rb1 * rb3;
// 1st and 2nd angle
b1mag2 = vb1x*vb1x + vb1y*vb1y + vb1z*vb1z;
b1mag2 = vb1x * vb1x + vb1y * vb1y + vb1z * vb1z;
b1mag = sqrt(b1mag2);
b2mag2 = vb2x*vb2x + vb2y*vb2y + vb2z*vb2z;
b2mag2 = vb2x * vb2x + vb2y * vb2y + vb2z * vb2z;
b2mag = sqrt(b2mag2);
b3mag2 = vb3x*vb3x + vb3y*vb3y + vb3z*vb3z;
b3mag2 = vb3x * vb3x + vb3y * vb3y + vb3z * vb3z;
b3mag = sqrt(b3mag2);
ctmp = vb1x*vb2x + vb1y*vb2y + vb1z*vb2z;
r12c1 = 1.0 / (b1mag*b2mag);
ctmp = vb1x * vb2x + vb1y * vb2y + vb1z * vb2z;
r12c1 = 1.0 / (b1mag * b2mag);
c1mag = ctmp * r12c1;
ctmp = vb2xm*vb3x + vb2ym*vb3y + vb2zm*vb3z;
r12c2 = 1.0 / (b2mag*b3mag);
ctmp = vb2xm * vb3x + vb2ym * vb3y + vb2zm * vb3z;
r12c2 = 1.0 / (b2mag * b3mag);
c2mag = ctmp * r12c2;
// cos and sin of 2 angles and final c
sin2 = MAX(1.0 - c1mag*c1mag,0.0);
sin2 = MAX(1.0 - c1mag * c1mag, 0.0);
sc1 = sqrt(sin2);
if (sc1 < SMALL) sc1 = SMALL;
sc1 = 1.0/sc1;
sc1 = 1.0 / sc1;
sin2 = MAX(1.0 - c2mag*c2mag,0.0);
sin2 = MAX(1.0 - c2mag * c2mag, 0.0);
sc2 = sqrt(sin2);
if (sc2 < SMALL) sc2 = SMALL;
sc2 = 1.0/sc2;
sc2 = 1.0 / sc2;
s1 = sc1 * sc1;
s2 = sc2 * sc2;
s12 = sc1 * sc2;
c = (c0 + c1mag*c2mag) * s12;
c = (c0 + c1mag * c2mag) * s12;
// error check
if (c > 1.0) c = 1.0;
@ -1074,36 +1082,35 @@ void ComputeStressMopProfile::compute_dihedrals()
double a = dudih;
c = c * a;
s12 = s12 * a;
a11 = c*sb1*s1;
a22 = -sb2 * (2.0*c0*s12 - c*(s1+s2));
a33 = c*sb3*s2;
a12 = -r12c1 * (c1mag*c*s1 + c2mag*s12);
a13 = -rb1*rb3*s12;
a23 = r12c2 * (c2mag*c*s2 + c1mag*s12);
a11 = c * sb1 * s1;
a22 = -sb2 * (2.0 * c0 * s12 - c * (s1 + s2));
a33 = c * sb3 * s2;
a12 = -r12c1 * (c1mag * c * s1 + c2mag * s12);
a13 = -rb1 * rb3 * s12;
a23 = r12c2 * (c2mag * c * s2 + c1mag * s12);
sx2 = a12*vb1x + a22*vb2x + a23*vb3x;
sy2 = a12*vb1y + a22*vb2y + a23*vb3y;
sz2 = a12*vb1z + a22*vb2z + a23*vb3z;
sx2 = a12 * vb1x + a22 * vb2x + a23 * vb3x;
sy2 = a12 * vb1y + a22 * vb2y + a23 * vb3y;
sz2 = a12 * vb1z + a22 * vb2z + a23 * vb3z;
f1[0] = a11*vb1x + a12*vb2x + a13*vb3x;
f1[1] = a11*vb1y + a12*vb2y + a13*vb3y;
f1[2] = a11*vb1z + a12*vb2z + a13*vb3z;
f1[0] = a11 * vb1x + a12 * vb2x + a13 * vb3x;
f1[1] = a11 * vb1y + a12 * vb2y + a13 * vb3y;
f1[2] = a11 * vb1z + a12 * vb2z + a13 * vb3z;
f2[0] = -sx2 - f1[0];
f2[1] = -sy2 - f1[1];
f2[2] = -sz2 - f1[2];
f4[0] = a13*vb1x + a23*vb2x + a33*vb3x;
f4[1] = a13*vb1y + a23*vb2y + a33*vb3y;
f4[2] = a13*vb1z + a23*vb2z + a33*vb3z;
f4[0] = a13 * vb1x + a23 * vb2x + a33 * vb3x;
f4[1] = a13 * vb1y + a23 * vb2y + a33 * vb3y;
f4[2] = a13 * vb1z + a23 * vb2z + a33 * vb3z;
f3[0] = sx2 - f4[0];
f3[1] = sy2 - f4[1];
f3[2] = sz2 - f4[2];
// only right bond crossing the plane
if (right_cross && !middle_cross && !left_cross)
{
if (right_cross && !middle_cross && !left_cross) {
double sgn = copysign(1.0, x_atom_1[dir] - pos);
df[0] = sgn * f1[0];
df[1] = sgn * f1[1];
@ -1111,8 +1118,7 @@ void ComputeStressMopProfile::compute_dihedrals()
}
// only middle bond crossing the plane
if (!right_cross && middle_cross && !left_cross)
{
if (!right_cross && middle_cross && !left_cross) {
double sgn = copysign(1.0, x_atom_2[dir] - pos);
df[0] = sgn * (f2[0] + f1[0]);
df[1] = sgn * (f2[1] + f1[1]);
@ -1120,8 +1126,7 @@ void ComputeStressMopProfile::compute_dihedrals()
}
// only left bond crossing the plane
if (!right_cross && !middle_cross && left_cross)
{
if (!right_cross && !middle_cross && left_cross) {
double sgn = copysign(1.0, x_atom_4[dir] - pos);
df[0] = sgn * f4[0];
df[1] = sgn * f4[1];
@ -1129,8 +1134,7 @@ void ComputeStressMopProfile::compute_dihedrals()
}
// only right & middle bonds crossing the plane
if (right_cross && middle_cross && !left_cross)
{
if (right_cross && middle_cross && !left_cross) {
double sgn = copysign(1.0, x_atom_2[dir] - pos);
df[0] = sgn * f2[0];
df[1] = sgn * f2[1];
@ -1138,8 +1142,7 @@ void ComputeStressMopProfile::compute_dihedrals()
}
// only right & left bonds crossing the plane
if (right_cross && !middle_cross && left_cross)
{
if (right_cross && !middle_cross && left_cross) {
double sgn = copysign(1.0, x_atom_1[dir] - pos);
df[0] = sgn * (f1[0] + f4[0]);
df[1] = sgn * (f1[1] + f4[1]);
@ -1147,8 +1150,7 @@ void ComputeStressMopProfile::compute_dihedrals()
}
// only middle & left bonds crossing the plane
if (!right_cross && middle_cross && left_cross)
{
if (!right_cross && middle_cross && left_cross) {
double sgn = copysign(1.0, x_atom_3[dir] - pos);
df[0] = sgn * f3[0];
df[1] = sgn * f3[1];
@ -1156,17 +1158,16 @@ void ComputeStressMopProfile::compute_dihedrals()
}
// all three bonds crossing the plane
if (right_cross && middle_cross && left_cross)
{
if (right_cross && middle_cross && left_cross) {
double sgn = copysign(1.0, x_atom_1[dir] - pos);
df[0] = sgn * (f1[0] + f3[0]);
df[1] = sgn * (f1[1] + f3[1]);
df[2] = sgn * (f1[2] + f3[2]);
}
local_contribution[ibin][0] += df[0]/area*nktv2p;
local_contribution[ibin][1] += df[1]/area*nktv2p;
local_contribution[ibin][2] += df[2]/area*nktv2p;
local_contribution[ibin][0] += df[0] / area * nktv2p;
local_contribution[ibin][1] += df[1] / area * nktv2p;
local_contribution[ibin][2] += df[2] / area * nktv2p;
}
}
}
@ -1177,13 +1178,12 @@ void ComputeStressMopProfile::compute_dihedrals()
if ((which[m] == CONF) || (which[m] == TOTAL) || (which[m] == DIHEDRAL)) {
for (int ibin = 0; ibin < nbins; ibin++) {
dihedral_local[ibin][m] = local_contribution[ibin][0];
dihedral_local[ibin][m+1] = local_contribution[ibin][1];
dihedral_local[ibin][m+2] = local_contribution[ibin][2];
dihedral_local[ibin][m + 1] = local_contribution[ibin][1];
dihedral_local[ibin][m + 2] = local_contribution[ibin][2];
}
}
m += 3;
}
}
/* ----------------------------------------------------------------------
@ -1203,11 +1203,11 @@ void ComputeStressMopProfile::setup_bins()
if ((origin > domain->boxhi[dir]) || (origin < domain->boxlo[dir]))
error->all(FLERR, "Origin of bins for compute stress/mop/profile is out of bounds");
n = static_cast<int> ((origin - boxlo[dir]) * invdelta);
lo = origin - n*delta;
n = static_cast<int>((origin - boxlo[dir]) * invdelta);
lo = origin - n * delta;
n = static_cast<int> ((boxhi[dir] - origin) * invdelta);
hi = origin + n*delta;
n = static_cast<int>((boxhi[dir] - origin) * invdelta);
hi = origin + n * delta;
offset = lo;
nbins = static_cast<int>((hi - lo) * invdelta + 1.5);
@ -1221,8 +1221,8 @@ void ComputeStressMopProfile::setup_bins()
memory->create(bond_global, nbins, nvalues, "stress/mop/profile:bond_global");
memory->create(angle_local, nbins, nvalues, "stress/mop/profile:angle_local");
memory->create(angle_global, nbins, nvalues, "stress/mop/profile:angle_global");
memory->create(dihedral_local,nbins,nvalues,"stress/mop/profile:dihedral_local");
memory->create(dihedral_global,nbins,nvalues,"stress/mop/profile:dihedral_global");
memory->create(dihedral_local, nbins, nvalues, "stress/mop/profile:dihedral_local");
memory->create(dihedral_global, nbins, nvalues, "stress/mop/profile:dihedral_global");
memory->create(local_contribution, nbins, 3, "stress/mop/profile:local_contribution");
// set bin coordinates

23
src/MC/fix_atom_swap.cpp
View File

@ -35,6 +35,7 @@
#include "modify.h"
#include "neighbor.h"
#include "pair.h"
#include "pair_hybrid.h"
#include "random_park.h"
#include "region.h"
#include "update.h"
@ -223,6 +224,28 @@ void FixAtomSwap::init()
error->all(FLERR, "Mu not allowed when not using semi-grand in fix atom/swap command");
}
// check if constraints for hybrid pair styles are fulfilled
if (utils::strmatch(force->pair_style, "^hybrid")) {
auto *hybrid = dynamic_cast<PairHybrid *>(force->pair);
if (hybrid) {
for (int i = 0; i < nswaptypes - 1; ++i) {
int type1 = type_list[i];
for (int j = i + 1; j < nswaptypes; ++j) {
int type2 = type_list[j];
if (hybrid->nmap[type1][type1] != hybrid->nmap[type2][type2])
error->all(FLERR, "Pair {} substyles for atom types {} and {} are not compatible",
force->pair_style, type1, type2);
for (int k = 0; k < hybrid->nmap[type1][type1]; ++k) {
if (hybrid->map[type1][type1][k] != hybrid->map[type2][type2][k])
error->all(FLERR, "Pair {} substyles for atom types {} and {} are not compatible",
force->pair_style, type1, type2);
}
}
}
}
}
// set index and check validity of region
if (idregion) {

156
src/compute_angle_local.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -13,38 +12,40 @@
------------------------------------------------------------------------- */
#include "compute_angle_local.h"
#include <cmath>
#include <cstring>
#include "angle.h"
#include "atom.h"
#include "atom_vec.h"
#include "molecule.h"
#include "update.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "angle.h"
#include "input.h"
#include "variable.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
#include "molecule.h"
#include "update.h"
#include "variable.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace MathConst;
static constexpr int DELTA = 10000;
enum{THETA,ENG,VARIABLE};
enum { THETA, ENG, VARIABLE };
/* ---------------------------------------------------------------------- */
ComputeAngleLocal::ComputeAngleLocal(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg),
bstyle(nullptr), vvar(nullptr), tstr(nullptr), vstr(nullptr), vlocal(nullptr), alocal(nullptr)
Compute(lmp, narg, arg), bstyle(nullptr), vvar(nullptr), tstr(nullptr), vstr(nullptr),
vlocal(nullptr), alocal(nullptr)
{
if (narg < 4) error->all(FLERR,"Illegal compute angle/local command");
if (narg < 4) error->all(FLERR, "Illegal compute angle/local command");
if (atom->avec->angles_allow == 0)
error->all(FLERR,"Compute angle/local used when angles are not allowed");
error->all(FLERR, "Compute angle/local used when angles are not allowed");
local_flag = 1;
@ -52,7 +53,7 @@ ComputeAngleLocal::ComputeAngleLocal(LAMMPS *lmp, int narg, char **arg) :
nvalues = narg - 3;
bstyle = new int[nvalues];
vstr = new char*[nvalues];
vstr = new char *[nvalues];
vvar = new int[nvalues];
nvalues = 0;
@ -61,16 +62,17 @@ ComputeAngleLocal::ComputeAngleLocal(LAMMPS *lmp, int narg, char **arg) :
int iarg;
for (iarg = 3; iarg < narg; iarg++) {
if (strcmp(arg[iarg],"theta") == 0) {
if (strcmp(arg[iarg], "theta") == 0) {
bstyle[nvalues++] = THETA;
tflag = 1;
} else if (strcmp(arg[iarg],"eng") == 0) {
} else if (strcmp(arg[iarg], "eng") == 0) {
bstyle[nvalues++] = ENG;
} else if (strncmp(arg[iarg],"v_",2) == 0) {
} else if (strncmp(arg[iarg], "v_", 2) == 0) {
bstyle[nvalues++] = VARIABLE;
vstr[nvar] = utils::strdup(&arg[iarg][2]);
nvar++;
} else break;
} else
break;
}
// optional args
@ -79,45 +81,46 @@ ComputeAngleLocal::ComputeAngleLocal(LAMMPS *lmp, int narg, char **arg) :
tstr = nullptr;
while (iarg < narg) {
if (strcmp(arg[iarg],"set") == 0) {
if (strcmp(arg[iarg], "set") == 0) {
setflag = 1;
if (iarg+3 > narg) error->all(FLERR,"Illegal compute angle/local command");
if (strcmp(arg[iarg+1],"theta") == 0) {
delete [] tstr;
tstr = utils::strdup(arg[iarg+2]);
if (iarg + 3 > narg) error->all(FLERR, "Illegal compute angle/local command");
if (strcmp(arg[iarg + 1], "theta") == 0) {
delete[] tstr;
tstr = utils::strdup(arg[iarg + 2]);
tflag = 1;
} else error->all(FLERR,"Illegal compute angle/local command");
} else
error->all(FLERR, "Illegal compute angle/local command");
iarg += 3;
} else error->all(FLERR,"Illegal compute angle/local command");
} else
error->all(FLERR, "Illegal compute angle/local command");
}
// error check
if (nvar) {
if (!setflag)
error->all(FLERR,"Compute angle/local variable requires a set variable");
if (!setflag) error->all(FLERR, "Compute angle/local variable requires a set variable");
for (int i = 0; i < nvar; i++) {
vvar[i] = input->variable->find(vstr[i]);
if (vvar[i] < 0)
error->all(FLERR,"Variable name for copute angle/local does not exist");
if (vvar[i] < 0) error->all(FLERR, "Variable name for compute angle/local does not exist");
if (!input->variable->equalstyle(vvar[i]))
error->all(FLERR,"Variable for compute angle/local is invalid style");
error->all(FLERR, "Variable for compute angle/local is invalid style");
}
if (tstr) {
tvar = input->variable->find(tstr);
if (tvar < 0)
error->all(FLERR,"Variable name for compute angle/local does not exist");
if (tvar < 0) error->all(FLERR, "Variable name for compute angle/local does not exist");
if (!input->variable->internalstyle(tvar))
error->all(FLERR,"Variable for compute angle/local is invalid style");
error->all(FLERR, "Variable for compute angle/local is invalid style");
}
} else if (setflag)
error->all(FLERR,"Compute angle/local set with no variable");
error->all(FLERR, "Compute angle/local set with no variable");
// initialize output
if (nvalues == 1) size_local_cols = 0;
else size_local_cols = nvalues;
if (nvalues == 1)
size_local_cols = 0;
else
size_local_cols = nvalues;
nmax = 0;
vlocal = nullptr;
@ -128,12 +131,12 @@ ComputeAngleLocal::ComputeAngleLocal(LAMMPS *lmp, int narg, char **arg) :
ComputeAngleLocal::~ComputeAngleLocal()
{
delete [] bstyle;
for (int i = 0; i < nvar; i++) delete [] vstr[i];
delete [] vstr;
delete [] vvar;
delete[] bstyle;
for (int i = 0; i < nvar; i++) delete[] vstr[i];
delete[] vstr;
delete[] vvar;
delete [] tstr;
delete[] tstr;
memory->destroy(vlocal);
memory->destroy(alocal);
@ -144,19 +147,17 @@ ComputeAngleLocal::~ComputeAngleLocal()
void ComputeAngleLocal::init()
{
if (force->angle == nullptr)
error->all(FLERR,"No angle style is defined for compute angle/local");
error->all(FLERR, "No angle style is defined for compute angle/local");
if (nvar) {
for (int i = 0; i < nvar; i++) {
vvar[i] = input->variable->find(vstr[i]);
if (vvar[i] < 0)
error->all(FLERR,"Variable name for compute angle/local does not exist");
if (vvar[i] < 0) error->all(FLERR, "Variable name for compute angle/local does not exist");
}
if (tstr) {
tvar = input->variable->find(tstr);
if (tvar < 0)
error->all(FLERR,"Variable name for compute angle/local does not exist");
if (tvar < 0) error->all(FLERR, "Variable name for compute angle/local does not exist");
}
}
@ -194,10 +195,10 @@ void ComputeAngleLocal::compute_local()
int ComputeAngleLocal::compute_angles(int flag)
{
int i,m,na,atom1,atom2,atom3,imol,iatom,atype,ivar;
int i, m, na, atom1, atom2, atom3, imol, iatom, atype, ivar;
tagint tagprev;
double delx1,dely1,delz1,delx2,dely2,delz2;
double rsq1,rsq2,r1,r2,c,theta;
double delx1, dely1, delz1, delx2, dely2, delz2;
double rsq1, rsq2, r1, r2, c, theta;
double *ptr;
double **x = atom->x;
@ -224,7 +225,8 @@ int ComputeAngleLocal::compute_angles(int flag)
for (atom2 = 0; atom2 < nlocal; atom2++) {
if (!(mask[atom2] & groupbit)) continue;
if (molecular == Atom::MOLECULAR) na = num_angle[atom2];
if (molecular == Atom::MOLECULAR)
na = num_angle[atom2];
else {
if (molindex[atom2] < 0) continue;
imol = molindex[atom2];
@ -242,8 +244,8 @@ int ComputeAngleLocal::compute_angles(int flag)
if (tag[atom2] != onemols[imol]->angle_atom2[atom2][i]) continue;
atype = onemols[imol]->angle_type[atom2][i];
tagprev = tag[atom2] - iatom - 1;
atom1 = atom->map(onemols[imol]->angle_atom1[atom2][i]+tagprev);
atom3 = atom->map(onemols[imol]->angle_atom3[atom2][i]+tagprev);
atom1 = atom->map(onemols[imol]->angle_atom1[atom2][i] + tagprev);
atom3 = atom->map(onemols[imol]->angle_atom3[atom2][i] + tagprev);
}
if (atom1 < 0 || !(mask[atom1] & groupbit)) continue;
@ -261,50 +263,54 @@ int ComputeAngleLocal::compute_angles(int flag)
delx1 = x[atom1][0] - x[atom2][0];
dely1 = x[atom1][1] - x[atom2][1];
delz1 = x[atom1][2] - x[atom2][2];
domain->minimum_image(delx1,dely1,delz1);
domain->minimum_image(delx1, dely1, delz1);
rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
rsq1 = delx1 * delx1 + dely1 * dely1 + delz1 * delz1;
r1 = sqrt(rsq1);
delx2 = x[atom3][0] - x[atom2][0];
dely2 = x[atom3][1] - x[atom2][1];
delz2 = x[atom3][2] - x[atom2][2];
domain->minimum_image(delx2,dely2,delz2);
domain->minimum_image(delx2, dely2, delz2);
rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
rsq2 = delx2 * delx2 + dely2 * dely2 + delz2 * delz2;
r2 = sqrt(rsq2);
// c = cosine of angle
// theta = angle in radians
c = delx1*delx2 + dely1*dely2 + delz1*delz2;
c /= r1*r2;
c = delx1 * delx2 + dely1 * dely2 + delz1 * delz2;
c /= r1 * r2;
if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0;
theta = acos(c);
}
if (nvalues == 1) ptr = &vlocal[m];
else ptr = alocal[m];
if (nvalues == 1)
ptr = &vlocal[m];
else
ptr = alocal[m];
if (nvar) {
ivar = 0;
if (tstr) input->variable->internal_set(tvar,theta);
if (tstr) input->variable->internal_set(tvar, theta);
}
for (int n = 0; n < nvalues; n++) {
switch (bstyle[n]) {
case THETA:
ptr[n] = 180.0*theta/MY_PI;
break;
case ENG:
if (atype > 0) ptr[n] = angle->single(atype,atom1,atom2,atom3);
else ptr[n] = 0.0;
break;
case VARIABLE:
ptr[n] = input->variable->compute_equal(vvar[ivar]);
ivar++;
break;
case THETA:
ptr[n] = 180.0 * theta / MY_PI;
break;
case ENG:
if (atype > 0)
ptr[n] = angle->single(atype, atom1, atom2, atom3);
else
ptr[n] = 0.0;
break;
case VARIABLE:
ptr[n] = input->variable->compute_equal(vvar[ivar]);
ivar++;
break;
}
}
@ -325,11 +331,11 @@ void ComputeAngleLocal::reallocate(int n)
if (nvalues == 1) {
memory->destroy(vlocal);
memory->create(vlocal,nmax,"angle/local:vector_local");
memory->create(vlocal, nmax, "angle/local:vector_local");
vector_local = vlocal;
} else {
memory->destroy(alocal);
memory->create(alocal,nmax,nvalues,"angle/local:array_local");
memory->create(alocal, nmax, nvalues, "angle/local:array_local");
array_local = alocal;
}
}
@ -340,6 +346,6 @@ void ComputeAngleLocal::reallocate(int n)
double ComputeAngleLocal::memory_usage()
{
double bytes = (double)nmax*nvalues * sizeof(double);
double bytes = (double) nmax * nvalues * sizeof(double);
return bytes;
}

344
src/compute_bond_local.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -35,18 +34,35 @@ using namespace LAMMPS_NS;
static constexpr int DELTA = 10000;
enum{DIST,DX,DY,DZ,VELVIB,OMEGA,ENGTRANS,ENGVIB,ENGROT,ENGPOT,FORCE,FX,FY,FZ,VARIABLE,BN};
enum {
DIST,
DX,
DY,
DZ,
VELVIB,
OMEGA,
ENGTRANS,
ENGVIB,
ENGROT,
ENGPOT,
FORCE,
FX,
FY,
FZ,
VARIABLE,
BN
};
/* ---------------------------------------------------------------------- */
ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg),
bstyle(nullptr), vvar(nullptr), dstr(nullptr), vstr(nullptr), vlocal(nullptr), alocal(nullptr)
Compute(lmp, narg, arg), bstyle(nullptr), vvar(nullptr), dstr(nullptr), vstr(nullptr),
vlocal(nullptr), alocal(nullptr)
{
if (narg < 4) error->all(FLERR,"Illegal compute bond/local command");
if (narg < 4) error->all(FLERR, "Illegal compute bond/local command");
if (atom->avec->bonds_allow == 0)
error->all(FLERR,"Compute bond/local used when bonds are not allowed");
error->all(FLERR, "Compute bond/local used when bonds are not allowed");
local_flag = 1;
comm_forward = 3;
@ -56,7 +72,7 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
nvalues = narg - 3;
bstyle = new int[nvalues];
bindex = new int[nvalues];
vstr = new char*[nvalues];
vstr = new char *[nvalues];
vvar = new int[nvalues];
nvalues = 0;
@ -64,30 +80,45 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
int iarg;
for (iarg = 3; iarg < narg; iarg++) {
if (strcmp(arg[iarg],"dist") == 0) bstyle[nvalues++] = DIST;
else if (strcmp(arg[iarg],"dx") == 0) bstyle[nvalues++] = DX;
else if (strcmp(arg[iarg],"dy") == 0) bstyle[nvalues++] = DY;
else if (strcmp(arg[iarg],"dz") == 0) bstyle[nvalues++] = DZ;
else if (strcmp(arg[iarg],"engpot") == 0) bstyle[nvalues++] = ENGPOT;
else if (strcmp(arg[iarg],"force") == 0) bstyle[nvalues++] = FORCE;
else if (strcmp(arg[iarg],"fx") == 0) bstyle[nvalues++] = FX;
else if (strcmp(arg[iarg],"fy") == 0) bstyle[nvalues++] = FY;
else if (strcmp(arg[iarg],"fz") == 0) bstyle[nvalues++] = FZ;
else if (strcmp(arg[iarg],"engvib") == 0) bstyle[nvalues++] = ENGVIB;
else if (strcmp(arg[iarg],"engrot") == 0) bstyle[nvalues++] = ENGROT;
else if (strcmp(arg[iarg],"engtrans") == 0) bstyle[nvalues++] = ENGTRANS;
else if (strcmp(arg[iarg],"omega") == 0) bstyle[nvalues++] = OMEGA;
else if (strcmp(arg[iarg],"velvib") == 0) bstyle[nvalues++] = VELVIB;
else if (strncmp(arg[iarg],"v_",2) == 0) {
if (strcmp(arg[iarg], "dist") == 0)
bstyle[nvalues++] = DIST;
else if (strcmp(arg[iarg], "dx") == 0)
bstyle[nvalues++] = DX;
else if (strcmp(arg[iarg], "dy") == 0)
bstyle[nvalues++] = DY;
else if (strcmp(arg[iarg], "dz") == 0)
bstyle[nvalues++] = DZ;
else if (strcmp(arg[iarg], "engpot") == 0)
bstyle[nvalues++] = ENGPOT;
else if (strcmp(arg[iarg], "force") == 0)
bstyle[nvalues++] = FORCE;
else if (strcmp(arg[iarg], "fx") == 0)
bstyle[nvalues++] = FX;
else if (strcmp(arg[iarg], "fy") == 0)
bstyle[nvalues++] = FY;
else if (strcmp(arg[iarg], "fz") == 0)
bstyle[nvalues++] = FZ;
else if (strcmp(arg[iarg], "engvib") == 0)
bstyle[nvalues++] = ENGVIB;
else if (strcmp(arg[iarg], "engrot") == 0)
bstyle[nvalues++] = ENGROT;
else if (strcmp(arg[iarg], "engtrans") == 0)
bstyle[nvalues++] = ENGTRANS;
else if (strcmp(arg[iarg], "omega") == 0)
bstyle[nvalues++] = OMEGA;
else if (strcmp(arg[iarg], "velvib") == 0)
bstyle[nvalues++] = VELVIB;
else if (strncmp(arg[iarg], "v_", 2) == 0) {
bstyle[nvalues++] = VARIABLE;
vstr[nvar] = utils::strdup(&arg[iarg][2]);
nvar++;
} else if (utils::strmatch(arg[iarg], "^b\\d+$")) { // b1, b2, b3, ... bN
} else if (utils::strmatch(arg[iarg], "^b\\d+$")) { // b1, b2, b3, ... bN
int n = std::stoi(&arg[iarg][1]);
if (n <= 0) error->all(FLERR, "Invalid keyword {} in compute bond/local command", arg[iarg]);
bstyle[nvalues] = BN;
bindex[nvalues++] = n - 1;
} else break;
} else
break;
}
// optional args
@ -96,40 +127,39 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
dstr = nullptr;
while (iarg < narg) {
if (strcmp(arg[iarg],"set") == 0) {
if (strcmp(arg[iarg], "set") == 0) {
setflag = 1;
if (iarg+3 > narg) utils::missing_cmd_args(FLERR,"compute bond/local set", error);
if (strcmp(arg[iarg+1],"dist") == 0) {
delete [] dstr;
dstr = utils::strdup(arg[iarg+2]);
} else error->all(FLERR,"Unknown compute bond/local set keyword: {}", arg[iarg+2]);
if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "compute bond/local set", error);
if (strcmp(arg[iarg + 1], "dist") == 0) {
delete[] dstr;
dstr = utils::strdup(arg[iarg + 2]);
} else
error->all(FLERR, "Unknown compute bond/local set keyword: {}", arg[iarg + 2]);
iarg += 3;
} else error->all(FLERR,"Unknown compute bond/local keyword: {}", arg[iarg]);
} else
error->all(FLERR, "Unknown compute bond/local keyword: {}", arg[iarg]);
}
// error check
if (nvar) {
if (!setflag)
error->all(FLERR,"Compute bond/local variable requires a set variable");
if (!setflag) error->all(FLERR, "Compute bond/local variable requires a set variable");
for (int i = 0; i < nvar; i++) {
vvar[i] = input->variable->find(vstr[i]);
if (vvar[i] < 0)
error->all(FLERR,"Variable name {} for copute bond/local does not exist", vstr[i]);
error->all(FLERR, "Variable name {} for compute bond/local does not exist", vstr[i]);
if (!input->variable->equalstyle(vvar[i]))
error->all(FLERR,"Variable {} for compute bond/local is invalid style", vstr[i]);
error->all(FLERR, "Variable {} for compute bond/local is invalid style", vstr[i]);
}
if (dstr) {
dvar = input->variable->find(dstr);
if (dvar < 0)
error->all(FLERR,"Variable name for compute bond/local does not exist");
if (dvar < 0) error->all(FLERR, "Variable name for compute bond/local does not exist");
if (!input->variable->internalstyle(dvar))
error->all(FLERR,"Variable for compute bond/local is invalid style");
error->all(FLERR, "Variable for compute bond/local is invalid style");
}
} else if (setflag)
error->all(FLERR,"Compute bond/local set used with without a variable");
error->all(FLERR, "Compute bond/local set used with without a variable");
// set singleflag if need to call bond->single()
// set velflag if compute any quantities based on velocities
@ -137,16 +167,20 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
singleflag = 0;
velflag = 0;
for (int i = 0; i < nvalues; i++) {
if (bstyle[i] == ENGPOT || bstyle[i] == FORCE || bstyle[i] == FX ||
bstyle[i] == FY || bstyle[i] == FZ || bstyle[i] == BN) singleflag = 1;
if (bstyle[i] == VELVIB || bstyle[i] == OMEGA || bstyle[i] == ENGTRANS ||
bstyle[i] == ENGVIB || bstyle[i] == ENGROT) velflag = 1;
if (bstyle[i] == ENGPOT || bstyle[i] == FORCE || bstyle[i] == FX || bstyle[i] == FY ||
bstyle[i] == FZ || bstyle[i] == BN)
singleflag = 1;
if (bstyle[i] == VELVIB || bstyle[i] == OMEGA || bstyle[i] == ENGTRANS || bstyle[i] == ENGVIB ||
bstyle[i] == ENGROT)
velflag = 1;
}
// initialize output
if (nvalues == 1) size_local_cols = 0;
else size_local_cols = nvalues;
if (nvalues == 1)
size_local_cols = 0;
else
size_local_cols = nvalues;
nmax = 0;
vlocal = nullptr;
@ -157,13 +191,13 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
ComputeBondLocal::~ComputeBondLocal()
{
delete [] bstyle;
delete [] bindex;
for (int i = 0; i < nvar; i++) delete [] vstr[i];
delete [] vstr;
delete [] vvar;
delete[] bstyle;
delete[] bindex;
for (int i = 0; i < nvar; i++) delete[] vstr[i];
delete[] vstr;
delete[] vvar;
delete [] dstr;
delete[] dstr;
memory->destroy(vlocal);
memory->destroy(alocal);
@ -173,8 +207,7 @@ ComputeBondLocal::~ComputeBondLocal()
void ComputeBondLocal::init()
{
if (force->bond == nullptr)
error->all(FLERR,"No bond style is defined for compute bond/local");
if (force->bond == nullptr) error->all(FLERR, "No bond style is defined for compute bond/local");
for (int i = 0; i < nvalues; i++)
if (bstyle[i] == BN && bindex[i] >= force->bond->single_extra)
@ -183,21 +216,21 @@ void ComputeBondLocal::init()
if (nvar) {
for (int i = 0; i < nvar; i++) {
vvar[i] = input->variable->find(vstr[i]);
if (vvar[i] < 0)
error->all(FLERR,"Variable name for compute bond/local does not exist");
if (vvar[i] < 0) error->all(FLERR, "Variable name for compute bond/local does not exist");
}
if (dstr) {
dvar = input->variable->find(dstr);
if (dvar < 0)
error->all(FLERR,"Variable name for compute bond/local does not exist");
if (dvar < 0) error->all(FLERR, "Variable name for compute bond/local does not exist");
}
}
// set ghostvelflag if need to acquire ghost atom velocities
if (velflag && !comm->ghost_velocity) ghostvelflag = 1;
else ghostvelflag = 0;
if (velflag && !comm->ghost_velocity)
ghostvelflag = 1;
else
ghostvelflag = 0;
// do initial memory allocation so that memory_usage() is correct
@ -236,17 +269,17 @@ void ComputeBondLocal::compute_local()
int ComputeBondLocal::compute_bonds(int flag)
{
int i,m,nb,atom1,atom2,imol,iatom,btype,ivar;
int i, m, nb, atom1, atom2, imol, iatom, btype, ivar;
tagint tagprev;
double dx,dy,dz,rsq;
double mass1,mass2,masstotal,invmasstotal;
double xcm[3],vcm[3];
double delr1[3],delr2[3],delv1[3],delv2[3];
double r12[3],vpar1,vpar2;
double vvib,vrotsq;
double inertia,omegasq;
double dx, dy, dz, rsq;
double mass1, mass2, masstotal, invmasstotal;
double xcm[3], vcm[3];
double delr1[3], delr2[3], delv1[3], delv2[3];
double r12[3], vpar1, vpar2;
double vvib, vrotsq;
double inertia, omegasq;
double mvv2e;
double engpot,engtrans,engvib,engrot,fbond;
double engpot, engtrans, engvib, engrot, fbond;
double *ptr;
double **x = atom->x;
@ -280,7 +313,8 @@ int ComputeBondLocal::compute_bonds(int flag)
for (atom1 = 0; atom1 < nlocal; atom1++) {
if (!(mask[atom1] & groupbit)) continue;
if (molecular == Atom::MOLECULAR) nb = num_bond[atom1];
if (molecular == Atom::MOLECULAR)
nb = num_bond[atom1];
else {
if (molindex[atom1] < 0) continue;
imol = molindex[atom1];
@ -295,7 +329,7 @@ int ComputeBondLocal::compute_bonds(int flag)
} else {
tagprev = tag[atom1] - iatom - 1;
btype = onemols[imol]->bond_type[iatom][i];
atom2 = atom->map(onemols[imol]->bond_atom[iatom][i]+tagprev);
atom2 = atom->map(onemols[imol]->bond_atom[iatom][i] + tagprev);
}
if (atom2 < 0 || !(mask[atom2] & groupbit)) continue;
@ -310,55 +344,56 @@ int ComputeBondLocal::compute_bonds(int flag)
dx = x[atom1][0] - x[atom2][0];
dy = x[atom1][1] - x[atom2][1];
dz = x[atom1][2] - x[atom2][2];
domain->minimum_image(dx,dy,dz);
rsq = dx*dx + dy*dy + dz*dz;
domain->minimum_image(dx, dy, dz);
rsq = dx * dx + dy * dy + dz * dz;
if (btype == 0) {
fbond = 0.0;
} else {
if (singleflag) engpot = bond->single(btype,rsq,atom1,atom2,fbond);
else fbond = engpot = 0.0;
if (singleflag)
engpot = bond->single(btype, rsq, atom1, atom2, fbond);
else
fbond = engpot = 0.0;
engvib = engrot = engtrans = omegasq = vvib = 0.0;
if (velflag) {
if (rmass) {
mass1 = rmass[atom1];
mass2 = rmass[atom2];
}
else {
} else {
mass1 = mass[type[atom1]];
mass2 = mass[type[atom2]];
}
masstotal = mass1+mass2;
masstotal = mass1 + mass2;
invmasstotal = 1.0 / (masstotal);
xcm[0] = (mass1*x[atom1][0] + mass2*x[atom2][0]) * invmasstotal;
xcm[1] = (mass1*x[atom1][1] + mass2*x[atom2][1]) * invmasstotal;
xcm[2] = (mass1*x[atom1][2] + mass2*x[atom2][2]) * invmasstotal;
vcm[0] = (mass1*v[atom1][0] + mass2*v[atom2][0]) * invmasstotal;
vcm[1] = (mass1*v[atom1][1] + mass2*v[atom2][1]) * invmasstotal;
vcm[2] = (mass1*v[atom1][2] + mass2*v[atom2][2]) * invmasstotal;
xcm[0] = (mass1 * x[atom1][0] + mass2 * x[atom2][0]) * invmasstotal;
xcm[1] = (mass1 * x[atom1][1] + mass2 * x[atom2][1]) * invmasstotal;
xcm[2] = (mass1 * x[atom1][2] + mass2 * x[atom2][2]) * invmasstotal;
vcm[0] = (mass1 * v[atom1][0] + mass2 * v[atom2][0]) * invmasstotal;
vcm[1] = (mass1 * v[atom1][1] + mass2 * v[atom2][1]) * invmasstotal;
vcm[2] = (mass1 * v[atom1][2] + mass2 * v[atom2][2]) * invmasstotal;
engtrans = 0.5 * masstotal * MathExtra::lensq3(vcm);
// r12 = unit bond vector from atom1 to atom2
MathExtra::sub3(x[atom2],x[atom1],r12);
MathExtra::sub3(x[atom2], x[atom1], r12);
MathExtra::norm3(r12);
// delr = vector from COM to each atom
// delv = velocity of each atom relative to COM
MathExtra::sub3(x[atom1],xcm,delr1);
MathExtra::sub3(x[atom2],xcm,delr2);
MathExtra::sub3(v[atom1],vcm,delv1);
MathExtra::sub3(v[atom2],vcm,delv2);
MathExtra::sub3(x[atom1], xcm, delr1);
MathExtra::sub3(x[atom2], xcm, delr2);
MathExtra::sub3(v[atom1], vcm, delv1);
MathExtra::sub3(v[atom2], vcm, delv2);
// vpar = component of delv parallel to bond vector
vpar1 = MathExtra::dot3(delv1,r12);
vpar2 = MathExtra::dot3(delv2,r12);
engvib = 0.5 * (mass1*vpar1*vpar1 + mass2*vpar2*vpar2);
vpar1 = MathExtra::dot3(delv1, r12);
vpar2 = MathExtra::dot3(delv2, r12);
engvib = 0.5 * (mass1 * vpar1 * vpar1 + mass2 * vpar2 * vpar2);
// vvib = relative velocity of 2 atoms along bond direction
// vvib < 0 for 2 atoms moving towards each other
@ -369,9 +404,8 @@ int ComputeBondLocal::compute_bonds(int flag)
// vrotsq = tangential speed squared of atom1 only
// omegasq = omega squared, and is the same for atom1 and atom2
inertia = mass1*MathExtra::lensq3(delr1) +
mass2*MathExtra::lensq3(delr2);
vrotsq = MathExtra::lensq3(delv1) - vpar1*vpar1;
inertia = mass1 * MathExtra::lensq3(delr1) + mass2 * MathExtra::lensq3(delr2);
vrotsq = MathExtra::lensq3(delv1) - vpar1 * vpar1;
omegasq = vrotsq / MathExtra::lensq3(delr1);
engrot = 0.5 * inertia * omegasq;
@ -384,12 +418,14 @@ int ComputeBondLocal::compute_bonds(int flag)
engrot *= mvv2e;
}
if (nvalues == 1) ptr = &vlocal[m];
else ptr = alocal[m];
if (nvalues == 1)
ptr = &vlocal[m];
else
ptr = alocal[m];
if (nvar) {
ivar = 0;
if (dstr) input->variable->internal_set(dvar,sqrt(rsq));
if (dstr) input->variable->internal_set(dvar, sqrt(rsq));
}
// to make sure dx, dy and dz are always from the lower to the higher id
@ -397,55 +433,55 @@ int ComputeBondLocal::compute_bonds(int flag)
for (int n = 0; n < nvalues; n++) {
switch (bstyle[n]) {
case DIST:
ptr[n] = sqrt(rsq);
break;
case DX:
ptr[n] = dx*directionCorrection;
break;
case DY:
ptr[n] = dy*directionCorrection;
break;
case DZ:
ptr[n] = dz*directionCorrection;
break;
case ENGPOT:
ptr[n] = engpot;
break;
case FORCE:
ptr[n] = sqrt(rsq)*fbond;
break;
case FX:
ptr[n] = dx*fbond;
break;
case FY:
ptr[n] = dy*fbond;
break;
case FZ:
ptr[n] = dz*fbond;
break;
case ENGVIB:
ptr[n] = engvib;
break;
case ENGROT:
ptr[n] = engrot;
break;
case ENGTRANS:
ptr[n] = engtrans;
break;
case OMEGA:
ptr[n] = sqrt(omegasq);
break;
case VELVIB:
ptr[n] = vvib;
break;
case VARIABLE:
ptr[n] = input->variable->compute_equal(vvar[ivar]);
ivar++;
break;
case BN:
ptr[n] = bond->svector[bindex[n]];
break;
case DIST:
ptr[n] = sqrt(rsq);
break;
case DX:
ptr[n] = dx * directionCorrection;
break;
case DY:
ptr[n] = dy * directionCorrection;
break;
case DZ:
ptr[n] = dz * directionCorrection;
break;
case ENGPOT:
ptr[n] = engpot;
break;
case FORCE:
ptr[n] = sqrt(rsq) * fbond;
break;
case FX:
ptr[n] = dx * fbond;
break;
case FY:
ptr[n] = dy * fbond;
break;
case FZ:
ptr[n] = dz * fbond;
break;
case ENGVIB:
ptr[n] = engvib;
break;
case ENGROT:
ptr[n] = engrot;
break;
case ENGTRANS:
ptr[n] = engtrans;
break;
case OMEGA:
ptr[n] = sqrt(omegasq);
break;
case VELVIB:
ptr[n] = vvib;
break;
case VARIABLE:
ptr[n] = input->variable->compute_equal(vvar[ivar]);
ivar++;
break;
case BN:
ptr[n] = bond->svector[bindex[n]];
break;
}
}
}
@ -459,10 +495,10 @@ int ComputeBondLocal::compute_bonds(int flag)
/* ---------------------------------------------------------------------- */
int ComputeBondLocal::pack_forward_comm(int n, int *list, double *buf,
int /*pbc_flag*/, int * /*pbc*/)
int ComputeBondLocal::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/,
int * /*pbc*/)
{
int i,j,m;
int i, j, m;
double **v = atom->v;
@ -481,7 +517,7 @@ int ComputeBondLocal::pack_forward_comm(int n, int *list, double *buf,
void ComputeBondLocal::unpack_forward_comm(int n, int first, double *buf)
{
int i,m,last;
int i, m, last;
double **v = atom->v;
@ -504,11 +540,11 @@ void ComputeBondLocal::reallocate(int n)
if (nvalues == 1) {
memory->destroy(vlocal);
memory->create(vlocal,nmax,"bond/local:vector_local");
memory->create(vlocal, nmax, "bond/local:vector_local");
vector_local = vlocal;
} else {
memory->destroy(alocal);
memory->create(alocal,nmax,nvalues,"bond/local:array_local");
memory->create(alocal, nmax, nvalues, "bond/local:array_local");
array_local = alocal;
}
}
@ -519,6 +555,6 @@ void ComputeBondLocal::reallocate(int n)
double ComputeBondLocal::memory_usage()
{
double bytes = (double)nmax*nvalues * sizeof(double);
double bytes = (double) nmax * nvalues * sizeof(double);
return bytes;
}

176
src/compute_dihedral_local.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -13,19 +12,21 @@
------------------------------------------------------------------------- */
#include "compute_dihedral_local.h"
#include <cmath>
#include <cstring>
#include "atom.h"
#include "atom_vec.h"
#include "molecule.h"
#include "update.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "input.h"
#include "variable.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
#include "molecule.h"
#include "update.h"
#include "variable.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace MathConst;
@ -37,14 +38,13 @@ enum { PHI, VARIABLE };
/* ---------------------------------------------------------------------- */
ComputeDihedralLocal::ComputeDihedralLocal(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg),
bstyle(nullptr), vvar(nullptr), pstr(nullptr), vstr(nullptr), vlocal(nullptr), alocal(nullptr)
Compute(lmp, narg, arg), bstyle(nullptr), vvar(nullptr), pstr(nullptr), vstr(nullptr),
vlocal(nullptr), alocal(nullptr)
{
if (narg < 4) error->all(FLERR,"Illegal compute dihedral/local command");
if (narg < 4) error->all(FLERR, "Illegal compute dihedral/local command");
if (atom->avec->dihedrals_allow == 0)
error->all(FLERR,
"Compute dihedral/local used when dihedrals are not allowed");
error->all(FLERR, "Compute dihedral/local used when dihedrals are not allowed");
local_flag = 1;
@ -52,7 +52,7 @@ ComputeDihedralLocal::ComputeDihedralLocal(LAMMPS *lmp, int narg, char **arg) :
nvalues = narg - 3;
bstyle = new int[nvalues];
vstr = new char*[nvalues];
vstr = new char *[nvalues];
vvar = new int[nvalues];
nvalues = 0;
@ -60,13 +60,14 @@ ComputeDihedralLocal::ComputeDihedralLocal(LAMMPS *lmp, int narg, char **arg) :
int iarg;
for (iarg = 3; iarg < narg; iarg++) {
if (strcmp(arg[iarg],"phi") == 0) {
if (strcmp(arg[iarg], "phi") == 0) {
bstyle[nvalues++] = PHI;
} else if (strncmp(arg[iarg],"v_",2) == 0) {
} else if (strncmp(arg[iarg], "v_", 2) == 0) {
bstyle[nvalues++] = VARIABLE;
vstr[nvar] = utils::strdup(&arg[iarg][2]);
nvar++;
} else break;
} else
break;
}
// optional args
@ -75,47 +76,45 @@ ComputeDihedralLocal::ComputeDihedralLocal(LAMMPS *lmp, int narg, char **arg) :
pstr = nullptr;
while (iarg < narg) {
if (strcmp(arg[iarg],"set") == 0) {
if (strcmp(arg[iarg], "set") == 0) {
setflag = 1;
if (iarg+3 > narg)
error->all(FLERR,"Illegal compute dihedral/local command");
if (strcmp(arg[iarg+1],"phi") == 0) {
delete [] pstr;
pstr = utils::strdup(arg[iarg+2]);
} else error->all(FLERR,"Illegal compute dihedral/local command");
if (iarg + 3 > narg) error->all(FLERR, "Illegal compute dihedral/local command");
if (strcmp(arg[iarg + 1], "phi") == 0) {
delete[] pstr;
pstr = utils::strdup(arg[iarg + 2]);
} else
error->all(FLERR, "Illegal compute dihedral/local command");
iarg += 3;
} else error->all(FLERR,"Illegal compute dihedral/local command");
} else
error->all(FLERR, "Illegal compute dihedral/local command");
}
// error check
if (nvar) {
if (!setflag)
error->all(FLERR,"Compute dihedral/local variable requires a set variable");
if (!setflag) error->all(FLERR, "Compute dihedral/local variable requires a set variable");
for (int i = 0; i < nvar; i++) {
vvar[i] = input->variable->find(vstr[i]);
if (vvar[i] < 0)
error->all(FLERR,
"Variable name for copute dihedral/local does not exist");
if (vvar[i] < 0) error->all(FLERR, "Variable name for compute dihedral/local does not exist");
if (!input->variable->equalstyle(vvar[i]))
error->all(FLERR,"Variable for compute dihedral/local is invalid style");
error->all(FLERR, "Variable for compute dihedral/local is invalid style");
}
if (pstr) {
pvar = input->variable->find(pstr);
if (pvar < 0)
error->all(FLERR,
"Variable name for compute dihedral/local does not exist");
if (pvar < 0) error->all(FLERR, "Variable name for compute dihedral/local does not exist");
if (!input->variable->internalstyle(pvar))
error->all(FLERR,"Variable for compute dihedral/local is invalid style");
error->all(FLERR, "Variable for compute dihedral/local is invalid style");
}
} else if (setflag)
error->all(FLERR,"Compute dihedral/local set with no variable");
error->all(FLERR, "Compute dihedral/local set with no variable");
// initialize output
if (nvalues == 1) size_local_cols = 0;
else size_local_cols = nvalues;
if (nvalues == 1)
size_local_cols = 0;
else
size_local_cols = nvalues;
nmax = 0;
vlocal = nullptr;
@ -126,12 +125,12 @@ ComputeDihedralLocal::ComputeDihedralLocal(LAMMPS *lmp, int narg, char **arg) :
ComputeDihedralLocal::~ComputeDihedralLocal()
{
delete [] bstyle;
for (int i = 0; i < nvar; i++) delete [] vstr[i];
delete [] vstr;
delete [] vvar;
delete[] bstyle;
for (int i = 0; i < nvar; i++) delete[] vstr[i];
delete[] vstr;
delete[] vvar;
delete [] pstr;
delete[] pstr;
memory->destroy(vlocal);
memory->destroy(alocal);
@ -142,21 +141,17 @@ ComputeDihedralLocal::~ComputeDihedralLocal()
void ComputeDihedralLocal::init()
{
if (force->dihedral == nullptr)
error->all(FLERR,"No dihedral style is defined for compute dihedral/local");
error->all(FLERR, "No dihedral style is defined for compute dihedral/local");
if (nvar) {
for (int i = 0; i < nvar; i++) {
vvar[i] = input->variable->find(vstr[i]);
if (vvar[i] < 0)
error->all(FLERR,
"Variable name for compute dihedral/local does not exist");
if (vvar[i] < 0) error->all(FLERR, "Variable name for compute dihedral/local does not exist");
}
if (pstr) {
pvar = input->variable->find(pstr);
if (pvar < 0)
error->all(FLERR,
"Variable name for compute dihedral/local does not exist");
if (pvar < 0) error->all(FLERR, "Variable name for compute dihedral/local does not exist");
}
}
@ -191,11 +186,11 @@ void ComputeDihedralLocal::compute_local()
int ComputeDihedralLocal::compute_dihedrals(int flag)
{
int i,m,nd,atom1,atom2,atom3,atom4,imol,iatom,ivar;
int i, m, nd, atom1, atom2, atom3, atom4, imol, iatom, ivar;
tagint tagprev;
double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm;
double ax,ay,az,bx,by,bz,rasq,rbsq,rgsq,rg,ra2inv,rb2inv,rabinv;
double s,c,phi;
double vb1x, vb1y, vb1z, vb2x, vb2y, vb2z, vb3x, vb3y, vb3z, vb2xm, vb2ym, vb2zm;
double ax, ay, az, bx, by, bz, rasq, rbsq, rgsq, rg, ra2inv, rb2inv, rabinv;
double s, c, phi;
double *ptr;
double **x = atom->x;
@ -220,7 +215,8 @@ int ComputeDihedralLocal::compute_dihedrals(int flag)
for (atom2 = 0; atom2 < nlocal; atom2++) {
if (!(mask[atom2] & groupbit)) continue;
if (molecular == Atom::MOLECULAR) nd = num_dihedral[atom2];
if (molecular == Atom::MOLECULAR)
nd = num_dihedral[atom2];
else {
if (molindex[atom2] < 0) continue;
imol = molindex[atom2];
@ -237,9 +233,9 @@ int ComputeDihedralLocal::compute_dihedrals(int flag)
} else {
if (tag[atom2] != onemols[imol]->dihedral_atom2[atom2][i]) continue;
tagprev = tag[atom2] - iatom - 1;
atom1 = atom->map(onemols[imol]->dihedral_atom1[atom2][i]+tagprev);
atom3 = atom->map(onemols[imol]->dihedral_atom3[atom2][i]+tagprev);
atom4 = atom->map(onemols[imol]->dihedral_atom4[atom2][i]+tagprev);
atom1 = atom->map(onemols[imol]->dihedral_atom1[atom2][i] + tagprev);
atom3 = atom->map(onemols[imol]->dihedral_atom3[atom2][i] + tagprev);
atom4 = atom->map(onemols[imol]->dihedral_atom4[atom2][i] + tagprev);
}
if (atom1 < 0 || !(mask[atom1] & groupbit)) continue;
@ -256,64 +252,66 @@ int ComputeDihedralLocal::compute_dihedrals(int flag)
vb1x = x[atom1][0] - x[atom2][0];
vb1y = x[atom1][1] - x[atom2][1];
vb1z = x[atom1][2] - x[atom2][2];
domain->minimum_image(vb1x,vb1y,vb1z);
domain->minimum_image(vb1x, vb1y, vb1z);
vb2x = x[atom3][0] - x[atom2][0];
vb2y = x[atom3][1] - x[atom2][1];
vb2z = x[atom3][2] - x[atom2][2];
domain->minimum_image(vb2x,vb2y,vb2z);
domain->minimum_image(vb2x, vb2y, vb2z);
vb2xm = -vb2x;
vb2ym = -vb2y;
vb2zm = -vb2z;
domain->minimum_image(vb2xm,vb2ym,vb2zm);
domain->minimum_image(vb2xm, vb2ym, vb2zm);
vb3x = x[atom4][0] - x[atom3][0];
vb3y = x[atom4][1] - x[atom3][1];
vb3z = x[atom4][2] - x[atom3][2];
domain->minimum_image(vb3x,vb3y,vb3z);
domain->minimum_image(vb3x, vb3y, vb3z);
ax = vb1y*vb2zm - vb1z*vb2ym;
ay = vb1z*vb2xm - vb1x*vb2zm;
az = vb1x*vb2ym - vb1y*vb2xm;
bx = vb3y*vb2zm - vb3z*vb2ym;
by = vb3z*vb2xm - vb3x*vb2zm;
bz = vb3x*vb2ym - vb3y*vb2xm;
ax = vb1y * vb2zm - vb1z * vb2ym;
ay = vb1z * vb2xm - vb1x * vb2zm;
az = vb1x * vb2ym - vb1y * vb2xm;
bx = vb3y * vb2zm - vb3z * vb2ym;
by = vb3z * vb2xm - vb3x * vb2zm;
bz = vb3x * vb2ym - vb3y * vb2xm;
rasq = ax*ax + ay*ay + az*az;
rbsq = bx*bx + by*by + bz*bz;
rgsq = vb2xm*vb2xm + vb2ym*vb2ym + vb2zm*vb2zm;
rasq = ax * ax + ay * ay + az * az;
rbsq = bx * bx + by * by + bz * bz;
rgsq = vb2xm * vb2xm + vb2ym * vb2ym + vb2zm * vb2zm;
rg = sqrt(rgsq);
ra2inv = rb2inv = 0.0;
if (rasq > 0) ra2inv = 1.0/rasq;
if (rbsq > 0) rb2inv = 1.0/rbsq;
rabinv = sqrt(ra2inv*rb2inv);
if (rasq > 0) ra2inv = 1.0 / rasq;
if (rbsq > 0) rb2inv = 1.0 / rbsq;
rabinv = sqrt(ra2inv * rb2inv);
c = (ax*bx + ay*by + az*bz)*rabinv;
s = rg*rabinv*(ax*vb3x + ay*vb3y + az*vb3z);
c = (ax * bx + ay * by + az * bz) * rabinv;
s = rg * rabinv * (ax * vb3x + ay * vb3y + az * vb3z);
if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0;
phi = atan2(s,c);
phi = atan2(s, c);
if (nvalues == 1) ptr = &vlocal[m];
else ptr = alocal[m];
if (nvalues == 1)
ptr = &vlocal[m];
else
ptr = alocal[m];
if (nvar) {
ivar = 0;
if (pstr) input->variable->internal_set(pvar,phi);
if (pstr) input->variable->internal_set(pvar, phi);
}
for (int n = 0; n < nvalues; n++) {
switch (bstyle[n]) {
case PHI:
ptr[n] = 180.0*phi/MY_PI;
break;
case VARIABLE:
ptr[n] = input->variable->compute_equal(vvar[ivar]);
ivar++;
break;
case PHI:
ptr[n] = 180.0 * phi / MY_PI;
break;
case VARIABLE:
ptr[n] = input->variable->compute_equal(vvar[ivar]);
ivar++;
break;
}
}
@ -334,11 +332,11 @@ void ComputeDihedralLocal::reallocate(int n)
if (nvalues == 1) {
memory->destroy(vlocal);
memory->create(vlocal,nmax,"dihedral/local:vector_local");
memory->create(vlocal, nmax, "dihedral/local:vector_local");
vector_local = vlocal;
} else {
memory->destroy(alocal);
memory->create(alocal,nmax,nvalues,"dihedral/local:array_local");
memory->create(alocal, nmax, nvalues, "dihedral/local:array_local");
array_local = alocal;
}
}
@ -349,6 +347,6 @@ void ComputeDihedralLocal::reallocate(int n)
double ComputeDihedralLocal::memory_usage()
{
double bytes = (double)nmax*nvalues * sizeof(double);
double bytes = (double) nmax * nvalues * sizeof(double);
return bytes;
}

89
src/compute_improper_local.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -13,36 +12,35 @@
------------------------------------------------------------------------- */
#include "compute_improper_local.h"
#include <cmath>
#include <cstring>
#include "atom.h"
#include "atom_vec.h"
#include "molecule.h"
#include "update.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
#include "molecule.h"
#include "update.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace MathConst;
static constexpr int DELTA = 10000;
static constexpr double SMALL = 0.001;
static constexpr double SMALL = 0.001;
/* ---------------------------------------------------------------------- */
ComputeImproperLocal::ComputeImproperLocal(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg),
vlocal(nullptr), alocal(nullptr)
Compute(lmp, narg, arg), vlocal(nullptr), alocal(nullptr)
{
if (narg < 4) error->all(FLERR,"Illegal compute improper/local command");
if (narg < 4) error->all(FLERR, "Illegal compute improper/local command");
if (atom->avec->impropers_allow == 0)
error->all(FLERR,
"Compute improper/local used when impropers are not allowed");
error->all(FLERR, "Compute improper/local used when impropers are not allowed");
local_flag = 1;
nvalues = narg - 3;
@ -50,12 +48,16 @@ ComputeImproperLocal::ComputeImproperLocal(LAMMPS *lmp, int narg, char **arg) :
nvalues = 0;
for (int iarg = 3; iarg < narg; iarg++) {
if (strcmp(arg[iarg],"chi") == 0) cflag = nvalues++;
else error->all(FLERR,"Invalid keyword in compute improper/local command");
if (strcmp(arg[iarg], "chi") == 0)
cflag = nvalues++;
else
error->all(FLERR, "Invalid keyword in compute improper/local command");
}
if (nvalues == 1) size_local_cols = 0;
else size_local_cols = nvalues;
if (nvalues == 1)
size_local_cols = 0;
else
size_local_cols = nvalues;
nmax = 0;
vlocal = nullptr;
@ -75,7 +77,7 @@ ComputeImproperLocal::~ComputeImproperLocal()
void ComputeImproperLocal::init()
{
if (force->improper == nullptr)
error->all(FLERR,"No improper style is defined for compute improper/local");
error->all(FLERR, "No improper style is defined for compute improper/local");
// do initial memory allocation so that memory_usage() is correct
@ -108,11 +110,11 @@ void ComputeImproperLocal::compute_local()
int ComputeImproperLocal::compute_impropers(int flag)
{
int i,m,n,ni,atom1,atom2,atom3,atom4,imol,iatom;
int i, m, n, ni, atom1, atom2, atom3, atom4, imol, iatom;
tagint tagprev;
double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z;
double ss1,ss2,ss3,r1,r2,r3,c0,c1,c2,s1,s2;
double s12,c;
double vb1x, vb1y, vb1z, vb2x, vb2y, vb2z, vb3x, vb3y, vb3z;
double ss1, ss2, ss3, r1, r2, r3, c0, c1, c2, s1, s2;
double s12, c;
double *cbuf;
double **x = atom->x;
@ -135,8 +137,10 @@ int ComputeImproperLocal::compute_impropers(int flag)
if (nvalues == 1) {
if (cflag >= 0) cbuf = vlocal;
} else {
if (cflag >= 0 && alocal) cbuf = &alocal[0][cflag];
else cbuf = nullptr;
if (cflag >= 0 && alocal)
cbuf = &alocal[0][cflag];
else
cbuf = nullptr;
}
}
@ -144,7 +148,8 @@ int ComputeImproperLocal::compute_impropers(int flag)
for (atom2 = 0; atom2 < nlocal; atom2++) {
if (!(mask[atom2] & groupbit)) continue;
if (molecular == Atom::MOLECULAR) ni = num_improper[atom2];
if (molecular == Atom::MOLECULAR)
ni = num_improper[atom2];
else {
if (molindex[atom2] < 0) continue;
imol = molindex[atom2];
@ -161,9 +166,9 @@ int ComputeImproperLocal::compute_impropers(int flag)
} else {
if (tag[atom2] != onemols[imol]->improper_atom2[atom2][i]) continue;
tagprev = tag[atom2] - iatom - 1;
atom1 = atom->map(onemols[imol]->improper_atom1[atom2][i]+tagprev);
atom3 = atom->map(onemols[imol]->improper_atom3[atom2][i]+tagprev);
atom4 = atom->map(onemols[imol]->improper_atom4[atom2][i]+tagprev);
atom1 = atom->map(onemols[imol]->improper_atom1[atom2][i] + tagprev);
atom3 = atom->map(onemols[imol]->improper_atom3[atom2][i] + tagprev);
atom4 = atom->map(onemols[imol]->improper_atom4[atom2][i] + tagprev);
}
if (atom1 < 0 || !(mask[atom1] & groupbit)) continue;
@ -178,21 +183,21 @@ int ComputeImproperLocal::compute_impropers(int flag)
vb1x = x[atom1][0] - x[atom2][0];
vb1y = x[atom1][1] - x[atom2][1];
vb1z = x[atom1][2] - x[atom2][2];
domain->minimum_image(vb1x,vb1y,vb1z);
domain->minimum_image(vb1x, vb1y, vb1z);
vb2x = x[atom3][0] - x[atom2][0];
vb2y = x[atom3][1] - x[atom2][1];
vb2z = x[atom3][2] - x[atom2][2];
domain->minimum_image(vb2x,vb2y,vb2z);
domain->minimum_image(vb2x, vb2y, vb2z);
vb3x = x[atom4][0] - x[atom3][0];
vb3y = x[atom4][1] - x[atom3][1];
vb3z = x[atom4][2] - x[atom3][2];
domain->minimum_image(vb3x,vb3y,vb3z);
domain->minimum_image(vb3x, vb3y, vb3z);
ss1 = 1.0 / (vb1x*vb1x + vb1y*vb1y + vb1z*vb1z);
ss2 = 1.0 / (vb2x*vb2x + vb2y*vb2y + vb2z*vb2z);
ss3 = 1.0 / (vb3x*vb3x + vb3y*vb3y + vb3z*vb3z);
ss1 = 1.0 / (vb1x * vb1x + vb1y * vb1y + vb1z * vb1z);
ss2 = 1.0 / (vb2x * vb2x + vb2y * vb2y + vb2z * vb2z);
ss3 = 1.0 / (vb3x * vb3x + vb3y * vb3y + vb3z * vb3z);
r1 = sqrt(ss1);
r2 = sqrt(ss2);
@ -202,20 +207,20 @@ int ComputeImproperLocal::compute_impropers(int flag)
c1 = (vb1x * vb2x + vb1y * vb2y + vb1z * vb2z) * r1 * r2;
c2 = -(vb3x * vb2x + vb3y * vb2y + vb3z * vb2z) * r3 * r2;
s1 = 1.0 - c1*c1;
s1 = 1.0 - c1 * c1;
if (s1 < SMALL) s1 = SMALL;
s1 = 1.0 / s1;
s2 = 1.0 - c2*c2;
s2 = 1.0 - c2 * c2;
if (s2 < SMALL) s2 = SMALL;
s2 = 1.0 / s2;
s12 = sqrt(s1*s2);
c = (c1*c2 + c0) * s12;
s12 = sqrt(s1 * s2);
c = (c1 * c2 + c0) * s12;
if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0;
cbuf[n] = 180.0*acos(c)/MY_PI;
cbuf[n] = 180.0 * acos(c) / MY_PI;
}
n += nvalues;
}
@ -237,11 +242,11 @@ void ComputeImproperLocal::reallocate(int n)
if (nvalues == 1) {
memory->destroy(vlocal);
memory->create(vlocal,nmax,"improper/local:vector_local");
memory->create(vlocal, nmax, "improper/local:vector_local");
vector_local = vlocal;
} else {
memory->destroy(alocal);
memory->create(alocal,nmax,nvalues,"improper/local:array_local");
memory->create(alocal, nmax, nvalues, "improper/local:array_local");
array_local = alocal;
}
}
@ -252,6 +257,6 @@ void ComputeImproperLocal::reallocate(int n)
double ComputeImproperLocal::memory_usage()
{
double bytes = (double)nmax*nvalues * sizeof(double);
double bytes = (double) nmax * nvalues * sizeof(double);
return bytes;
}

4
src/compute_pair_local.cpp
View File

@ -38,7 +38,7 @@ enum { TYPE, RADIUS };
ComputePairLocal::ComputePairLocal(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), pstyle(nullptr), pindex(nullptr), vlocal(nullptr), alocal(nullptr)
{
if (narg < 4) error->all(FLERR, "Illegal compute pair/local command");
if (narg < 4) utils::missing_cmd_args(FLERR, "compute pair/local", error);
local_flag = 1;
nvalues = narg - 3;
@ -97,7 +97,7 @@ ComputePairLocal::ComputePairLocal(LAMMPS *lmp, int narg, char **arg) :
// error check
if (cutstyle == RADIUS && !atom->radius_flag)
error->all(FLERR, "This compute pair/local requires atom attribute radius");
error->all(FLERR, "Compute pair/local with ID {} requires atom attribute radius", id);
// set singleflag if need to call pair->single()

16
src/compute_property_local.cpp
View File

@ -39,7 +39,7 @@ ComputePropertyLocal::ComputePropertyLocal(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), vlocal(nullptr), alocal(nullptr), indices(nullptr),
pack_choice(nullptr)
{
if (narg < 4) error->all(FLERR, "Illegal compute property/local command");
if (narg < 4) utils::missing_cmd_args(FLERR, "compute property/local", error);
local_flag = 1;
nvalues = narg - 3;
@ -202,25 +202,23 @@ ComputePropertyLocal::ComputePropertyLocal(LAMMPS *lmp, int narg, char **arg) :
while (iarg < narg) {
if (strcmp(arg[iarg], "cutoff") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute property/local command");
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "compute property/local cutoff", error);
if (strcmp(arg[iarg + 1], "type") == 0)
cutstyle = TYPE;
else if (strcmp(arg[iarg + 1], "radius") == 0)
cutstyle = RADIUS;
else
error->all(FLERR, "Illegal compute property/local command");
error->all(FLERR, "Unknown compute property/local cutoff keyword: {}", arg[iarg + 1]);
iarg += 2;
} else
error->all(FLERR, "Illegal compute property/local command");
error->all(FLERR, "Unknown compute property/local keyword: {}", arg[iarg]);
}
// error check
if (atom->molecular == 2 &&
(kindflag == BOND || kindflag == ANGLE || kindflag == DIHEDRAL || kindflag == IMPROPER))
error->all(FLERR,
"Compute property/local does not (yet) work "
"with atom_style template");
error->all(FLERR, "Compute property/local does not (yet) work with atom_style template");
if (kindflag == BOND && atom->avec->bonds_allow == 0)
error->all(FLERR, "Compute property/local for property that isn't allocated");
@ -231,7 +229,7 @@ ComputePropertyLocal::ComputePropertyLocal(LAMMPS *lmp, int narg, char **arg) :
if (kindflag == IMPROPER && atom->avec->impropers_allow == 0)
error->all(FLERR, "Compute property/local for property that isn't allocated");
if (cutstyle == RADIUS && !atom->radius_flag)
error->all(FLERR, "Compute property/local requires atom attribute radius");
error->all(FLERR, "Compute property/local with ID {} requires atom attribute radius", id);
nmax = 0;
vlocal = nullptr;
@ -255,8 +253,6 @@ void ComputePropertyLocal::init()
if (kindflag == NEIGH || kindflag == PAIR) {
if (force->pair == nullptr)
error->all(FLERR, "No pair style is defined for compute property/local");
if (force->pair->single_enable == 0)
error->all(FLERR, "Pair style does not support compute property/local");
}
// for NEIGH/PAIR need an occasional half neighbor list

7
src/fix_nve.cpp
View File

@ -15,6 +15,7 @@
#include "fix_nve.h"
#include "atom.h"
#include "error.h"
#include "force.h"
#include "respa.h"
#include "update.h"
@ -28,7 +29,11 @@ FixNVE::FixNVE(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (!utils::strmatch(style,"^nve/sphere") && narg < 3)
utils::missing_cmd_args(FLERR, "fix nve", error);
utils::missing_cmd_args(FLERR, fmt::format("fix {}", style), error);
auto plain_style = utils::strip_style_suffix(style, lmp);
if (utils::strmatch(plain_style, "^nve$") && narg > 3)
error->all(FLERR, "Unsupported additional arguments for fix {}", style);
dynamic_group_allow = 1;
time_integrate = 1;

1
src/pair_hybrid.h
View File

@ -28,6 +28,7 @@ namespace LAMMPS_NS {
class PairHybrid : public Pair {
friend class AtomVecDielectric;
friend class ComputeSpin;
friend class FixAtomSwap;
friend class FixGPU;
friend class FixIntel;
friend class FixNVESpin;

21
tools/lammps-gui/chartviewer.cpp
View File

@ -423,27 +423,6 @@ void ChartViewer::add_data(int step, double data)
{
if (last_step < step) {
last_step = step;
// do not add data that deviates by more than 4 sigma from the average
// over the last 5 to 20 data items. this is a hack to work around
// getting corrupted data from lammps_get_last_thermo()
const auto &points = series->points();
const auto count = points.count();
if (count > 4) {
double ysum = 0.0;
double ysumsq = 0.0;
int first = count - 20;
if (first < 0) first = 0;
for (int i = first; i < count; ++i) {
double val = points[i].y();
ysum += val;
ysumsq += val * val;
}
const double num = count - first;
const double avg = ysum / num;
const double avgsq = ysumsq / num;
if (fabs(data - avg) > (5.0 * sqrt(avgsq - (avg * avg)))) return;
}
series->append(step, data);
QSettings settings;

5
tools/lammps-gui/lammpsgui.cpp
View File

@ -988,6 +988,7 @@ void LammpsGui::logupdate()
if (logwindow) {
const auto text = capturer->GetChunk();
if (text.size() > 0) {
logwindow->moveCursor(QTextCursor::End);
logwindow->insertPlainText(text.c_str());
logwindow->moveCursor(QTextCursor::End);
logwindow->textCursor().deleteChar();
@ -2016,8 +2017,8 @@ bool LammpsGui::eventFilter(QObject *watched, QEvent *event)
}
// LAMMPS geturl command with current location of the input and solution files on the web
static const QString geturl = "geturl https://raw.githubusercontent.com/akohlmey/"
"lammps-tutorials-inputs/main/tutorial%1/%2 output %2 verify no";
static const QString geturl = "geturl https://raw.githubusercontent.com/lammpstutorials/"
"lammpstutorials-article/refs/heads/main/files/tutorial%1/%2 output %2 verify no";
void LammpsGui::setup_tutorial(int tutno, const QString &dir, bool purgedir, bool getsolution)
{