Merge branch 'develop' into collected-small-changes
# Conflicts: # python/lammps/pylammps.py
This commit is contained in:
Axel Kohlmeyer
@ -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
80
doc/src/angle_cosine_squared_restricted.rst
Normal 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).
|
||||
@ -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
71
doc/src/dihedral_cosine_squared_restricted.rst
Normal 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).
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
4
src/.gitignore
vendored
4
src/.gitignore
vendored
@ -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
298
src/EXTRA-MOLECULE/angle_cosine_squared_restricted.cpp
Normal 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
49
src/EXTRA-MOLECULE/angle_cosine_squared_restricted.h
Normal 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
395
src/EXTRA-MOLECULE/dihedral_cosine_squared_restricted.cpp
Normal 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
47
src/EXTRA-MOLECULE/dihedral_cosine_squared_restricted.h
Normal 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
170
src/OPENMP/angle_cosine_squared_restricted_omp.cpp
Normal 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
46
src/OPENMP/angle_cosine_squared_restricted_omp.h
Normal 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
|
||||
@ -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;
|
||||
}
|
||||
}
|
||||
|
||||
@ -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"
|
||||
|
||||
@ -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
|
||||
...
|
||||
@ -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
|
||||
...
|
||||
Reference in New Issue
Block a user