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

Merge branch 'develop' into general-triclinic

Browse Source
This commit is contained in:
Axel Kohlmeyer
2024-01-25 09:04:46 -05:00
parent 8cf2705d68 54794a45de
commit 9cae812c51
38 changed files with 2271 additions and 177 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
doc/src/Commands_bond.rst
View File

@ -124,7 +124,7 @@ OPT.
* *
* *
* :doc:`charmm (iko) <dihedral_charmm>` * :doc:`charmm (iko) <dihedral_charmm>`
* :doc:`charmmfsw <dihedral_charmm>` * :doc:`charmmfsw (k) <dihedral_charmm>`
* :doc:`class2 (ko) <dihedral_class2>` * :doc:`class2 (ko) <dihedral_class2>`
* :doc:`cosine/shift/exp (o) <dihedral_cosine_shift_exp>` * :doc:`cosine/shift/exp (o) <dihedral_cosine_shift_exp>`
* :doc:`fourier (io) <dihedral_fourier>` * :doc:`fourier (io) <dihedral_fourier>`

2
doc/src/Commands_pair.rst
View File

@ -146,7 +146,7 @@ OPT.
* :doc:`lj/charmm/coul/long/soft (o) <pair_fep_soft>` * :doc:`lj/charmm/coul/long/soft (o) <pair_fep_soft>`
* :doc:`lj/charmm/coul/msm (o) <pair_charmm>` * :doc:`lj/charmm/coul/msm (o) <pair_charmm>`
* :doc:`lj/charmmfsw/coul/charmmfsh <pair_charmm>` * :doc:`lj/charmmfsw/coul/charmmfsh <pair_charmm>`
* :doc:`lj/charmmfsw/coul/long <pair_charmm>` * :doc:`lj/charmmfsw/coul/long (k) <pair_charmm>`
* :doc:`lj/class2 (gko) <pair_class2>` * :doc:`lj/class2 (gko) <pair_class2>`
* :doc:`lj/class2/coul/cut (ko) <pair_class2>` * :doc:`lj/class2/coul/cut (ko) <pair_class2>`
* :doc:`lj/class2/coul/cut/soft <pair_fep_soft>` * :doc:`lj/class2/coul/cut/soft <pair_fep_soft>`

4
doc/src/angle_charmm.rst
View File

@ -70,7 +70,9 @@ for more info.
Related commands Related commands
"""""""""""""""" """"""""""""""""
:doc:`angle_coeff <angle_coeff>` :doc:`angle_coeff <angle_coeff>`, :doc:`pair_style lj/charmm variants <pair_charmm>`,
:doc:`dihedral_style charmm <dihedral_charmm>`,
:doc:`dihedral_style charmmfsw <dihedral_charmm>`, :doc:`fix cmap <fix_cmap>`
Default Default
""""""" """""""

19
doc/src/angle_lepton.rst
View File

@ -11,7 +11,16 @@ Syntax
.. code-block:: LAMMPS .. code-block:: LAMMPS
angle_style lepton angle_style style args
* style = *lepton*
* args = optional arguments
.. parsed-literal::
args = *auto_offset* or *no_offset*
*auto_offset* = offset the potential energy so that the value at theta0 is 0.0 (default)
*no_offset* = do not offset the potential energy
Examples Examples
"""""""" """"""""
@ -19,6 +28,7 @@ Examples
.. code-block:: LAMMPS .. code-block:: LAMMPS
angle_style lepton angle_style lepton
angle_style lepton no_offset
angle_coeff 1 120.0 "k*theta^2; k=250.0" angle_coeff 1 120.0 "k*theta^2; k=250.0"
angle_coeff 2 90.0 "k2*theta^2 + k3*theta^3 + k4*theta^4; k2=300.0; k3=-100.0; k4=50.0" angle_coeff 2 90.0 "k2*theta^2 + k3*theta^3 + k4*theta^4; k2=300.0; k3=-100.0; k4=50.0"
@ -41,6 +51,13 @@ angle coefficient. For example `"200.0*theta^2"` represents a
U_{angle,i} = K (\theta_i - \theta_0)^2 = K \theta^2 \qquad \theta = \theta_i - \theta_0 U_{angle,i} = K (\theta_i - \theta_0)^2 = K \theta^2 \qquad \theta = \theta_i - \theta_0
.. versionchanged:: TBD
By default the potential energy U is shifted so that the value U is 0.0
for $theta = theta_0$. This is equivalent to using the optional keyword
*auto_offset*. When using the keyword *no_offset* instead, the
potential energy is not shifted.
The `Lepton library <https://simtk.org/projects/lepton>`_, that the The `Lepton library <https://simtk.org/projects/lepton>`_, that the
*lepton* angle style interfaces with, evaluates this expression string *lepton* angle style interfaces with, evaluates this expression string
at run time to compute the pairwise energy. It also creates an at run time to compute the pairwise energy. It also creates an

19
doc/src/bond_lepton.rst
View File

@ -11,7 +11,16 @@ Syntax
.. code-block:: LAMMPS .. code-block:: LAMMPS
bond_style lepton bond_style style args
* style = *lepton*
* args = optional arguments
.. parsed-literal::
args = *auto_offset* or *no_offset*
*auto_offset* = offset the potential energy so that the value at r0 is 0.0 (default)
*no_offset* = do not offset the potential energy
Examples Examples
"""""""" """"""""
@ -19,6 +28,7 @@ Examples
.. code-block:: LAMMPS .. code-block:: LAMMPS
bond_style lepton bond_style lepton
bond_style lepton no_offset
bond_coeff 1 1.5 "k*r^2; k=250.0" bond_coeff 1 1.5 "k*r^2; k=250.0"
bond_coeff 2 1.1 "k2*r^2 + k3*r^3 + k4*r^4; k2=300.0; k3=-100.0; k4=50.0" bond_coeff 2 1.1 "k2*r^2 + k3*r^3 + k4*r^4; k2=300.0; k3=-100.0; k4=50.0"
@ -40,6 +50,13 @@ constant *K* of 200.0 energy units:
U_{bond,i} = K (r_i - r_0)^2 = K r^2 \qquad r = r_i - r_0 U_{bond,i} = K (r_i - r_0)^2 = K r^2 \qquad r = r_i - r_0
.. versionchanged:: TBD
By default the potential energy U is shifted so that he value U is 0.0
for $r = r_0$. This is equivalent to using the optional keyword
*auto_offset*. When using the keyword *no_offset* instead, the
potential energy is not shifted.
The `Lepton library <https://simtk.org/projects/lepton>`_, that the The `Lepton library <https://simtk.org/projects/lepton>`_, that the
*lepton* bond style interfaces with, evaluates this expression string at *lepton* bond style interfaces with, evaluates this expression string at
run time to compute the pairwise energy. It also creates an analytical run time to compute the pairwise energy. It also creates an analytical

7
doc/src/dihedral_charmm.rst
View File

@ -3,6 +3,7 @@
.. index:: dihedral_style charmm/kk .. index:: dihedral_style charmm/kk
.. index:: dihedral_style charmm/omp .. index:: dihedral_style charmm/omp
.. index:: dihedral_style charmmfsw .. index:: dihedral_style charmmfsw
.. index:: dihedral_style charmmfsw/kk
dihedral_style charmm command dihedral_style charmm command
============================= =============================
@ -12,6 +13,8 @@ Accelerator Variants: *charmm/intel*, *charmm/kk*, *charmm/omp*
dihedral_style charmmfsw command dihedral_style charmmfsw command
================================ ================================
Accelerator Variants: *charmmfsw/kk*
Syntax Syntax
"""""" """"""
@ -144,7 +147,9 @@ for more info.
Related commands Related commands
"""""""""""""""" """"""""""""""""
:doc:`dihedral_coeff <dihedral_coeff>` :doc:`dihedral_coeff <dihedral_coeff>`,
:doc:`pair_style lj/charmm variants <pair_charmm>`,
:doc:`angle_style charmm <angle_charmm>`, :doc:`fix cmap <fix_cmap>`
Default Default
""""""" """""""

2
doc/src/molecule.rst
View File

@ -376,7 +376,7 @@ not listed, the default diameter of each atom in the molecule is 1.0.
---------- ----------
..versionadded:: TBD .. versionadded:: TBD
*Dipoles* section: *Dipoles* section:

7
doc/src/pair_charmm.rst
View File

@ -16,6 +16,7 @@
.. index:: pair_style lj/charmm/coul/msm/omp .. index:: pair_style lj/charmm/coul/msm/omp
.. index:: pair_style lj/charmmfsw/coul/charmmfsh .. index:: pair_style lj/charmmfsw/coul/charmmfsh
.. index:: pair_style lj/charmmfsw/coul/long .. index:: pair_style lj/charmmfsw/coul/long
.. index:: pair_style lj/charmmfsw/coul/long/kk
pair_style lj/charmm/coul/charmm command pair_style lj/charmm/coul/charmm command
======================================== ========================================
@ -43,6 +44,8 @@ pair_style lj/charmmfsw/coul/charmmfsh command
pair_style lj/charmmfsw/coul/long command pair_style lj/charmmfsw/coul/long command
========================================= =========================================
Accelerator Variants: *lj/charmmfsw/coul/long/kk*
Syntax Syntax
"""""" """"""
@ -281,7 +284,9 @@ page for more info.
Related commands Related commands
"""""""""""""""" """"""""""""""""
:doc:`pair_coeff <pair_coeff>` :doc:`pair_coeff <pair_coeff>`, :doc:`angle_style charmm <angle_charmm>`,
:doc:`dihedral_style charmm <dihedral_charmm>`,
:doc:`dihedral_style charmmfsw <dihedral_charmm>`, :doc:`fix cmap <fix_cmap>`
Default Default
""""""" """""""

6
doc/src/pair_lepton.rst
View File

@ -72,7 +72,7 @@ interactions between particles which depend on the distance and have a
cutoff. The potential function must be provided as an expression string cutoff. The potential function must be provided as an expression string
using "r" as the distance variable. With pair style *lepton/coul* one using "r" as the distance variable. With pair style *lepton/coul* one
may additionally reference the charges of the two atoms of the pair with may additionally reference the charges of the two atoms of the pair with
"qi" and "qj", respectively. With pair style *lepton/coul* one may "qi" and "qj", respectively. With pair style *lepton/sphere* one may
instead reference the radii of the two atoms of the pair with "radi" and instead reference the radii of the two atoms of the pair with "radi" and
"radj", respectively; this is half of the diameter that can be set in "radj", respectively; this is half of the diameter that can be set in
:doc:`data files <read_data>` or the :doc:`set command <set>`. :doc:`data files <read_data>` or the :doc:`set command <set>`.
@ -166,8 +166,8 @@ mixing. Thus, expressions for *all* I,J pairs must be specified
explicitly. explicitly.
Only pair style *lepton* supports the :doc:`pair_modify shift <pair_modify>` Only pair style *lepton* supports the :doc:`pair_modify shift <pair_modify>`
option for shifting the energy of the pair interaction so that it is option for shifting the potential energy of the pair interaction so that
0 at the cutoff, pair styles *lepton/coul* and *lepton/sphere* do *not*. it is 0 at the cutoff, pair styles *lepton/coul* and *lepton/sphere* do *not*.
The :doc:`pair_modify table <pair_modify>` options are not relevant for The :doc:`pair_modify table <pair_modify>` options are not relevant for
the these pair styles. the these pair styles.

4
src/KOKKOS/Install.sh
View File

@ -106,6 +106,8 @@ action compute_temp_kokkos.cpp
action compute_temp_kokkos.h action compute_temp_kokkos.h
action dihedral_charmm_kokkos.cpp dihedral_charmm.cpp action dihedral_charmm_kokkos.cpp dihedral_charmm.cpp
action dihedral_charmm_kokkos.h dihedral_charmm.h action dihedral_charmm_kokkos.h dihedral_charmm.h
action dihedral_charmmfsw_kokkos.cpp dihedral_charmmfsw.cpp
action dihedral_charmmfsw_kokkos.h dihedral_charmmfsw.h
action dihedral_class2_kokkos.cpp dihedral_class2.cpp action dihedral_class2_kokkos.cpp dihedral_class2.cpp
action dihedral_class2_kokkos.h dihedral_class2.h action dihedral_class2_kokkos.h dihedral_class2.h
action dihedral_harmonic_kokkos.cpp dihedral_harmonic.cpp action dihedral_harmonic_kokkos.cpp dihedral_harmonic.cpp
@ -310,6 +312,8 @@ action pair_lj_charmm_coul_charmm_kokkos.cpp pair_lj_charmm_coul_charmm.cpp
action pair_lj_charmm_coul_charmm_kokkos.h pair_lj_charmm_coul_charmm.h action pair_lj_charmm_coul_charmm_kokkos.h pair_lj_charmm_coul_charmm.h
action pair_lj_charmm_coul_long_kokkos.cpp pair_lj_charmm_coul_long.cpp action pair_lj_charmm_coul_long_kokkos.cpp pair_lj_charmm_coul_long.cpp
action pair_lj_charmm_coul_long_kokkos.h pair_lj_charmm_coul_long.h action pair_lj_charmm_coul_long_kokkos.h pair_lj_charmm_coul_long.h
action pair_lj_charmmfsw_coul_long_kokkos.cpp pair_lj_charmmfsw_coul_long.cpp
action pair_lj_charmmfsw_coul_long_kokkos.h pair_lj_charmmfsw_coul_long.h
action pair_lj_class2_coul_cut_kokkos.cpp pair_lj_class2_coul_cut.cpp action pair_lj_class2_coul_cut_kokkos.cpp pair_lj_class2_coul_cut.cpp
action pair_lj_class2_coul_cut_kokkos.h pair_lj_class2_coul_cut.h action pair_lj_class2_coul_cut_kokkos.h pair_lj_class2_coul_cut.h
action pair_lj_class2_coul_long_kokkos.cpp pair_lj_class2_coul_long.cpp action pair_lj_class2_coul_long_kokkos.cpp pair_lj_class2_coul_long.cpp

1
src/KOKKOS/atom_map_kokkos.cpp
View File

@ -273,6 +273,7 @@ void AtomKokkos::map_set()
error->one(FLERR,"Failed to insert into Kokkos hash atom map"); error->one(FLERR,"Failed to insert into Kokkos hash atom map");
k_sametag.modify_device(); k_sametag.modify_device();
k_sametag.sync_host();
if (map_style == MAP_ARRAY) if (map_style == MAP_ARRAY)
k_map_array.modify_device(); k_map_array.modify_device();

815
src/KOKKOS/dihedral_charmmfsw_kokkos.cpp Normal file
View File

