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

Merge branch 'develop' into collected-small-changes

Browse Source
This commit is contained in:
Axel Kohlmeyer
2022-12-16 18:25:12 -05:00
parent a51f31fa6d 3b267682b5
commit 111faa758d
23 changed files with 933 additions and 7 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
doc/src/Commands_fix.rst
View File

@ -65,7 +65,7 @@ OPT.
* :doc:`drude <fix_drude>`
* :doc:`drude/transform/direct <fix_drude_transform>`
* :doc:`drude/transform/inverse <fix_drude_transform>`
* :doc:`dt/reset <fix_dt_reset>`
* :doc:`dt/reset (k) <fix_dt_reset>`
* :doc:`edpd/source <fix_dpd_source>`
* :doc:`efield <fix_efield>`
* :doc:`ehex <fix_ehex>`
@ -250,7 +250,7 @@ OPT.
* :doc:`tune/kspace <fix_tune_kspace>`
* :doc:`vector <fix_vector>`
* :doc:`viscosity <fix_viscosity>`
* :doc:`viscous <fix_viscous>`
* :doc:`viscous (k) <fix_viscous>`
* :doc:`viscous/sphere <fix_viscous_sphere>`
* :doc:`wall/body/polygon <fix_wall_body_polygon>`
* :doc:`wall/body/polyhedron <fix_wall_body_polyhedron>`

20
doc/src/Tools.rst
View File

