diff --git a/doc/src/Commands_bond.rst b/doc/src/Commands_bond.rst index ef36b6b7c4..0c389df399 100644 --- a/doc/src/Commands_bond.rst +++ b/doc/src/Commands_bond.rst @@ -89,6 +89,7 @@ OPT. * :doc:`cosine/shift (o) ` * :doc:`cosine/shift/exp (o) ` * :doc:`cosine/squared (o) ` + * :doc:`cosine/squared/restricted (o) ` * :doc:`cross ` * :doc:`dipole (o) ` * :doc:`fourier (o) ` @@ -127,6 +128,7 @@ OPT. * :doc:`charmmfsw (k) ` * :doc:`class2 (ko) ` * :doc:`cosine/shift/exp (o) ` + * :doc:`cosine/squared/restricted ` * :doc:`fourier (io) ` * :doc:`harmonic (iko) ` * :doc:`helix (o) ` diff --git a/doc/src/angle_cosine_squared_restricted.rst b/doc/src/angle_cosine_squared_restricted.rst new file mode 100644 index 0000000000..c80a9142ff --- /dev/null +++ b/doc/src/angle_cosine_squared_restricted.rst @@ -0,0 +1,80 @@ +.. index:: angle_style cosine/squared/restricted +.. index:: angle_style cosine/squared/restricted/omp + +angle_style cosine/squared/restricted command +============================================= + +Accelerator Variants: *cosine/squared/restricted/omp* + +Syntax +"""""" + +.. code-block:: LAMMPS + + angle_style cosine/squared/restricted + +Examples +"""""""" + +.. code-block:: LAMMPS + + angle_style cosine/squared/restricted + angle_coeff 2*4 75.0 100.0 + +Description +""""""""""" + +.. versionadded:: TBD + +The *cosine/squared/restricted* angle style uses the potential + +.. math:: + + E = K [\cos(\theta) - \cos(\theta_0)]^2 / \sin^2(\theta) + +, which is commonly used in the MARTINI force field, +where :math:`\theta_0` is the equilibrium value of the angle, and :math:`K` +is a prefactor. Note that the usual 1/2 factor is included in :math:`K`. + +See :ref:`(Bulacu) ` for a description of the restricted angle for the MARTINI force field. + +The following coefficients must be defined for each angle type via the +:doc:`angle_coeff ` command as in the example above, or in +the data file or restart files read by the :doc:`read_data ` +or :doc:`read_restart ` commands: + +* :math:`K` (energy) +* :math:`\theta_0` (degrees) + +:math:`\theta_0` is specified in degrees, but LAMMPS converts it to radians +internally. + +---------- + +.. include:: accel_styles.rst + +---------- + +Restrictions +"""""""""""" + +This angle style can only be used if LAMMPS was built with the +EXTRA-MOLECULE package. See the :doc:`Build package ` doc page +for more info. + +Related commands +"""""""""""""""" + +:doc:`angle_coeff ` + +Default +""""""" + +none + +---------- + +.. _restricted-Bulacu: + +**(Bulacu)** Bulacu, Goga, Zhao, Rossi, Monticelli, Periole, Tieleman, Marrink, J Chem Theory Comput, 9, 3282-3292 +(2013). diff --git a/doc/src/angle_style.rst b/doc/src/angle_style.rst index 1f1ae72647..eeb174e772 100644 --- a/doc/src/angle_style.rst +++ b/doc/src/angle_style.rst @@ -10,7 +10,7 @@ Syntax angle_style style -* style = *none* or *zero* or *hybrid* or *amoeba* or *charmm* or *class2* or *class2/p6* or *cosine* or *cosine/buck6d* or *cosine/delta* or *cosine/periodic* or *cosine/shift* or *cosine/shift/exp* or *cosine/squared* or *cross* or *dipole* or *fourier* or *fourier/simple* or *gaussian* or *harmonic* or *lepton* or *mm3* or *quartic* or *spica* or *table* +* style = *none* or *zero* or *hybrid* or *amoeba* or *charmm* or *class2* or *class2/p6* or *cosine* or *cosine/buck6d* or *cosine/delta* or *cosine/periodic* or *cosine/shift* or *cosine/shift/exp* or *cosine/squared* or *cosine/squared/restricted* or *cross* or *dipole* or *fourier* or *fourier/simple* or *gaussian* or *harmonic* or *lepton* or *mm3* or *quartic* or *spica* or *table* Examples """""""" @@ -84,6 +84,7 @@ of (g,i,k,o,t) to indicate which accelerated styles exist. * :doc:`cosine/shift ` - angle cosine with a shift * :doc:`cosine/shift/exp ` - cosine with shift and exponential term in spring constant * :doc:`cosine/squared ` - angle with cosine squared term +* :doc:`cosine/squared/restricted ` - angle with restricted cosine squared term * :doc:`cross ` - cross term coupling angle and bond lengths * :doc:`dipole ` - angle that controls orientation of a point dipole * :doc:`fourier ` - angle with multiple cosine terms diff --git a/doc/src/dihedral_cosine_squared_restricted.rst b/doc/src/dihedral_cosine_squared_restricted.rst new file mode 100644 index 0000000000..cba846b965 --- /dev/null +++ b/doc/src/dihedral_cosine_squared_restricted.rst @@ -0,0 +1,71 @@ +.. index:: dihedral_style cosine/squared/restricted + +dihedral_style cosine/squared/restricted command +================================================ + + +Syntax +"""""" + +.. code-block:: LAMMPS + + dihedral_style cosine/squared/restricted + +Examples +"""""""" + +.. code-block:: LAMMPS + + dihedral_style cosine/squared/restricted + dihedral_coeff 1 10.0 120 + +Description +""""""""""" + +.. versionadded:: TBD + +The *cosine/squared/restricted* dihedral style uses the potential + +.. math:: + + E = K [\cos(\phi) - \cos(\phi_0)]^2 / \sin^2(\phi) + +, which is commonly used in the MARTINI force field. + +See :ref:`(Bulacu) ` for a description of the restricted dihedral for the MARTINI force field. + +The following coefficients must be defined for each dihedral type via the +:doc:`dihedral_coeff ` command as in the example above, or in +the data file or restart files read by the :doc:`read_data ` +or :doc:`read_restart ` commands: + +* :math:`K` (energy) +* :math:`\phi_0` (degrees) + +:math:`\phi_0` is specified in degrees, but LAMMPS converts it to radians internally. + +---------- + +Restrictions +"""""""""""" + +This dihedral style can only be used if LAMMPS was built with the +EXTRA-MOLECULE package. See the :doc:`Build package ` doc page +for more info. + +Related commands +"""""""""""""""" + +:doc:`dihedral_coeff ` + +Default +""""""" + +none + +---------- + +.. _restricted-Bul: + +**(Bulacu)** Bulacu, Goga, Zhao, Rossi, Monticelli, Periole, Tieleman, Marrink, J Chem Theory Comput, 9, 3282-3292 +(2013). diff --git a/doc/src/dihedral_style.rst b/doc/src/dihedral_style.rst index 45dd66e750..56b8ff5226 100644 --- a/doc/src/dihedral_style.rst +++ b/doc/src/dihedral_style.rst @@ -10,7 +10,7 @@ Syntax dihedral_style style -* style = *none* or *zero* or *hybrid* or *charmm* or *charmmfsw* or *class2* or *cosine/shift/exp* or *fourier* or *harmonic* or *helix* or *lepton* or *multi/harmonic* or *nharmonic* or *opls* or *spherical* or *table* or *table/cut* +* style = *none* or *zero* or *hybrid* or *charmm* or *charmmfsw* or *class2* or *cosine/shift/exp* or *cosine/squared/restricted* or *fourier* or *harmonic* or *helix* or *lepton* or *multi/harmonic* or *nharmonic* or *opls* or *spherical* or *table* or *table/cut* Examples """""""" @@ -105,6 +105,7 @@ exist. * :doc:`charmmfsw ` - CHARMM dihedral with force switching * :doc:`class2 ` - COMPASS (class 2) dihedral * :doc:`cosine/shift/exp ` - dihedral with exponential in spring constant +* :doc:`cosine/squared/restricted ` - squared cosine dihedral with restricted term * :doc:`fourier ` - dihedral with multiple cosine terms * :doc:`harmonic ` - harmonic dihedral * :doc:`helix ` - helix dihedral diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index cf4e48feff..b531b7bb51 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -380,6 +380,7 @@ btype buckPlusAttr buf builtin +Bulacu Bulatov Bureekaew burlywood @@ -1303,6 +1304,7 @@ gmres gname gneb GNEB +Goga Goldfarb Gompper Gonzalez-Melchor @@ -2262,6 +2264,7 @@ Montalenti Monterey Montero Monti +Monticelli Mora Morefoo Morfill @@ -2790,6 +2793,7 @@ peridynamic Peridynamic peridynamics Peridynamics +Periole perl permittivity perp @@ -3205,6 +3209,7 @@ Ronevich Rosati Rosato Rosenberger +Rossi Rossky rosybrown rotationally @@ -3648,6 +3653,7 @@ thrid ThunderX thylakoid THz +Tieleman Tigran Tij Tildesley diff --git a/python/lammps/pylammps.py b/python/lammps/pylammps.py index b959cf7d65..d8f175c5b8 100644 --- a/python/lammps/pylammps.py +++ b/python/lammps/pylammps.py @@ -544,6 +544,18 @@ class PyLammps(object): """ self._cmd_history = [] + + def append_cmd_history(self, cmd): + """ + Commands will be added to the command history but not executed. + + Add `commands` only to the command history, but do not execute them, so that you can + conveniently create Lammps input files, using + :py:meth:`PyLammps.write_script()`. + """ + self._cmd_history.append(cmd) + + def command(self, cmd): """ Execute LAMMPS command diff --git a/src/.gitignore b/src/.gitignore index 02f2a6a6ea..b145f81159 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -463,6 +463,8 @@ /angle_cosine_shift_exp.h /angle_cosine_squared.cpp /angle_cosine_squared.h +/angle_cosine_squared_restricted.cpp +/angle_cosine_squared_restricted.h /angle_cross.cpp /angle_cross.h /angle_dipole.cpp @@ -681,6 +683,8 @@ /dihedral_class2.h /dihedral_cosine_shift_exp.cpp /dihedral_cosine_shift_exp.h +/dihedral_cosine_squared_restricted.cpp +/dihedral_cosine_squared_restricted.h /dihedral_fourier.cpp /dihedral_fourier.h /dihedral_harmonic.cpp diff --git a/src/EXTRA-MOLECULE/angle_cosine_squared_restricted.cpp b/src/EXTRA-MOLECULE/angle_cosine_squared_restricted.cpp new file mode 100644 index 0000000000..2da31ef893 --- /dev/null +++ b/src/EXTRA-MOLECULE/angle_cosine_squared_restricted.cpp @@ -0,0 +1,298 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "angle_cosine_squared_restricted.h" + +#include "atom.h" +#include "comm.h" +#include "domain.h" +#include "error.h" +#include "force.h" +#include "math_const.h" +#include "memory.h" +#include "neighbor.h" + +#include + +using namespace LAMMPS_NS; +using MathConst::DEG2RAD; +using MathConst::RAD2DEG; + +/* ---------------------------------------------------------------------- */ + +AngleCosineSquaredRestricted::AngleCosineSquaredRestricted(LAMMPS *_lmp) : Angle(_lmp) +{ + k = nullptr; + theta0 = nullptr; + born_matrix_enable = 1; +} + +/* ---------------------------------------------------------------------- */ + +AngleCosineSquaredRestricted::~AngleCosineSquaredRestricted() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(k); + memory->destroy(theta0); + } +} + +/* ---------------------------------------------------------------------- */ + +void AngleCosineSquaredRestricted::compute(int eflag, int vflag) +{ + int i1, i2, i3, n, type; + double delx1, dely1, delz1, delx2, dely2, delz2; + double eangle, f1[3], f3[3]; + double tk, rsq1, rsq2, r1, r2, c, a, a11, a12, a22; + + eangle = 0.0; + ev_init(eflag, vflag); + + double **x = atom->x; + double **f = atom->f; + int **anglelist = neighbor->anglelist; + int nanglelist = neighbor->nanglelist; + int nlocal = atom->nlocal; + int newton_bond = force->newton_bond; + + for (n = 0; n < nanglelist; n++) { + i1 = anglelist[n][0]; + i2 = anglelist[n][1]; + i3 = anglelist[n][2]; + type = anglelist[n][3]; + + // 1st bond + + delx1 = x[i1][0] - x[i2][0]; + dely1 = x[i1][1] - x[i2][1]; + delz1 = x[i1][2] - x[i2][2]; + + rsq1 = delx1 * delx1 + dely1 * dely1 + delz1 * delz1; + r1 = sqrt(rsq1); + + // 2nd bond + + delx2 = x[i3][0] - x[i2][0]; + dely2 = x[i3][1] - x[i2][1]; + delz2 = x[i3][2] - x[i2][2]; + + rsq2 = delx2 * delx2 + dely2 * dely2 + delz2 * delz2; + r2 = sqrt(rsq2); + + // angle (cos and sin) + + c = delx1 * delx2 + dely1 * dely2 + delz1 * delz2; + c /= r1 * r2; + + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + + // force & energy + + double sq_sin = 1.0 - c * c; + double c0 = cos(theta0[type]); + + tk = k[type] * (c - c0) * (1.0 - c * c0) / (sq_sin * sq_sin); + + if (eflag) eangle = k[type] * (c - c0) * (c - c0) / sq_sin; + + a = 2.0 * tk; + a11 = a * c / rsq1; + a12 = -a / (r1 * r2); + a22 = a * c / rsq2; + + f1[0] = a11 * delx1 + a12 * delx2; + f1[1] = a11 * dely1 + a12 * dely2; + f1[2] = a11 * delz1 + a12 * delz2; + f3[0] = a22 * delx2 + a12 * delx1; + f3[1] = a22 * dely2 + a12 * dely1; + f3[2] = a22 * delz2 + a12 * delz1; + + // apply force to each of 3 atoms + + if (newton_bond || i1 < nlocal) { + f[i1][0] += f1[0]; + f[i1][1] += f1[1]; + f[i1][2] += f1[2]; + } + + if (newton_bond || i2 < nlocal) { + f[i2][0] -= f1[0] + f3[0]; + f[i2][1] -= f1[1] + f3[1]; + f[i2][2] -= f1[2] + f3[2]; + } + + if (newton_bond || i3 < nlocal) { + f[i3][0] += f3[0]; + f[i3][1] += f3[1]; + f[i3][2] += f3[2]; + } + + if (evflag) + ev_tally(i1, i2, i3, nlocal, newton_bond, eangle, f1, f3, delx1, dely1, delz1, delx2, dely2, + delz2); + } +} + +/* ---------------------------------------------------------------------- */ + +void AngleCosineSquaredRestricted::allocate() +{ + allocated = 1; + const int np1 = atom->nangletypes + 1; + + memory->create(k, np1, "angle:k"); + memory->create(theta0, np1, "angle:theta0"); + + memory->create(setflag, np1, "angle:setflag"); + for (int i = 1; i < np1; i++) setflag[i] = 0; +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more types +------------------------------------------------------------------------- */ + +void AngleCosineSquaredRestricted::coeff(int narg, char **arg) +{ + if (narg != 3) error->all(FLERR, "Incorrect args for angle coefficients"); + if (!allocated) allocate(); + + int ilo, ihi; + utils::bounds(FLERR, arg[0], 1, atom->nangletypes, ilo, ihi, error); + + double k_one = utils::numeric(FLERR, arg[1], false, lmp); + double theta0_one = utils::numeric(FLERR, arg[2], false, lmp); + + // convert theta0 from degrees to radians + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + k[i] = k_one; + theta0[i] = DEG2RAD * theta0_one; + setflag[i] = 1; + count++; + } + + if (count == 0) error->all(FLERR, "Incorrect args for angle coefficients"); +} + +/* ---------------------------------------------------------------------- */ + +double AngleCosineSquaredRestricted::equilibrium_angle(int i) +{ + return theta0[i]; +} + +/* ---------------------------------------------------------------------- + proc 0 writes out coeffs to restart file +------------------------------------------------------------------------- */ + +void AngleCosineSquaredRestricted::write_restart(FILE *fp) +{ + fwrite(&k[1], sizeof(double), atom->nangletypes, fp); + fwrite(&theta0[1], sizeof(double), atom->nangletypes, fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads coeffs from restart file, bcasts them +------------------------------------------------------------------------- */ + +void AngleCosineSquaredRestricted::read_restart(FILE *fp) +{ + allocate(); + + if (comm->me == 0) { + utils::sfread(FLERR, &k[1], sizeof(double), atom->nangletypes, fp, nullptr, error); + utils::sfread(FLERR, &theta0[1], sizeof(double), atom->nangletypes, fp, nullptr, error); + } + MPI_Bcast(&k[1], atom->nangletypes, MPI_DOUBLE, 0, world); + MPI_Bcast(&theta0[1], atom->nangletypes, MPI_DOUBLE, 0, world); + + for (int i = 1; i <= atom->nangletypes; i++) setflag[i] = 1; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to data file +------------------------------------------------------------------------- */ + +void AngleCosineSquaredRestricted::write_data(FILE *fp) +{ + for (int i = 1; i <= atom->nangletypes; i++) + fprintf(fp, "%d %g %g\n", i, k[i], RAD2DEG * theta0[i]); +} + +/* ---------------------------------------------------------------------- */ + +double AngleCosineSquaredRestricted::single(int type, int i1, int i2, int i3) +{ + double **x = atom->x; + + double delx1 = x[i1][0] - x[i2][0]; + double dely1 = x[i1][1] - x[i2][1]; + double delz1 = x[i1][2] - x[i2][2]; + domain->minimum_image(delx1, dely1, delz1); + double r1 = sqrt(delx1 * delx1 + dely1 * dely1 + delz1 * delz1); + + double delx2 = x[i3][0] - x[i2][0]; + double dely2 = x[i3][1] - x[i2][1]; + double delz2 = x[i3][2] - x[i2][2]; + domain->minimum_image(delx2, dely2, delz2); + double r2 = sqrt(delx2 * delx2 + dely2 * dely2 + delz2 * delz2); + + double c = delx1 * delx2 + dely1 * dely2 + delz1 * delz2; + c /= r1 * r2; + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + + double sq_sin = 1.0 - c * c; + double dcostheta = c - cos(theta0[type]); + double tk = k[type] * dcostheta / sq_sin; + return tk * dcostheta; +} + +/* ---------------------------------------------------------------------- */ + +void AngleCosineSquaredRestricted::born_matrix(int type, int i1, int i2, int i3, double &du, + double &du2) +{ + double **x = atom->x; + + double delx1 = x[i1][0] - x[i2][0]; + double dely1 = x[i1][1] - x[i2][1]; + double delz1 = x[i1][2] - x[i2][2]; + domain->minimum_image(delx1, dely1, delz1); + double r1 = sqrt(delx1 * delx1 + dely1 * dely1 + delz1 * delz1); + + double delx2 = x[i3][0] - x[i2][0]; + double dely2 = x[i3][1] - x[i2][1]; + double delz2 = x[i3][2] - x[i2][2]; + domain->minimum_image(delx2, dely2, delz2); + double r2 = sqrt(delx2 * delx2 + dely2 * dely2 + delz2 * delz2); + + double c = delx1 * delx2 + dely1 * dely2 + delz1 * delz2; + c /= r1 * r2; + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + + double c0 = cos(theta0[type]); + double sq_sin = 1.0 - c * c; + + du = 2 * k[type] * (c - c0) * (1.0 - c * c0) / (sq_sin * sq_sin); + + double numerator = 2.0 * c0 * c * c * c - 3.0 * (c0 * c0 + 1) * c * c + 6 * c0 * c - c0 * c0 - 1; + double denominator = sq_sin * sq_sin * sq_sin; + + du2 = 2 * k[type] * numerator / denominator; +} diff --git a/src/EXTRA-MOLECULE/angle_cosine_squared_restricted.h b/src/EXTRA-MOLECULE/angle_cosine_squared_restricted.h new file mode 100644 index 0000000000..674252b7d0 --- /dev/null +++ b/src/EXTRA-MOLECULE/angle_cosine_squared_restricted.h @@ -0,0 +1,49 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef ANGLE_CLASS +// clang-format off +AngleStyle(cosine/squared/restricted,AngleCosineSquaredRestricted); +// clang-format on +#else + +#ifndef LMP_ANGLE_COSINE_SQUARED_RESTRICTED_H +#define LMP_ANGLE_COSINE_SQUARED_RESTRICTED_H + +#include "angle.h" + +namespace LAMMPS_NS { + +class AngleCosineSquaredRestricted : public Angle { + public: + AngleCosineSquaredRestricted(class LAMMPS *); + ~AngleCosineSquaredRestricted() override; + void compute(int, int) override; + void coeff(int, char **) override; + double equilibrium_angle(int) override; + void write_restart(FILE *) override; + void read_restart(FILE *) override; + void write_data(FILE *) override; + double single(int, int, int, int) override; + void born_matrix(int type, int i1, int i2, int i3, double &du, double &du2) override; + + protected: + double *k, *theta0; + + void allocate(); +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/EXTRA-MOLECULE/dihedral_cosine_squared_restricted.cpp b/src/EXTRA-MOLECULE/dihedral_cosine_squared_restricted.cpp new file mode 100644 index 0000000000..5b2feb6897 --- /dev/null +++ b/src/EXTRA-MOLECULE/dihedral_cosine_squared_restricted.cpp @@ -0,0 +1,395 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "dihedral_cosine_squared_restricted.h" + +#include "atom.h" +#include "comm.h" +#include "error.h" +#include "force.h" +#include "math_const.h" +#include "memory.h" +#include "neighbor.h" + +#include + +using namespace LAMMPS_NS; +using MathConst::DEG2RAD; +using MathConst::RAD2DEG; + +static constexpr double TOLERANCE = 0.05; +static constexpr double SMALL = 0.001; + +/* ---------------------------------------------------------------------- */ + +DihedralCosineSquaredRestricted::DihedralCosineSquaredRestricted(LAMMPS *_lmp) : + Dihedral(_lmp), k(nullptr), phi0(nullptr) +{ + writedata = 1; + born_matrix_enable = 1; +} + +/* ---------------------------------------------------------------------- */ + +DihedralCosineSquaredRestricted::~DihedralCosineSquaredRestricted() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(k); + memory->destroy(phi0); + } +} + +/* ---------------------------------------------------------------------- */ + +void DihedralCosineSquaredRestricted::compute(int eflag, int vflag) +{ + int i1, i2, i3, i4, n, type; + double vb1x, vb1y, vb1z, vb2x, vb2y, vb2z, vb3x, vb3y, vb3z, vb2xm, vb2ym, vb2zm; + double edihedral, f1[3], f2[3], f3[3], f4[3]; + double sb1, sb2, sb3, rb1, rb3, c0, b1mag2, b1mag, b2mag2; + double b2mag, b3mag2, b3mag, ctmp, r12c1, c1mag, r12c2; + double c2mag, sc1, sc2, s1, s12, c, pd, a, a11, a22; + double a33, a12, a13, a23, sx2, sy2, sz2; + double s2, sin2; + + edihedral = 0.0; + ev_init(eflag, vflag); + + double **x = atom->x; + double **f = atom->f; + int **dihedrallist = neighbor->dihedrallist; + int ndihedrallist = neighbor->ndihedrallist; + int nlocal = atom->nlocal; + int newton_bond = force->newton_bond; + + for (n = 0; n < ndihedrallist; n++) { + i1 = dihedrallist[n][0]; + i2 = dihedrallist[n][1]; + i3 = dihedrallist[n][2]; + i4 = dihedrallist[n][3]; + type = dihedrallist[n][4]; + + // 1st bond + + vb1x = x[i1][0] - x[i2][0]; + vb1y = x[i1][1] - x[i2][1]; + vb1z = x[i1][2] - x[i2][2]; + + // 2nd bond + + vb2x = x[i3][0] - x[i2][0]; + vb2y = x[i3][1] - x[i2][1]; + vb2z = x[i3][2] - x[i2][2]; + + vb2xm = -vb2x; + vb2ym = -vb2y; + vb2zm = -vb2z; + + // 3rd bond + + vb3x = x[i4][0] - x[i3][0]; + vb3y = x[i4][1] - x[i3][1]; + vb3z = x[i4][2] - x[i3][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); + + rb1 = sqrt(sb1); + rb3 = sqrt(sb3); + + c0 = (vb1x * vb3x + vb1y * vb3y + vb1z * vb3z) * rb1 * rb3; + + // 1st and 2nd angle + + b1mag2 = vb1x * vb1x + vb1y * vb1y + vb1z * vb1z; + b1mag = sqrt(b1mag2); + b2mag2 = vb2x * vb2x + vb2y * vb2y + vb2z * vb2z; + b2mag = sqrt(b2mag2); + b3mag2 = vb3x * vb3x + vb3y * vb3y + vb3z * vb3z; + b3mag = sqrt(b3mag2); + + 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); + c2mag = ctmp * r12c2; + + // cos and sin of 2 angles and final c + + sin2 = MAX(1.0 - c1mag * c1mag, 0.0); + sc1 = sqrt(sin2); + if (sc1 < SMALL) sc1 = SMALL; + sc1 = 1.0 / sc1; + + sin2 = MAX(1.0 - c2mag * c2mag, 0.0); + sc2 = sqrt(sin2); + if (sc2 < SMALL) sc2 = SMALL; + sc2 = 1.0 / sc2; + + s1 = sc1 * sc1; + s2 = sc2 * sc2; + s12 = sc1 * sc2; + c = (c0 + c1mag * c2mag) * s12; + + // error check + + if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) problem(FLERR, i1, i2, i3, i4); + + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + + // force & energy + + double p0 = cos(phi0[type]); + double sq_sin = 1.0 - c * c; + + pd = 2 * k[type] * (c - p0) * (1.0 - c * p0) / (sq_sin * sq_sin); + + if (eflag) edihedral = k[type] * (c - p0) * (c - p0) / sq_sin; + + a = pd; + 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); + + 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; + + 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; + + f3[0] = sx2 - f4[0]; + f3[1] = sy2 - f4[1]; + f3[2] = sz2 - f4[2]; + + // apply force to each of 4 atoms + + if (newton_bond || i1 < nlocal) { + f[i1][0] += f1[0]; + f[i1][1] += f1[1]; + f[i1][2] += f1[2]; + } + + if (newton_bond || i2 < nlocal) { + f[i2][0] += f2[0]; + f[i2][1] += f2[1]; + f[i2][2] += f2[2]; + } + + if (newton_bond || i3 < nlocal) { + f[i3][0] += f3[0]; + f[i3][1] += f3[1]; + f[i3][2] += f3[2]; + } + + if (newton_bond || i4 < nlocal) { + f[i4][0] += f4[0]; + f[i4][1] += f4[1]; + f[i4][2] += f4[2]; + } + + if (evflag) + ev_tally(i1, i2, i3, i4, nlocal, newton_bond, edihedral, f1, f3, f4, vb1x, vb1y, vb1z, vb2x, + vb2y, vb2z, vb3x, vb3y, vb3z); + } +} + +/* ---------------------------------------------------------------------- */ + +void DihedralCosineSquaredRestricted::allocate() +{ + allocated = 1; + const int np1 = atom->ndihedraltypes + 1; + + memory->create(k, np1, "dihedral:k"); + memory->create(phi0, np1, "dihedral:phi0"); + + memory->create(setflag, np1, "dihedral:setflag"); + for (int i = 1; i < np1; i++) setflag[i] = 0; +} + +/* ---------------------------------------------------------------------- + set coeffs for one type +------------------------------------------------------------------------- */ + +void DihedralCosineSquaredRestricted::coeff(int narg, char **arg) +{ + if (narg != 3) error->all(FLERR, "Incorrect args for dihedral coefficients"); + if (!allocated) allocate(); + + int ilo, ihi; + utils::bounds(FLERR, arg[0], 1, atom->ndihedraltypes, ilo, ihi, error); + + double k_one = utils::numeric(FLERR, arg[1], false, lmp); + double phi0_one = utils::numeric(FLERR, arg[2], false, lmp); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + k[i] = k_one; + phi0[i] = DEG2RAD * phi0_one; + setflag[i] = 1; + count++; + } + + if (count == 0) error->all(FLERR, "Incorrect args for dihedral coefficients"); +} + +/* ---------------------------------------------------------------------- + proc 0 writes out coeffs to restart file +------------------------------------------------------------------------- */ + +void DihedralCosineSquaredRestricted::write_restart(FILE *fp) +{ + fwrite(&k[1], sizeof(double), atom->ndihedraltypes, fp); + fwrite(&phi0[1], sizeof(double), atom->ndihedraltypes, fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads coeffs from restart file, bcasts them +------------------------------------------------------------------------- */ + +void DihedralCosineSquaredRestricted::read_restart(FILE *fp) +{ + allocate(); + + if (comm->me == 0) { + utils::sfread(FLERR, &k[1], sizeof(double), atom->ndihedraltypes, fp, nullptr, error); + utils::sfread(FLERR, &phi0[1], sizeof(double), atom->ndihedraltypes, fp, nullptr, error); + } + MPI_Bcast(&k[1], atom->ndihedraltypes, MPI_DOUBLE, 0, world); + MPI_Bcast(&phi0[1], atom->ndihedraltypes, MPI_DOUBLE, 0, world); + + for (int i = 1; i <= atom->ndihedraltypes; i++) setflag[i] = 1; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to data file +------------------------------------------------------------------------- */ + +void DihedralCosineSquaredRestricted::write_data(FILE *fp) +{ + for (int i = 1; i <= atom->ndihedraltypes; i++) + fprintf(fp, "%d %g %g\n", i, k[i], RAD2DEG * phi0[i]); +} + +/* ---------------------------------------------------------------------- */ + +void DihedralCosineSquaredRestricted::born_matrix(int nd, int i1, int i2, int i3, int i4, + double &du, double &du2) +{ + double vb1x, vb1y, vb1z, vb2x, vb2y, vb2z, vb3x, vb3y, vb3z, vb2xm, vb2ym, vb2zm; + double sb1, sb3, rb1, rb3, c0, b1mag2, b1mag, b2mag2; + double b2mag, b3mag2, b3mag, ctmp, r12c1, c1mag, r12c2; + double c2mag, sc1, sc2, s12, c; + double sin2; + + double **x = atom->x; + int **dihedrallist = neighbor->dihedrallist; + + int type = dihedrallist[nd][4]; + + // 1st bond + vb1x = x[i1][0] - x[i2][0]; + vb1y = x[i1][1] - x[i2][1]; + vb1z = x[i1][2] - x[i2][2]; + + // 2nd bond + vb2x = x[i3][0] - x[i2][0]; + vb2y = x[i3][1] - x[i2][1]; + vb2z = x[i3][2] - x[i2][2]; + + vb2xm = -vb2x; + vb2ym = -vb2y; + vb2zm = -vb2z; + + // 3rd bond + vb3x = x[i4][0] - x[i3][0]; + vb3y = x[i4][1] - x[i3][1]; + vb3z = x[i4][2] - x[i3][2]; + + // c0 calculation + sb1 = 1.0 / (vb1x * vb1x + vb1y * vb1y + vb1z * vb1z); + sb3 = 1.0 / (vb3x * vb3x + vb3y * vb3y + vb3z * vb3z); + + rb1 = sqrt(sb1); + rb3 = sqrt(sb3); + + c0 = (vb1x * vb3x + vb1y * vb3y + vb1z * vb3z) * rb1 * rb3; + + // 1st and 2nd angle + b1mag2 = vb1x * vb1x + vb1y * vb1y + vb1z * vb1z; + b1mag = sqrt(b1mag2); + b2mag2 = vb2x * vb2x + vb2y * vb2y + vb2z * vb2z; + b2mag = sqrt(b2mag2); + b3mag2 = vb3x * vb3x + vb3y * vb3y + vb3z * vb3z; + b3mag = sqrt(b3mag2); + + 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); + c2mag = ctmp * r12c2; + + // cos and sin of 2 angles and final c + sin2 = MAX(1.0 - c1mag * c1mag, 0.0); + sc1 = sqrt(sin2); + if (sc1 < SMALL) sc1 = SMALL; + sc1 = 1.0 / sc1; + + sin2 = MAX(1.0 - c2mag * c2mag, 0.0); + sc2 = sqrt(sin2); + if (sc2 < SMALL) sc2 = SMALL; + sc2 = 1.0 / sc2; + + s12 = sc1 * sc2; + c = (c0 + c1mag * c2mag) * s12; + + // error check + if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) problem(FLERR, i1, i2, i3, i4); + + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + + double p0 = cos(phi0[type]); + double sq_sin = 1.0 - c * c; + + du = 2 * k[type] * (c - p0) * (1.0 - c * p0) / (sq_sin * sq_sin); + + double numerator = 2.0 * p0 * c * c * c - 3.0 * (p0 * p0 + 1) * c * c + 6 * p0 * c - p0 * p0 - 1; + double denominator = sq_sin * sq_sin * sq_sin; + + du2 = 2 * k[type] * numerator / denominator; +} diff --git a/src/EXTRA-MOLECULE/dihedral_cosine_squared_restricted.h b/src/EXTRA-MOLECULE/dihedral_cosine_squared_restricted.h new file mode 100644 index 0000000000..b9e2e1d9d8 --- /dev/null +++ b/src/EXTRA-MOLECULE/dihedral_cosine_squared_restricted.h @@ -0,0 +1,47 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef DIHEDRAL_CLASS +// clang-format off +DihedralStyle(cosine/squared/restricted,DihedralCosineSquaredRestricted); +// clang-format on +#else + +#ifndef LMP_DIHEDRAL_COSINE_SQUARED_RESTRICTED_H +#define LMP_DIHEDRAL_COSINE_SQUARED_RESTRICTED_H + +#include "dihedral.h" + +namespace LAMMPS_NS { + +class DihedralCosineSquaredRestricted : public Dihedral { + public: + DihedralCosineSquaredRestricted(class LAMMPS *); + ~DihedralCosineSquaredRestricted() override; + void compute(int, int) override; + void coeff(int, char **) override; + void write_restart(FILE *) override; + void read_restart(FILE *) override; + void write_data(FILE *) override; + void born_matrix(int, int, int, int, int, double &, double &) override; + + protected: + double *k, *phi0; + + void allocate(); +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/OPENMP/angle_cosine_squared_restricted_omp.cpp b/src/OPENMP/angle_cosine_squared_restricted_omp.cpp new file mode 100644 index 0000000000..80960653b4 --- /dev/null +++ b/src/OPENMP/angle_cosine_squared_restricted_omp.cpp @@ -0,0 +1,170 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Axel Kohlmeyer (Temple U) +------------------------------------------------------------------------- */ + +#include "angle_cosine_squared_restricted_omp.h" + +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "neighbor.h" + +#include + +#include "omp_compat.h" +#include "suffix.h" +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +AngleCosineSquaredRestrictedOMP::AngleCosineSquaredRestrictedOMP(class LAMMPS *lmp) + : AngleCosineSquaredRestricted(lmp), ThrOMP(lmp,THR_ANGLE) +{ + suffix_flag |= Suffix::OMP; +} + +/* ---------------------------------------------------------------------- */ + +void AngleCosineSquaredRestrictedOMP::compute(int eflag, int vflag) +{ + ev_init(eflag,vflag); + + const int nall = atom->nlocal + atom->nghost; + const int nthreads = comm->nthreads; + const int inum = neighbor->nanglelist; + +#if defined(_OPENMP) +#pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(eflag,vflag) +#endif + { + int ifrom, ito, tid; + + loop_setup_thr(ifrom, ito, tid, inum, nthreads); + ThrData *thr = fix->get_thr(tid); + thr->timer(Timer::START); + ev_setup_thr(eflag, vflag, nall, eatom, vatom, cvatom, thr); + + if (inum > 0) { + if (evflag) { + if (eflag) { + if (force->newton_bond) eval<1,1,1>(ifrom, ito, thr); + else eval<1,1,0>(ifrom, ito, thr); + } else { + if (force->newton_bond) eval<1,0,1>(ifrom, ito, thr); + else eval<1,0,0>(ifrom, ito, thr); + } + } else { + if (force->newton_bond) eval<0,0,1>(ifrom, ito, thr); + else eval<0,0,0>(ifrom, ito, thr); + } + } + thr->timer(Timer::BOND); + reduce_thr(this, eflag, vflag, thr); + } // end of omp parallel region +} + +template +void AngleCosineSquaredRestrictedOMP::eval(int nfrom, int nto, ThrData * const thr) +{ + int i1,i2,i3,n,type; + double delx1,dely1,delz1,delx2,dely2,delz2; + double eangle,f1[3],f3[3]; + double tk,rsq1,rsq2,r1,r2,c,a,a11,a12,a22; + + const auto * _noalias const x = (dbl3_t *) atom->x[0]; + auto * _noalias const f = (dbl3_t *) thr->get_f()[0]; + const int4_t * _noalias const anglelist = (int4_t *) neighbor->anglelist[0]; + const int nlocal = atom->nlocal; + eangle = 0.0; + + for (n = nfrom; n < nto; n++) { + i1 = anglelist[n].a; + i2 = anglelist[n].b; + i3 = anglelist[n].c; + type = anglelist[n].t; + + // 1st bond + + delx1 = x[i1].x - x[i2].x; + dely1 = x[i1].y - x[i2].y; + delz1 = x[i1].z - x[i2].z; + + rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; + r1 = sqrt(rsq1); + + // 2nd bond + + delx2 = x[i3].x - x[i2].x; + dely2 = x[i3].y - x[i2].y; + delz2 = x[i3].z - x[i2].z; + + rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; + r2 = sqrt(rsq2); + + // angle (cos and sin) + + c = delx1*delx2 + dely1*dely2 + delz1*delz2; + c /= r1*r2; + + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + + // force & energy + + double sq_sin = 1.0 - c * c; + double c0 = cos(theta0[type]); + + tk = k[type] * (c - c0) * (1.0 - c * c0) / (sq_sin * sq_sin); + + if (EFLAG) eangle = k[type] * (c - c0) * (c - c0) / sq_sin; + + a = 2.0 * tk; + a11 = a*c / rsq1; + a12 = -a / (r1*r2); + a22 = a*c / rsq2; + + f1[0] = a11*delx1 + a12*delx2; + f1[1] = a11*dely1 + a12*dely2; + f1[2] = a11*delz1 + a12*delz2; + f3[0] = a22*delx2 + a12*delx1; + f3[1] = a22*dely2 + a12*dely1; + f3[2] = a22*delz2 + a12*delz1; + + // apply force to each of 3 atoms + + if (NEWTON_BOND || i1 < nlocal) { + f[i1].x += f1[0]; + f[i1].y += f1[1]; + f[i1].z += f1[2]; + } + + if (NEWTON_BOND || i2 < nlocal) { + f[i2].x -= f1[0] + f3[0]; + f[i2].y -= f1[1] + f3[1]; + f[i2].z -= f1[2] + f3[2]; + } + + if (NEWTON_BOND || i3 < nlocal) { + f[i3].x += f3[0]; + f[i3].y += f3[1]; + f[i3].z += f3[2]; + } + + if (EVFLAG) ev_tally_thr(this,i1,i2,i3,nlocal,NEWTON_BOND,eangle,f1,f3, + delx1,dely1,delz1,delx2,dely2,delz2,thr); + } +} diff --git a/src/OPENMP/angle_cosine_squared_restricted_omp.h b/src/OPENMP/angle_cosine_squared_restricted_omp.h new file mode 100644 index 0000000000..b2ad545230 --- /dev/null +++ b/src/OPENMP/angle_cosine_squared_restricted_omp.h @@ -0,0 +1,46 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Axel Kohlmeyer (Temple U) +------------------------------------------------------------------------- */ + +#ifdef ANGLE_CLASS +// clang-format off +AngleStyle(cosine/squared/restricted/omp,AngleCosineSquaredRestrictedOMP); +// clang-format on +#else + +#ifndef LMP_ANGLE_COSINE_SQUARED_RESTRICTED_OMP_H +#define LMP_ANGLE_COSINE_SQUARED_RESTRICTED_OMP_H + +#include "angle_cosine_squared_restricted.h" +#include "thr_omp.h" + +namespace LAMMPS_NS { + +class AngleCosineSquaredRestrictedOMP : public AngleCosineSquaredRestricted, public ThrOMP { + + public: + AngleCosineSquaredRestrictedOMP(class LAMMPS *lmp); + void compute(int, int) override; + + private: + template + void eval(int ifrom, int ito, ThrData *const thr); +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/thermo.cpp b/src/thermo.cpp index 5ef5eb59b8..128f8573cf 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -1515,15 +1515,16 @@ void Thermo::compute_fix() int m = field2index[ifield]; Fix *fix = fixes[m]; - // check for out-of-range access if vector/array is variable length - if (argindex1[ifield] == 0) { dvalue = fix->compute_scalar(); if (normflag && fix->extscalar) dvalue /= natoms; } else if (argindex2[ifield] == 0) { - if (fix->size_vector_variable && argindex1[ifield] > fix->size_vector) - error->all(FLERR, "Thermo fix vector is accessed out-of-range"); - dvalue = fix->compute_vector(argindex1[ifield] - 1); + + // if index exceeds variable vector length, use a zero value + // this can be useful if vector length is not known a priori + + if (fix->size_vector_variable && argindex1[ifield] > fix->size_vector) dvalue = 0.0; + else dvalue = fix->compute_vector(argindex1[ifield] - 1); if (normflag) { if (fix->extvector == 0) return; @@ -1533,9 +1534,12 @@ void Thermo::compute_fix() dvalue /= natoms; } } else { - if (fix->size_array_rows_variable && argindex1[ifield] > fix->size_array_rows) - error->all(FLERR, "Thermo fix array is accessed out-of-range"); - dvalue = fix->compute_array(argindex1[ifield] - 1, argindex2[ifield] - 1); + + // if index exceeds variable array rows, use a zero value + // this can be useful if array size is not known a priori + + if (fix->size_array_rows_variable && argindex1[ifield] > fix->size_array_rows) dvalue = 0.0; + else dvalue = fix->compute_array(argindex1[ifield] - 1, argindex2[ifield] - 1); if (normflag && fix->extarray) dvalue /= natoms; } } diff --git a/src/variable.cpp b/src/variable.cpp index 89d15d38f5..be9239c027 100644 --- a/src/variable.cpp +++ b/src/variable.cpp @@ -1809,12 +1809,16 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) if (!fix->vector_flag) print_var_error(FLERR,"Mismatched fix in variable formula",ivar); - if (index1 > fix->size_vector) + if (index1 > fix->size_vector && fix->size_vector_variable == 0) print_var_error(FLERR,"Variable formula fix vector is accessed out-of-range",ivar,0); if (update->whichflag > 0 && update->ntimestep % fix->global_freq) print_var_error(FLERR,"Fix in variable not computed at a compatible time",ivar); - value1 = fix->compute_vector(index1-1); + // if index exceeds variable vector length, use a zero value + // this can be useful if vector length is not known a priori + + if (fix->size_vector_variable && index1 > fix->size_vector) value1 = 0.0; + else value1 = fix->compute_vector(index1-1); argstack[nargstack++] = value1; // f_ID[i][j] = scalar from global array @@ -1823,14 +1827,18 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) if (!fix->array_flag) print_var_error(FLERR,"Mismatched fix in variable formula",ivar); - if (index1 > fix->size_array_rows) + if (index1 > fix->size_array_rows && fix->size_array_rows_variable == 0) print_var_error(FLERR,"Variable formula fix array is accessed out-of-range",ivar,0); if (index2 > fix->size_array_cols) print_var_error(FLERR,"Variable formula fix array is accessed out-of-range",ivar,0); if (update->whichflag > 0 && update->ntimestep % fix->global_freq) print_var_error(FLERR,"Fix in variable not computed at a compatible time",ivar); - value1 = fix->compute_array(index1-1,index2-1); + // if index exceeds variable array rows, use a zero value + // this can be useful if array size is not known a priori + + if (fix->size_array_rows_variable && index1 > fix->size_array_rows) value1 = 0.0; + else value1 = fix->compute_array(index1-1,index2-1); argstack[nargstack++] = value1; // F_ID[i] = scalar element of per-atom vector, note uppercase "F" diff --git a/unittest/force-styles/tests/angle-cosine_squared_restricted.yaml b/unittest/force-styles/tests/angle-cosine_squared_restricted.yaml new file mode 100644 index 0000000000..400babb3c0 --- /dev/null +++ b/unittest/force-styles/tests/angle-cosine_squared_restricted.yaml @@ -0,0 +1,88 @@ +--- +lammps_version: 7 Feb 2024 +tags: +date_generated: Fri Mar 29 21:14:05 2024 +epsilon: 2.5e-13 +skip_tests: +prerequisites: ! | + atom full + angle cosine/squared/restricted +pre_commands: ! "" +post_commands: ! "" +input_file: in.fourmol +angle_style: cosine/squared/restricted +angle_coeff: ! | + 1 75.0 110.1 + 2 45.0 111.0 + 3 50.0 120.0 + 4 100.0 108.5 +equilibrium: 4 1.9216075064457567 1.9373154697137058 2.0943951023931953 1.8936822384138476 +extract: ! "" +natoms: 29 +init_energy: 43.16721849625078 +init_stress: ! |2- + 9.1896481583329006e+01 -9.9073798939810871e+01 7.1773173564818959e+00 1.0147556441011812e+02 -3.1614702670042060e+01 6.3320341995277474e-01 +init_forces: ! |2 + 1 5.4616646260609812e+01 6.4912761997628223e+00 -3.3960574186374622e+01 + 2 -1.1520735029225428e+00 -9.3280585589864415e+00 -7.5009000871307308e+00 + 3 -1.3676356546913162e+01 4.9291191335900074e+01 5.1495824934074705e+01 + 4 -1.7548953236977070e+01 -2.4649435647138144e+01 1.4994156277002995e+01 + 5 -4.7191706109161828e+01 -4.7108395361010651e+01 -2.3419810899443089e+01 + 6 5.4456046768792298e+01 -1.0628221895648938e+01 -4.4205117451573471e+01 + 7 -1.5479402058953555e+01 1.3800535101019676e+01 -3.3710945974014352e+00 + 8 -4.5010195114618039e+00 -1.9041823240486167e+01 1.0790072547105123e+02 + 9 -1.5817751557854660e+01 1.6558437937477478e+01 -4.5571617839569933e-01 + 10 4.9597370008417769e+00 4.0344934206635401e+01 -7.9652184287503914e+01 + 11 -1.9904147169399383e+01 -5.8247456312686090e+00 8.4471337004560691e+00 + 12 -1.0750141644440905e+01 -1.0959289684407434e+01 -9.7375245407417310e+00 + 13 3.0110873732763217e+00 7.6125328583508187e+00 -1.8728278765280493e-02 + 14 1.4804534890961024e+01 -9.4162987076082505e+00 -6.4475703726793121e+00 + 15 1.2405121083152171e+01 -6.1967835746148943e+00 1.8294821492807777e+01 + 16 3.1352949560138343e-01 1.1035513814862394e+01 6.5673781465816949e+00 + 17 1.4548484648501341e+00 -1.9813691528391315e+00 1.0691808580347990e+00 + 18 1.5557496930974857e-01 1.6257536594477804e+00 -6.5978344782819773e+00 + 19 -1.9049016835509454e+00 -2.3731208324249899e+00 2.8255706424396045e+00 + 20 1.7493267142411968e+00 7.4736717297720956e-01 3.7722638358423728e+00 + 21 3.3279741492378183e+00 4.1660737555188891e+00 -1.2450454543863602e+01 + 22 -5.9261551336076748e+00 -4.0404494135493199e+00 4.4014087371258892e+00 + 23 2.5981809843698564e+00 -1.2562434196956951e-01 8.0490458067377126e+00 + 24 -1.2888190495929852e+00 5.0097160214605641e+00 -2.9825927599219058e+00 + 25 -1.0338174598680057e+00 -3.4636797875161851e+00 5.6646346770669309e-01 + 26 2.3226365094609909e+00 -1.5460362339443785e+00 2.4161292922152127e+00 + 27 -2.8563709770508883e-01 2.3448568102582543e+00 -7.9873605851552354e-01 + 28 -6.2456251656312733e-01 -1.4020243449848773e+00 3.2568906896433880e-02 + 29 9.1019961426821616e-01 -9.4283246527337705e-01 7.6616715161908966e-01 +run_energy: 42.9022488368839 +run_stress: ! |2- + 9.0388148558231080e+01 -9.8306784914711429e+01 7.9186363564803592e+00 9.9981546403593825e+01 -3.2129110268269365e+01 -3.6431126570529515e-01 +run_forces: ! |2 + 1 5.3981820216746733e+01 6.9247252423091306e+00 -3.3185401895138106e+01 + 2 -1.1944338217682868e+00 -9.7346413974706021e+00 -7.8615437717021557e+00 + 3 -1.2688555593426315e+01 4.8769286528132106e+01 5.0839055814558762e+01 + 4 -1.7617326624105857e+01 -2.4269207832769677e+01 1.5260301967540304e+01 + 5 -4.7323547490196731e+01 -4.6894886975434055e+01 -2.3440133100826486e+01 + 6 5.4382423772065998e+01 -1.0696866841386006e+01 -4.4201170087829055e+01 + 7 -1.5492082837122805e+01 1.3818527237625043e+01 -3.3812483934830673e+00 + 8 -4.4634623862144078e+00 -1.9093779476665727e+01 1.0771575923879119e+02 + 9 -1.5843321387394678e+01 1.6589965876806524e+01 -4.4411228983080475e-01 + 10 5.0768945646127097e+00 4.0040817044035947e+01 -7.9450255761393905e+01 + 11 -1.9901406164815324e+01 -5.7840333581333585e+00 8.4316140350776880e+00 + 12 -1.0895919045638994e+01 -1.1165043741641941e+01 -9.2828865105939240e+00 + 13 2.9272358235905029e+00 7.4354786772387680e+00 1.0142661372278194e-01 + 14 1.4792360939950459e+01 -8.9353318022993182e+00 -6.3910772582380133e+00 + 15 1.2418288954456656e+01 -6.0153396231171570e+00 1.7743757855168997e+01 + 16 4.0781740822800217e-01 1.0963878023768416e+01 6.4917802461985135e+00 + 17 1.4332136710323162e+00 -1.9535475809980913e+00 1.0541332979772831e+00 + 18 1.4023773325248556e-01 1.4086440965467837e+00 -5.6940807148906991e+00 + 19 -1.6400859684012328e+00 -2.0475288381692676e+00 2.4445067235235367e+00 + 20 1.4998482351487472e+00 6.3888474162248376e-01 3.2495739913671629e+00 + 21 3.1432629198624089e+00 3.9072243542777101e+00 -1.1745189375890282e+01 + 22 -5.5777858442915553e+00 -3.7895069817029867e+00 4.1606823590513633e+00 + 23 2.4345229244291464e+00 -1.1771737257472337e-01 7.5845070168389190e+00 + 24 -1.1529613091971900e+00 4.4757566436680412e+00 -2.6619956178204864e+00 + 25 -9.1659948782954015e-01 -3.0936731063721679e+00 5.0626909884369264e-01 + 26 2.0695607970267300e+00 -1.3820835372958733e+00 2.1557265189767940e+00 + 27 -2.6669219302298908e-01 2.1720439184593126e+00 -7.3558544897318312e-01 + 28 -5.7634641435086742e-01 -1.2987031699370619e+00 2.8222187310339070e-02 + 29 8.4303860737385650e-01 -8.7334074852225052e-01 7.0736326166284402e-01 +... diff --git a/unittest/force-styles/tests/dihedral-cosine_squared_restricted.yaml b/unittest/force-styles/tests/dihedral-cosine_squared_restricted.yaml new file mode 100644 index 0000000000..f67a093017 --- /dev/null +++ b/unittest/force-styles/tests/dihedral-cosine_squared_restricted.yaml @@ -0,0 +1,88 @@ +--- +lammps_version: 7 Feb 2024 +tags: +date_generated: Sat Apr 13 11:41:16 2024 +epsilon: 2.5e-13 +skip_tests: +prerequisites: ! | + atom full + dihedral cosine/squared/restricted +pre_commands: ! "" +post_commands: ! "" +input_file: in.fourmol +dihedral_style: cosine/squared/restricted +dihedral_coeff: ! | + 1 17.0 10.0 + 2 14.0 20.0 + 3 15.0 -10.0 + 4 12.0 0.0 + 5 11.0 45.0 +extract: ! "" +natoms: 29 +init_energy: 10643.96352037142 +init_stress: ! |2- + 8.2655807554971598e+03 6.2324906154519340e+03 -1.4498071370948830e+04 1.0953420376656040e+04 1.1057094124317168e+03 -8.3313158421382595e+03 +init_forces: ! |2 + 1 -2.9998360221989606e+03 1.5614240759467791e+03 -1.4993484270227148e+03 + 2 2.3464901141409141e+03 -1.2117018270409203e+03 1.1464619463537790e+03 + 3 -8.4367569419635300e+04 -8.8610214478341717e+04 1.3641091254798606e+04 + 4 -7.7728947837396186e+02 -1.1834661410404278e+03 6.2790819006201627e+02 + 5 1.9271341641256136e+03 -1.1638057759737167e+03 2.7159431719235641e+01 + 6 4.4134274351326589e+04 5.5926488430371595e+04 2.3644507487725932e+03 + 7 5.6253242167919998e+03 6.0326849037643433e+03 -9.8937413117578353e+02 + 8 1.5398571438042008e+05 1.4826543321068442e+05 -4.6786535127649200e+04 + 9 5.2872438366384667e+03 3.9367087924112984e+03 -2.1528982310842534e+03 + 10 -1.2076914217812811e+05 -1.1928499037053883e+05 3.1759537647021025e+04 + 11 -2.7869283077129630e+03 -2.7624462832092408e+03 2.3692905638286411e+03 + 12 -1.7137278451679694e+02 -4.0643861446907772e+02 -8.4482327773661689e+02 + 13 1.8033469195161777e+02 4.1452681735267817e+02 8.7002410666914204e+02 + 14 1.0785063285146067e+02 -1.5007679697175007e+02 -5.3363703034369962e+01 + 15 -1.2637123349071543e+02 4.2783038317466065e+01 -1.2928289731925594e+02 + 16 -1.3974766748247812e+03 -1.2440022591573602e+03 -3.1834363647575998e+02 + 17 -1.9838028936516241e+02 -1.6290672210557804e+02 -3.1954457727099452e+01 + 18 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 19 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 20 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 21 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 22 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 23 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 24 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 25 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 26 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 27 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 28 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 29 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 +run_energy: 4829.55350497352 +run_stress: ! |2- + 3.3195623927709175e+03 3.0229558326855845e+03 -6.3425182254565516e+03 4.3274493033758445e+03 3.3064220676415192e+02 -1.1283617668722684e+03 +run_forces: ! |2 + 1 -1.6686851537334303e+03 9.1157319328820222e+02 -8.8504748962571853e+02 + 2 1.2605943337473971e+03 -6.8570419879992460e+02 6.5988055106439663e+02 + 3 -2.4907715351961961e+03 1.5938154791899667e+03 -4.0018442338176410e+03 + 4 -3.9608971224931270e+03 -6.6732751207320443e+03 2.0299175931026020e+03 + 5 6.9745717717678440e+02 -1.1129862611363344e+03 -1.9316966111163129e+02 + 6 -9.8181046153612078e+03 -1.0177815528067989e+04 6.1127410299948697e+03 + 7 5.2268917470447668e+03 6.8715264070828716e+03 -1.6421183728373735e+03 + 8 2.4089311767862888e+04 2.3281334463282841e+04 -2.9014860835053078e+03 + 9 1.8251257887454533e+03 5.6771778492275894e+02 -8.4486828193101837e+02 + 10 -8.7189544342799436e+03 -8.9637840944091586e+03 2.1674670163757610e+03 + 11 -6.1742160434459652e+02 -5.4300765195316922e+02 8.4050088549381860e+02 + 12 -1.6333025955329197e+02 -2.0286364513967663e+02 -5.7218852607845849e+02 + 13 1.1676829617621111e+02 2.8391630256639763e+02 6.2310026378133114e+02 + 14 1.2782423911823358e+02 -1.7234482742516084e+02 -6.2457926948088009e+01 + 15 -1.0593900405480906e+02 3.6959553351961858e+01 -1.1139986642128649e+02 + 16 -5.7436221922065306e+03 -4.9695083643645858e+03 -1.2100596896535424e+03 + 17 -5.6247428648599033e+01 -4.5553491656954620e+01 -8.9672078827135397e+00 + 18 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 19 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 20 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 21 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 22 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 23 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 24 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 25 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 26 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 27 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 28 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 29 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 +...