@ -0,0 +1,815 @@
// 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: Mitch Murphy (alphataubio)
Based on serial dihedral_charmmfsw.cpp lj-fsw sections (force-switched)
provided by Robert Meissner and Lucio Colombi Ciacchi of Bremen
University, Germany, with additional assistance from
Robert A. Latour, Clemson University.
------------------------------------------------------------------------- */
#include "dihedral_charmmfsw_kokkos.h"
#include "atom_kokkos.h"
#include "atom_masks.h"
#include "error.h"
#include "force.h"
#include "kokkos.h"
#include "math_const.h"
#include "memory_kokkos.h"
#include "neighbor_kokkos.h"
#include "pair.h"
#include <cmath>
using namespace LAMMPS_NS;
using namespace MathConst;
#define TOLERANCE 0.05
/* ---------------------------------------------------------------------- */
template<class DeviceType>
DihedralCharmmfswKokkos<DeviceType>::DihedralCharmmfswKokkos(LAMMPS *lmp) : DihedralCharmmfsw(lmp)
{
atomKK = (AtomKokkos *) atom;
neighborKK = (NeighborKokkos *) neighbor;
execution_space = ExecutionSpaceFromDevice<DeviceType>::space;
datamask_read = X_MASK | F_MASK | Q_MASK | ENERGY_MASK | VIRIAL_MASK | TYPE_MASK;
datamask_modify = F_MASK | ENERGY_MASK | VIRIAL_MASK;
k_warning_flag = Kokkos::DualView<int,DeviceType>("Dihedral:warning_flag");
d_warning_flag = k_warning_flag.template view<DeviceType>();
h_warning_flag = k_warning_flag.h_view;
centroidstressflag = CENTROID_NOTAVAIL;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
DihedralCharmmfswKokkos<DeviceType>::~DihedralCharmmfswKokkos()
{
if (!copymode) {
memoryKK->destroy_kokkos(k_eatom,eatom);
memoryKK->destroy_kokkos(k_vatom,vatom);
}
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void DihedralCharmmfswKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
{
eflag = eflag_in;
vflag = vflag_in;
if (lmp->kokkos->neighflag == FULL)
error->all(FLERR,"Dihedral_style charmm/kk requires half neighbor list");
ev_init(eflag,vflag,0);
// ensure pair->ev_tally() will use 1-4 virial contribution
if (weightflag && vflag_global == VIRIAL_FDOTR)
force->pair->vflag_either = force->pair->vflag_global = 1;
// reallocate per-atom arrays if necessary
if (eflag_atom) {
//if(k_eatom.extent(0)<maxeatom) { // won't work without adding zero functor
memoryKK->destroy_kokkos(k_eatom,eatom);
memoryKK->create_kokkos(k_eatom,eatom,maxeatom,"dihedral:eatom");
d_eatom = k_eatom.template view<KKDeviceType>();
k_eatom_pair = Kokkos::DualView<E_FLOAT*,Kokkos::LayoutRight,KKDeviceType>("dihedral:eatom_pair",maxeatom);
d_eatom_pair = k_eatom_pair.template view<KKDeviceType>();
//}
}
if (vflag_atom) {
//if(k_vatom.extent(0)<maxvatom) { // won't work without adding zero functor
memoryKK->destroy_kokkos(k_vatom,vatom);
memoryKK->create_kokkos(k_vatom,vatom,maxvatom,"dihedral:vatom");
d_vatom = k_vatom.template view<KKDeviceType>();
k_vatom_pair = Kokkos::DualView<F_FLOAT*[6],Kokkos::LayoutRight,KKDeviceType>("dihedral:vatom_pair",maxvatom);
d_vatom_pair = k_vatom_pair.template view<KKDeviceType>();
//}
}
x = atomKK->k_x.view<DeviceType>();
f = atomKK->k_f.view<DeviceType>();
q = atomKK->k_q.view<DeviceType>();
atomtype = atomKK->k_type.view<DeviceType>();
neighborKK->k_dihedrallist.template sync<DeviceType>();
dihedrallist = neighborKK->k_dihedrallist.view<DeviceType>();
int ndihedrallist = neighborKK->ndihedrallist;
nlocal = atom->nlocal;
newton_bond = force->newton_bond;
qqrd2e = force->qqrd2e;
h_warning_flag() = 0;
k_warning_flag.template modify<LMPHostType>();
k_warning_flag.template sync<DeviceType>();
copymode = 1;
// loop over neighbors of my atoms
EVM_FLOAT evm;
if (evflag) {
if (newton_bond) {
Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagDihedralCharmmfswCompute<1,1> >(0,ndihedrallist),*this,evm);
} else {
Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagDihedralCharmmfswCompute<0,1> >(0,ndihedrallist),*this,evm);
}
} else {
if (newton_bond) {
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagDihedralCharmmfswCompute<1,0> >(0,ndihedrallist),*this);
} else {
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagDihedralCharmmfswCompute<0,0> >(0,ndihedrallist),*this);
}
}
// error check
k_warning_flag.template modify<DeviceType>();
k_warning_flag.template sync<LMPHostType>();
if (h_warning_flag())
error->warning(FLERR,"Dihedral problem");
if (eflag_global) {
energy += evm.emol;
force->pair->eng_vdwl += evm.evdwl;
force->pair->eng_coul += evm.ecoul;
}
if (vflag_global) {
virial[0] += evm.v[0];
virial[1] += evm.v[1];
virial[2] += evm.v[2];
virial[3] += evm.v[3];
virial[4] += evm.v[4];
virial[5] += evm.v[5];
force->pair->virial[0] += evm.vp[0];
force->pair->virial[1] += evm.vp[1];
force->pair->virial[2] += evm.vp[2];
force->pair->virial[3] += evm.vp[3];
force->pair->virial[4] += evm.vp[4];
force->pair->virial[5] += evm.vp[5];
}
// don't yet have dualviews for eatom and vatom in pair_kokkos,
// so need to manually copy these to pair style
int n = nlocal;
if (newton_bond) n += atom->nghost;
if (eflag_atom) {
k_eatom.template modify<DeviceType>();
k_eatom.template sync<LMPHostType>();
k_eatom_pair.template modify<DeviceType>();
k_eatom_pair.template sync<LMPHostType>();
for (int i = 0; i < n; i++)
force->pair->eatom[i] += k_eatom_pair.h_view(i);
}
if (vflag_atom) {
k_vatom.template modify<DeviceType>();
k_vatom.template sync<LMPHostType>();
k_vatom_pair.template modify<DeviceType>();
k_vatom_pair.template sync<LMPHostType>();
for (int i = 0; i < n; i++) {
force->pair->vatom[i][0] += k_vatom_pair.h_view(i,0);
force->pair->vatom[i][1] += k_vatom_pair.h_view(i,1);
force->pair->vatom[i][2] += k_vatom_pair.h_view(i,2);
force->pair->vatom[i][3] += k_vatom_pair.h_view(i,3);
force->pair->vatom[i][4] += k_vatom_pair.h_view(i,4);
force->pair->vatom[i][5] += k_vatom_pair.h_view(i,5);
}
}
copymode = 0;
}
template<class DeviceType>
template<int NEWTON_BOND, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void DihedralCharmmfswKokkos<DeviceType>::operator()(TagDihedralCharmmfswCompute<NEWTON_BOND,EVFLAG>, const int &n, EVM_FLOAT& evm) const {
// The f array is atomic
Kokkos::View<F_FLOAT*[3], typename DAT::t_f_array::array_layout,typename KKDevice<DeviceType>::value,Kokkos::MemoryTraits<Kokkos::Atomic|Kokkos::Unmanaged> > a_f = f;
const int i1 = dihedrallist(n,0);
const int i2 = dihedrallist(n,1);
const int i3 = dihedrallist(n,2);
const int i4 = dihedrallist(n,3);
const int type = dihedrallist(n,4);
// 1st bond
const F_FLOAT vb1x = x(i1,0) - x(i2,0);
const F_FLOAT vb1y = x(i1,1) - x(i2,1);
const F_FLOAT vb1z = x(i1,2) - x(i2,2);
// 2nd bond
const F_FLOAT vb2x = x(i3,0) - x(i2,0);
const F_FLOAT vb2y = x(i3,1) - x(i2,1);
const F_FLOAT vb2z = x(i3,2) - x(i2,2);
const F_FLOAT vb2xm = -vb2x;
const F_FLOAT vb2ym = -vb2y;
const F_FLOAT vb2zm = -vb2z;
// 3rd bond
const F_FLOAT vb3x = x(i4,0) - x(i3,0);
const F_FLOAT vb3y = x(i4,1) - x(i3,1);
const F_FLOAT vb3z = x(i4,2) - x(i3,2);
const F_FLOAT ax = vb1y*vb2zm - vb1z*vb2ym;
const F_FLOAT ay = vb1z*vb2xm - vb1x*vb2zm;
const F_FLOAT az = vb1x*vb2ym - vb1y*vb2xm;
const F_FLOAT bx = vb3y*vb2zm - vb3z*vb2ym;
const F_FLOAT by = vb3z*vb2xm - vb3x*vb2zm;
const F_FLOAT bz = vb3x*vb2ym - vb3y*vb2xm;
const F_FLOAT rasq = ax*ax + ay*ay + az*az;
const F_FLOAT rbsq = bx*bx + by*by + bz*bz;
const F_FLOAT rgsq = vb2xm*vb2xm + vb2ym*vb2ym + vb2zm*vb2zm;
const F_FLOAT rg = sqrt(rgsq);
F_FLOAT rginv,ra2inv,rb2inv;
rginv = ra2inv = rb2inv = 0.0;
if (rg > 0) rginv = 1.0/rg;
if (rasq > 0) ra2inv = 1.0/rasq;
if (rbsq > 0) rb2inv = 1.0/rbsq;
const F_FLOAT rabinv = sqrt(ra2inv*rb2inv);
F_FLOAT c = (ax*bx + ay*by + az*bz)*rabinv;
F_FLOAT s = rg*rabinv*(ax*vb3x + ay*vb3y + az*vb3z);
// error check
if ((c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) && !d_warning_flag())
d_warning_flag() = 1;
if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0;
const int m = d_multiplicity[type];
F_FLOAT p = 1.0;
F_FLOAT ddf1,df1;
ddf1 = df1 = 0.0;
for (int i = 0; i < m; i++) {
ddf1 = p*c - df1*s;
df1 = p*s + df1*c;
p = ddf1;
}
p = p*d_cos_shift[type] + df1*d_sin_shift[type];
df1 = df1*d_cos_shift[type] - ddf1*d_sin_shift[type];
df1 *= -m;
p += 1.0;
if (m == 0) {
p = 1.0 + d_cos_shift[type];
df1 = 0.0;
}
E_FLOAT edihedral = 0.0;
if (eflag) edihedral = d_k[type] * p;
const F_FLOAT fg = vb1x*vb2xm + vb1y*vb2ym + vb1z*vb2zm;
const F_FLOAT hg = vb3x*vb2xm + vb3y*vb2ym + vb3z*vb2zm;
const F_FLOAT fga = fg*ra2inv*rginv;
const F_FLOAT hgb = hg*rb2inv*rginv;
const F_FLOAT gaa = -ra2inv*rg;
const F_FLOAT gbb = rb2inv*rg;
const F_FLOAT dtfx = gaa*ax;
const F_FLOAT dtfy = gaa*ay;
const F_FLOAT dtfz = gaa*az;
const F_FLOAT dtgx = fga*ax - hgb*bx;
const F_FLOAT dtgy = fga*ay - hgb*by;
const F_FLOAT dtgz = fga*az - hgb*bz;
const F_FLOAT dthx = gbb*bx;
const F_FLOAT dthy = gbb*by;
const F_FLOAT dthz = gbb*bz;
const F_FLOAT df = -d_k[type] * df1;
const F_FLOAT sx2 = df*dtgx;
const F_FLOAT sy2 = df*dtgy;
const F_FLOAT sz2 = df*dtgz;
F_FLOAT f1[3],f2[3],f3[3],f4[3];
f1[0] = df*dtfx;
f1[1] = df*dtfy;
f1[2] = df*dtfz;
f2[0] = sx2 - f1[0];
f2[1] = sy2 - f1[1];
f2[2] = sz2 - f1[2];
f4[0] = df*dthx;
f4[1] = df*dthy;
f4[2] = df*dthz;
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) {
a_f(i1,0) += f1[0];
a_f(i1,1) += f1[1];
a_f(i1,2) += f1[2];
}
if (NEWTON_BOND || i2 < nlocal) {
a_f(i2,0) += f2[0];
a_f(i2,1) += f2[1];
a_f(i2,2) += f2[2];
}
if (NEWTON_BOND || i3 < nlocal) {
a_f(i3,0) += f3[0];
a_f(i3,1) += f3[1];
a_f(i3,2) += f3[2];
}
if (NEWTON_BOND || i4 < nlocal) {
a_f(i4,0) += f4[0];
a_f(i4,1) += f4[1];
a_f(i4,2) += f4[2];
}
if (EVFLAG)
ev_tally(evm,i1,i2,i3,i4,edihedral,f1,f3,f4,
vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z);
// 1-4 LJ and Coulomb interactions
// tally energy/virial in pair, using newton_bond as newton flag
if (d_weight[type] > 0.0) {
const int itype = atomtype[i1];
const int jtype = atomtype[i4];
const F_FLOAT delx = x(i1,0) - x(i4,0);
const F_FLOAT dely = x(i1,1) - x(i4,1);
const F_FLOAT delz = x(i1,2) - x(i4,2);
const F_FLOAT rsq = delx*delx + dely*dely + delz*delz;
const F_FLOAT r2inv = 1.0/rsq;
const F_FLOAT r6inv = r2inv*r2inv*r2inv;
F_FLOAT forcecoul;
if (implicit) forcecoul = qqrd2e * q[i1]*q[i4]*r2inv;
else forcecoul = qqrd2e * q[i1]*q[i4]*sqrt(r2inv);
const F_FLOAT forcelj = r6inv * (d_lj14_1(itype,jtype)*r6inv - d_lj14_2(itype,jtype));
const F_FLOAT fpair = d_weight[type] * (forcelj+forcecoul)*r2inv;
const F_FLOAT r = sqrt(rsq);
F_FLOAT ecoul = 0.0;
F_FLOAT evdwl = 0.0;
F_FLOAT evdwl14_12, evdwl14_6;
if (eflag) {
if (dihedflag)
ecoul = d_weight[type] * forcecoul;
else
ecoul = d_weight[type] * qqrd2e * q[i1] * q[i4] *
(sqrt(r2inv) + r * cut_coulinv14 * cut_coulinv14 - 2.0 * cut_coulinv14);
evdwl14_12 = r6inv * d_lj14_3(itype,jtype) * r6inv -
d_lj14_3(itype,jtype) * cut_lj_inner6inv * cut_lj6inv;
evdwl14_6 =
-d_lj14_4(itype,jtype) * r6inv + d_lj14_4(itype,jtype) * cut_lj_inner3inv * cut_lj3inv;
evdwl = evdwl14_12 + evdwl14_6;
evdwl *= d_weight[type];
}
if (newton_bond || i1 < nlocal) {
a_f(i1,0) += delx*fpair;
a_f(i1,1) += dely*fpair;
a_f(i1,2) += delz*fpair;
}
if (newton_bond || i4 < nlocal) {
a_f(i4,0) -= delx*fpair;
a_f(i4,1) -= dely*fpair;
a_f(i4,2) -= delz*fpair;
}
if (EVFLAG) ev_tally(evm,i1,i4,evdwl,ecoul,fpair,delx,dely,delz);
}
}
template<class DeviceType>
template<int NEWTON_BOND, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void DihedralCharmmfswKokkos<DeviceType>::operator()(TagDihedralCharmmfswCompute<NEWTON_BOND,EVFLAG>, const int &n) const {
EVM_FLOAT evm;
this->template operator()<NEWTON_BOND,EVFLAG>(TagDihedralCharmmfswCompute<NEWTON_BOND,EVFLAG>(), n, evm);
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void DihedralCharmmfswKokkos<DeviceType>::allocate()
{
DihedralCharmmfsw::allocate();
}
/* ----------------------------------------------------------------------
set coeffs for one or more types
------------------------------------------------------------------------- */
template<class DeviceType>
void DihedralCharmmfswKokkos<DeviceType>::coeff(int narg, char **arg)
{
DihedralCharmmfsw::coeff(narg, arg);
int nd = atom->ndihedraltypes;
typename AT::tdual_ffloat_1d k_k("DihedralCharmmfsw::k",nd+1);
typename AT::tdual_ffloat_1d k_multiplicity("DihedralCharmmfsw::multiplicity",nd+1);
typename AT::tdual_ffloat_1d k_shift("DihedralCharmmfsw::shift",nd+1);
typename AT::tdual_ffloat_1d k_cos_shift("DihedralCharmmfsw::cos_shift",nd+1);
typename AT::tdual_ffloat_1d k_sin_shift("DihedralCharmmfsw::sin_shift",nd+1);
typename AT::tdual_ffloat_1d k_weight("DihedralCharmmfsw::weight",nd+1);
d_k = k_k.template view<DeviceType>();
d_multiplicity = k_multiplicity.template view<DeviceType>();
d_shift = k_shift.template view<DeviceType>();
d_cos_shift = k_cos_shift.template view<DeviceType>();
d_sin_shift = k_sin_shift.template view<DeviceType>();
d_weight = k_weight.template view<DeviceType>();
int n = atom->ndihedraltypes;
for (int i = 1; i <= n; i++) {
k_k.h_view[i] = k[i];
k_multiplicity.h_view[i] = multiplicity[i];
k_shift.h_view[i] = shift[i];
k_cos_shift.h_view[i] = cos_shift[i];
k_sin_shift.h_view[i] = sin_shift[i];
k_weight.h_view[i] = weight[i];
}
k_k.template modify<LMPHostType>();
k_multiplicity.template modify<LMPHostType>();
k_shift.template modify<LMPHostType>();
k_cos_shift.template modify<LMPHostType>();
k_sin_shift.template modify<LMPHostType>();
k_weight.template modify<LMPHostType>();
k_k.template sync<DeviceType>();
k_multiplicity.template sync<DeviceType>();
k_shift.template sync<DeviceType>();
k_cos_shift.template sync<DeviceType>();
k_sin_shift.template sync<DeviceType>();
k_weight.template sync<DeviceType>();
}
/* ----------------------------------------------------------------------
error check and initialize all values needed for force computation
------------------------------------------------------------------------- */
template<class DeviceType>
void DihedralCharmmfswKokkos<DeviceType>::init_style()
{
DihedralCharmmfsw::init_style();
int n = atom->ntypes;
DAT::tdual_ffloat_2d k_lj14_1("DihedralCharmmfsw:lj14_1",n+1,n+1);
DAT::tdual_ffloat_2d k_lj14_2("DihedralCharmmfsw:lj14_2",n+1,n+1);
DAT::tdual_ffloat_2d k_lj14_3("DihedralCharmmfsw:lj14_3",n+1,n+1);
DAT::tdual_ffloat_2d k_lj14_4("DihedralCharmmfsw:lj14_4",n+1,n+1);
d_lj14_1 = k_lj14_1.template view<DeviceType>();
d_lj14_2 = k_lj14_2.template view<DeviceType>();
d_lj14_3 = k_lj14_3.template view<DeviceType>();
d_lj14_4 = k_lj14_4.template view<DeviceType>();
if (weightflag) {
int n = atom->ntypes;
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= n; j++) {
k_lj14_1.h_view(i,j) = lj14_1[i][j];
k_lj14_2.h_view(i,j) = lj14_2[i][j];
k_lj14_3.h_view(i,j) = lj14_3[i][j];
k_lj14_4.h_view(i,j) = lj14_4[i][j];
}
}
}
k_lj14_1.template modify<LMPHostType>();
k_lj14_2.template modify<LMPHostType>();
k_lj14_3.template modify<LMPHostType>();
k_lj14_4.template modify<LMPHostType>();
k_lj14_1.template sync<DeviceType>();
k_lj14_2.template sync<DeviceType>();
k_lj14_3.template sync<DeviceType>();
k_lj14_4.template sync<DeviceType>();
}
/* ----------------------------------------------------------------------
proc 0 reads coeffs from restart file, bcasts them
------------------------------------------------------------------------- */
template<class DeviceType>
void DihedralCharmmfswKokkos<DeviceType>::read_restart(FILE *fp)
{
DihedralCharmmfsw::read_restart(fp);
int nd = atom->ndihedraltypes;
typename AT::tdual_ffloat_1d k_k("DihedralCharmmfsw::k",nd+1);
typename AT::tdual_ffloat_1d k_multiplicity("DihedralCharmmfsw::multiplicity",nd+1);
typename AT::tdual_ffloat_1d k_shift("DihedralCharmmfsw::shift",nd+1);
typename AT::tdual_ffloat_1d k_cos_shift("DihedralCharmmfsw::cos_shift",nd+1);
typename AT::tdual_ffloat_1d k_sin_shift("DihedralCharmmfsw::sin_shift",nd+1);
typename AT::tdual_ffloat_1d k_weight("DihedralCharmmfsw::weight",nd+1);
d_k = k_k.template view<DeviceType>();
d_multiplicity = k_multiplicity.template view<DeviceType>();
d_shift = k_shift.template view<DeviceType>();
d_cos_shift = k_cos_shift.template view<DeviceType>();
d_sin_shift = k_sin_shift.template view<DeviceType>();
d_weight = k_weight.template view<DeviceType>();
int n = atom->ndihedraltypes;
for (int i = 1; i <= n; i++) {
k_k.h_view[i] = k[i];
k_multiplicity.h_view[i] = multiplicity[i];
k_shift.h_view[i] = shift[i];
k_cos_shift.h_view[i] = cos_shift[i];
k_sin_shift.h_view[i] = sin_shift[i];
k_weight.h_view[i] = weight[i];
}
k_k.template modify<LMPHostType>();
k_multiplicity.template modify<LMPHostType>();
k_shift.template modify<LMPHostType>();
k_cos_shift.template modify<LMPHostType>();
k_sin_shift.template modify<LMPHostType>();
k_weight.template modify<LMPHostType>();
k_k.template sync<DeviceType>();
k_multiplicity.template sync<DeviceType>();
k_shift.template sync<DeviceType>();
k_cos_shift.template sync<DeviceType>();
k_sin_shift.template sync<DeviceType>();
k_weight.template sync<DeviceType>();
}
/* ----------------------------------------------------------------------
tally energy and virial into global and per-atom accumulators
virial = r1F1 + r2F2 + r3F3 + r4F4 = (r1-r2) F1 + (r3-r2) F3 + (r4-r2) F4
= (r1-r2) F1 + (r3-r2) F3 + (r4-r3 + r3-r2) F4
= vb1*f1 + vb2*f3 + (vb3+vb2)*f4
------------------------------------------------------------------------- */
template<class DeviceType>
//template<int NEWTON_BOND>
KOKKOS_INLINE_FUNCTION
void DihedralCharmmfswKokkos<DeviceType>::ev_tally(EVM_FLOAT &evm, const int i1, const int i2, const int i3, const int i4,
F_FLOAT &edihedral, F_FLOAT *f1, F_FLOAT *f3, F_FLOAT *f4,
const F_FLOAT &vb1x, const F_FLOAT &vb1y, const F_FLOAT &vb1z,
const F_FLOAT &vb2x, const F_FLOAT &vb2y, const F_FLOAT &vb2z,
const F_FLOAT &vb3x, const F_FLOAT &vb3y, const F_FLOAT &vb3z) const
{
E_FLOAT edihedralquarter;
F_FLOAT v[6];
if (eflag_either) {
if (eflag_global) {
if (newton_bond) evm.emol += edihedral;
else {
edihedralquarter = 0.25*edihedral;
if (i1 < nlocal) evm.emol += edihedralquarter;
if (i2 < nlocal) evm.emol += edihedralquarter;
if (i3 < nlocal) evm.emol += edihedralquarter;
if (i4 < nlocal) evm.emol += edihedralquarter;
}
}
if (eflag_atom) {
edihedralquarter = 0.25*edihedral;
if (newton_bond || i1 < nlocal) d_eatom[i1] += edihedralquarter;
if (newton_bond || i2 < nlocal) d_eatom[i2] += edihedralquarter;
if (newton_bond || i3 < nlocal) d_eatom[i3] += edihedralquarter;
if (newton_bond || i4 < nlocal) d_eatom[i4] += edihedralquarter;
}
}
if (vflag_either) {
v[0] = vb1x*f1[0] + vb2x*f3[0] + (vb3x+vb2x)*f4[0];
v[1] = vb1y*f1[1] + vb2y*f3[1] + (vb3y+vb2y)*f4[1];
v[2] = vb1z*f1[2] + vb2z*f3[2] + (vb3z+vb2z)*f4[2];
v[3] = vb1x*f1[1] + vb2x*f3[1] + (vb3x+vb2x)*f4[1];
v[4] = vb1x*f1[2] + vb2x*f3[2] + (vb3x+vb2x)*f4[2];
v[5] = vb1y*f1[2] + vb2y*f3[2] + (vb3y+vb2y)*f4[2];
if (vflag_global) {
if (newton_bond) {
evm.v[0] += v[0];
evm.v[1] += v[1];
evm.v[2] += v[2];
evm.v[3] += v[3];
evm.v[4] += v[4];
evm.v[5] += v[5];
} else {
if (i1 < nlocal) {
evm.v[0] += 0.25*v[0];
evm.v[1] += 0.25*v[1];
evm.v[2] += 0.25*v[2];
evm.v[3] += 0.25*v[3];
evm.v[4] += 0.25*v[4];
evm.v[5] += 0.25*v[5];
}
if (i2 < nlocal) {
evm.v[0] += 0.25*v[0];
evm.v[1] += 0.25*v[1];
evm.v[2] += 0.25*v[2];
evm.v[3] += 0.25*v[3];
evm.v[4] += 0.25*v[4];
evm.v[5] += 0.25*v[5];
}
if (i3 < nlocal) {
evm.v[0] += 0.25*v[0];
evm.v[1] += 0.25*v[1];
evm.v[2] += 0.25*v[2];
evm.v[3] += 0.25*v[3];
evm.v[4] += 0.25*v[4];
evm.v[5] += 0.25*v[5];
}
if (i4 < nlocal) {
evm.v[0] += 0.25*v[0];
evm.v[1] += 0.25*v[1];
evm.v[2] += 0.25*v[2];
evm.v[3] += 0.25*v[3];
evm.v[4] += 0.25*v[4];
evm.v[5] += 0.25*v[5];
}
}
}
if (vflag_atom) {
if (newton_bond || i1 < nlocal) {
d_vatom(i1,0) += 0.25*v[0];
d_vatom(i1,1) += 0.25*v[1];
d_vatom(i1,2) += 0.25*v[2];
d_vatom(i1,3) += 0.25*v[3];
d_vatom(i1,4) += 0.25*v[4];
d_vatom(i1,5) += 0.25*v[5];
}
if (newton_bond || i2 < nlocal) {
d_vatom(i2,0) += 0.25*v[0];
d_vatom(i2,1) += 0.25*v[1];
d_vatom(i2,2) += 0.25*v[2];
d_vatom(i2,3) += 0.25*v[3];
d_vatom(i2,4) += 0.25*v[4];
d_vatom(i2,5) += 0.25*v[5];
}
if (newton_bond || i3 < nlocal) {
d_vatom(i3,0) += 0.25*v[0];
d_vatom(i3,1) += 0.25*v[1];
d_vatom(i3,2) += 0.25*v[2];
d_vatom(i3,3) += 0.25*v[3];
d_vatom(i3,4) += 0.25*v[4];
d_vatom(i3,5) += 0.25*v[5];
}
if (newton_bond || i4 < nlocal) {
d_vatom(i4,0) += 0.25*v[0];
d_vatom(i4,1) += 0.25*v[1];
d_vatom(i4,2) += 0.25*v[2];
d_vatom(i4,3) += 0.25*v[3];
d_vatom(i4,4) += 0.25*v[4];
d_vatom(i4,5) += 0.25*v[5];
}
}
}
}
/* ----------------------------------------------------------------------
tally eng_vdwl and virial into global and per-atom accumulators
need i < nlocal test since called by bond_quartic and dihedral_charmm
------------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void DihedralCharmmfswKokkos<DeviceType>::ev_tally(EVM_FLOAT &evm, const int i, const int j,
const F_FLOAT &evdwl, const F_FLOAT &ecoul, const F_FLOAT &fpair, const F_FLOAT &delx,
const F_FLOAT &dely, const F_FLOAT &delz) const
{
E_FLOAT evdwlhalf,ecoulhalf,epairhalf;
F_FLOAT v[6];
if (eflag_either) {
if (eflag_global) {
if (newton_bond) {
evm.evdwl += evdwl;
evm.ecoul += ecoul;
} else {
evdwlhalf = 0.5*evdwl;
ecoulhalf = 0.5*ecoul;
if (i < nlocal) {
evm.evdwl += evdwlhalf;
evm.ecoul += ecoulhalf;
}
if (j < nlocal) {
evm.evdwl += evdwlhalf;
evm.ecoul += ecoulhalf;
}
}
}
if (eflag_atom) {
epairhalf = 0.5 * (evdwl + ecoul);
if (newton_bond || i < nlocal) d_eatom_pair[i] += epairhalf;
if (newton_bond || j < nlocal) d_eatom_pair[j] += epairhalf;
}
}
if (vflag_either) {
v[0] = delx*delx*fpair;
v[1] = dely*dely*fpair;
v[2] = delz*delz*fpair;
v[3] = delx*dely*fpair;
v[4] = delx*delz*fpair;
v[5] = dely*delz*fpair;
if (vflag_global) {
if (newton_bond) {
evm.vp[0] += v[0];
evm.vp[1] += v[1];
evm.vp[2] += v[2];
evm.vp[3] += v[3];
evm.vp[4] += v[4];
evm.vp[5] += v[5];
} else {
if (i < nlocal) {
evm.vp[0] += 0.5*v[0];
evm.vp[1] += 0.5*v[1];
evm.vp[2] += 0.5*v[2];
evm.vp[3] += 0.5*v[3];
evm.vp[4] += 0.5*v[4];
evm.vp[5] += 0.5*v[5];
}
if (j < nlocal) {
evm.vp[0] += 0.5*v[0];
evm.vp[1] += 0.5*v[1];
evm.vp[2] += 0.5*v[2];
evm.vp[3] += 0.5*v[3];
evm.vp[4] += 0.5*v[4];
evm.vp[5] += 0.5*v[5];
}
}
}
if (vflag_atom) {
if (newton_bond || i < nlocal) {
d_vatom_pair(i,0) += 0.5*v[0];
d_vatom_pair(i,1) += 0.5*v[1];
d_vatom_pair(i,2) += 0.5*v[2];
d_vatom_pair(i,3) += 0.5*v[3];
d_vatom_pair(i,4) += 0.5*v[4];
d_vatom_pair(i,5) += 0.5*v[5];
}
if (newton_bond || j < nlocal) {
d_vatom_pair(j,0) += 0.5*v[0];
d_vatom_pair(j,1) += 0.5*v[1];
d_vatom_pair(j,2) += 0.5*v[2];
d_vatom_pair(j,3) += 0.5*v[3];
d_vatom_pair(j,4) += 0.5*v[4];
d_vatom_pair(j,5) += 0.5*v[5];
}
}
}
}
/* ---------------------------------------------------------------------- */
namespace LAMMPS_NS {
template class DihedralCharmmfswKokkos<LMPDeviceType>;
#ifdef LMP_KOKKOS_GPU
template class DihedralCharmmfswKokkos<LMPHostType>;
#endif
}

