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

Merge pull request #4113 from evoyiatzis/restricted-angle

Browse Source 
Restricted squared cosine angle and dihedral potentials
This commit is contained in:
Axel Kohlmeyer
2024-04-15 13:13:11 -04:00
committed by GitHub
parent 1346be4168 5af3626cf7
commit d34a9bd33e
15 changed files with 1348 additions and 2 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
doc/src/Commands_bond.rst
View File

@ -89,6 +89,7 @@ OPT.
* :doc:`cosine/shift (o) <angle_cosine_shift>`
* :doc:`cosine/shift/exp (o) <angle_cosine_shift_exp>`
* :doc:`cosine/squared (o) <angle_cosine_squared>`
* :doc:`cosine/squared/restricted (o) <angle_cosine_squared_restricted>`
* :doc:`cross <angle_cross>`
* :doc:`dipole (o) <angle_dipole>`
* :doc:`fourier (o) <angle_fourier>`
@ -127,6 +128,7 @@ OPT.
* :doc:`charmmfsw (k) <dihedral_charmm>`
* :doc:`class2 (ko) <dihedral_class2>`
* :doc:`cosine/shift/exp (o) <dihedral_cosine_shift_exp>`
* :doc:`cosine/squared/restricted <dihedral_cosine_squared_restricted>`
* :doc:`fourier (io) <dihedral_fourier>`
* :doc:`harmonic (iko) <dihedral_harmonic>`
* :doc:`helix (o) <dihedral_helix>`

80
doc/src/angle_cosine_squared_restricted.rst Normal file
View File

@ -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) <restricted-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 <angle_coeff>` command as in the example above, or in
the data file or restart files read by the :doc:`read_data <read_data>`
or :doc:`read_restart <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 <Build_package>` doc page
for more info.
Related commands
""""""""""""""""
:doc:`angle_coeff <angle_coeff>`
Default
"""""""
none
----------
.. _restricted-Bulacu:
**(Bulacu)** Bulacu, Goga, Zhao, Rossi, Monticelli, Periole, Tieleman, Marrink, J Chem Theory Comput, 9, 3282-3292
(2013).

3
doc/src/angle_style.rst
View File

@ -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_shift>` - angle cosine with a shift
* :doc:`cosine/shift/exp <angle_cosine_shift_exp>` - cosine with shift and exponential term in spring constant
* :doc:`cosine/squared <angle_cosine_squared>` - angle with cosine squared term
* :doc:`cosine/squared/restricted <angle_cosine_squared_restricted>` - angle with restricted cosine squared term
* :doc:`cross <angle_cross>` - cross term coupling angle and bond lengths
* :doc:`dipole <angle_dipole>` - angle that controls orientation of a point dipole
* :doc:`fourier <angle_fourier>` - angle with multiple cosine terms

71
doc/src/dihedral_cosine_squared_restricted.rst Normal file
View File

@ -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) <restricted-Bul>` 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 <dihedral_coeff>` command as in the example above, or in
the data file or restart files read by the :doc:`read_data <read_data>`
or :doc:`read_restart <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 <Build_package>` doc page
for more info.
Related commands
""""""""""""""""
:doc:`dihedral_coeff <dihedral_coeff>`
Default
"""""""
none
----------
.. _restricted-Bul:
**(Bulacu)** Bulacu, Goga, Zhao, Rossi, Monticelli, Periole, Tieleman, Marrink, J Chem Theory Comput, 9, 3282-3292
(2013).

3
doc/src/dihedral_style.rst
View File

@ -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 <dihedral_charmm>` - CHARMM dihedral with force switching
* :doc:`class2 <dihedral_class2>` - COMPASS (class 2) dihedral
* :doc:`cosine/shift/exp <dihedral_cosine_shift_exp>` - dihedral with exponential in spring constant
* :doc:`cosine/squared/restricted <dihedral_cosine_squared_restricted>` - squared cosine dihedral with restricted term
* :doc:`fourier <dihedral_fourier>` - dihedral with multiple cosine terms
* :doc:`harmonic <dihedral_harmonic>` - harmonic dihedral
* :doc:`helix <dihedral_helix>` - helix dihedral

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

@ -380,6 +380,7 @@ btype
buckPlusAttr
buf
builtin
Bulacu
Bulatov
Bureekaew
burlywood
@ -1302,6 +1303,7 @@ gmres
gname
gneb
GNEB
Goga
Goldfarb
Gompper
Gonzalez-Melchor
@ -2261,6 +2263,7 @@ Montalenti
Monterey
Montero
Monti
Monticelli
Mora
Morefoo
Morfill
@ -2789,6 +2792,7 @@ peridynamic
Peridynamic
peridynamics
Peridynamics
Periole
perl
permittivity
perp
@ -3204,6 +3208,7 @@ Ronevich
Rosati
Rosato
Rosenberger
Rossi
Rossky
rosybrown
rotationally
@ -3647,6 +3652,7 @@ thrid
ThunderX
thylakoid
THz
Tieleman
Tigran
Tij
Tildesley

4
src/.gitignore vendored
View File

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

298
src/EXTRA-MOLECULE/angle_cosine_squared_restricted.cpp Normal file
View File

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

49
src/EXTRA-MOLECULE/angle_cosine_squared_restricted.h Normal file
View File

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

395
src/EXTRA-MOLECULE/dihedral_cosine_squared_restricted.cpp Normal file
View File

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

47
src/EXTRA-MOLECULE/dihedral_cosine_squared_restricted.h Normal file
View File

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

170
src/OPENMP/angle_cosine_squared_restricted_omp.cpp Normal file
View File

@ -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 <cmath>
#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 <int EVFLAG, int EFLAG, int NEWTON_BOND>
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);
}
}

46
src/OPENMP/angle_cosine_squared_restricted_omp.h Normal file
View File

@ -0,0 +1,46 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
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 <int EVFLAG, int EFLAG, int NEWTON_BOND>
void eval(int ifrom, int ito, ThrData *const thr);
};
} // namespace LAMMPS_NS
#endif
#endif

88
unittest/force-styles/tests/angle-cosine_squared_restricted.yaml Normal file
View File

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

88
unittest/force-styles/tests/dihedral-cosine_squared_restricted.yaml Normal file
View File

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