@ -57,7 +57,7 @@ Pre-processing tools
* :ref:`msi2lmp <msi>`
* :ref:`polybond <polybond>`
* :ref:`stl_bin2txt <stlconvert>`
* :ref:`tabulate <tabulate>`
Post-processing tools
=====================
@ -1159,13 +1159,27 @@ For illustration purposes below is a part of the Tcl example script.
----------
.. _tabulate:
tabulate tool
--------------
.. versionadded:: TBD
The ``tabulate`` folder contains Python scripts scripts to generate tabulated
potential files for LAMMPS. The bulk of the code is in the ``tabulate`` module
in the ``tabulate.py`` file. Some example files demonstrating its use are
included. See the README file for more information.
----------
.. _vim:
vim tool
------------------
The files in the tools/vim directory are add-ons to the VIM editor
that allow easier editing of LAMMPS input scripts. See the README.txt
The files in the ``tools/vim`` directory are add-ons to the VIM editor
that allow easier editing of LAMMPS input scripts. See the ``README.txt``
file for details.
These files were provided by Gerolf Ziegenhain (gerolf at

4
doc/src/angle_table.rst
View File

@ -59,6 +59,10 @@ format of this file is described below.
----------
Suitable tables for use with this angle style can be created using the
Python code in the ``tools/tabulate`` folder of the LAMMPS source code
distribution.
The format of a tabulated file is as follows (without the
parenthesized comments):

10
doc/src/bond_table.rst
View File

@ -59,6 +59,13 @@ this file is described below.
----------
Suitable tables for use with this bond style can be created by LAMMPS
itself from existing bond styles using the :doc:`bond_write
<bond_write>` command. This can be useful to have a template file for
testing the bond style settings and to build a compatible custom file.
Another option to generate tables is the Python code in the
``tools/tabulate`` folder of the LAMMPS source code distribution.
The format of a tabulated file is as follows (without the
parenthesized comments):
@ -149,7 +156,8 @@ info.
Related commands
""""""""""""""""
:doc:`bond_coeff <bond_coeff>`, :doc:`delete_bonds <delete_bonds>`
:doc:`bond_coeff <bond_coeff>`, :doc:`delete_bonds <delete_bonds>`,
:doc:`bond_write <bond_write>`
Default
"""""""

4
doc/src/dihedral_table.rst
View File

@ -114,6 +114,10 @@ below.
----------
Suitable tables for use with this dihedral style can be created using
the Python code in the ``tools/tabulate`` folder of the LAMMPS source
code distribution.
The format of a tabulated file is as follows (without the
parenthesized comments). It can begin with one or more comment
or blank lines.

6
doc/src/fix_dt_reset.rst
View File

@ -3,6 +3,8 @@
fix dt/reset command
====================
Accelerator Variants: *dt/reset/kk*
Syntax
""""""
@ -86,6 +88,10 @@ allows dump files to be written at intervals specified by simulation
time, rather than by timesteps. Simulation time is in time units;
see the :doc:`units <units>` doc page for details.
----------
.. include:: accel_styles.rst
Restart, fix_modify, output, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""

6
doc/src/fix_viscous.rst
View File

@ -3,6 +3,8 @@
fix viscous command
===================
Accelerator Variants: *viscous/kk*
Syntax
""""""
@ -84,6 +86,10 @@ more easily be used as a thermostat.
----------
.. include:: accel_styles.rst
----------
Restart, fix_modify, output, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""

7
doc/src/pair_table.rst
View File

@ -120,6 +120,13 @@ best effect:
----------
Suitable tables in the correct format for use with these pair styles can
be created by LAMMPS itself using the :doc:`pair_write <pair_write>`
command. In combination with the :doc:`pair style python <pair_python>`
this can be a powerful mechanism to implement and test tables for use
with LAMMPS. Another option to generate tables is the Python code in
the ``tools/tabulate`` folder of the LAMMPS source code distribution.
The format of a tabulated file has an (optional) header followed by a
series of one or more sections, defined as follows (without the
parenthesized comments). The header must start with a `#` character

4
src/KOKKOS/Install.sh
View File

@ -117,6 +117,8 @@ action fix_acks2_reaxff_kokkos.cpp fix_acks2_reaxff.cpp
action fix_acks2_reaxff_kokkos.h fix_acks2_reaxff.h
action fix_deform_kokkos.cpp
action fix_deform_kokkos.h
action fix_dt_reset_kokkos.cpp
action fix_dt_reset_kokkos.h
action fix_enforce2d_kokkos.cpp
action fix_enforce2d_kokkos.h
action fix_eos_table_rx_kokkos.cpp fix_eos_table_rx.cpp
@ -165,6 +167,8 @@ action fix_wall_lj93_kokkos.cpp
action fix_wall_lj93_kokkos.h
action fix_wall_reflect_kokkos.cpp
action fix_wall_reflect_kokkos.h
action fix_viscous_kokkos.cpp
action fix_viscous_kokkos.h
action fix_dpd_energy_kokkos.cpp fix_dpd_energy.cpp
action fix_dpd_energy_kokkos.h fix_dpd_energy.h
action fix_rx_kokkos.cpp fix_rx.cpp

195
src/KOKKOS/fix_dt_reset_kokkos.cpp Normal file
View File

@ -0,0 +1,195 @@
// 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.
------------------------------------------------------------------------- */
#include "fix_dt_reset_kokkos.h"
#include "atom_kokkos.h"
#include "atom_masks.h"
#include "error.h"
#include "force.h"
#include "input.h"
#include "integrate.h"
#include "kokkos_base.h"
#include "memory_kokkos.h"
#include "modify.h"
#include "output.h"
#include "pair.h"
#include "update.h"
using namespace LAMMPS_NS;
using namespace FixConst;
#define BIG 1.0e20
/* ---------------------------------------------------------------------- */
template<class DeviceType>
FixDtResetKokkos<DeviceType>::FixDtResetKokkos(LAMMPS *lmp, int narg, char **arg) :
FixDtReset(lmp, narg, arg)
{
kokkosable = 1;
atomKK = (AtomKokkos *) atom;
execution_space = ExecutionSpaceFromDevice<DeviceType>::space;
datamask_read = EMPTY_MASK;
datamask_modify = EMPTY_MASK;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void FixDtResetKokkos<DeviceType>::init()
{
FixDtReset::init();
k_emax = Kokkos::DualView<double*, Kokkos::LayoutRight, DeviceType>("FixDtResetKokkos:gamma", 1);
k_emax.h_view(0) = emax;
k_emax.template modify<LMPHostType>();
k_emax.template sync<DeviceType>();
if (utils::strmatch(update->integrate_style,"^respa"))
error->all(FLERR,"Cannot (yet) use respa with Kokkos");
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void FixDtResetKokkos<DeviceType>::end_of_step()
{
atomKK->sync(execution_space, V_MASK | F_MASK | MASK_MASK | TYPE_MASK | RMASS_MASK);
v = atomKK->k_v.view<DeviceType>();
f = atomKK->k_f.view<DeviceType>();
mask = atomKK->k_mask.view<DeviceType>();
type = atomKK->k_type.view<DeviceType>();
if (atomKK->rmass)
rmass = atomKK->k_rmass.view<DeviceType>();
else
mass = atomKK->k_mass.view<DeviceType>();
int nlocal = atom->nlocal;
double dt;
copymode = 1;
if (atomKK->rmass)
Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagFixDtResetRMass>(0,nlocal), *this, Kokkos::Min<double>(dt));
else
Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagFixDtResetMass>(0,nlocal), *this, Kokkos::Min<double>(dt));
copymode = 0;
MPI_Allreduce(MPI_IN_PLACE, &dt, 1, MPI_DOUBLE, MPI_MIN, world);
atomKK->modified(execution_space, F_MASK);
if (minbound) dt = MAX(dt, tmin);
if (maxbound) dt = MIN(dt, tmax);
// if timestep didn't change, just return
// else reset update->dt and other classes that depend on it
// rRESPA, pair style, fixes
if (dt == update->dt) return;
laststep = update->ntimestep;
// calls to other classes that need to know timestep size changed
// similar logic is in Input::timestep()
update->update_time();
update->dt = dt;
update->dt_default = 0;
if (force->pair) force->pair->reset_dt();
for (int i = 0; i < modify->nfix; i++) modify->fix[i]->reset_dt();
output->reset_dt();
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void FixDtResetKokkos<DeviceType>::operator()(TagFixDtResetMass, const int &i, double &k_dt) const {
double dtv, dtf, dte, dtsq;
double vsq, fsq, massinv;
double delx, dely, delz, delr;
double emax = k_emax.d_view(0);
if (mask[i] & groupbit) {
massinv = 1.0 / mass[type[i]];
vsq = v(i,0) * v(i,0) + v(i,1) * v(i,1) + v(i,2) * v(i,2);
fsq = f(i,0) * f(i,0) + f(i,1) * f(i,1) + f(i,2) * f(i,2);
dtv = dtf = dte = BIG;
if (vsq > 0.0) dtv = xmax / sqrt(vsq);
if (fsq > 0.0) dtf = sqrt(2.0 * xmax / (ftm2v * sqrt(fsq) * massinv));
k_dt = MIN(dtv, dtf);
if ((emax > 0.0) && (fsq * vsq > 0.0)) {
dte = emax / sqrt(fsq * vsq) / sqrt(ftm2v * mvv2e);
k_dt = MIN(dt, dte);
}
dtsq = k_dt * k_dt;
delx = k_dt * v(i,0) + 0.5 * dtsq * massinv * f(i,0) * ftm2v;
dely = k_dt * v(i,1) + 0.5 * dtsq * massinv * f(i,1) * ftm2v;
delz = k_dt * v(i,2) + 0.5 * dtsq * massinv * f(i,2) * ftm2v;
delr = sqrt(delx * delx + dely * dely + delz * delz);
if (delr > xmax) k_dt *= xmax / delr;
}
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void FixDtResetKokkos<DeviceType>::operator()(TagFixDtResetRMass, const int &i, double &k_dt) const {
double dtv, dtf, dte, dtsq;
double vsq, fsq, massinv;
double delx, dely, delz, delr;
double emax = k_emax.d_view(0);
if (mask[i] & groupbit) {
massinv = 1.0 / rmass[i];
vsq = v(i,0) * v(i,0) + v(i,1) * v(i,1) + v(i,2) * v(i,2);
fsq = f(i,0) * f(i,0) + f(i,1) * f(i,1) + f(i,2) * f(i,2);
dtv = dtf = dte = BIG;
if (vsq > 0.0) dtv = xmax / sqrt(vsq);
if (fsq > 0.0) dtf = sqrt(2.0 * xmax / (ftm2v * sqrt(fsq) * massinv));
k_dt = MIN(dtv, dtf);
if ((emax > 0.0) && (fsq * vsq > 0.0)) {
dte = emax / sqrt(fsq * vsq) / sqrt(ftm2v * mvv2e);
k_dt = MIN(dt, dte);
}
dtsq = k_dt * k_dt;
delx = k_dt * v(i,0) + 0.5 * dtsq * massinv * f(i,0) * ftm2v;
dely = k_dt * v(i,1) + 0.5 * dtsq * massinv * f(i,1) * ftm2v;
delz = k_dt * v(i,2) + 0.5 * dtsq * massinv * f(i,2) * ftm2v;
delr = sqrt(delx * delx + dely * dely + delz * delz);
if (delr > xmax) k_dt *= xmax / delr;
}
}
namespace LAMMPS_NS {
template class FixDtResetKokkos<LMPDeviceType>;
#ifdef LMP_KOKKOS_GPU
template class FixDtResetKokkos<LMPHostType>;
#endif
}

66
src/KOKKOS/fix_dt_reset_kokkos.h Normal file
View File

@ -0,0 +1,66 @@
/* -*- 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 FIX_CLASS
// clang-format off
FixStyle(dt/reset/kk,FixDtResetKokkos<LMPDeviceType>);
FixStyle(dt/reset/kk/device,FixDtResetKokkos<LMPDeviceType>);
FixStyle(dt/reset/kk/host,FixDtResetKokkos<LMPHostType>);
// clang-format on
#else
// clang-format off
#ifndef LMP_FIX_DT_RESET_KOKKOS_H
#define LMP_FIX_DT_RESET_KOKKOS_H
#include "fix_dt_reset.h"
#include "kokkos_type.h"
namespace LAMMPS_NS {
struct TagFixDtResetMass{};
struct TagFixDtResetRMass{};
template<class DeviceType>
class FixDtResetKokkos : public FixDtReset {
public:
typedef DeviceType device_type;
typedef ArrayTypes<DeviceType> AT;
FixDtResetKokkos(class LAMMPS *, int, char **);
// ~FixDtResetKokkos() override;
void init() override;
void end_of_step() override;
KOKKOS_INLINE_FUNCTION
void operator()(TagFixDtResetMass, const int&, double&) const;
KOKKOS_INLINE_FUNCTION
void operator()(TagFixDtResetRMass, const int&, double&) const;
private:
typename AT::t_v_array v;
typename AT::t_f_array f;
typename AT::t_int_1d_randomread mask;
typename AT::t_int_1d_randomread type;
typename ArrayTypes<DeviceType>::t_float_1d_randomread rmass;
typename ArrayTypes<DeviceType>::t_float_1d_randomread mass;
Kokkos::DualView<double*, Kokkos::LayoutRight, DeviceType> k_emax;
};
}
#endif
#endif

106
src/KOKKOS/fix_viscous_kokkos.cpp Normal file
View File

@ -0,0 +1,106 @@
// 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.
------------------------------------------------------------------------- */
#include "fix_viscous_kokkos.h"
#include "atom_kokkos.h"
#include "update.h"
#include "modify.h"
#include "input.h"
#include "memory_kokkos.h"
#include "error.h"
#include "atom_masks.h"
#include "kokkos_base.h"
using namespace LAMMPS_NS;
using namespace FixConst;
/* ---------------------------------------------------------------------- */
template<class DeviceType>
FixViscousKokkos<DeviceType>::FixViscousKokkos(LAMMPS *lmp, int narg, char **arg) :
FixViscous(lmp, narg, arg)
{
kokkosable = 1;
atomKK = (AtomKokkos *) atom;
execution_space = ExecutionSpaceFromDevice<DeviceType>::space;
datamask_read = EMPTY_MASK;
datamask_modify = EMPTY_MASK;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
FixViscousKokkos<DeviceType>::~FixViscousKokkos()
{
if (copymode) return;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void FixViscousKokkos<DeviceType>::init()
{
FixViscous::init();
k_gamma = Kokkos::DualView<double*, Kokkos::LayoutRight, DeviceType>("FixViscousKokkos:gamma",atom->ntypes+1);
for (int i = 1; i <= atom->ntypes; i++) k_gamma.h_view(i) = gamma[i];
k_gamma.template modify<LMPHostType>();
k_gamma.template sync<DeviceType>();
if (utils::strmatch(update->integrate_style,"^respa"))
error->all(FLERR,"Cannot (yet) use respa with Kokkos");
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void FixViscousKokkos<DeviceType>::post_force(int /*vflag*/)
{
atomKK->sync(execution_space, V_MASK | F_MASK | MASK_MASK | TYPE_MASK);
v = atomKK->k_v.view<DeviceType>();
f = atomKK->k_f.view<DeviceType>();
mask = atomKK->k_mask.view<DeviceType>();
type = atomKK->k_type.view<DeviceType>();
int nlocal = atom->nlocal;
copymode = 1;
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagFixViscous>(0,nlocal),*this);
copymode = 0;
atomKK->modified(execution_space, F_MASK);
}
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void FixViscousKokkos<DeviceType>::operator()(TagFixViscous, const int &i) const {
if (mask[i] & groupbit) {
double drag = k_gamma.d_view(type[i]);
f(i,0) -= drag*v(i,0);
f(i,1) -= drag*v(i,1);
f(i,2) -= drag*v(i,2);
}
}
namespace LAMMPS_NS {
template class FixViscousKokkos<LMPDeviceType>;
#ifdef LMP_KOKKOS_GPU
template class FixViscousKokkos<LMPHostType>;
#endif
}

60
src/KOKKOS/fix_viscous_kokkos.h Normal file
View File

@ -0,0 +1,60 @@
/* -*- 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 FIX_CLASS
// clang-format off
FixStyle(viscous/kk,FixViscousKokkos<LMPDeviceType>);
FixStyle(viscous/kk/device,FixViscousKokkos<LMPDeviceType>);
FixStyle(viscous/kk/host,FixViscousKokkos<LMPHostType>);
// clang-format on
#else
// clang-format off
#ifndef LMP_FIX_VISCOUS_KOKKOS_H
#define LMP_FIX_VISCOUS_KOKKOS_H
#include "fix_viscous.h"
#include "kokkos_type.h"
namespace LAMMPS_NS {
struct TagFixViscous{};
template<class DeviceType>
class FixViscousKokkos : public FixViscous {
public:
typedef DeviceType device_type;
typedef ArrayTypes<DeviceType> AT;
FixViscousKokkos(class LAMMPS *, int, char **);
~FixViscousKokkos() override;
void init() override;
void post_force(int) override;
KOKKOS_INLINE_FUNCTION
void operator()(TagFixViscous, const int&) const;
private:
typename AT::t_v_array v;
typename AT::t_f_array f;
typename AT::t_int_1d_randomread mask;
typename AT::t_int_1d_randomread type;
Kokkos::DualView<double*, Kokkos::LayoutRight, DeviceType> k_gamma;
};
}
#endif
#endif

2
src/fix_dt_reset.h
View File

@ -34,7 +34,7 @@ class FixDtReset : public Fix {
void end_of_step() override;
double compute_scalar() override;
private:
protected:
bigint laststep;
int minbound, maxbound;
double tmin, tmax, xmax, emax;

2
src/fix_viscous.cpp
View File

@ -61,6 +61,8 @@ FixViscous::FixViscous(LAMMPS *lmp, int narg, char **arg) :
FixViscous::~FixViscous()
{
if (copymode) return;
delete [] gamma;
}

1
tools/README
View File

@ -48,6 +48,7 @@ smd convert Smooth Mach Dynamics triangles to VTK
spin perform a cubic polynomial interpolation of a GNEB MEP
stl_bin2txt convert binary STL files to ASCII
swig Interface file and demo scripts for SWIG wrappers for the LAMMPS C library interface
tabulate Python scripts to generate tabulated potential files for LAMMPS
tinker2lmp.py convert Tinker input files to LAMMPS input files
valgrind suppression files for use with valgrind's memcheck tool
vim add-ons to VIM editor for editing LAMMPS input scripts

39
tools/tabulate/README.md Normal file
View File

@ -0,0 +1,39 @@
# Python scripts to generate tabulated potential files for LAMMPS
This directory contains a Python module 'tabulate' that can be used to
generate tabulated potential files for pair style table, bond style
table, and angle style table
To create tables, you need to define your energy and - optionally -
force functions and then an instance of either the PairTabulate(),
BondTabulate(), AngleTabulate(), or DihedralTabulate() class from the
tabulate module and call its run() method to generate the table. Most
of the settings (number of points, minimum, maximum etc.) are provided
as command line flags. The run() method may be called multiple times to
generate multiple tables, for instance after changing parameters of the
energy/force functions. If the force function is not provided, the
force will be determined through numerical differentiation.
Please see the individual tabulation scripts in this folder for examples:
| ------------------------------|----------------------------------------------------------------------------|
| File | Description |
| ------------------------------|----------------------------------------------------------------------------|
| pair_lj_tabulate.py | creates two Lennard-Jones pair potential tables with different parameters |
| bond_morse_tabulate.py | creates a table for a Morse bond potential table |
| angle_harmonic_tabulate.py | creates a table for a harmonic angle potential table |
| dihedral_harmonic_tabulate.py | creates a table for a harmonic dihedral potential table |
| pair_hybrid_tabulate.py | creates a Morse/Lennard-Jones hybrid potential table with smooth switching |
| ------------------------------|----------------------------------------------------------------------------|
Common command line flags:
```
options:
-h, --help show this help message and exit
--num-points NUM, -n NUM Number of tabulated points (default: 1000)
--filename FILENAME, -f FILENAME Name of output file (default: -)
--diff-num, -d Differentiate energy function numerically
--inner XMIN, -i XMIN Inner cutoff of table (required for pair)
--outer XMAX, -o XMAX Outer cutoff of table (required)
```

21
tools/tabulate/angle_harmonic_tabulate.py Executable file
View File

@ -0,0 +1,21 @@
#!/usr/bin/env python3
from tabulate import AngleTabulate
################################################################################
import math
def harmonic_energy(theta):
k = 100.0
thetazero = 120.0
# the force constant in LAMMPS is in energy per radians^2 so convert from degrees to radians
deg2rad = math.pi / 180.0
t = (theta - thetazero) * deg2rad
f = k * t * t
return f
################################################################################
if __name__ == "__main__":
atable = AngleTabulate(harmonic_energy, units='real')
atable.run('HARM')

28
tools/tabulate/bond_morse_tabulate.py Executable file
View File

@ -0,0 +1,28 @@
#!/usr/bin/env python3
from tabulate import BondTabulate
################################################################################
import math
def morse_energy(r):
depth = 1.0
width = 2.0
rzero = 1.2
ralpha = math.exp(-width*(r-rzero))
f = depth * (-1.0 + (1.0 - ralpha) * (1.0 - ralpha))
return f
def morse_force(r):
depth = 1.0
width = 2.0
rzero = 1.2
ralpha = math.exp(-width*(r-rzero))
f = -2.0 * depth * width * (1.0 -ralpha) * ralpha
return f
################################################################################
if __name__ == "__main__":
btable = BondTabulate(morse_energy, morse_force, units='lj')
btable.run('MORSE')

19
tools/tabulate/dihedral_harmonic_tabulate.py Executable file
View File

@ -0,0 +1,19 @@
#!/usr/bin/env python3
from tabulate import DihedralTabulate
################################################################################
import math
def harmonic_energy(theta):
k = 100.0
# the force constant in LAMMPS is in energy per radians^2 so convert from degrees to radians
deg2rad = math.pi / 180.0
f = k * (1 - math.cos(2.0 * deg2rad * theta))
return f
################################################################################
if __name__ == "__main__":
dtable = DihedralTabulate(harmonic_energy, units='metal')
dtable.run('HARM')

31
tools/tabulate/pair_hybrid_tabulate.py Executable file
View File

@ -0,0 +1,31 @@
#!/usr/bin/env python3
from tabulate import PairTabulate
################################################################################
import math
# hybrid potential using Morse for repulsive and LJ for attractive.
# using tanh() as smooth switching function between the two.
def hybrid_energy(r):
depth = 1.0
width = 5.0
epsilon = 1.0
sigma = 1.0
rzero = 1.2
f1 = 4.0*epsilon*(math.pow(sigma/r,12.0) - math.pow(sigma/r,6.0))
ralpha = math.exp(-width*(r-rzero))
f2 = depth * (-1.0 + (1.0 - ralpha) * (1.0 - ralpha))
switch = 0.5*math.tanh(10.0*(r - rzero)) + 0.5
f = switch*f1 + (1.0 - switch)*f2
return f
################################################################################
if __name__ == "__main__":
ptable = PairTabulate(hybrid_energy, units='lj', comment='Morse repulsion + LJ attraction')
ptable.run('MORSE_LJ')

26
tools/tabulate/pair_lj_tabulate.py Executable file
View File

@ -0,0 +1,26 @@
#!/usr/bin/env python
from tabulate import PairTabulate
################################################################################
import math
epsilon = 1.0
sigma = 1.0
def lj_energy(r):
f = 4.0*epsilon*(math.pow(sigma/r,12.0) - math.pow(sigma/r,6.0))
return f
def lj_force(r):
epsilon = 1.0
sigma = 1.0
f = -4.0*epsilon*(-12.0*math.pow(sigma/r,12.0)/r + 6.0*math.pow(sigma/r,6.0)/r)
return f
################################################################################
if __name__ == "__main__":
ptable = PairTabulate(lj_energy, lj_force)
ptable.run('LJ_11')
epsilon = 1.0
sigma = 1.5
ptable.run('LJ_12')

279
tools/tabulate/tabulate.py Executable file
View File

@ -0,0 +1,279 @@
# ----------------------------------------------------------------------
# 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.
# -------------------------------------------------------------------------
"""
Python classes and support functions to generate tabulated potential files for LAMMPS.
"""
# for python2/3 compatibility
from __future__ import print_function
import sys
import argparse
from os import path
from datetime import datetime
########################################################################
def numdiff(x, func):
""" Get the value of the first derivative of a function 'func(x)' at 'x'
from numerical differentiation."""
# optimal delta x value for 2-point numerical differentiation of two floating point numbers
epsilon = 6.05504e-6
fval1 = func(x - epsilon)
fval2 = func(x + epsilon)
return 0.5 * (fval2-fval1) / epsilon
########################################################################
def mktable(tstyle, label, num, xmin, xmax, efunc, diff=False, ffunc=None):
""" Do the tabulation of the provided energy function. Compute force from
numerical differentiation if no force function is provided. Also detect
minimum for use to determine potential shifting in bonded potentials."""
# must use numerical differentiation if no force function provided
if not ffunc:
diff = True
print("# Creating %s table %s with %d points from %g to %g" % (tstyle, label, num, xmin, xmax))
table = []
delx = (xmax - xmin) / (float(num) - 1.0)
emin = 999999999.0
xzero = 0.0
for i in range(0,num):
x = xmin + i*delx
energy = efunc(x)
if energy < emin:
emin = energy
xzero = x
if diff:
force = -numdiff(x, efunc)
else:
force = ffunc(x)
table.append([i+1, x, energy, force])
return table, xzero
########################################################################
# base class with shared functionality
class Tabulate(object):
"""Base tabulation class. Contains all shared functionality: common argument parsing,
output file handling, table writing"""
def __init__(self, style, efunc, ffunc=None, units=None, comment=None):
self.fp = sys.stdout
self.tstyle = style
self.energyfunc = efunc
self.forcefunc = ffunc
self.units = units
self.comment = comment
self.eshift = 0.0
self.args = None
self.diff = True
self.parser = argparse.ArgumentParser(description='Tool to generate tabulated '
+ self.tstyle + ' potential files for LAMMPS')
self.parser.add_argument('--num-points', '-n', dest='num', default=1000, type=int,
help="Number of tabulated points")
self.parser.add_argument('--filename', '-f', dest='filename', default='-',
help="Name of output file")
self.parser.add_argument('--diff-num', '-d', dest='diff',default=False,
action='store_true',
help="Differentiate energy function numerically")
self.parser.add_argument('--inner', '-i', dest='xmin', required=True, type=float,
help="Inner cutoff of table")
self.parser.add_argument('--outer', '-o', dest='xmax', required=True, type=float,
help="Outer cutoff of table")
def openfile(self, label):
"""Open table file, if needed and print label for new table entry"""
if self.args and self.args.filename != '-':
try:
if path.isfile(self.args.filename):
self.fp = open(self.args.filename, 'a')
print("# Appending table to file " + self.args.filename)
else:
self.fp = open(self.args.filename, 'w')
print("# Writing table to new file " + self.args.filename)
self.fp.write('# DATE: ' + datetime.now().date().isoformat())
if self.units:
self.fp.write(' UNITS: ' + str(self.units))
if self.comment:
self.fp.write(' COMMENT: ' + str(self.comment) + '\n')
except IOError:
sys.exit("Cannot open file %s for writing table data" % self.args.filename)
self.fp.write('\n' + label + '\n')
def writetable(self, table, offset):
""" Formatted output tabulated data with 4 columns"""
for i,r,energy,force in table:
self.fp.write("%8d %- 22.15g %- 22.15g %- 22.15g\n" % (i, r, energy - offset, force))
def helpexit(self, text):
""" Convenience function to exit program with error and help message"""
sys.exit('\n' + text + '\n\n' + self.parser.format_help())
################################################################################
# create tabulation for pair styles
class PairTabulate(Tabulate):
def __init__(self, efunc, ffunc=None, units=None, comment=None):
super(PairTabulate, self).__init__('pair', efunc, ffunc, units, comment)
self.parser.add_argument('--eshift', '-e', dest='eshift', default=False,
action='store_true',
help="Shift potential energy to be zero at outer cutoff")
try:
self.args = self.parser.parse_args()
except argparse.ArgumentError:
sys.exit()
def run(self, label):
# sanity checks
if self.args.num < 2:
self.helpexit('Expect 2 or more points in table for tabulation')
if self.args.xmin <= 0.0:
self.helpexit('Inner tabulation cutoff must be > 0 for pair style table')
if self.args.xmax <= self.args.xmin:
self.helpexit('Outer cutoff must be larger than inner cutoff')
self.diff = self.args.diff
if not self.forcefunc:
self.diff = True
offset = 0.0
if self.args.eshift:
offset=self.energyfunc(self.args.xmax)
table, dummy = mktable(self.tstyle, label, self.args.num, self.args.xmin, self.args.xmax,
self.energyfunc, self.args.diff, self.forcefunc)
# open table file
self.openfile(label)
# write pair style specific header
if self.forcefunc:
diffmin = -numdiff(self.args.xmin, self.forcefunc)
diffmax = -numdiff(self.args.xmax, self.forcefunc)
self.fp.write("N %d R %g %g FPRIME %- 22.15g %- 22.15g\n\n"
% (self.args.num, self.args.xmin, self.args.xmax, diffmin, diffmax))
else:
self.fp.write("N %d R %g %g\n\n" % (self.args.num, self.args.xmin, self.args.xmax))
self.writetable(table, offset)
if self.args.filename != '-':
self.fp.close()
################################################################################
# shared functionality to create tabulation for bond or angle styles
class BondAngleTabulate(Tabulate):
def __init__(self, style, efunc, ffunc=None, units=None, comment=None):
super(BondAngleTabulate, self).__init__(style, efunc, ffunc, units, comment)
self.parser.add_argument('--eshift', '-e', dest='eshift', default=False,
action='store_true',
help="Shift potential energy to be zero at minimum")
idx = [a.dest for a in self.parser._actions].index('xmin')
self.parser._actions[idx].required=False
self.parser._actions[idx].default=0.0
if style == 'angle':
idx = [a.dest for a in self.parser._actions].index('xmax')
self.parser._actions[idx].required=False
self.parser._actions[idx].default=180.0
try:
self.args = self.parser.parse_args()
except argparse.ArgumentError:
sys.exit()
def run(self, label):
# sanity checks
if self.args.num < 2:
self.helpexit('Expect 2 or more points in table for tabulation')
if self.args.xmin < 0.0:
self.helpexit('Inner cutoff must not be negative')
if self.tstyle == 'angle' and self.args.xmax > 180.0:
self.helpexit('Outer cutoff must not be larger than 180.0 degrees')
self.diff = self.args.diff
if not self.forcefunc:
self.diff = True
table, xzero = mktable(self.tstyle, label, self.args.num, self.args.xmin, self.args.xmax,
self.energyfunc, self.args.diff, self.forcefunc)
print("# Minimum energy of tabulated potential is at %g" % xzero)
offset = 0.0
if self.args.eshift:
offset=self.energyfunc(xzero)
self.openfile(label)
if self.forcefunc:
diffmin = -numdiff(self.args.xmin, self.forcefunc)
diffmax = -numdiff(self.args.xmax, self.forcefunc)
self.fp.write("N %d FP %- 22.15g %- 22.15g EQ %g\n\n" %
(self.args.num, diffmin, diffmax, xzero))
else:
self.fp.write("N %d EQ %g\n\n" % (self.args.num, xzero))
self.writetable(table, offset)
if self.args.filename != '-':
self.fp.close()
################################################################################
class BondTabulate(BondAngleTabulate):
def __init__(self, efunc, ffunc=None, units=None, comment=None):
super(BondTabulate, self).__init__('bond', efunc, ffunc, units, comment)
################################################################################
class AngleTabulate(BondAngleTabulate):
def __init__(self, efunc, ffunc=None, units=None, comment=None):
super(AngleTabulate, self).__init__('angle', efunc, ffunc, units, comment)
################################################################################
# create tabulation for dihdedral
class DihedralTabulate(Tabulate):
def __init__(self, efunc, ffunc=None, units=None, comment=None):
super(DihedralTabulate, self).__init__('dihedral', efunc, ffunc, units, comment)
idx = [a.dest for a in self.parser._actions].index('xmin')
self.parser._actions[idx].required=False
self.parser._actions[idx].default=-180.0
idx = [a.dest for a in self.parser._actions].index('xmax')
self.parser._actions[idx].required=False
self.parser._actions[idx].default=179.999999
try:
self.args = self.parser.parse_args()
except argparse.ArgumentError:
sys.exit()
def run(self, label):
# sanity checks
if self.args.num < 2:
self.helpexit('Expect 2 or more points in table for tabulation')
if self.args.xmin < -180 or self.args.xmin > 0.0:
self.helpexit('Inner cutoff must be within -180.0 and 0.0 degrees')
if self.args.xmax < 0.0 or self.args.xmin > 360.0:
self.helpexit('Outer cutoff must be within 0.0 and 360.0 degrees')
if (self.args.xmax - self.args.xmin) >= 360.0:
self.helpexit('Inner and outer cutoff range must be less than 360.0 degrees')
self.diff = self.args.diff
if not self.forcefunc:
self.diff = True
table, dummy = mktable(self.tstyle, label, self.args.num, self.args.xmin, self.args.xmax,
self.energyfunc, self.args.diff, self.forcefunc)
self.openfile(label)
self.fp.write("N %d DEGREES \n\n" % (self.args.num))
self.writetable(table, 0.0)
if self.args.filename != '-':
self.fp.close()