118
src/KOKKOS/dihedral_charmmfsw_kokkos.h Normal file
View File

@ -0,0 +1,118 @@
/* -*- 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(charmmfsw/kk,DihedralCharmmfswKokkos<LMPDeviceType>);
DihedralStyle(charmmfsw/kk/device,DihedralCharmmfswKokkos<LMPDeviceType>);
DihedralStyle(charmmfsw/kk/host,DihedralCharmmfswKokkos<LMPHostType>);
// clang-format on
#else
// clang-format off
#ifndef LMP_DIHEDRAL_CHARMMFSW_KOKKOS_H
#define LMP_DIHEDRAL_CHARMMFSW_KOKKOS_H
#include "dihedral_charmmfsw.h"
#include "kokkos_type.h"
#include "dihedral_charmm_kokkos.h" // needed for s_EVM_FLOAT
namespace LAMMPS_NS {
template<int NEWTON_BOND, int EVFLAG>
struct TagDihedralCharmmfswCompute{};
template<class DeviceType>
class DihedralCharmmfswKokkos : public DihedralCharmmfsw {
public:
typedef DeviceType device_type;
typedef EVM_FLOAT value_type;
typedef ArrayTypes<DeviceType> AT;
DihedralCharmmfswKokkos(class LAMMPS *);
~DihedralCharmmfswKokkos() override;
void compute(int, int) override;
void coeff(int, char **) override;
void init_style() override;
void read_restart(FILE *) override;
template<int NEWTON_BOND, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void operator()(TagDihedralCharmmfswCompute<NEWTON_BOND,EVFLAG>, const int&, EVM_FLOAT&) const;
template<int NEWTON_BOND, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void operator()(TagDihedralCharmmfswCompute<NEWTON_BOND,EVFLAG>, const int&) const;
//template<int NEWTON_BOND>
KOKKOS_INLINE_FUNCTION
void ev_tally(EVM_FLOAT &evm, const int i1, const int i2, const int i3, const int i4,
F_FLOAT &edihedral, F_FLOAT *f1, F_FLOAT *f3, F_FLOAT *f4,
const F_FLOAT &vb1x, const F_FLOAT &vb1y, const F_FLOAT &vb1z,
const F_FLOAT &vb2x, const F_FLOAT &vb2y, const F_FLOAT &vb2z,
const F_FLOAT &vb3x, const F_FLOAT &vb3y, const F_FLOAT &vb3z) const;
KOKKOS_INLINE_FUNCTION
void ev_tally(EVM_FLOAT &evm, const int i, const int j,
const F_FLOAT &evdwl, const F_FLOAT &ecoul, const F_FLOAT &fpair, const F_FLOAT &delx,
const F_FLOAT &dely, const F_FLOAT &delz) const;
protected:
class NeighborKokkos *neighborKK;
typename AT::t_x_array_randomread x;
typename AT::t_int_1d_randomread atomtype;
typename AT::t_ffloat_1d_randomread q;
typename AT::t_f_array f;
typename AT::t_int_2d dihedrallist;
typedef typename KKDevice<DeviceType>::value KKDeviceType;
Kokkos::DualView<E_FLOAT*,Kokkos::LayoutRight,KKDeviceType> k_eatom;
Kokkos::DualView<F_FLOAT*[6],Kokkos::LayoutRight,KKDeviceType> k_vatom;
Kokkos::View<E_FLOAT*,Kokkos::LayoutRight,KKDeviceType,Kokkos::MemoryTraits<Kokkos::Atomic> > d_eatom;
Kokkos::View<F_FLOAT*[6],Kokkos::LayoutRight,KKDeviceType,Kokkos::MemoryTraits<Kokkos::Atomic> > d_vatom;
Kokkos::DualView<E_FLOAT*,Kokkos::LayoutRight,KKDeviceType> k_eatom_pair;
Kokkos::DualView<F_FLOAT*[6],Kokkos::LayoutRight,KKDeviceType> k_vatom_pair;
Kokkos::View<E_FLOAT*,Kokkos::LayoutRight,KKDeviceType,Kokkos::MemoryTraits<Kokkos::Atomic> > d_eatom_pair;
Kokkos::View<F_FLOAT*[6],Kokkos::LayoutRight,KKDeviceType,Kokkos::MemoryTraits<Kokkos::Atomic> > d_vatom_pair;
int nlocal,newton_bond;
int eflag,vflag;
double qqrd2e;
Kokkos::DualView<int,DeviceType> k_warning_flag;
typename Kokkos::DualView<int,DeviceType>::t_dev d_warning_flag;
typename Kokkos::DualView<int,DeviceType>::t_host h_warning_flag;
typename AT::t_ffloat_2d d_lj14_1;
typename AT::t_ffloat_2d d_lj14_2;
typename AT::t_ffloat_2d d_lj14_3;
typename AT::t_ffloat_2d d_lj14_4;
typename AT::t_ffloat_1d d_k;
typename AT::t_ffloat_1d d_multiplicity;
typename AT::t_ffloat_1d d_shift;
typename AT::t_ffloat_1d d_sin_shift;
typename AT::t_ffloat_1d d_cos_shift;
typename AT::t_ffloat_1d d_weight;
void allocate() override;
};
}
#endif
#endif

10
src/KOKKOS/pair_kokkos.h
View File

@ -627,7 +627,7 @@ struct PairComputeFunctor {
const int itype = c.type(i); const int itype = c.type(i);
const F_FLOAT qtmp = c.q(i); const F_FLOAT qtmp = c.q(i);
if (ZEROFLAG) { if (NEIGHFLAG == FULL && ZEROFLAG) {
Kokkos::single(Kokkos::PerThread(team), [&] (){ Kokkos::single(Kokkos::PerThread(team), [&] (){
f(i,0) = 0.0; f(i,0) = 0.0;
f(i,1) = 0.0; f(i,1) = 0.0;
@ -674,7 +674,7 @@ struct PairComputeFunctor {
const int J_CONTRIB = ((NEIGHFLAG == HALF || NEIGHFLAG == HALFTHREAD) && j < c.nlocal); const int J_CONTRIB = ((NEIGHFLAG == HALF || NEIGHFLAG == HALFTHREAD) && j < c.nlocal);
const E_FLOAT factor = J_CONTRIB?1.0:0.5; const E_FLOAT factor = J_CONTRIB?1.0:0.5;
if ((NEIGHFLAG == HALF || NEIGHFLAG == HALFTHREAD) && j < c.nlocal) { if (J_CONTRIB) {
a_f(j,0) -= fx; a_f(j,0) -= fx;
a_f(j,1) -= fy; a_f(j,1) -= fy;
a_f(j,2) -= fz; a_f(j,2) -= fz;
@ -746,8 +746,10 @@ struct PairComputeFunctor {
a_f(i,1) += fev.f[1]; a_f(i,1) += fev.f[1];
a_f(i,2) += fev.f[2]; a_f(i,2) += fev.f[2];
if (c.eflag_global) if (c.eflag_global) {
ev.evdwl += fev.evdwl; ev.evdwl += fev.evdwl;
ev.ecoul += fev.ecoul;
}
if (c.vflag_global) { if (c.vflag_global) {
ev.v[0] += fev.v[0]; ev.v[0] += fev.v[0];
@ -761,7 +763,7 @@ struct PairComputeFunctor {
if (NEIGHFLAG == FULL) { if (NEIGHFLAG == FULL) {
if (c.eflag_atom) if (c.eflag_atom)
a_eatom(i) += fev.evdwl; a_eatom(i) += fev.evdwl + fev.ecoul;
if (c.vflag_atom) { if (c.vflag_atom) {
a_vatom(i,0) += fev.v[0]; a_vatom(i,0) += fev.v[0];

3
src/KOKKOS/pair_lj_charmm_coul_long_kokkos.cpp
View File

@ -214,9 +214,7 @@ compute_evdwl(const F_FLOAT& rsq, const int& /*i*/, const int& /*j*/,
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj; (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
englj *= switch1; englj *= switch1;
} }
return englj; return englj;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -488,4 +486,3 @@ template class PairLJCharmmCoulLongKokkos<LMPDeviceType>;
template class PairLJCharmmCoulLongKokkos<LMPHostType>; template class PairLJCharmmCoulLongKokkos<LMPHostType>;
#endif #endif
} }

497
src/KOKKOS/pair_lj_charmmfsw_coul_long_kokkos.cpp Normal file
View File

@ -0,0 +1,497 @@
// 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: Mitch Murphy (alphataubio)
Based on serial kspace lj-fsw sections (force-switched) provided by
Robert Meissner and Lucio Colombi Ciacchi of Bremen University, Germany,
with additional assistance from Robert A. Latour, Clemson University
------------------------------------------------------------------------- */
#include "pair_lj_charmmfsw_coul_long_kokkos.h"
#include "atom_kokkos.h"
#include "atom_masks.h"
#include "error.h"
#include "force.h"
#include "kokkos.h"
#include "memory_kokkos.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "neighbor.h"
#include "respa.h"
#include "update.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
#define EWALD_F 1.12837917
#define EWALD_P 0.3275911
#define A1 0.254829592
#define A2 -0.284496736
#define A3 1.421413741
#define A4 -1.453152027
#define A5 1.061405429
/* ---------------------------------------------------------------------- */
template<class DeviceType>
PairLJCharmmfswCoulLongKokkos<DeviceType>::PairLJCharmmfswCoulLongKokkos(LAMMPS *lmp):PairLJCharmmfswCoulLong(lmp)
{
respa_enable = 0;
kokkosable = 1;
atomKK = (AtomKokkos *) atom;
execution_space = ExecutionSpaceFromDevice<DeviceType>::space;
datamask_read = X_MASK | F_MASK | TYPE_MASK | Q_MASK | ENERGY_MASK | VIRIAL_MASK;
datamask_modify = F_MASK | ENERGY_MASK | VIRIAL_MASK;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
PairLJCharmmfswCoulLongKokkos<DeviceType>::~PairLJCharmmfswCoulLongKokkos()
{
if (copymode) return;
if (allocated) {
memoryKK->destroy_kokkos(k_eatom,eatom);
memoryKK->destroy_kokkos(k_vatom,vatom);
memoryKK->destroy_kokkos(k_cutsq,cutsq);
}
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void PairLJCharmmfswCoulLongKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
{
eflag = eflag_in;
vflag = vflag_in;
if (neighflag == FULL) no_virial_fdotr_compute = 1;
ev_init(eflag,vflag,0);
// reallocate per-atom arrays if necessary
if (eflag_atom) {
memoryKK->destroy_kokkos(k_eatom,eatom);
memoryKK->create_kokkos(k_eatom,eatom,maxeatom,"pair:eatom");
d_eatom = k_eatom.view<DeviceType>();
}
if (vflag_atom) {
memoryKK->destroy_kokkos(k_vatom,vatom);
memoryKK->create_kokkos(k_vatom,vatom,maxvatom,"pair:vatom");
d_vatom = k_vatom.view<DeviceType>();
}
atomKK->sync(execution_space,datamask_read);
k_cutsq.template sync<DeviceType>();
k_params.template sync<DeviceType>();
if (eflag || vflag) atomKK->modified(execution_space,datamask_modify);
else atomKK->modified(execution_space,F_MASK);
x = atomKK->k_x.view<DeviceType>();
c_x = atomKK->k_x.view<DeviceType>();
f = atomKK->k_f.view<DeviceType>();
q = atomKK->k_q.view<DeviceType>();
type = atomKK->k_type.view<DeviceType>();
nlocal = atom->nlocal;
nall = atom->nlocal + atom->nghost;
special_lj[0] = force->special_lj[0];
special_lj[1] = force->special_lj[1];
special_lj[2] = force->special_lj[2];
special_lj[3] = force->special_lj[3];
special_coul[0] = force->special_coul[0];
special_coul[1] = force->special_coul[1];
special_coul[2] = force->special_coul[2];
special_coul[3] = force->special_coul[3];
qqrd2e = force->qqrd2e;
newton_pair = force->newton_pair;
// loop over neighbors of my atoms
copymode = 1;
EV_FLOAT ev;
if (ncoultablebits)
ev = pair_compute<PairLJCharmmfswCoulLongKokkos<DeviceType>,CoulLongTable<1> >
(this,(NeighListKokkos<DeviceType>*)list);
else
ev = pair_compute<PairLJCharmmfswCoulLongKokkos<DeviceType>,CoulLongTable<0> >
(this,(NeighListKokkos<DeviceType>*)list);
if (eflag) {
eng_vdwl += ev.evdwl;
eng_coul += ev.ecoul;
}
if (vflag_global) {
virial[0] += ev.v[0];
virial[1] += ev.v[1];
virial[2] += ev.v[2];
virial[3] += ev.v[3];
virial[4] += ev.v[4];
virial[5] += ev.v[5];
}
if (eflag_atom) {
k_eatom.template modify<DeviceType>();
k_eatom.template sync<LMPHostType>();
}
if (vflag_atom) {
k_vatom.template modify<DeviceType>();
k_vatom.template sync<LMPHostType>();
}
if (vflag_fdotr) pair_virial_fdotr_compute(this);
copymode = 0;
}
/* ----------------------------------------------------------------------
compute LJ CHARMM pair force between atoms i and j
---------------------------------------------------------------------- */
template<class DeviceType>
template<bool STACKPARAMS, class Specialisation>
KOKKOS_INLINE_FUNCTION
F_FLOAT PairLJCharmmfswCoulLongKokkos<DeviceType>::
compute_fpair(const F_FLOAT& rsq, const int& /*i*/, const int& /*j*/,
const int& itype, const int& jtype) const {
const F_FLOAT r2inv = 1.0/rsq;
const F_FLOAT r6inv = r2inv*r2inv*r2inv;
F_FLOAT forcelj, switch1;
forcelj = r6inv *
((STACKPARAMS?m_params[itype][jtype].lj1:params(itype,jtype).lj1)*r6inv -
(STACKPARAMS?m_params[itype][jtype].lj2:params(itype,jtype).lj2));
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
forcelj = forcelj*switch1;
}
return forcelj*r2inv;
}
/* ----------------------------------------------------------------------
compute LJ CHARMM pair potential energy between atoms i and j
---------------------------------------------------------------------- */
template<class DeviceType>
template<bool STACKPARAMS, class Specialisation>
KOKKOS_INLINE_FUNCTION
F_FLOAT PairLJCharmmfswCoulLongKokkos<DeviceType>::
compute_evdwl(const F_FLOAT& rsq, const int& /*i*/, const int& /*j*/,
const int& itype, const int& jtype) const {
const F_FLOAT r2inv = 1.0/rsq;
const F_FLOAT r6inv = r2inv*r2inv*r2inv;
const F_FLOAT r = sqrt(rsq);
const F_FLOAT rinv = 1.0/r;
const F_FLOAT r3inv = rinv*rinv*rinv;
F_FLOAT englj, englj12, englj6;
if (rsq > cut_lj_innersq) {
englj12 = (STACKPARAMS?m_params[itype][jtype].lj3:params(itype,jtype).lj3)*cut_lj6*
denom_lj12 * (r6inv - cut_lj6inv)*(r6inv - cut_lj6inv);
englj6 = -(STACKPARAMS?m_params[itype][jtype].lj4:params(itype,jtype).lj4)*
cut_lj3*denom_lj6 * (r3inv - cut_lj3inv)*(r3inv - cut_lj3inv);
englj = englj12 + englj6;
} else {
englj12 = r6inv*(STACKPARAMS?m_params[itype][jtype].lj3:params(itype,jtype).lj3)*r6inv -
(STACKPARAMS?m_params[itype][jtype].lj3:params(itype,jtype).lj3)*cut_lj_inner6inv*cut_lj6inv;
englj6 = -(STACKPARAMS?m_params[itype][jtype].lj4:params(itype,jtype).lj4)*r6inv +
(STACKPARAMS?m_params[itype][jtype].lj4:params(itype,jtype).lj4)*
cut_lj_inner3inv*cut_lj3inv;
englj = englj12 + englj6;
}
return englj;
}
/* ----------------------------------------------------------------------
compute coulomb pair force between atoms i and j
---------------------------------------------------------------------- */
template<class DeviceType>
template<bool STACKPARAMS, class Specialisation>
KOKKOS_INLINE_FUNCTION
F_FLOAT PairLJCharmmfswCoulLongKokkos<DeviceType>::
compute_fcoul(const F_FLOAT& rsq, const int& /*i*/, const int&j,
const int& /*itype*/, const int& /*jtype*/,
const F_FLOAT& factor_coul, const F_FLOAT& qtmp) const {
if (Specialisation::DoTable && rsq > tabinnersq) {
union_int_float_t rsq_lookup;
rsq_lookup.f = rsq;
const int itable = (rsq_lookup.i & ncoulmask) >> ncoulshiftbits;
const F_FLOAT fraction = (rsq_lookup.f - d_rtable[itable]) * d_drtable[itable];
const F_FLOAT table = d_ftable[itable] + fraction*d_dftable[itable];
F_FLOAT forcecoul = qtmp*q[j] * table;
if (factor_coul < 1.0) {
const F_FLOAT table = d_ctable[itable] + fraction*d_dctable[itable];
const F_FLOAT prefactor = qtmp*q[j] * table;
forcecoul -= (1.0-factor_coul)*prefactor;
}
return forcecoul/rsq;
} else {
const F_FLOAT r = sqrt(rsq);
const F_FLOAT grij = g_ewald * r;
const F_FLOAT expm2 = exp(-grij*grij);
const F_FLOAT t = 1.0 / (1.0 + EWALD_P*grij);
const F_FLOAT rinv = 1.0/r;
const F_FLOAT erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
const F_FLOAT prefactor = qqrd2e * qtmp*q[j]*rinv;
F_FLOAT forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
return forcecoul*rinv*rinv;
}
}
/* ----------------------------------------------------------------------
compute coulomb pair potential energy between atoms i and j
---------------------------------------------------------------------- */
template<class DeviceType>
template<bool STACKPARAMS, class Specialisation>
KOKKOS_INLINE_FUNCTION
F_FLOAT PairLJCharmmfswCoulLongKokkos<DeviceType>::
compute_ecoul(const F_FLOAT& rsq, const int& /*i*/, const int&j,
const int& /*itype*/, const int& /*jtype*/, const F_FLOAT& factor_coul, const F_FLOAT& qtmp) const {
if (Specialisation::DoTable && rsq > tabinnersq) {
union_int_float_t rsq_lookup;
rsq_lookup.f = rsq;
const int itable = (rsq_lookup.i & ncoulmask) >> ncoulshiftbits;
const F_FLOAT fraction = (rsq_lookup.f - d_rtable[itable]) * d_drtable[itable];
const F_FLOAT table = d_etable[itable] + fraction*d_detable[itable];
F_FLOAT ecoul = qtmp*q[j] * table;
if (factor_coul < 1.0) {
const F_FLOAT table = d_ctable[itable] + fraction*d_dctable[itable];
const F_FLOAT prefactor = qtmp*q[j] * table;
ecoul -= (1.0-factor_coul)*prefactor;
}
return ecoul;
} else {
const F_FLOAT r = sqrt(rsq);
const F_FLOAT grij = g_ewald * r;
const F_FLOAT expm2 = exp(-grij*grij);
const F_FLOAT t = 1.0 / (1.0 + EWALD_P*grij);
const F_FLOAT erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
const F_FLOAT prefactor = qqrd2e * qtmp*q[j]/r;
F_FLOAT ecoul = prefactor * erfc;
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
return ecoul;
}
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
template<class DeviceType>
void PairLJCharmmfswCoulLongKokkos<DeviceType>::allocate()
{
PairLJCharmmfswCoulLong::allocate();
int n = atom->ntypes;
memory->destroy(cutsq);
memoryKK->create_kokkos(k_cutsq,cutsq,n+1,n+1,"pair:cutsq");
d_cutsq = k_cutsq.template view<DeviceType>();
d_cut_ljsq = typename AT::t_ffloat_2d("pair:cut_ljsq",n+1,n+1);
d_cut_coulsq = typename AT::t_ffloat_2d("pair:cut_coulsq",n+1,n+1);
k_params = Kokkos::DualView<params_lj_coul**,Kokkos::LayoutRight,DeviceType>("PairLJCharmmfswCoulLong::params",n+1,n+1);
params = k_params.template view<DeviceType>();
}
template<class DeviceType>
void PairLJCharmmfswCoulLongKokkos<DeviceType>::init_tables(double cut_coul, double *cut_respa)
{
Pair::init_tables(cut_coul,cut_respa);
typedef typename ArrayTypes<DeviceType>::t_ffloat_1d table_type;
typedef typename ArrayTypes<LMPHostType>::t_ffloat_1d host_table_type;
int ntable = 1;
for (int i = 0; i < ncoultablebits; i++) ntable *= 2;
// Copy rtable and drtable
{
host_table_type h_table("HostTable",ntable);
table_type d_table("DeviceTable",ntable);
for (int i = 0; i < ntable; i++) {
h_table(i) = rtable[i];
}
Kokkos::deep_copy(d_table,h_table);
d_rtable = d_table;
}
{
host_table_type h_table("HostTable",ntable);
table_type d_table("DeviceTable",ntable);
for (int i = 0; i < ntable; i++) {
h_table(i) = drtable[i];
}
Kokkos::deep_copy(d_table,h_table);
d_drtable = d_table;
}
{
host_table_type h_table("HostTable",ntable);
table_type d_table("DeviceTable",ntable);
// Copy ftable and dftable
for (int i = 0; i < ntable; i++) {
h_table(i) = ftable[i];
}
Kokkos::deep_copy(d_table,h_table);
d_ftable = d_table;
}
{
host_table_type h_table("HostTable",ntable);
table_type d_table("DeviceTable",ntable);
for (int i = 0; i < ntable; i++) {
h_table(i) = dftable[i];
}
Kokkos::deep_copy(d_table,h_table);
d_dftable = d_table;
}
{
host_table_type h_table("HostTable",ntable);
table_type d_table("DeviceTable",ntable);
// Copy ctable and dctable
for (int i = 0; i < ntable; i++) {
h_table(i) = ctable[i];
}
Kokkos::deep_copy(d_table,h_table);
d_ctable = d_table;
}
{
host_table_type h_table("HostTable",ntable);
table_type d_table("DeviceTable",ntable);
for (int i = 0; i < ntable; i++) {
h_table(i) = dctable[i];
}
Kokkos::deep_copy(d_table,h_table);
d_dctable = d_table;
}
{
host_table_type h_table("HostTable",ntable);
table_type d_table("DeviceTable",ntable);
// Copy etable and detable
for (int i = 0; i < ntable; i++) {
h_table(i) = etable[i];
}
Kokkos::deep_copy(d_table,h_table);
d_etable = d_table;
}
{
host_table_type h_table("HostTable",ntable);
table_type d_table("DeviceTable",ntable);
for (int i = 0; i < ntable; i++) {
h_table(i) = detable[i];
}
Kokkos::deep_copy(d_table,h_table);
d_detable = d_table;
}
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
template<class DeviceType>
void PairLJCharmmfswCoulLongKokkos<DeviceType>::init_style()
{
PairLJCharmmfswCoulLong::init_style();
Kokkos::deep_copy(d_cut_ljsq,cut_ljsq);
Kokkos::deep_copy(d_cut_coulsq,cut_coulsq);
// error if rRESPA with inner levels
if (update->whichflag == 1 && utils::strmatch(update->integrate_style,"^respa")) {
int respa = 0;
if (((Respa *) update->integrate)->level_inner >= 0) respa = 1;
if (((Respa *) update->integrate)->level_middle >= 0) respa = 2;
if (respa)
error->all(FLERR,"Cannot use Kokkos pair style with rRESPA inner/middle");
}
// adjust neighbor list request for KOKKOS
neighflag = lmp->kokkos->neighflag;
auto request = neighbor->find_request(this);
request->set_kokkos_host(std::is_same_v<DeviceType,LMPHostType> &&
!std::is_same_v<DeviceType,LMPDeviceType>);
request->set_kokkos_device(std::is_same_v<DeviceType,LMPDeviceType>);
if (neighflag == FULL) request->enable_full();
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
template<class DeviceType>
double PairLJCharmmfswCoulLongKokkos<DeviceType>::init_one(int i, int j)
{
double cutone = PairLJCharmmfswCoulLong::init_one(i,j);
k_params.h_view(i,j).lj1 = lj1[i][j];
k_params.h_view(i,j).lj2 = lj2[i][j];
k_params.h_view(i,j).lj3 = lj3[i][j];
k_params.h_view(i,j).lj4 = lj4[i][j];
//k_params.h_view(i,j).offset = offset[i][j];
k_params.h_view(i,j).cut_ljsq = cut_ljsq;
k_params.h_view(i,j).cut_coulsq = cut_coulsq;
k_params.h_view(j,i) = k_params.h_view(i,j);
if (i<MAX_TYPES_STACKPARAMS+1 && j<MAX_TYPES_STACKPARAMS+1) {
m_params[i][j] = m_params[j][i] = k_params.h_view(i,j);
m_cutsq[j][i] = m_cutsq[i][j] = cutone*cutone;
m_cut_ljsq[j][i] = m_cut_ljsq[i][j] = cut_ljsq;
m_cut_coulsq[j][i] = m_cut_coulsq[i][j] = cut_coulsq;
}
k_cutsq.h_view(i,j) = k_cutsq.h_view(j,i) = cutone*cutone;
k_cutsq.template modify<LMPHostType>();
k_params.template modify<LMPHostType>();
return cutone;
}
namespace LAMMPS_NS {
template class PairLJCharmmfswCoulLongKokkos<LMPDeviceType>;
#ifdef LMP_KOKKOS_GPU
template class PairLJCharmmfswCoulLongKokkos<LMPHostType>;
#endif
}

145
src/KOKKOS/pair_lj_charmmfsw_coul_long_kokkos.h Normal file
View File

@ -0,0 +1,145 @@
/* -*- 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 PAIR_CLASS
// clang-format off
PairStyle(lj/charmmfsw/coul/long/kk,PairLJCharmmfswCoulLongKokkos<LMPDeviceType>);
PairStyle(lj/charmmfsw/coul/long/kk/device,PairLJCharmmfswCoulLongKokkos<LMPDeviceType>);
PairStyle(lj/charmmfsw/coul/long/kk/host,PairLJCharmmfswCoulLongKokkos<LMPHostType>);
// clang-format on
#else
// clang-format off
#ifndef LMP_PAIR_LJ_CHARMMFSW_COUL_LONG_KOKKOS_H
#define LMP_PAIR_LJ_CHARMMFSW_COUL_LONG_KOKKOS_H
#include "pair_kokkos.h"
#include "pair_lj_charmmfsw_coul_long.h"
#include "neigh_list_kokkos.h"
namespace LAMMPS_NS {
template<class DeviceType>
class PairLJCharmmfswCoulLongKokkos : public PairLJCharmmfswCoulLong {
public:
enum {EnabledNeighFlags=FULL|HALFTHREAD|HALF};
enum {COUL_FLAG=1};
typedef DeviceType device_type;
typedef ArrayTypes<DeviceType> AT;
PairLJCharmmfswCoulLongKokkos(class LAMMPS *);
~PairLJCharmmfswCoulLongKokkos() override;
void compute(int, int) override;
void init_tables(double cut_coul, double *cut_respa) override;
void init_style() override;
double init_one(int, int) override;
protected:
template<bool STACKPARAMS, class Specialisation>
KOKKOS_INLINE_FUNCTION
F_FLOAT compute_fpair(const F_FLOAT& rsq, const int& i, const int&j,
const int& itype, const int& jtype) const;
template<bool STACKPARAMS, class Specialisation>
KOKKOS_INLINE_FUNCTION
F_FLOAT compute_fcoul(const F_FLOAT& rsq, const int& i, const int&j, const int& itype,
const int& jtype, const F_FLOAT& factor_coul, const F_FLOAT& qtmp) const;
template<bool STACKPARAMS, class Specialisation>
KOKKOS_INLINE_FUNCTION
F_FLOAT compute_evdwl(const F_FLOAT& rsq, const int& i, const int&j,
const int& itype, const int& jtype) const;
template<bool STACKPARAMS, class Specialisation>
KOKKOS_INLINE_FUNCTION
F_FLOAT compute_ecoul(const F_FLOAT& rsq, const int& i, const int&j,
const int& itype, const int& jtype, const F_FLOAT& factor_coul, const F_FLOAT& qtmp) const;
Kokkos::DualView<params_lj_coul**,Kokkos::LayoutRight,DeviceType> k_params;
typename Kokkos::DualView<params_lj_coul**,
Kokkos::LayoutRight,DeviceType>::t_dev_const_um params;
// hardwired to space for 12 atom types
params_lj_coul m_params[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
F_FLOAT m_cutsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
F_FLOAT m_cut_ljsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
F_FLOAT m_cut_coulsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
typename AT::t_x_array_randomread x;
typename AT::t_x_array c_x;
typename AT::t_f_array f;
typename AT::t_int_1d_randomread type;
typename AT::t_float_1d_randomread q;
DAT::tdual_efloat_1d k_eatom;
DAT::tdual_virial_array k_vatom;
typename AT::t_efloat_1d d_eatom;
typename AT::t_virial_array d_vatom;
int newton_pair;
typename AT::tdual_ffloat_2d k_cutsq;
typename AT::t_ffloat_2d d_cutsq;
typename AT::t_ffloat_2d d_cut_ljsq;
typename AT::t_ffloat_2d d_cut_coulsq;
typename AT::t_ffloat_1d_randomread
d_rtable, d_drtable, d_ftable, d_dftable,
d_ctable, d_dctable, d_etable, d_detable;
int neighflag;
int nlocal,nall,eflag,vflag;
double special_coul[4];
double special_lj[4];
double qqrd2e;
void allocate() override;
friend struct PairComputeFunctor<PairLJCharmmfswCoulLongKokkos,FULL,true,0,CoulLongTable<1>>;
friend struct PairComputeFunctor<PairLJCharmmfswCoulLongKokkos,FULL,true,1,CoulLongTable<1>>;
friend struct PairComputeFunctor<PairLJCharmmfswCoulLongKokkos,HALF,true,0,CoulLongTable<1>>;
friend struct PairComputeFunctor<PairLJCharmmfswCoulLongKokkos,HALFTHREAD,true,0,CoulLongTable<1>>;
friend struct PairComputeFunctor<PairLJCharmmfswCoulLongKokkos,FULL,false,0,CoulLongTable<1>>;
friend struct PairComputeFunctor<PairLJCharmmfswCoulLongKokkos,FULL,false,1,CoulLongTable<1>>;
friend struct PairComputeFunctor<PairLJCharmmfswCoulLongKokkos,HALF,false,0,CoulLongTable<1>>;
friend struct PairComputeFunctor<PairLJCharmmfswCoulLongKokkos,HALFTHREAD,false,0,CoulLongTable<1>>;
friend EV_FLOAT pair_compute_neighlist<PairLJCharmmfswCoulLongKokkos,FULL,0,CoulLongTable<1>>(PairLJCharmmfswCoulLongKokkos*,NeighListKokkos<DeviceType>*);
friend EV_FLOAT pair_compute_neighlist<PairLJCharmmfswCoulLongKokkos,FULL,1,CoulLongTable<1>>(PairLJCharmmfswCoulLongKokkos*,NeighListKokkos<DeviceType>*);
friend EV_FLOAT pair_compute_neighlist<PairLJCharmmfswCoulLongKokkos,HALF,0,CoulLongTable<1>>(PairLJCharmmfswCoulLongKokkos*,NeighListKokkos<DeviceType>*);
friend EV_FLOAT pair_compute_neighlist<PairLJCharmmfswCoulLongKokkos,HALFTHREAD,0,CoulLongTable<1>>(PairLJCharmmfswCoulLongKokkos*,NeighListKokkos<DeviceType>*);
friend EV_FLOAT pair_compute<PairLJCharmmfswCoulLongKokkos,CoulLongTable<1>>(PairLJCharmmfswCoulLongKokkos*,
NeighListKokkos<DeviceType>*);
friend struct PairComputeFunctor<PairLJCharmmfswCoulLongKokkos,FULL,true,0,CoulLongTable<0>>;
friend struct PairComputeFunctor<PairLJCharmmfswCoulLongKokkos,FULL,true,1,CoulLongTable<0>>;
friend struct PairComputeFunctor<PairLJCharmmfswCoulLongKokkos,HALF,true,0,CoulLongTable<0>>;
friend struct PairComputeFunctor<PairLJCharmmfswCoulLongKokkos,HALFTHREAD,true,0,CoulLongTable<0>>;
friend struct PairComputeFunctor<PairLJCharmmfswCoulLongKokkos,FULL,false,0,CoulLongTable<0>>;
friend struct PairComputeFunctor<PairLJCharmmfswCoulLongKokkos,FULL,false,1,CoulLongTable<0>>;
friend struct PairComputeFunctor<PairLJCharmmfswCoulLongKokkos,HALF,false,0,CoulLongTable<0>>;
friend struct PairComputeFunctor<PairLJCharmmfswCoulLongKokkos,HALFTHREAD,false,0,CoulLongTable<0>>;
friend EV_FLOAT pair_compute_neighlist<PairLJCharmmfswCoulLongKokkos,FULL,0,CoulLongTable<0>>(PairLJCharmmfswCoulLongKokkos*,NeighListKokkos<DeviceType>*);
friend EV_FLOAT pair_compute_neighlist<PairLJCharmmfswCoulLongKokkos,FULL,1,CoulLongTable<0>>(PairLJCharmmfswCoulLongKokkos*,NeighListKokkos<DeviceType>*);
friend EV_FLOAT pair_compute_neighlist<PairLJCharmmfswCoulLongKokkos,HALF,0,CoulLongTable<0>>(PairLJCharmmfswCoulLongKokkos*,NeighListKokkos<DeviceType>*);
friend EV_FLOAT pair_compute_neighlist<PairLJCharmmfswCoulLongKokkos,HALFTHREAD,0,CoulLongTable<0>>(PairLJCharmmfswCoulLongKokkos*,NeighListKokkos<DeviceType>*);
friend EV_FLOAT pair_compute<PairLJCharmmfswCoulLongKokkos,CoulLongTable<0>>(PairLJCharmmfswCoulLongKokkos*,
NeighListKokkos<DeviceType>*);
friend void pair_virial_fdotr_compute<PairLJCharmmfswCoulLongKokkos>(PairLJCharmmfswCoulLongKokkos*);
};
}
#endif
#endif

4
src/KSPACE/pair_lj_charmmfsw_coul_long.cpp
View File

@ -76,6 +76,8 @@ PairLJCharmmfswCoulLong::PairLJCharmmfswCoulLong(LAMMPS *lmp) : Pair(lmp)
PairLJCharmmfswCoulLong::~PairLJCharmmfswCoulLong() PairLJCharmmfswCoulLong::~PairLJCharmmfswCoulLong()
{ {
if (copymode) return;
// switch qqr2e back from CHARMM value to LAMMPS value // switch qqr2e back from CHARMM value to LAMMPS value
if (update && strcmp(update->unit_style,"real") == 0) { if (update && strcmp(update->unit_style,"real") == 0) {
@ -85,8 +87,6 @@ PairLJCharmmfswCoulLong::~PairLJCharmmfswCoulLong()
force->qqr2e = force->qqr2e_lammps_real; force->qqr2e = force->qqr2e_lammps_real;
} }
if (copymode) return;
if (allocated) { if (allocated) {
memory->destroy(setflag); memory->destroy(setflag);
memory->destroy(cutsq); memory->destroy(cutsq);

74
src/LEPTON/angle_lepton.cpp
View File

@ -44,6 +44,7 @@ AngleLepton::AngleLepton(LAMMPS *_lmp) :
{ {
writedata = 1; writedata = 1;
reinitflag = 0; reinitflag = 0;
auto_offset = 1;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -90,10 +91,21 @@ template <int EVFLAG, int EFLAG, int NEWTON_BOND> void AngleLepton::eval()
{ {
std::vector<Lepton::CompiledExpression> angleforce; std::vector<Lepton::CompiledExpression> angleforce;
std::vector<Lepton::CompiledExpression> anglepot; std::vector<Lepton::CompiledExpression> anglepot;
for (const auto &expr : expressions) { std::vector<bool> has_ref;
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp)); try {
angleforce.emplace_back(parsed.differentiate("theta").createCompiledExpression()); for (const auto &expr : expressions) {
if (EFLAG) anglepot.emplace_back(parsed.createCompiledExpression()); auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp));
angleforce.emplace_back(parsed.differentiate("theta").createCompiledExpression());
has_ref.push_back(true);
try {
angleforce.back().getVariableReference("theta");
} catch (Lepton::Exception &) {
has_ref.back() = false;
}
if (EFLAG) anglepot.emplace_back(parsed.createCompiledExpression());
}
} catch (std::exception &e) {
error->all(FLERR, e.what());
} }
const double *const *const x = atom->x; const double *const *const x = atom->x;
@ -142,8 +154,7 @@ template <int EVFLAG, int EFLAG, int NEWTON_BOND> void AngleLepton::eval()
const double dtheta = acos(c) - theta0[type]; const double dtheta = acos(c) - theta0[type];
const int idx = type2expression[type]; const int idx = type2expression[type];
angleforce[idx].getVariableReference("theta") = dtheta; if (has_ref[idx]) angleforce[idx].getVariableReference("theta") = dtheta;
const double a = -angleforce[idx].evaluate() * s; const double a = -angleforce[idx].evaluate() * s;
const double a11 = a * c / rsq1; const double a11 = a * c / rsq1;
const double a12 = -a / (r1 * r2); const double a12 = -a / (r1 * r2);
@ -179,7 +190,11 @@ template <int EVFLAG, int EFLAG, int NEWTON_BOND> void AngleLepton::eval()
double eangle = 0.0; double eangle = 0.0;
if (EFLAG) { if (EFLAG) {
anglepot[idx].getVariableReference("theta") = dtheta; try {
anglepot[idx].getVariableReference("theta") = dtheta;
} catch (Lepton::Exception &) {
; // ignore -> constant force
}
eangle = anglepot[idx].evaluate() - offset[type]; eangle = anglepot[idx].evaluate() - offset[type];
} }
if (EVFLAG) if (EVFLAG)
@ -202,6 +217,24 @@ void AngleLepton::allocate()
for (int i = 1; i < np1; i++) setflag[i] = 0; for (int i = 1; i < np1; i++) setflag[i] = 0;
} }
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void AngleLepton::settings(int narg, char **arg)
{
auto_offset = 1;
if (narg > 0) {
if (strcmp(arg[0],"auto_offset") == 0) {
auto_offset = 1;
} else if (strcmp(arg[0],"no_offset") == 0) {
auto_offset = 0;
} else {
error->all(FLERR, "Unknown angle style lepton setting {}", arg[0]);
}
}
}
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
set coeffs for one or more types set coeffs for one or more types
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
@ -224,9 +257,20 @@ void AngleLepton::coeff(int narg, char **arg)
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(exp_one, lmp)); auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(exp_one, lmp));
auto anglepot = parsed.createCompiledExpression(); auto anglepot = parsed.createCompiledExpression();
auto angleforce = parsed.differentiate("theta").createCompiledExpression(); auto angleforce = parsed.differentiate("theta").createCompiledExpression();
anglepot.getVariableReference("theta") = 0.0; try {
angleforce.getVariableReference("theta") = 0.0; anglepot.getVariableReference("theta") = 0.0;
offset_one = anglepot.evaluate(); } catch (Lepton::Exception &) {
if (comm->me == 0)
error->warning(FLERR, "Lepton potential expression {} does not depend on 'theta'", exp_one);
}
try {
angleforce.getVariableReference("theta") = 0.0;
} catch (Lepton::Exception &) {
if (comm->me == 0)
error->warning(FLERR, "Force from Lepton expression {} does not depend on 'theta'",
exp_one);
}
if (auto_offset) offset_one = anglepot.evaluate();
angleforce.evaluate(); angleforce.evaluate();
} catch (std::exception &e) { } catch (std::exception &e) {
error->all(FLERR, e.what()); error->all(FLERR, e.what());
@ -284,6 +328,7 @@ void AngleLepton::write_restart(FILE *fp)
fwrite(&n, sizeof(int), 1, fp); fwrite(&n, sizeof(int), 1, fp);
fwrite(exp.c_str(), sizeof(char), n, fp); fwrite(exp.c_str(), sizeof(char), n, fp);
} }
fwrite(&auto_offset, sizeof(int), 1, fp);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -323,6 +368,9 @@ void AngleLepton::read_restart(FILE *fp)
expressions.emplace_back(buf); expressions.emplace_back(buf);
} }
if (comm->me == 0) utils::sfread(FLERR, &auto_offset, sizeof(int), 1, fp, nullptr, error);
MPI_Bcast(&auto_offset, 1, MPI_INT, 0, world);
delete[] buf; delete[] buf;
} }
@ -363,7 +411,11 @@ double AngleLepton::single(int type, int i1, int i2, int i3)
const auto &expr = expressions[type2expression[type]]; const auto &expr = expressions[type2expression[type]];
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp)); auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp));
auto anglepot = parsed.createCompiledExpression(); auto anglepot = parsed.createCompiledExpression();
anglepot.getVariableReference("theta") = dtheta; try {
anglepot.getVariableReference("theta") = dtheta;
} catch (Lepton::Exception &) {
; // ignore -> constant potential
}
return anglepot.evaluate() - offset[type]; return anglepot.evaluate() - offset[type];
} }

2
src/LEPTON/angle_lepton.h
View File

@ -29,6 +29,7 @@ class AngleLepton : public Angle {
AngleLepton(class LAMMPS *); AngleLepton(class LAMMPS *);
~AngleLepton() override; ~AngleLepton() override;
void compute(int, int) override; void compute(int, int) override;
void settings(int, char **) override;
void coeff(int, char **) override; void coeff(int, char **) override;
double equilibrium_angle(int) override; double equilibrium_angle(int) override;
void write_restart(FILE *) override; void write_restart(FILE *) override;
@ -42,6 +43,7 @@ class AngleLepton : public Angle {
double *theta0; double *theta0;
int *type2expression; int *type2expression;
double *offset; double *offset;
int auto_offset;
virtual void allocate(); virtual void allocate();

62
src/LEPTON/bond_lepton.cpp
View File

@ -37,6 +37,7 @@ BondLepton::BondLepton(LAMMPS *_lmp) :
{ {
writedata = 1; writedata = 1;
reinitflag = 0; reinitflag = 0;
auto_offset = 1;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -82,10 +83,17 @@ template <int EVFLAG, int EFLAG, int NEWTON_BOND> void BondLepton::eval()
{ {
std::vector<Lepton::CompiledExpression> bondforce; std::vector<Lepton::CompiledExpression> bondforce;
std::vector<Lepton::CompiledExpression> bondpot; std::vector<Lepton::CompiledExpression> bondpot;
std::vector<bool> has_ref;
try { try {
for (const auto &expr : expressions) { for (const auto &expr : expressions) {
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp)); auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp));
bondforce.emplace_back(parsed.differentiate("r").createCompiledExpression()); bondforce.emplace_back(parsed.differentiate("r").createCompiledExpression());
has_ref.push_back(true);
try {
bondforce.back().getVariableReference("r");
} catch (Lepton::Exception &) {
has_ref.back() = false;
}
if (EFLAG) bondpot.emplace_back(parsed.createCompiledExpression()); if (EFLAG) bondpot.emplace_back(parsed.createCompiledExpression());
} }
} catch (std::exception &e) { } catch (std::exception &e) {
@ -116,7 +124,7 @@ template <int EVFLAG, int EFLAG, int NEWTON_BOND> void BondLepton::eval()
double fbond = 0.0; double fbond = 0.0;
if (r > 0.0) { if (r > 0.0) {
bondforce[idx].getVariableReference("r") = dr; if (has_ref[idx]) bondforce[idx].getVariableReference("r") = dr;
fbond = -bondforce[idx].evaluate() / r; fbond = -bondforce[idx].evaluate() / r;
} }
@ -136,7 +144,11 @@ template <int EVFLAG, int EFLAG, int NEWTON_BOND> void BondLepton::eval()
double ebond = 0.0; double ebond = 0.0;
if (EFLAG) { if (EFLAG) {
bondpot[idx].getVariableReference("r") = dr; try {
bondpot[idx].getVariableReference("r") = dr;
} catch (Lepton::Exception &) {
; // ignore -> constant potential
}
ebond = bondpot[idx].evaluate() - offset[type]; ebond = bondpot[idx].evaluate() - offset[type];
} }
if (EVFLAG) ev_tally(i1, i2, nlocal, NEWTON_BOND, ebond, fbond, delx, dely, delz); if (EVFLAG) ev_tally(i1, i2, nlocal, NEWTON_BOND, ebond, fbond, delx, dely, delz);
@ -157,6 +169,24 @@ void BondLepton::allocate()
for (int i = 1; i < np1; i++) setflag[i] = 0; for (int i = 1; i < np1; i++) setflag[i] = 0;
} }
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void BondLepton::settings(int narg, char **arg)
{
auto_offset = 1;
if (narg > 0) {
if (strcmp(arg[0],"auto_offset") == 0) {
auto_offset = 1;
} else if (strcmp(arg[0],"no_offset") == 0) {
auto_offset = 0;
} else {
error->all(FLERR, "Unknown bond style lepton setting {}", arg[0]);
}
}
}
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
set coeffs for one or more types set coeffs for one or more types
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
@ -179,9 +209,19 @@ void BondLepton::coeff(int narg, char **arg)
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(exp_one, lmp)); auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(exp_one, lmp));
auto bondpot = parsed.createCompiledExpression(); auto bondpot = parsed.createCompiledExpression();
auto bondforce = parsed.differentiate("r").createCompiledExpression(); auto bondforce = parsed.differentiate("r").createCompiledExpression();
bondpot.getVariableReference("r") = 0.0; try {
bondforce.getVariableReference("r") = 0.0; bondpot.getVariableReference("r") = 0.0;
offset_one = bondpot.evaluate(); } catch (Lepton::Exception &e) {
if (comm->me == 0)
error->warning(FLERR, "Lepton potential expression {} does not depend on 'r'", exp_one);
}
try {
bondforce.getVariableReference("r") = 0.0;
} catch (Lepton::Exception &e) {
if (comm->me == 0)
error->warning(FLERR, "Force from Lepton expression {} does not depend on 'r'", exp_one);
}
if (auto_offset) offset_one = bondpot.evaluate();
bondforce.evaluate(); bondforce.evaluate();
} catch (std::exception &e) { } catch (std::exception &e) {
error->all(FLERR, e.what()); error->all(FLERR, e.what());
@ -239,6 +279,7 @@ void BondLepton::write_restart(FILE *fp)
fwrite(&n, sizeof(int), 1, fp); fwrite(&n, sizeof(int), 1, fp);
fwrite(exp.c_str(), sizeof(char), n, fp); fwrite(exp.c_str(), sizeof(char), n, fp);
} }
fwrite(&auto_offset, sizeof(int), 1, fp);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -278,6 +319,9 @@ void BondLepton::read_restart(FILE *fp)
expressions.emplace_back(buf); expressions.emplace_back(buf);
} }
if (comm->me == 0) utils::sfread(FLERR, &auto_offset, sizeof(int), 1, fp, nullptr, error);
MPI_Bcast(&auto_offset, 1, MPI_INT, 0, world);
delete[] buf; delete[] buf;
} }
@ -302,8 +346,12 @@ double BondLepton::single(int type, double rsq, int /*i*/, int /*j*/, double &ff
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp)); auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp));
auto bondpot = parsed.createCompiledExpression(); auto bondpot = parsed.createCompiledExpression();
auto bondforce = parsed.differentiate("r").createCompiledExpression(); auto bondforce = parsed.differentiate("r").createCompiledExpression();
bondforce.getVariableReference("r") = dr; try {
bondpot.getVariableReference("r") = dr; bondpot.getVariableReference("r") = dr;
bondforce.getVariableReference("r") = dr;
} catch (Lepton::Exception &) {
; // ignore -> constant potential or force
}
// force and energy // force and energy

2
src/LEPTON/bond_lepton.h
View File

@ -29,6 +29,7 @@ class BondLepton : public Bond {
BondLepton(class LAMMPS *); BondLepton(class LAMMPS *);
~BondLepton() override; ~BondLepton() override;
void compute(int, int) override; void compute(int, int) override;
void settings(int, char **) override;
void coeff(int, char **) override; void coeff(int, char **) override;
double equilibrium_distance(int) override; double equilibrium_distance(int) override;
void write_restart(FILE *) override; void write_restart(FILE *) override;
@ -42,6 +43,7 @@ class BondLepton : public Bond {
double *r0; double *r0;
int *type2expression; int *type2expression;
double *offset; double *offset;
int auto_offset;
virtual void allocate(); virtual void allocate();

29
src/LEPTON/dihedral_lepton.cpp
View File

@ -92,10 +92,17 @@ template <int EVFLAG, int EFLAG, int NEWTON_BOND> void DihedralLepton::eval()
{ {
std::vector<Lepton::CompiledExpression> dihedralforce; std::vector<Lepton::CompiledExpression> dihedralforce;
std::vector<Lepton::CompiledExpression> dihedralpot; std::vector<Lepton::CompiledExpression> dihedralpot;
std::vector<bool> has_ref;
try { try {
for (const auto &expr : expressions) { for (const auto &expr : expressions) {
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp)); auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp));
dihedralforce.emplace_back(parsed.differentiate("phi").createCompiledExpression()); dihedralforce.emplace_back(parsed.differentiate("phi").createCompiledExpression());
has_ref.push_back(true);
try {
dihedralforce.back().getVariableReference("phi");
} catch (Lepton::Exception &) {
has_ref.back() = false;
}
if (EFLAG) dihedralpot.emplace_back(parsed.createCompiledExpression()); if (EFLAG) dihedralpot.emplace_back(parsed.createCompiledExpression());
} }
} catch (std::exception &e) { } catch (std::exception &e) {
@ -278,7 +285,7 @@ template <int EVFLAG, int EFLAG, int NEWTON_BOND> void DihedralLepton::eval()
} }
const int idx = type2expression[type]; const int idx = type2expression[type];
dihedralforce[idx].getVariableReference("phi") = phi; if (has_ref[idx]) dihedralforce[idx].getVariableReference("phi") = phi;
double m_du_dphi = -dihedralforce[idx].evaluate(); double m_du_dphi = -dihedralforce[idx].evaluate();
// ----- Step 4: Calculate the force direction in real space ----- // ----- Step 4: Calculate the force direction in real space -----
@ -322,7 +329,11 @@ template <int EVFLAG, int EFLAG, int NEWTON_BOND> void DihedralLepton::eval()
double edihedral = 0.0; double edihedral = 0.0;
if (EFLAG) { if (EFLAG) {
dihedralpot[idx].getVariableReference("phi") = phi; try {
dihedralpot[idx].getVariableReference("phi") = phi;
} catch (Lepton::Exception &) {
; // ignore -> constant potential
}
edihedral = dihedralpot[idx].evaluate(); edihedral = dihedralpot[idx].evaluate();
} }
if (EVFLAG) if (EVFLAG)
@ -362,8 +373,18 @@ void DihedralLepton::coeff(int narg, char **arg)
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(exp_one, lmp)); auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(exp_one, lmp));
auto dihedralpot = parsed.createCompiledExpression(); auto dihedralpot = parsed.createCompiledExpression();
auto dihedralforce = parsed.differentiate("phi").createCompiledExpression(); auto dihedralforce = parsed.differentiate("phi").createCompiledExpression();
dihedralpot.getVariableReference("phi") = 0.0; try {
dihedralforce.getVariableReference("phi") = 0.0; dihedralpot.getVariableReference("phi") = 0.0;
} catch (Lepton::Exception &) {
if (comm->me == 0)
error->warning(FLERR, "Lepton potential expression {} does not depend on 'phi'", exp_one);
}
try {
dihedralforce.getVariableReference("phi") = 0.0;
} catch (Lepton::Exception &) {
if (comm->me == 0)
error->warning(FLERR, "Force from Lepton expression {} does not depend on 'phi'", exp_one);
}
dihedralforce.evaluate(); dihedralforce.evaluate();
} catch (std::exception &e) { } catch (std::exception &e) {
error->all(FLERR, e.what()); error->all(FLERR, e.what());

15
src/LEPTON/fix_wall_lepton.cpp
View File

@ -13,6 +13,7 @@
#include "fix_wall_lepton.h" #include "fix_wall_lepton.h"
#include "atom.h" #include "atom.h"
#include "comm.h"
#include "error.h" #include "error.h"
#include "Lepton.h" #include "Lepton.h"
@ -41,8 +42,18 @@ void FixWallLepton::post_constructor()
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(exp_one, lmp)); auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(exp_one, lmp));
auto wallpot = parsed.createCompiledExpression(); auto wallpot = parsed.createCompiledExpression();
auto wallforce = parsed.differentiate("r").createCompiledExpression(); auto wallforce = parsed.differentiate("r").createCompiledExpression();
wallpot.getVariableReference("r") = 0.0; try {
wallforce.getVariableReference("r") = 0.0; wallpot.getVariableReference("r") = 0.0;
} catch (Lepton::Exception &) {
if (comm->me == 0)
error->warning(FLERR, "Lepton potential expression {} does not depend on 'r'", exp_one);
}
try {
wallforce.getVariableReference("r") = 0.0;
} catch (Lepton::Exception &) {
if (comm->me == 0)
error->warning(FLERR, "Force from Lepton expression {} does not depend on 'r'", exp_one);
}
wallpot.evaluate(); wallpot.evaluate();
wallforce.evaluate(); wallforce.evaluate();
} catch (std::exception &e) { } catch (std::exception &e) {

41
src/LEPTON/pair_lepton.cpp
View File

@ -27,6 +27,7 @@
#include "Lepton.h" #include "Lepton.h"
#include "lepton_utils.h" #include "lepton_utils.h"
#include <array>
#include <cmath> #include <cmath>
#include <map> #include <map>
@ -105,11 +106,17 @@ template <int EVFLAG, int EFLAG, int NEWTON_PAIR> void PairLepton::eval()
std::vector<Lepton::CompiledExpression> pairforce; std::vector<Lepton::CompiledExpression> pairforce;
std::vector<Lepton::CompiledExpression> pairpot; std::vector<Lepton::CompiledExpression> pairpot;
std::vector<bool> has_ref;
try { try {
for (const auto &expr : expressions) { for (const auto &expr : expressions) {
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp), functions); auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp), functions);
pairforce.emplace_back(parsed.differentiate("r").createCompiledExpression()); pairforce.emplace_back(parsed.differentiate("r").createCompiledExpression());
pairforce.back().getVariableReference("r"); has_ref.push_back(true);
try {
pairforce.back().getVariableReference("r");
} catch (Lepton::Exception &) {
has_ref.back() = false;
}
if (EFLAG) pairpot.emplace_back(parsed.createCompiledExpression()); if (EFLAG) pairpot.emplace_back(parsed.createCompiledExpression());
} }
} catch (std::exception &e) { } catch (std::exception &e) {
@ -142,8 +149,7 @@ template <int EVFLAG, int EFLAG, int NEWTON_PAIR> void PairLepton::eval()
if (rsq < cutsq[itype][jtype]) { if (rsq < cutsq[itype][jtype]) {
const double r = sqrt(rsq); const double r = sqrt(rsq);
const int idx = type2expression[itype][jtype]; const int idx = type2expression[itype][jtype];
double &r_for = pairforce[idx].getVariableReference("r"); if (has_ref[idx]) pairforce[idx].getVariableReference("r") = r;
r_for = r;
const double fpair = -pairforce[idx].evaluate() / r * factor_lj; const double fpair = -pairforce[idx].evaluate() / r * factor_lj;
fxtmp += delx * fpair; fxtmp += delx * fpair;
@ -157,7 +163,11 @@ template <int EVFLAG, int EFLAG, int NEWTON_PAIR> void PairLepton::eval()
double evdwl = 0.0; double evdwl = 0.0;
if (EFLAG) { if (EFLAG) {
pairpot[idx].getVariableReference("r") = r; try {
pairpot[idx].getVariableReference("r") = r;
} catch (Lepton::Exception &) {
; // ignore -> constant potential
}
evdwl = pairpot[idx].evaluate() - offset[itype][jtype]; evdwl = pairpot[idx].evaluate() - offset[itype][jtype];
evdwl *= factor_lj; evdwl *= factor_lj;
} }
@ -229,8 +239,12 @@ void PairLepton::coeff(int narg, char **arg)
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(exp_one, lmp), functions); auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(exp_one, lmp), functions);
auto pairforce = parsed.differentiate("r").createCompiledExpression(); auto pairforce = parsed.differentiate("r").createCompiledExpression();
auto pairpot = parsed.createCompiledExpression(); auto pairpot = parsed.createCompiledExpression();
pairpot.getVariableReference("r") = 1.0; try {
pairforce.getVariableReference("r") = 1.0; pairpot.getVariableReference("r") = 1.0;
pairforce.getVariableReference("r") = 1.0;
} catch (Lepton::Exception &) {
; // ignore -> constant potential or force
}
pairpot.evaluate(); pairpot.evaluate();
pairforce.evaluate(); pairforce.evaluate();
} catch (std::exception &e) { } catch (std::exception &e) {
@ -270,7 +284,11 @@ double PairLepton::init_one(int i, int j)
try { try {
auto expr = LeptonUtils::substitute(expressions[type2expression[i][j]], lmp); auto expr = LeptonUtils::substitute(expressions[type2expression[i][j]], lmp);
auto pairpot = Lepton::Parser::parse(expr, functions).createCompiledExpression(); auto pairpot = Lepton::Parser::parse(expr, functions).createCompiledExpression();
pairpot.getVariableReference("r") = cut[i][j]; try {
pairpot.getVariableReference("r") = cut[i][j];
} catch (Lepton::Exception &) {
; // ignore -> constant potential
}
offset[i][j] = pairpot.evaluate(); offset[i][j] = pairpot.evaluate();
} catch (std::exception &) { } catch (std::exception &) {
} }
@ -429,9 +447,12 @@ double PairLepton::single(int /* i */, int /* j */, int itype, int jtype, double
auto pairforce = parsed.differentiate("r").createCompiledExpression(); auto pairforce = parsed.differentiate("r").createCompiledExpression();
const double r = sqrt(rsq); const double r = sqrt(rsq);
pairpot.getVariableReference("r") = r; try {
pairforce.getVariableReference("r") = r; pairpot.getVariableReference("r") = r;
pairforce.getVariableReference("r") = r;
} catch (Lepton::Exception &) {
; // ignore -> constant potential or force
}
fforce = -pairforce.evaluate() / r * factor_lj; fforce = -pairforce.evaluate() / r * factor_lj;
return (pairpot.evaluate() - offset[itype][jtype]) * factor_lj; return (pairpot.evaluate() - offset[itype][jtype]) * factor_lj;
} }

50
src/LEPTON/pair_lepton_coul.cpp
View File

@ -28,6 +28,8 @@
#include "Lepton.h" #include "Lepton.h"
#include "lepton_utils.h" #include "lepton_utils.h"
#include <array>
#include <cmath> #include <cmath>
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
@ -79,25 +81,30 @@ template <int EVFLAG, int EFLAG, int NEWTON_PAIR> void PairLeptonCoul::eval()
std::vector<Lepton::CompiledExpression> pairforce; std::vector<Lepton::CompiledExpression> pairforce;
std::vector<Lepton::CompiledExpression> pairpot; std::vector<Lepton::CompiledExpression> pairpot;
std::vector<std::pair<bool, bool>> have_q; std::vector<std::array<bool, 3>> has_ref;
try { try {
for (const auto &expr : expressions) { for (const auto &expr : expressions) {
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp), functions); auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp), functions);
pairforce.emplace_back(parsed.differentiate("r").createCompiledExpression()); pairforce.emplace_back(parsed.differentiate("r").createCompiledExpression());
has_ref.push_back({true, true, true});
try {
pairforce.back().getVariableReference("r");
} catch (Lepton::Exception &) {
has_ref.back()[0] = false;
}
if (EFLAG) pairpot.emplace_back(parsed.createCompiledExpression()); if (EFLAG) pairpot.emplace_back(parsed.createCompiledExpression());
pairforce.back().getVariableReference("r");
have_q.emplace_back(true, true);
// check if there are references to charges // check if there are references to charges
try { try {
pairforce.back().getVariableReference("qi"); pairforce.back().getVariableReference("qi");
} catch (std::exception &) { } catch (Lepton::Exception &) {
have_q.back().first = false; has_ref.back()[1] = false;
} }
try { try {
pairforce.back().getVariableReference("qj"); pairforce.back().getVariableReference("qj");
} catch (std::exception &) { } catch (Lepton::Exception &) {
have_q.back().second = false; has_ref.back()[2] = false;
} }
} }
} catch (std::exception &e) { } catch (std::exception &e) {
@ -130,9 +137,9 @@ template <int EVFLAG, int EFLAG, int NEWTON_PAIR> void PairLeptonCoul::eval()
if (rsq < cutsq[itype][jtype]) { if (rsq < cutsq[itype][jtype]) {
const double r = sqrt(rsq); const double r = sqrt(rsq);
const int idx = type2expression[itype][jtype]; const int idx = type2expression[itype][jtype];
pairforce[idx].getVariableReference("r") = r; if (has_ref[idx][0]) pairforce[idx].getVariableReference("r") = r;
if (have_q[idx].first) pairforce[idx].getVariableReference("qi") = q2e * q[i]; if (has_ref[idx][1]) pairforce[idx].getVariableReference("qi") = q2e * q[i];
if (have_q[idx].second) pairforce[idx].getVariableReference("qj") = q2e * q[j]; if (has_ref[idx][2]) pairforce[idx].getVariableReference("qj") = q2e * q[j];
const double fpair = -pairforce[idx].evaluate() / r * factor_coul; const double fpair = -pairforce[idx].evaluate() / r * factor_coul;
fxtmp += delx * fpair; fxtmp += delx * fpair;
@ -146,9 +153,14 @@ template <int EVFLAG, int EFLAG, int NEWTON_PAIR> void PairLeptonCoul::eval()
double ecoul = 0.0; double ecoul = 0.0;
if (EFLAG) { if (EFLAG) {
pairpot[idx].getVariableReference("r") = r; try {
if (have_q[idx].first) pairpot[idx].getVariableReference("qi") = q2e * q[i]; pairpot[idx].getVariableReference("r") = r;
if (have_q[idx].second) pairpot[idx].getVariableReference("qj") = q2e * q[j]; } catch (Lepton::Exception &) {
; // ignore -> constant potential
}
if (has_ref[idx][1]) pairpot[idx].getVariableReference("qi") = q2e * q[i];
if (has_ref[idx][2]) pairpot[idx].getVariableReference("qj") = q2e * q[j];
ecoul = pairpot[idx].evaluate(); ecoul = pairpot[idx].evaluate();
ecoul *= factor_coul; ecoul *= factor_coul;
} }
@ -249,18 +261,22 @@ double PairLeptonCoul::single(int i, int j, int itype, int jtype, double rsq, do
const double r = sqrt(rsq); const double r = sqrt(rsq);
const double q2e = sqrt(force->qqrd2e); const double q2e = sqrt(force->qqrd2e);
pairpot.getVariableReference("r") = r; try {
pairforce.getVariableReference("r") = r; pairpot.getVariableReference("r") = r;
pairforce.getVariableReference("r") = r;
} catch (Lepton::Exception &) {
; // ignore -> constant potential or force
}
try { try {
pairpot.getVariableReference("qi") = q2e * atom->q[i]; pairpot.getVariableReference("qi") = q2e * atom->q[i];
pairforce.getVariableReference("qi") = q2e * atom->q[i]; pairforce.getVariableReference("qi") = q2e * atom->q[i];
} catch (std::exception &) { } catch (Lepton::Exception &) {
/* ignore */ /* ignore */
} }
try { try {
pairpot.getVariableReference("qj") = q2e * atom->q[j]; pairpot.getVariableReference("qj") = q2e * atom->q[j];
pairforce.getVariableReference("qj") = q2e * atom->q[j]; pairforce.getVariableReference("qj") = q2e * atom->q[j];
} catch (std::exception &) { } catch (Lepton::Exception &) {
/* ignore */ /* ignore */
} }

55
src/LEPTON/pair_lepton_sphere.cpp
View File

@ -28,6 +28,7 @@
#include "Lepton.h" #include "Lepton.h"
#include "lepton_utils.h" #include "lepton_utils.h"
#include <array>
#include <cmath> #include <cmath>
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
@ -77,25 +78,30 @@ template <int EVFLAG, int EFLAG, int NEWTON_PAIR> void PairLeptonSphere::eval()
std::vector<Lepton::CompiledExpression> pairforce; std::vector<Lepton::CompiledExpression> pairforce;
std::vector<Lepton::CompiledExpression> pairpot; std::vector<Lepton::CompiledExpression> pairpot;
std::vector<std::pair<bool, bool>> have_rad; std::vector<std::array<bool, 3>> has_ref;
try { try {
for (const auto &expr : expressions) { for (const auto &expr : expressions) {
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp), functions); auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp), functions);
pairforce.emplace_back(parsed.differentiate("r").createCompiledExpression()); pairforce.emplace_back(parsed.differentiate("r").createCompiledExpression());
has_ref.push_back({true, true, true});
try {
pairforce.back().getVariableReference("r");
} catch (Lepton::Exception &) {
has_ref.back()[0] = false;
}
if (EFLAG) pairpot.emplace_back(parsed.createCompiledExpression()); if (EFLAG) pairpot.emplace_back(parsed.createCompiledExpression());
pairforce.back().getVariableReference("r");
have_rad.emplace_back(true, true);
// check if there are references to charges // check if there are references to radii
try { try {
pairforce.back().getVariableReference("radi"); pairforce.back().getVariableReference("radi");
} catch (std::exception &) { } catch (Lepton::Exception &) {
have_rad.back().first = false; has_ref.back()[1] = false;
} }
try { try {
pairforce.back().getVariableReference("radj"); pairforce.back().getVariableReference("radj");
} catch (std::exception &) { } catch (Lepton::Exception &) {
have_rad.back().second = false; has_ref.back()[2] = false;
} }
} }
} catch (std::exception &e) { } catch (std::exception &e) {
@ -128,9 +134,9 @@ template <int EVFLAG, int EFLAG, int NEWTON_PAIR> void PairLeptonSphere::eval()
if (rsq < cutsq[itype][jtype]) { if (rsq < cutsq[itype][jtype]) {
const double r = sqrt(rsq); const double r = sqrt(rsq);
const int idx = type2expression[itype][jtype]; const int idx = type2expression[itype][jtype];
pairforce[idx].getVariableReference("r") = r; if (has_ref[idx][0]) pairforce[idx].getVariableReference("r") = r;
if (have_rad[idx].first) pairforce[idx].getVariableReference("radi") = radius[i]; if (has_ref[idx][1]) pairforce[idx].getVariableReference("radi") = radius[i];
if (have_rad[idx].second) pairforce[idx].getVariableReference("radj") = radius[j]; if (has_ref[idx][2]) pairforce[idx].getVariableReference("radj") = radius[j];
const double fpair = -pairforce[idx].evaluate() / r * factor_lj; const double fpair = -pairforce[idx].evaluate() / r * factor_lj;
fxtmp += delx * fpair; fxtmp += delx * fpair;
@ -144,9 +150,14 @@ template <int EVFLAG, int EFLAG, int NEWTON_PAIR> void PairLeptonSphere::eval()
double evdwl = 0.0; double evdwl = 0.0;
if (EFLAG) { if (EFLAG) {
pairpot[idx].getVariableReference("r") = r; try {
if (have_rad[idx].first) pairpot[idx].getVariableReference("radi") = radius[i]; pairpot[idx].getVariableReference("r") = r;
if (have_rad[idx].second) pairpot[idx].getVariableReference("radj") = radius[j]; } catch (Lepton::Exception &) {
; // ignore -> constant potential
}
if (has_ref[idx][1]) pairpot[idx].getVariableReference("radi") = radius[i];
if (has_ref[idx][2]) pairpot[idx].getVariableReference("radj") = radius[j];
evdwl = pairpot[idx].evaluate(); evdwl = pairpot[idx].evaluate();
evdwl *= factor_lj; evdwl *= factor_lj;
} }
@ -211,19 +222,23 @@ double PairLeptonSphere::single(int i, int j, int itype, int jtype, double rsq,
auto pairforce = parsed.differentiate("r").createCompiledExpression(); auto pairforce = parsed.differentiate("r").createCompiledExpression();
const double r = sqrt(rsq); const double r = sqrt(rsq);
pairpot.getVariableReference("r") = r; try {
pairforce.getVariableReference("r") = r; pairpot.getVariableReference("r") = r;
pairforce.getVariableReference("r") = r;
} catch (Lepton::Exception &) {
; // ignore -> constant potential or force
}
try { try {
pairpot.getVariableReference("radi") = atom->radius[i]; pairpot.getVariableReference("radi") = atom->radius[i];
pairforce.getVariableReference("radi") = atom->radius[i]; pairforce.getVariableReference("radi") = atom->radius[i];
} catch (std::exception &) { } catch (Lepton::Exception &) {
/* ignore */ ; // ignore
} }
try { try {
pairpot.getVariableReference("radj") = atom->radius[j]; pairpot.getVariableReference("radj") = atom->radius[j];
pairforce.getVariableReference("radj") = atom->radius[j]; pairforce.getVariableReference("radj") = atom->radius[j];
} catch (std::exception &) { } catch (Lepton::Exception &) {
/* ignore */ ; // ignore
} }
fforce = -pairforce.evaluate() / r * factor_lj; fforce = -pairforce.evaluate() / r * factor_lj;

16
src/OPENMP/angle_lepton_omp.cpp
View File

@ -91,10 +91,17 @@ void AngleLeptonOMP::eval(int nfrom, int nto, ThrData *const thr)
{ {
std::vector<Lepton::CompiledExpression> angleforce; std::vector<Lepton::CompiledExpression> angleforce;
std::vector<Lepton::CompiledExpression> anglepot; std::vector<Lepton::CompiledExpression> anglepot;
std::vector<bool> has_ref;
try { try {
for (const auto &expr : expressions) { for (const auto &expr : expressions) {
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, Pointers::lmp)); auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, Pointers::lmp));
angleforce.emplace_back(parsed.differentiate("theta").createCompiledExpression()); angleforce.emplace_back(parsed.differentiate("theta").createCompiledExpression());
has_ref.push_back(true);
try {
angleforce.back().getVariableReference("theta");
} catch (Lepton::Exception &) {
has_ref.back() = false;
}
if (EFLAG) anglepot.emplace_back(parsed.createCompiledExpression()); if (EFLAG) anglepot.emplace_back(parsed.createCompiledExpression());
} }
} catch (std::exception &e) { } catch (std::exception &e) {
@ -146,8 +153,7 @@ void AngleLeptonOMP::eval(int nfrom, int nto, ThrData *const thr)
const double dtheta = acos(c) - theta0[type]; const double dtheta = acos(c) - theta0[type];
const int idx = type2expression[type]; const int idx = type2expression[type];
angleforce[idx].getVariableReference("theta") = dtheta; if (has_ref[idx]) angleforce[idx].getVariableReference("theta") = dtheta;
const double a = -angleforce[idx].evaluate() * s; const double a = -angleforce[idx].evaluate() * s;
const double a11 = a * c / rsq1; const double a11 = a * c / rsq1;
const double a12 = -a / (r1 * r2); const double a12 = -a / (r1 * r2);
@ -183,7 +189,11 @@ void AngleLeptonOMP::eval(int nfrom, int nto, ThrData *const thr)
double eangle = 0.0; double eangle = 0.0;
if (EFLAG) { if (EFLAG) {
anglepot[idx].getVariableReference("theta") = dtheta; try {
anglepot[idx].getVariableReference("theta") = dtheta;
} catch (Lepton::Exception &) {
; // ignore -> constant force
}
eangle = anglepot[idx].evaluate() - offset[type]; eangle = anglepot[idx].evaluate() - offset[type];
} }
if (EVFLAG) if (EVFLAG)

15
src/OPENMP/bond_lepton_omp.cpp
View File

@ -89,10 +89,17 @@ void BondLeptonOMP::eval(int nfrom, int nto, ThrData *const thr)
{ {
std::vector<Lepton::CompiledExpression> bondforce; std::vector<Lepton::CompiledExpression> bondforce;
std::vector<Lepton::CompiledExpression> bondpot; std::vector<Lepton::CompiledExpression> bondpot;
std::vector<bool> has_ref;
try { try {
for (const auto &expr : expressions) { for (const auto &expr : expressions) {
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, Pointers::lmp)); auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, Pointers::lmp));
bondforce.emplace_back(parsed.differentiate("r").createCompiledExpression()); bondforce.emplace_back(parsed.differentiate("r").createCompiledExpression());
has_ref.push_back(true);
try {
bondforce.back().getVariableReference("r");
} catch (Lepton::Exception &) {
has_ref.back() = false;
}
if (EFLAG) bondpot.emplace_back(parsed.createCompiledExpression()); if (EFLAG) bondpot.emplace_back(parsed.createCompiledExpression());
} }
} catch (std::exception &e) { } catch (std::exception &e) {
@ -122,7 +129,7 @@ void BondLeptonOMP::eval(int nfrom, int nto, ThrData *const thr)
double fbond = 0.0; double fbond = 0.0;
if (r > 0.0) { if (r > 0.0) {
bondforce[idx].getVariableReference("r") = dr; if (has_ref[idx]) bondforce[idx].getVariableReference("r") = dr;
fbond = -bondforce[idx].evaluate() / r; fbond = -bondforce[idx].evaluate() / r;
} }
@ -142,7 +149,11 @@ void BondLeptonOMP::eval(int nfrom, int nto, ThrData *const thr)
double ebond = 0.0; double ebond = 0.0;
if (EFLAG) { if (EFLAG) {
bondpot[idx].getVariableReference("r") = dr; try {
bondpot[idx].getVariableReference("r") = dr;
} catch (Lepton::Exception &) {
; // ignore -> constant potential
}
ebond = bondpot[idx].evaluate() - offset[type]; ebond = bondpot[idx].evaluate() - offset[type];
} }
if (EVFLAG) if (EVFLAG)

19
src/OPENMP/dihedral_lepton_omp.cpp
View File

@ -19,9 +19,9 @@
#include "atom.h" #include "atom.h"
#include "comm.h" #include "comm.h"
#include "force.h" #include "force.h"
#include "math_extra.h"
#include "neighbor.h" #include "neighbor.h"
#include "suffix.h" #include "suffix.h"
#include "math_extra.h"
#include <cmath> #include <cmath>
@ -94,10 +94,17 @@ void DihedralLeptonOMP::eval(int nfrom, int nto, ThrData *const thr)
{ {
std::vector<Lepton::CompiledExpression> dihedralforce; std::vector<Lepton::CompiledExpression> dihedralforce;
std::vector<Lepton::CompiledExpression> dihedralpot; std::vector<Lepton::CompiledExpression> dihedralpot;
std::vector<bool> has_ref;
try { try {
for (const auto &expr : expressions) { for (const auto &expr : expressions) {
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, Pointers::lmp)); auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, Pointers::lmp));
dihedralforce.emplace_back(parsed.differentiate("phi").createCompiledExpression()); dihedralforce.emplace_back(parsed.differentiate("phi").createCompiledExpression());
has_ref.push_back(true);
try {
dihedralforce.back().getVariableReference("phi");
} catch (Lepton::Exception &) {
has_ref.back() = false;
}
if (EFLAG) dihedralpot.emplace_back(parsed.createCompiledExpression()); if (EFLAG) dihedralpot.emplace_back(parsed.createCompiledExpression());
} }
} catch (std::exception &e) { } catch (std::exception &e) {
@ -106,7 +113,7 @@ void DihedralLeptonOMP::eval(int nfrom, int nto, ThrData *const thr)
const double *const *const x = atom->x; const double *const *const x = atom->x;
auto *_noalias const f = (dbl3_t *) thr->get_f()[0]; auto *_noalias const f = (dbl3_t *) thr->get_f()[0];
const int * const * const dihedrallist = neighbor->dihedrallist; const int *const *const dihedrallist = neighbor->dihedrallist;
const int nlocal = atom->nlocal; const int nlocal = atom->nlocal;
// The dihedral angle "phi" is the angle between n123 and n234 // The dihedral angle "phi" is the angle between n123 and n234
@ -279,7 +286,7 @@ void DihedralLeptonOMP::eval(int nfrom, int nto, ThrData *const thr)
} }
const int idx = type2expression[type]; const int idx = type2expression[type];
dihedralforce[idx].getVariableReference("phi") = phi; if (has_ref[idx]) dihedralforce[idx].getVariableReference("phi") = phi;
double m_du_dphi = -dihedralforce[idx].evaluate(); double m_du_dphi = -dihedralforce[idx].evaluate();
// ----- Step 4: Calculate the force direction in real space ----- // ----- Step 4: Calculate the force direction in real space -----
@ -323,7 +330,11 @@ void DihedralLeptonOMP::eval(int nfrom, int nto, ThrData *const thr)
double edihedral = 0.0; double edihedral = 0.0;
if (EFLAG) { if (EFLAG) {
dihedralpot[idx].getVariableReference("phi") = phi; try {
dihedralpot[idx].getVariableReference("phi") = phi;
} catch (Lepton::Exception &) {
; // ignore -> constant potential
}
edihedral = dihedralpot[idx].evaluate(); edihedral = dihedralpot[idx].evaluate();
} }
if (EVFLAG) if (EVFLAG)

42
src/OPENMP/pair_lepton_coul_omp.cpp
View File

@ -20,11 +20,13 @@
#include "neigh_list.h" #include "neigh_list.h"
#include "suffix.h" #include "suffix.h"
#include <cmath>
#include "Lepton.h" #include "Lepton.h"
#include "lepton_utils.h" #include "lepton_utils.h"
#include "omp_compat.h" #include "omp_compat.h"
#include <array>
#include <cmath>
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -101,25 +103,30 @@ void PairLeptonCoulOMP::eval(int iifrom, int iito, ThrData *const thr)
std::vector<Lepton::CompiledExpression> pairforce; std::vector<Lepton::CompiledExpression> pairforce;
std::vector<Lepton::CompiledExpression> pairpot; std::vector<Lepton::CompiledExpression> pairpot;
std::vector<std::pair<bool, bool>> have_q; std::vector<std::array<bool, 3>> has_ref;
try { try {
for (const auto &expr : expressions) { for (const auto &expr : expressions) {
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, Pointers::lmp), functions); auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, Pointers::lmp), functions);
pairforce.emplace_back(parsed.differentiate("r").createCompiledExpression()); pairforce.emplace_back(parsed.differentiate("r").createCompiledExpression());
has_ref.push_back({true, true, true});
try {
pairforce.back().getVariableReference("r");
} catch (Lepton::Exception &) {
has_ref.back()[0] = false;
}
if (EFLAG) pairpot.emplace_back(parsed.createCompiledExpression()); if (EFLAG) pairpot.emplace_back(parsed.createCompiledExpression());
pairforce.back().getVariableReference("r");
have_q.emplace_back(true, true);
// check if there are references to charges // check if there are references to charges
try { try {
pairforce.back().getVariableReference("qi"); pairforce.back().getVariableReference("qi");
} catch (std::exception &) { } catch (Lepton::Exception &) {
have_q.back().first = false; has_ref.back()[1] = false;
} }
try { try {
pairforce.back().getVariableReference("qj"); pairforce.back().getVariableReference("qj");
} catch (std::exception &) { } catch (Lepton::Exception &) {
have_q.back().second = false; has_ref.back()[2] = false;
} }
} }
} catch (std::exception &e) { } catch (std::exception &e) {
@ -152,9 +159,9 @@ void PairLeptonCoulOMP::eval(int iifrom, int iito, ThrData *const thr)
if (rsq < cutsq[itype][jtype]) { if (rsq < cutsq[itype][jtype]) {
const double r = sqrt(rsq); const double r = sqrt(rsq);
const int idx = type2expression[itype][jtype]; const int idx = type2expression[itype][jtype];
pairforce[idx].getVariableReference("r") = r; if (has_ref[idx][0]) pairforce[idx].getVariableReference("r") = r;
if (have_q[idx].first) pairforce[idx].getVariableReference("qi") = q2e * q[i]; if (has_ref[idx][1]) pairforce[idx].getVariableReference("qi") = q2e * q[i];
if (have_q[idx].second) pairforce[idx].getVariableReference("qj") = q2e * q[j]; if (has_ref[idx][2]) pairforce[idx].getVariableReference("qj") = q2e * q[j];
const double fpair = -pairforce[idx].evaluate() / r * factor_coul; const double fpair = -pairforce[idx].evaluate() / r * factor_coul;
fxtmp += delx * fpair; fxtmp += delx * fpair;
@ -168,9 +175,14 @@ void PairLeptonCoulOMP::eval(int iifrom, int iito, ThrData *const thr)
double ecoul = 0.0; double ecoul = 0.0;
if (EFLAG) { if (EFLAG) {
pairpot[idx].getVariableReference("r") = r; try {
if (have_q[idx].first) pairpot[idx].getVariableReference("qi") = q2e * q[i]; pairpot[idx].getVariableReference("r") = r;
if (have_q[idx].second) pairpot[idx].getVariableReference("qj") = q2e * q[j]; } catch (Lepton::Exception &) {
; // ignore -> constant potential
}
if (has_ref[idx][1]) pairpot[idx].getVariableReference("qi") = q2e * q[i];
if (has_ref[idx][2]) pairpot[idx].getVariableReference("qj") = q2e * q[j];
ecoul = pairpot[idx].evaluate(); ecoul = pairpot[idx].evaluate();
ecoul *= factor_coul; ecoul *= factor_coul;
} }

20
src/OPENMP/pair_lepton_omp.cpp
View File

@ -20,11 +20,12 @@
#include "neigh_list.h" #include "neigh_list.h"
#include "suffix.h" #include "suffix.h"
#include <cmath>
#include "Lepton.h" #include "Lepton.h"
#include "lepton_utils.h" #include "lepton_utils.h"
#include "omp_compat.h" #include "omp_compat.h"
#include <array>
#include <cmath>
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -96,10 +97,17 @@ void PairLeptonOMP::eval(int iifrom, int iito, ThrData *const thr)
std::vector<Lepton::CompiledExpression> pairforce; std::vector<Lepton::CompiledExpression> pairforce;
std::vector<Lepton::CompiledExpression> pairpot; std::vector<Lepton::CompiledExpression> pairpot;
std::vector<bool> have_ref;
try { try {
for (const auto &expr : expressions) { for (const auto &expr : expressions) {
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, Pointers::lmp), functions); auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, Pointers::lmp), functions);
pairforce.emplace_back(parsed.differentiate("r").createCompiledExpression()); pairforce.emplace_back(parsed.differentiate("r").createCompiledExpression());
have_ref.push_back(true);
try {
pairforce.back().getVariableReference("r");
} catch (Lepton::Exception &) {
have_ref.back() = false;
}
if (EFLAG) pairpot.emplace_back(parsed.createCompiledExpression()); if (EFLAG) pairpot.emplace_back(parsed.createCompiledExpression());
} }
} catch (std::exception &e) { } catch (std::exception &e) {
@ -132,7 +140,7 @@ void PairLeptonOMP::eval(int iifrom, int iito, ThrData *const thr)
if (rsq < cutsq[itype][jtype]) { if (rsq < cutsq[itype][jtype]) {
const double r = sqrt(rsq); const double r = sqrt(rsq);
const int idx = type2expression[itype][jtype]; const int idx = type2expression[itype][jtype];
pairforce[idx].getVariableReference("r") = r; if (have_ref[idx]) pairforce[idx].getVariableReference("r") = r;
const double fpair = -pairforce[idx].evaluate() / r * factor_lj; const double fpair = -pairforce[idx].evaluate() / r * factor_lj;
fxtmp += delx * fpair; fxtmp += delx * fpair;
@ -146,7 +154,11 @@ void PairLeptonOMP::eval(int iifrom, int iito, ThrData *const thr)
double evdwl = 0.0; double evdwl = 0.0;
if (EFLAG) { if (EFLAG) {
pairpot[idx].getVariableReference("r") = r; try {
pairpot[idx].getVariableReference("r") = r;
} catch (Lepton::Exception &) {
; // ignore -> constant potential
}
evdwl = pairpot[idx].evaluate() - offset[itype][jtype]; evdwl = pairpot[idx].evaluate() - offset[itype][jtype];
evdwl *= factor_lj; evdwl *= factor_lj;
} }

44
src/OPENMP/pair_lepton_sphere_omp.cpp
View File

@ -20,11 +20,13 @@
#include "neigh_list.h" #include "neigh_list.h"
#include "suffix.h" #include "suffix.h"
#include <cmath>
#include "Lepton.h" #include "Lepton.h"
#include "lepton_utils.h" #include "lepton_utils.h"
#include "omp_compat.h" #include "omp_compat.h"
#include <array>
#include <cmath>
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -99,25 +101,30 @@ void PairLeptonSphereOMP::eval(int iifrom, int iito, ThrData *const thr)
std::vector<Lepton::CompiledExpression> pairforce; std::vector<Lepton::CompiledExpression> pairforce;
std::vector<Lepton::CompiledExpression> pairpot; std::vector<Lepton::CompiledExpression> pairpot;
std::vector<std::pair<bool, bool>> have_rad; std::vector<std::array<bool, 3>> has_ref;
try { try {
for (const auto &expr : expressions) { for (const auto &expr : expressions) {
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, Pointers::lmp), functions); auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, Pointers::lmp), functions);
pairforce.emplace_back(parsed.differentiate("r").createCompiledExpression()); pairforce.emplace_back(parsed.differentiate("r").createCompiledExpression());
has_ref.push_back({true, true, true});
try {
pairforce.back().getVariableReference("r");
} catch (Lepton::Exception &) {
has_ref.back()[0] = false;
}
if (EFLAG) pairpot.emplace_back(parsed.createCompiledExpression()); if (EFLAG) pairpot.emplace_back(parsed.createCompiledExpression());
pairforce.back().getVariableReference("r");
have_rad.emplace_back(true, true);
// check if there are references to charges // check if there are references to radii
try { try {
pairforce.back().getVariableReference("radi"); pairforce.back().getVariableReference("radi");
} catch (std::exception &) { } catch (Lepton::Exception &) {
have_rad.back().first = false; has_ref.back()[1] = false;
} }
try { try {
pairforce.back().getVariableReference("radj"); pairforce.back().getVariableReference("radj");
} catch (std::exception &) { } catch (Lepton::Exception &) {
have_rad.back().second = false; has_ref.back()[2] = false;
} }
} }
} catch (std::exception &e) { } catch (std::exception &e) {
@ -150,9 +157,9 @@ void PairLeptonSphereOMP::eval(int iifrom, int iito, ThrData *const thr)
if (rsq < cutsq[itype][jtype]) { if (rsq < cutsq[itype][jtype]) {
const double r = sqrt(rsq); const double r = sqrt(rsq);
const int idx = type2expression[itype][jtype]; const int idx = type2expression[itype][jtype];
pairforce[idx].getVariableReference("r") = r; if (has_ref[idx][0]) pairforce[idx].getVariableReference("r") = r;
if (have_rad[idx].first) pairforce[idx].getVariableReference("radi") = radius[i]; if (has_ref[idx][1]) pairforce[idx].getVariableReference("radi") = radius[i];
if (have_rad[idx].second) pairforce[idx].getVariableReference("radj") = radius[j]; if (has_ref[idx][2]) pairforce[idx].getVariableReference("radj") = radius[j];
const double fpair = -pairforce[idx].evaluate() / r * factor_lj; const double fpair = -pairforce[idx].evaluate() / r * factor_lj;
fxtmp += delx * fpair; fxtmp += delx * fpair;
@ -166,9 +173,14 @@ void PairLeptonSphereOMP::eval(int iifrom, int iito, ThrData *const thr)
double evdwl = 0.0; double evdwl = 0.0;
if (EFLAG) { if (EFLAG) {
pairpot[idx].getVariableReference("r") = r; try {
if (have_rad[idx].first) pairpot[idx].getVariableReference("radi") = radius[i]; pairpot[idx].getVariableReference("r") = r;
if (have_rad[idx].second) pairpot[idx].getVariableReference("radj") = radius[j]; } catch (Lepton::Exception &) {
; // ignore -> constant potential
}
if (has_ref[idx][1]) pairpot[idx].getVariableReference("radi") = radius[i];
if (has_ref[idx][2]) pairpot[idx].getVariableReference("radj") = radius[j];
evdwl = pairpot[idx].evaluate(); evdwl = pairpot[idx].evaluate();
evdwl *= factor_lj; evdwl *= factor_lj;
} }

2
unittest/force-styles/tests/angle-lepton.yaml
View File

@ -9,7 +9,7 @@ prerequisites: ! |
pre_commands: ! "" pre_commands: ! ""
post_commands: ! "" post_commands: ! ""
input_file: in.fourmol input_file: in.fourmol
angle_style: lepton angle_style: lepton auto_offset
angle_coeff: ! | angle_coeff: ! |
1 110.1 "k*theta^2; k=75.0" 1 110.1 "k*theta^2; k=75.0"
2 111.0 "k*theta^2; k=45.0" 2 111.0 "k*theta^2; k=45.0"

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

@ -0,0 +1,88 @@
---
lammps_version: 22 Dec 2022
date_generated: Fri Dec 23 15:10:29 2022
epsilon: 7.5e-13
skip_tests:
prerequisites: ! |
atom full
angle lepton
pre_commands: ! ""
post_commands: ! ""
input_file: in.fourmol
angle_style: lepton no_offset
angle_coeff: ! |
1 110.1 "k*theta^2; k=75.0"
2 111.0 "k*theta^2; k=45.0"
3 120.0 "k*theta^2; k=50.0"
4 108.5 "k*theta^2; k=100.0"
equilibrium: 4 1.9216075064457567 1.9373154697137058 2.0943951023931953 1.8936822384138476
extract: ! |
theta0 1
natoms: 29
init_energy: 41.53081789649104
init_stress: ! |2-
8.9723357320869297e+01 -8.7188643750026529e+01 -2.5347135708427655e+00 9.2043419883119782e+01 -2.8187238090404904e+01 -1.5291148024926793e+00
init_forces: ! |2
1 4.7865489310693540e+01 7.8760925902181516e+00 -3.2694525514709866e+01
2 -1.1124882516177341e+00 -9.0075464203887403e+00 -7.2431691227364459e+00
3 -5.9057050592859328e+00 5.3263619873546261e+01 5.2353380124691469e+01
4 -1.6032230038990633e+01 -2.4560529343731403e+01 1.2891625920422307e+01
5 -4.4802331573497639e+01 -4.8300919461089379e+01 -2.3310767889219324e+01
6 4.7083124388174824e+01 -9.5212933434476312e+00 -3.2526392870546800e+01
7 -1.6208182775476303e+01 1.4458587960739102e+01 -3.5314745459502710e+00
8 -6.5664612141881040e+00 -2.5126850154274202e+01 8.2187944731423329e+01
9 -1.5504395262358301e+01 1.6121044185227817e+01 -4.2007069622477866e-01
10 9.9863759179365275e+00 4.1873540105704549e+01 -6.6085640966037403e+01
11 -2.0441876158908627e+01 -6.5186824168985984e+00 9.0023620309811072e+00
12 -1.0772126658369565e+01 -1.0807367300158219e+01 -9.6049647456797871e+00
13 2.8847886813946291e+00 7.2973241014859198e+00 -1.0414233993842981e-01
14 1.5267407478336393e+01 -9.4754911480231776e+00 -6.6307012925544200e+00
15 1.2402914209534773e+01 -6.2644630791613967e+00 1.8484576795819933e+01
16 3.8927757686508357e-01 1.0690061587911176e+01 6.1542759189377696e+00
17 1.4664194297570785e+00 -1.9971277376602425e+00 1.0776844613215999e+00
18 1.5785371874873322e-01 1.6495665212200166e+00 -6.6944747776990434e+00
19 -1.9328033033421670e+00 -2.4078805870919706e+00 2.8669575541313534e+00
20 1.7749495845934338e+00 7.5831406587195394e-01 3.8275172235676900e+00
21 3.4186149299343742e+00 4.2795410364249484e+00 -1.2789555411020650e+01
22 -6.0875600315279677e+00 -4.1504951869796605e+00 4.5212856070195766e+00
23 2.6689451015935934e+00 -1.2904584944528752e-01 8.2682698040010738e+00
24 -1.3053945393770587e+00 5.0741459325183271e+00 -3.0209518576073018e+00
25 -1.0471133765834284e+00 -3.5082261409793856e+00 5.7374874908501228e-01
26 2.3525079159604871e+00 -1.5659197915389413e+00 2.4472031085222894e+00
27 -2.8720725187343754e-01 2.3577465459557132e+00 -8.0312673032168869e-01
28 -6.2799575211500369e-01 -1.4097313073755862e+00 3.2747938980616453e-02
29 9.1520300398844123e-01 -9.4801523858012704e-01 7.7037879134107223e-01
run_energy: 41.28323739029462
run_stress: ! |2-
8.8236221596506681e+01 -8.6492260623309562e+01 -1.7439609731970940e+00 9.0601855980531312e+01 -2.8735005690484968e+01 -2.6097632235197477e+00
run_forces: ! |2
1 4.7316793853445830e+01 8.2815577813110188e+00 -3.2021703111755464e+01
2 -1.1508196824491330e+00 -9.3814982172707460e+00 -7.5761211707510139e+00
3 -5.1083163691832576e+00 5.2667553294971619e+01 5.1784852458007592e+01
4 -1.6078177452605999e+01 -2.4156048365236213e+01 1.3140924677013103e+01
5 -4.4915734474022280e+01 -4.8095168640411821e+01 -2.3331149037574161e+01
6 4.7077916942842350e+01 -9.5906213020090156e+00 -3.2570331503075487e+01
7 -1.6228599672412471e+01 1.4485102617342370e+01 -3.5441153194985300e+00
8 -6.5097893981550730e+00 -2.5117582302614530e+01 8.2131369512416001e+01
9 -1.5527440970965937e+01 1.6147270375910470e+01 -4.0812004993325646e-01
10 1.0070812216240984e+01 4.1571532807578805e+01 -6.5968810328796337e+01
11 -2.0431584971707451e+01 -6.4817395192247664e+00 8.9879981618991636e+00
12 -1.0884695976714678e+01 -1.1067390190389006e+01 -9.1551242768940568e+00
13 2.8052913970098801e+00 7.1296301666594912e+00 1.3173039168682621e-02
14 1.5254877537873529e+01 -8.9700095533297350e+00 -6.5719846903613162e+00
15 1.2392009100170984e+01 -6.0827695435257292e+00 1.7929674392339596e+01
16 4.7158712437377481e-01 1.0631038523396533e+01 6.0960085687560355e+00
17 1.4458707962589659e+00 -1.9708579331587350e+00 1.0634586790394520e+00
18 1.4201882413835909e-01 1.4265339757773337e+00 -5.7663956896747992e+00
19 -1.6609130686729365e+00 -2.0735307593211125e+00 2.4755525101127143e+00
20 1.5188942445345774e+00 6.4699678354377899e-01 3.2908431795620849e+00
21 3.2242729509516406e+00 4.0079233768386153e+00 -1.2047892238650988e+01
22 -5.7215184687399772e+00 -3.8871624402883409e+00 4.2679223469272234e+00
23 2.4972455177883366e+00 -1.2076093655027398e-01 7.7799698917237645e+00
24 -1.1661978296905471e+00 4.5271404898674854e+00 -2.6925565853370195e+00
25 -9.2712094527152167e-01 -3.1291890525017125e+00 5.1208215565053827e-01
26 2.0933187749620688e+00 -1.3979514373657731e+00 2.1804744296864813e+00
27 -2.6804542538020537e-01 2.1830651328698103e+00 -7.3931790038945400e-01
28 -5.7927072943128310e-01 -1.3052929090347909e+00 2.8365455885795865e-02
29 8.4731615481148848e-01 -8.7777222383501941e-01 7.1095244450365813e-01
...

89
unittest/force-styles/tests/bond-lepton_nooffset.yaml Normal file
View File

@ -0,0 +1,89 @@
---
lammps_version: 21 Nov 2023
date_generated: Thu Jan 18 10:15:41 2024
epsilon: 2.5e-13
skip_tests:
prerequisites: ! |
atom full
bond lepton
pre_commands: ! ""
post_commands: ! ""
input_file: in.fourmol
bond_style: lepton no_offset
bond_coeff: ! |
1 1.5 "k*r^2; k=250.0"
2 1.1 "k2*r^2 + k3*r^3 + k4*r^4; k2=300.0; k3=-100.0; k4=50.0"
3 1.3 "k*r^2; k=350.0"
4 1.2 "k*(r-0.2)^2; k=500.0"
5 1.0 "k*r^2; k=450.0"
equilibrium: 5 1.5 1.1 1.3 1.2 1
extract: ! |
r0 1
natoms: 29
init_energy: 38.295825321689215
init_stress: ! |-
-4.7778964706834920e+01 -9.3066674567350432e+01 3.4789470658440035e+02 -3.0023920169312170e+01 -8.0421418879842847e+01 5.8592449335969732e+01
init_forces: ! |2
1 -5.9149914305071416e+00 -3.7728809612345245e+01 -2.7769433362963369e+01
2 -9.4281609567839944e+00 -7.7586487054273015e+00 1.1096676787527940e+01
3 3.2211742366572125e+01 2.7682361264425523e+01 -7.0109911672970497e+00
4 4.9260777576375503e+00 -1.3809750102765932e+00 3.4951785613141868e+00
5 -1.2606902198593501e+00 -1.9373397933007170e+00 6.4372463095041841e+00
6 -3.8858476307965482e+01 6.8567296300319640e+01 1.9889888806614337e+02
7 7.5297927100028144e+00 -3.8622600737556944e+01 -1.9268793182212875e+02
8 1.3018665172824681e+01 -1.2902789438539877e+01 3.2406676637830003e+00
9 7.4343536239661590e-01 8.0072549738604493e-01 3.2899591078538779e+00
10 6.1558871886113291e+00 -2.2419470219698296e+00 1.0080175092279852e+01
11 -3.7020922615305768e-01 -9.1704102274126453e-01 -1.5046795827370363e+00
12 5.2437190958790678e+00 3.4225915524442998e+00 -2.5523597276998897e+00
13 -1.1277007635800260e+01 4.4610677459696646e+00 2.1195215396108269e-01
14 2.9813926585641828e+00 -6.0667387499775116e-01 7.7317115100728788e+00
15 2.5872825164662799e-01 -9.9415365173790704e+00 -3.5428115826174169e+00
16 5.2775953236493464e+01 -3.1855535724919463e+01 -1.6524229620195118e+02
17 -5.8735858023559175e+01 4.0959855098908882e+01 1.5582804819495431e+02
18 -9.0963607969319646e+00 -4.3343406270234155e+00 -1.7623055551859267e+01
19 1.2597490501067170e+01 8.0591915019111742e+00 1.5261489294231819e+01
20 -3.5011297041352050e+00 -3.7248508748877587e+00 2.3615662576274494e+00
21 -1.5332952658285048e+00 5.9630208068632040e-01 -7.4967230017303281e+00
22 4.2380253233105529e+00 1.0270453290850614e+00 6.6489894421385651e+00
23 -2.7047300574820481e+00 -1.6233474097713818e+00 8.4773355959176278e-01
24 -6.6588083188726532e+00 3.5110922792825918e+00 -6.5625174267043489e+00
25 7.9844426562464141e+00 -1.2853795683286129e+00 6.7123710742192300e+00
26 -1.3256343373737607e+00 -2.2257127109539789e+00 -1.4985364751488087e-01
27 6.6999960289138851e+00 6.3808952243186141e+00 2.0100808779497248e+00
28 -8.8466157439236681e-01 3.8018717064230995e-01 -5.9857060538593476e-01
29 -5.8153344545215182e+00 -6.7610823949609244e+00 -1.4115102725637900e+00
run_energy: 37.78424389351509
run_stress: ! |-
-4.6127506998693484e+01 -9.2129732247211749e+01 3.4548310342284810e+02 -2.9841348469661163e+01 -7.8434962689387717e+01 5.9253167412123155e+01
run_forces: ! |2
1 -5.8451208652159004e+00 -3.7483084455000643e+01 -2.7706576989352534e+01
2 -9.4646964278974774e+00 -7.8058897724822449e+00 1.1098831256058579e+01
3 3.1827086102630346e+01 2.7573911030624821e+01 -6.9576662575837211e+00
4 5.1502169659901655e+00 -1.4367546726785101e+00 3.6631301025186187e+00
5 -1.2208420775139264e+00 -1.8781699435112362e+00 6.2332639085051911e+00
6 -3.8491523409043303e+01 6.8063273218541468e+01 1.9723141045830272e+02
7 7.4838209349394775e+00 -3.8394258853636330e+01 -1.9092625515909930e+02
8 1.2676329319901857e+01 -1.2475162287097550e+01 3.3659783337736577e+00
9 6.8845241565874460e-01 7.3814593866184031e-01 3.0434095400342533e+00
10 6.2545583994797553e+00 -2.9600470917047201e+00 9.4247125735981765e+00
11 -1.9554747834212524e-01 -4.8434314068172696e-01 -7.9452259566032057e-01
12 5.2092795750960841e+00 3.1431929551776721e+00 -3.1346654851373348e+00
13 -1.1496483840617872e+01 4.5245217971580018e+00 2.1348220240918236e-01
14 3.1913399826660909e+00 -6.3760720126489068e-01 8.2740980433927742e+00
15 2.7338564489784484e-01 -9.7206665011069671e+00 -3.4841809697094543e+00
16 5.2461611410918316e+01 -3.1639255494702798e+01 -1.6483607587596811e+02
17 -5.8501866653548078e+01 4.0872194473703807e+01 1.5529162691391761e+02
18 -7.0990354207248405e+00 -2.4743922643289666e+00 -1.7824398936159682e+01
19 1.2019842510974870e+01 7.7105128268768715e+00 1.4523712108141252e+01
20 -4.9208070902500296e+00 -5.2361205625479048e+00 3.3006868280184283e+00
21 -1.8548628650934149e+00 2.7467524264262122e-01 -6.7601469408617412e+00
22 3.9136757840663186e+00 9.5561415744904055e-01 6.1181929861632272e+00
23 -2.0588129189729036e+00 -1.2302894000916618e+00 6.4195395469851357e-01
24 -5.7681973234153086e+00 2.0209144998436366e+00 -5.2864044021513967e+00
25 6.3696975292216704e+00 -1.0109756418053095e+00 5.3564043759405795e+00
26 -6.0150020580636188e-01 -1.0099388580383271e+00 -6.9999973789182365e-02
27 6.8467535469188450e+00 5.7500299184200578e+00 2.2775780974490298e+00
28 -1.3929430925479587e+00 5.9772788540443345e-01 -9.4056106886485980e-01
29 -5.4538104543708865e+00 -6.3477578038244911e+00 -1.3370170285841700e+00
...

83
unittest/force-styles/tests/mol-pair-lepton.yaml
View File

@ -1,6 +1,6 @@
--- ---
lammps_version: 22 Dec 2022 lammps_version: 21 Nov 2023
date_generated: Thu Dec 22 09:57:30 2022 date_generated: Thu Jan 18 11:01:50 2024
epsilon: 5e-14 epsilon: 5e-14
skip_tests: intel skip_tests: intel
prerequisites: ! | prerequisites: ! |
@ -23,23 +23,24 @@ pair_coeff: ! |
2 4 "4.0*eps*((sig/r)^12-(sig/r)^6);eps=0.005;sig=0.5" 2 4 "4.0*eps*((sig/r)^12-(sig/r)^6);eps=0.005;sig=0.5"
2 5 "4.0*eps*((sig/r)^12-(sig/r)^6);eps=0.00866025;sig=2.05" 2 5 "4.0*eps*((sig/r)^12-(sig/r)^6);eps=0.00866025;sig=2.05"
3 3 "4.0*eps*((sig/r)^12-(sig/r)^6);eps=0.02;sig=3.2" 3 3 "4.0*eps*((sig/r)^12-(sig/r)^6);eps=0.02;sig=3.2"
3 4 "4.0*eps*((sig/r)^12-(sig/r)^6);eps=0.0173205;sig=3.15" 3 4 "-eps*r;eps=0.0173205;sig=3.15"
3 5 "4.0*eps*((sig/r)^12-(sig/r)^6);eps=0.0173205;sig=3.15" 3 5 "4.0*eps*((sig/r)^12-(sig/r)^6);eps=0.0173205;sig=3.15"
4 4 "10.0"
extract: ! "" extract: ! ""
natoms: 29 natoms: 29
init_vdwl: 749.2468149791969 init_vdwl: 746.1575578155301
init_coul: 0 init_coul: 0
init_stress: ! |2- init_stress: ! |2-
2.1793853434038242e+03 2.1988955172192768e+03 4.6653977523326257e+03 -7.5956547636050584e+02 2.4751536734032861e+01 6.6652028436400667e+02 2.1723526811665593e+03 2.1959162890293533e+03 4.6328064825512138e+03 -7.5509180369489252e+02 9.4506578600439983e+00 6.7585028859193505e+02
init_forces: ! |2 init_forces: ! |2
1 -2.3333390280895912e+01 2.6994567613322641e+02 3.3272827850356805e+02 1 -2.3359983837422618e+01 2.6996030011590727e+02 3.3274783233743295e+02
2 1.5828554630414899e+02 1.3025008843535872e+02 -1.8629682358935722e+02 2 1.5828554630414899e+02 1.3025008843535872e+02 -1.8629682358935722e+02
3 -1.3528903738169066e+02 -3.8704313358319990e+02 -1.4568978437133106e+02 3 -1.3528903738169066e+02 -3.8704313358319990e+02 -1.4568978437133106e+02
4 -7.8711096705893366e+00 2.1350518625373538e+00 -5.5954532185548134e+00 4 -7.8711096705893366e+00 2.1350518625373538e+00 -5.5954532185548134e+00
5 -2.5176757268228540e+00 -4.0521510681020239e+00 1.2152704057877019e+01 5 -2.5176757268228540e+00 -4.0521510681020239e+00 1.2152704057877019e+01
6 -8.3190662465252137e+02 9.6394149462625603e+02 1.1509093566509248e+03 6 -8.3190662465252137e+02 9.6394149462625603e+02 1.1509093566509248e+03
7 5.8203388932513583e+01 -3.3608997951626793e+02 -1.7179617996573040e+03 7 6.6340523101244187e+01 -3.4078810185436379e+02 -1.7003039516942540e+03
8 1.4451392284291535e+02 -1.0927475861088995e+02 3.9990593492420442e+02 8 1.3674478037618434e+02 -1.0517874373121482e+02 3.8291074246191346e+02
9 7.9156945283097443e+01 8.5273009783986538e+01 3.5032175698445189e+02 9 7.9156945283097443e+01 8.5273009783986538e+01 3.5032175698445189e+02
10 5.3118875219105360e+02 -6.1040990859419412e+02 -1.8355872642619292e+02 10 5.3118875219105360e+02 -6.1040990859419412e+02 -1.8355872642619292e+02
11 -2.3530157267965532e+00 -5.9077640073819717e+00 -9.6590723955414290e+00 11 -2.3530157267965532e+00 -5.9077640073819717e+00 -9.6590723955414290e+00
@ -48,8 +49,8 @@ init_forces: ! |2
14 -3.3852721292265153e+00 6.8636181241903649e-01 -8.7507190862499868e+00 14 -3.3852721292265153e+00 6.8636181241903649e-01 -8.7507190862499868e+00
15 -2.0454999188605300e-01 8.4846165523049883e+00 3.0131615419406712e+00 15 -2.0454999188605300e-01 8.4846165523049883e+00 3.0131615419406712e+00
16 4.6326310311812108e+02 -3.3087715736498188e+02 -1.1893024561782554e+03 16 4.6326310311812108e+02 -3.3087715736498188e+02 -1.1893024561782554e+03
17 -4.5334300923766727e+02 3.1554283255882569e+02 1.2058417793481203e+03 17 -4.5371128972368928e+02 3.1609940794953951e+02 1.2052011419527653e+03
18 -1.8862623280672661e-02 -3.3402010907951661e-02 3.1000479299095260e-02 18 8.0197172683943874e-03 -2.4939258820032362e-03 -1.0571459969936936e-02
19 3.1843079640570047e-04 -2.3918627818763426e-04 1.7427252638513439e-03 19 3.1843079640570047e-04 -2.3918627818763426e-04 1.7427252638513439e-03
20 -9.9760831209706009e-04 -1.0209184826753090e-03 3.6910972636601454e-04 20 -9.9760831209706009e-04 -1.0209184826753090e-03 3.6910972636601454e-04
21 -7.1566125273265186e+01 -8.1615678329920655e+01 2.2589561408339890e+02 21 -7.1566125273265186e+01 -8.1615678329920655e+01 2.2589561408339890e+02
@ -61,38 +62,38 @@ init_forces: ! |2
27 5.1810388677546001e+01 -2.2705458321213797e+02 9.0849111082069669e+01 27 5.1810388677546001e+01 -2.2705458321213797e+02 9.0849111082069669e+01
28 -1.8041307121444069e+02 7.7534042932772905e+01 -1.2206956760706598e+02 28 -1.8041307121444069e+02 7.7534042932772905e+01 -1.2206956760706598e+02
29 1.2861057254925012e+02 1.4952711274394568e+02 3.1216025556267880e+01 29 1.2861057254925012e+02 1.4952711274394568e+02 3.1216025556267880e+01
run_vdwl: 719.4530651193046 run_vdwl: 716.5213000416621
run_coul: 0 run_coul: 0
run_stress: ! |2- run_stress: ! |2-
2.1330153957371017e+03 2.1547728168285516e+03 4.3976497417710125e+03 -7.3873328448298525e+02 4.1743821105370067e+01 6.2788012209191027e+02 2.1263870112744726e+03 2.1520080341389726e+03 4.3663519512361027e+03 -7.3456213833770062e+02 2.6927285459244832e+01 6.3691834104928068e+02
run_forces: ! |2 run_forces: ! |2
1 -2.0299419751359164e+01 2.6686193378823020e+02 3.2358785870694015e+02 1 -2.0326040164905073e+01 2.6687684422507328e+02 3.2360752654223910e+02
2 1.5298617928491225e+02 1.2596516341409203e+02 -1.7961292655338619e+02 2 1.5298608857690186e+02 1.2596506573447739e+02 -1.7961281277841888e+02
3 -1.3353630652439830e+02 -3.7923748696131315e+02 -1.4291839793625817e+02 3 -1.3353631293077220e+02 -3.7923732277833739e+02 -1.4291833260989750e+02
4 -7.8374717836161762e+00 2.1276610789823409e+00 -5.5845014473820616e+00 4 -7.8374717116975035e+00 2.1276610267113969e+00 -5.5845014524498486e+00
5 -2.5014258630866735e+00 -4.0250131424704412e+00 1.2103512372025639e+01 5 -2.5014258756924157e+00 -4.0250131713717776e+00 1.2103512280982228e+01
6 -8.0681462887292457e+02 9.2165637136761688e+02 1.0270795806932783e+03 6 -8.0714971444536457e+02 9.2203068890526424e+02 1.0274502514782534e+03
7 5.5780279349903523e+01 -3.1117530951561656e+02 -1.5746991292869018e+03 7 6.3722543724608350e+01 -3.1586173092061807e+02 -1.5580743968587681e+03
8 1.3452983055535049e+02 -1.0064659350255911e+02 3.8851791558207651e+02 8 1.2737293861904031e+02 -9.6945064279519002e+01 3.7231518354375891e+02
9 7.6746213883425980e+01 8.2501469877402130e+01 3.3944351200617882e+02 9 7.6709940036396304e+01 8.2451980339096536e+01 3.3926849385746954e+02
10 5.2128033527695595e+02 -5.9920098848285863e+02 -1.8126029815043339e+02 10 5.2123408713149831e+02 -5.9914309504622599e+02 -1.8121478407355445e+02
11 -2.3573118090915246e+00 -5.8616944550888359e+00 -9.6049808811326205e+00 11 -2.3573086824741427e+00 -5.8616969504300931e+00 -9.6049799947287671e+00
12 1.7503975847822900e+01 1.0626930310560814e+01 -8.0603160272054968e+00 12 1.7504108236707797e+01 1.0626901299509713e+01 -8.0602444903747301e+00
13 8.0530313322973104e+00 -3.1756495170399117e+00 -1.4618315664740528e-01 13 8.0530313558451159e+00 -3.1756495145404533e+00 -1.4618321144421534e-01
14 -3.3416065168069773e+00 6.6492606336082150e-01 -8.6345131440469700e+00 14 -3.3416062225209915e+00 6.6492609500227240e-01 -8.6345136470911594e+00
15 -2.2253843262374914e-01 8.5025661635348779e+00 3.0369735873081622e+00 15 -2.2253820242887132e-01 8.5025660110994483e+00 3.0369741645942137e+00
16 4.3476311264989465e+02 -3.1171086735551415e+02 -1.1135217194927448e+03 16 4.3476708820318731e+02 -3.1171425443331651e+02 -1.1135289618967258e+03
17 -4.2469846140777133e+02 2.9615411776780593e+02 1.1302573488400669e+03 17 -4.2507048343681140e+02 2.9671384825884064e+02 1.1296230654445915e+03
18 -1.8849981672825908e-02 -3.3371636477421307e-02 3.0986293443778727e-02 18 8.0130752607770750e-03 -2.4895867517657545e-03 -1.0574351684568857e-02
19 3.0940277774414027e-04 -2.4634536455373044e-04 1.7433360008861016e-03 19 3.0939970262803125e-04 -2.4635874092791046e-04 1.7433490521479268e-03
20 -9.8648131277150790e-04 -1.0112587134526946e-03 3.6932948773965417e-04 20 -9.8648319666298735e-04 -1.0112621691758337e-03 3.6933139856766442e-04
21 -7.0490745283106378e+01 -7.9749153581142139e+01 2.2171003384646431e+02 21 -7.0490745298133859e+01 -7.9749153568373742e+01 2.2171003384665224e+02
22 -1.0638717908920071e+02 -2.5949502163177968e+01 -1.6645589526812276e+02 22 -1.0638717908973166e+02 -2.5949502162671845e+01 -1.6645589526807785e+02
23 1.7686797710735027e+02 1.0571018898885514e+02 -5.5243337084099387e+01 23 1.7686797710711278e+02 1.0571018898899243e+02 -5.5243337084327727e+01
24 3.8206017656281375e+01 -2.1022820141992960e+02 1.1260711266189014e+02 24 3.8206017659583978e+01 -2.1022820135505594e+02 1.1260711269986750e+02
25 -1.4918881473530880e+02 2.3762151395876508e+01 -1.2549188139143085e+02 25 -1.4918881473631544e+02 2.3762151403215309e+01 -1.2549188138812220e+02
26 1.1097059498808308e+02 1.8645503634228518e+02 1.2861559677865248e+01 26 1.1097059498835199e+02 1.8645503634383900e+02 1.2861559678659969e+01
27 5.0800844984832125e+01 -2.2296588090685469e+02 8.8607367716323253e+01 27 5.0800844960383969e+01 -2.2296588092255456e+02 8.8607367714616288e+01
28 -1.7694190504288886e+02 7.6029945485182026e+01 -1.1950518150242071e+02 28 -1.7694190504410764e+02 7.6029945484553380e+01 -1.1950518150262033e+02
29 1.2614894925528141e+02 1.4694250820033548e+02 3.0893386672863034e+01 29 1.2614894924957088e+02 1.4694250819500266e+02 3.0893386676150566e+01
... ...

35
unittest/utils/test_lepton.cpp
View File

@ -542,6 +542,41 @@ TEST(Lepton, Optimize)
out.str(""); out.str("");
} }
TEST(Lepton, Exception)
{
Lepton::CompiledExpression function, derivative;
auto parsed = Lepton::Parser::parse("x*x");
function = parsed.createCompiledExpression();
derivative = parsed.differentiate("x").createCompiledExpression();
double x = 1.5;
EXPECT_NO_THROW(function.getVariableReference("x") = x;);
EXPECT_NO_THROW(derivative.getVariableReference("x") = x;);
EXPECT_DOUBLE_EQ(function.evaluate(), 2.25);
EXPECT_DOUBLE_EQ(derivative.evaluate(), 3.0);
parsed = Lepton::Parser::parse("x");
function = parsed.createCompiledExpression();
derivative = parsed.differentiate("x").createCompiledExpression();
x = 2.5;
EXPECT_NO_THROW(function.getVariableReference("x") = x;);
EXPECT_THROW(derivative.getVariableReference("x") = x;, Lepton::Exception);
EXPECT_DOUBLE_EQ(function.evaluate(), 2.5);
EXPECT_DOUBLE_EQ(derivative.evaluate(), 1.0);
parsed = Lepton::Parser::parse("1.0");
function = parsed.createCompiledExpression();
derivative = parsed.differentiate("x").createCompiledExpression();
x = 0.5;
EXPECT_THROW(function.getVariableReference("x") = x;, Lepton::Exception);
EXPECT_THROW(derivative.getVariableReference("x") = x;, Lepton::Exception);
EXPECT_DOUBLE_EQ(function.evaluate(), 1.0);
EXPECT_DOUBLE_EQ(derivative.evaluate(), 0.0);
}
int main(int argc, char **argv) int main(int argc, char **argv)
{ {
MPI_Init(&argc, &argv); MPI_Init(&argc, &argv);