tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'develop' of github.com:lammps/lammps into hagerty_issue4140_kokkos_collective

Browse Source
This commit is contained in:
Stan Moore
2024-07-01 12:33:17 -06:00
parent 3defe567df a74500f416
commit abcc9ed08d
408 changed files with 69505 additions and 9307 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
.github/CODEOWNERS vendored
View File

@ -38,6 +38,7 @@ src/ML-HDNNP/* @singraber
src/ML-IAP/* @athomps
src/ML-PACE/* @yury-lysogorskiy
src/ML-POD/* @exapde
src/ML-UF3/* @monk-04
src/MOFFF/* @hheenen
src/MOLFILE/* @akohlmey
src/NETCDF/* @pastewka
@ -58,7 +59,8 @@ src/VTK/* @rbberger
# individual files in packages
src/GPU/pair_vashishta_gpu.* @andeplane
src/KOKKOS/pair_vashishta_kokkos.* @andeplane
src/KOKKOS/pair_vashishta_kokkos.* @andeplane @stanmoore1
src/KOSSOS/pair_pod_kokkos.* @exapde @stanmoore1
src/MANYBODY/pair_vashishta_table.* @andeplane
src/MANYBODY/pair_atm.* @sergeylishchuk
src/MANYBODY/pair_nb3b_screened.* @flodesani
@ -72,6 +74,8 @@ src/MC/fix_sgcmc.* @athomps
src/REAXFF/compute_reaxff_atom.* @rbberger
src/KOKKOS/compute_reaxff_atom_kokkos.* @rbberger
src/REPLICA/fix_pimd_langevin.* @Yi-FanLi
src/DPD-BASIC/pair_dpd_coul_slater_long.* @Eddy-Barraud
src/GPU/pair_dpd_coul_slater_long.* @Eddy-Barraud
# core LAMMPS classes
src/lammps.* @sjplimp

11
.gitignore vendored
View File

@ -43,12 +43,12 @@ Thumbs.db
#cmake
/build*
/CMakeCache.txt
/CMakeFiles/
/Testing
CMakeCache.txt
CMakeFiles
/Makefile
/Testing
/cmake_install.cmake
Testing
Temporary
cmake_install.cmake
/lmp
out/Debug
out/RelWithDebInfo
@ -60,3 +60,4 @@ src/Makefile.package.settings-e
/cmake/build/x64-Debug-Clang
/install/x64-GUI-MSVC
/install
.Rhistory

9
cmake/CMakeLists.txt
View File

@ -23,6 +23,7 @@ project(lammps CXX)
set(SOVERSION 0)
get_property(BUILD_IS_MULTI_CONFIG GLOBAL PROPERTY GENERATOR_IS_MULTI_CONFIG)
include(GNUInstallDirs)
get_filename_component(LAMMPS_DIR ${CMAKE_CURRENT_SOURCE_DIR}/.. ABSOLUTE)
get_filename_component(LAMMPS_LIB_BINARY_DIR ${CMAKE_BINARY_DIR}/lib ABSOLUTE)
# collect all executables and shared libs in the top level build folder
@ -208,7 +209,7 @@ else()
unset(CMAKE_CXX_CLANG_TIDY CACHE)
endif()
include(GNUInstallDirs)
file(GLOB ALL_SOURCES CONFIGURE_DEPENDS ${LAMMPS_SOURCE_DIR}/[^.]*.cpp)
file(GLOB MAIN_SOURCES CONFIGURE_DEPENDS ${LAMMPS_SOURCE_DIR}/main.cpp)
list(REMOVE_ITEM ALL_SOURCES ${MAIN_SOURCES})
@ -256,6 +257,7 @@ set(STANDARD_PACKAGES
DRUDE
EFF
ELECTRODE
EXTRA-COMMAND
EXTRA-COMPUTE
EXTRA-DUMP
EXTRA-FIX
@ -281,10 +283,11 @@ set(STANDARD_PACKAGES
ML-HDNNP
ML-IAP
ML-PACE
ML-POD
ML-QUIP
ML-RANN
ML-SNAP
ML-POD
ML-UF3
MOFFF
MOLECULE
MOLFILE
@ -689,7 +692,7 @@ endif()
# packages which selectively include variants based on enabled styles
# e.g. accelerator packages
######################################################################
foreach(PKG_WITH_INCL CORESHELL DPD-SMOOTH MC MISC PHONON QEQ OPENMP KOKKOS OPT INTEL GPU)
foreach(PKG_WITH_INCL CORESHELL DPD-BASIC DPD-SMOOTH MC MISC PHONON QEQ OPENMP KOKKOS OPT INTEL GPU)
if(PKG_${PKG_WITH_INCL})
include(Packages/${PKG_WITH_INCL})
endif()

9
cmake/Modules/Packages/DPD-BASIC.cmake Normal file
View File

@ -0,0 +1,9 @@
# pair style dpd/coul/slater/long may only be installed if also KSPACE is installed
if(NOT PKG_KSPACE)
get_property(LAMMPS_PAIR_HEADERS GLOBAL PROPERTY PAIR)
list(REMOVE_ITEM LAMMPS_PAIR_HEADERS ${LAMMPS_SOURCE_DIR}/DPD-BASIC/pair_dpd_coul_slater_long.h)
set_property(GLOBAL PROPERTY PAIR "${LAMMPS_PAIR_HEADERS}")
get_target_property(LAMMPS_SOURCES lammps SOURCES)
list(REMOVE_ITEM LAMMPS_SOURCES ${LAMMPS_SOURCE_DIR}/DPD-BASIC/pair_dpd_coul_slater_long.cpp)
set_property(TARGET lammps PROPERTY SOURCES "${LAMMPS_SOURCES}")
endif()

20
cmake/Modules/Packages/PLUMED.cmake
View File

@ -1,5 +1,9 @@
# Plumed2 support for PLUMED package
# for supporting multiple concurrent plumed2 installations for debugging and testing
set(PLUMED_SUFFIX "" CACHE STRING "Suffix for Plumed2 library")
mark_as_advanced(PLUMED_SUFFIX)
if(BUILD_MPI)
set(PLUMED_CONFIG_MPI "--enable-mpi")
set(PLUMED_CONFIG_CC ${CMAKE_MPI_C_COMPILER})
@ -21,9 +25,11 @@ else()
set(PLUMED_CONFIG_OMP "--disable-openmp")
endif()
set(PLUMED_URL "https://github.com/plumed/plumed2/releases/download/v2.8.3/plumed-src-2.8.3.tgz"
# Note: must also adjust check for supported API versions in
# fix_plumed.cpp when version changes from v2.n.x to v2.n+1.y
set(PLUMED_URL "https://github.com/plumed/plumed2/releases/download/v2.9.1/plumed-src-2.9.1.tgz"
CACHE STRING "URL for PLUMED tarball")
set(PLUMED_MD5 "76d23cd394eba9e6530316ed1184e219" CACHE STRING "MD5 checksum of PLUMED tarball")
set(PLUMED_MD5 "c3b2d31479c1e9ce211719d40e9efbd7" CACHE STRING "MD5 checksum of PLUMED tarball")
mark_as_advanced(PLUMED_URL)
mark_as_advanced(PLUMED_MD5)
@ -151,15 +157,15 @@ else()
file(MAKE_DIRECTORY ${INSTALL_DIR}/include)
else()
find_package(PkgConfig REQUIRED)
pkg_check_modules(PLUMED REQUIRED plumed)
pkg_check_modules(PLUMED REQUIRED plumed${PLUMED_SUFFIX})
add_library(LAMMPS::PLUMED INTERFACE IMPORTED)
if(PLUMED_MODE STREQUAL "STATIC")
include(${PLUMED_LIBDIR}/plumed/src/lib/Plumed.cmake.static)
include(${PLUMED_LIBDIR}/plumed${PLUMED_SUFFIX}/src/lib/Plumed.cmake.static)
elseif(PLUMED_MODE STREQUAL "SHARED")
include(${PLUMED_LIBDIR}/plumed/src/lib/Plumed.cmake.shared)
include(${PLUMED_LIBDIR}/plumed${PLUMED_SUFFIX}/src/lib/Plumed.cmake.shared)
elseif(PLUMED_MODE STREQUAL "RUNTIME")
set_target_properties(LAMMPS::PLUMED PROPERTIES INTERFACE_COMPILE_DEFINITIONS "__PLUMED_DEFAULT_KERNEL=${PLUMED_LIBDIR}/${CMAKE_SHARED_LIBRARY_PREFIX}plumedKernel${CMAKE_SHARED_LIBRARY_SUFFIX}")
include(${PLUMED_LIBDIR}/plumed/src/lib/Plumed.cmake.runtime)
set_target_properties(LAMMPS::PLUMED PROPERTIES INTERFACE_COMPILE_DEFINITIONS "__PLUMED_DEFAULT_KERNEL=${PLUMED_LIBDIR}/${CMAKE_SHARED_LIBRARY_PREFIX}plumed${PLUMED_SUFFIX}Kernel${CMAKE_SHARED_LIBRARY_SUFFIX}")
include(${PLUMED_LIBDIR}/plumed${PLUMED_SUFFIX}/src/lib/Plumed.cmake.runtime)
endif()
set_target_properties(LAMMPS::PLUMED PROPERTIES INTERFACE_LINK_LIBRARIES "${PLUMED_LOAD}")
set_target_properties(LAMMPS::PLUMED PROPERTIES INTERFACE_INCLUDE_DIRECTORIES "${PLUMED_INCLUDE_DIRS}")

2
cmake/packaging/build_linux_tgz.sh
View File

@ -59,12 +59,14 @@ done
echo "Set up wrapper script"
MYDIR=$(dirname "$0")
cp ${MYDIR}/xdg-open ${DESTDIR}/bin
cp ${MYDIR}/linux_wrapper.sh ${DESTDIR}/bin
for s in ${DESTDIR}/bin/*
do \
EXE=$(basename $s)
test ${EXE} = linux_wrapper.sh && continue
test ${EXE} = qt.conf && continue
test ${EXE} = xdg-open && continue
ln -s bin/linux_wrapper.sh ${DESTDIR}/${EXE}
done

14
cmake/packaging/linux_wrapper.sh
View File

@ -4,15 +4,17 @@
# reset locale to avoid problems with decimal numbers
export LC_ALL=C
BASEDIR=$(dirname "$0")
EXENAME=$(basename "$0")
BASEDIR="$(dirname "$0")"
EXENAME="$(basename "$0")"
PATH="${BASEDIR}/bin:${PATH}"
# append to LD_LIBRARY_PATH to prefer local (newer) libs
LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${BASEDIR}/lib
LD_LIBRARY_PATH="${LD_LIBRARY_PATH}:${BASEDIR}/lib"
# set some environment variables for LAMMPS etc.
LAMMPS_POTENTIALS=${BASEDIR}/share/lammps/potentials
MSI2LMP_LIBRARY=${BASEDIR}/share/lammps/frc_files
export LD_LIBRARY_PATH LAMMPS_POTENTIALS MSI2LMP_LIBRARY
LAMMPS_POTENTIALS="${BASEDIR}/share/lammps/potentials"
MSI2LMP_LIBRARY="${BASEDIR}/share/lammps/frc_files"
export LD_LIBRARY_PATH LAMMPS_POTENTIALS MSI2LMP_LIBRARY PATH
exec "${BASEDIR}/bin/${EXENAME}" "$@"

1074
cmake/packaging/xdg-open Executable file
View File

File diff suppressed because it is too large Load Diff

4
cmake/presets/all_off.cmake
View File

@ -26,8 +26,9 @@ set(ALL_PACKAGES
DPD-REACT
DPD-SMOOTH
DRUDE
ELECTRODE
EFF
ELECTRODE
EXTRA-COMMAND
EXTRA-COMPUTE
EXTRA-DUMP
EXTRA-FIX
@ -60,6 +61,7 @@ set(ALL_PACKAGES
ML-QUIP
ML-RANN
ML-SNAP
ML-UF3
MOFFF
MOLECULE
MOLFILE

4
cmake/presets/all_on.cmake
View File

@ -28,8 +28,9 @@ set(ALL_PACKAGES
DPD-REACT
DPD-SMOOTH
DRUDE
ELECTRODE
EFF
ELECTRODE
EXTRA-COMMAND
EXTRA-COMPUTE
EXTRA-DUMP
EXTRA-FIX
@ -62,6 +63,7 @@ set(ALL_PACKAGES
ML-QUIP
ML-RANN
ML-SNAP
ML-UF3
MOFFF
MOLECULE
MOLFILE

4
cmake/presets/mingw-cross.cmake
View File

@ -22,8 +22,9 @@ set(WIN_PACKAGES
DPD-REACT
DPD-SMOOTH
DRUDE
ELECTRODE
EFF
ELECTRODE
EXTRA-COMMAND
EXTRA-COMPUTE
EXTRA-DUMP
EXTRA-FIX
@ -50,6 +51,7 @@ set(WIN_PACKAGES
ML-POD
ML-RANN
ML-SNAP
ML-UF3
MOFFF
MOLECULE
MOLFILE

2
cmake/presets/most.cmake
View File

@ -26,6 +26,7 @@ set(ALL_PACKAGES
DRUDE
EFF
ELECTRODE
EXTRA-COMMAND
EXTRA-COMPUTE
EXTRA-DUMP
EXTRA-FIX
@ -45,6 +46,7 @@ set(ALL_PACKAGES
ML-IAP
ML-POD
ML-SNAP
ML-UF3
MOFFF
MOLECULE
OPENMP

4
cmake/presets/windows.cmake
View File

@ -22,6 +22,7 @@ set(WIN_PACKAGES
DRUDE
EFF
ELECTRODE
EXTRA-COMMAND
EXTRA-COMPUTE
EXTRA-DUMP
EXTRA-FIX
@ -42,6 +43,7 @@ set(WIN_PACKAGES
ML-IAP
ML-POD
ML-SNAP
ML-UF3
MOFFF
MOLECULE
MOLFILE
@ -50,8 +52,8 @@ set(WIN_PACKAGES
ORIENT
PERI
PHONON
POEMS
PLUGIN
POEMS
PTM
QEQ
QTB

4
doc/lammps.1
View File

@ -1,7 +1,7 @@
.TH LAMMPS "1" "17 April 2024" "2024-04-17"
.TH LAMMPS "1" "27 June 2024" "2024-06-27"
.SH NAME
.B LAMMPS
\- Molecular Dynamics Simulator. Version 17 April 2024
\- Molecular Dynamics Simulator. Version 27 June 2024
.SH SYNOPSIS
.B lmp

4
doc/src/Commands_bond.rst
View File

@ -27,7 +27,7 @@ OPT.
* :doc:`none <bond_none>`
* :doc:`zero <bond_zero>`
* :doc:`hybrid <bond_hybrid>`
* :doc:`hybrid (k) <bond_hybrid>`
*
*
*
@ -72,7 +72,7 @@ OPT.
* :doc:`none <angle_none>`
* :doc:`zero <angle_zero>`
* :doc:`hybrid (k) <angle_hybrid>`
* :doc:`hybrid <angle_hybrid>`
*
*
*

4
doc/src/Commands_compute.rst
View File

@ -108,6 +108,10 @@ KOKKOS, o = OPENMP, t = OPT.
* :doc:`pe/mol/tally <compute_tally>`
* :doc:`pe/tally <compute_tally>`
* :doc:`plasticity/atom <compute_plasticity_atom>`
* :doc:`pod/atom <compute_pod_atom>`
* :doc:`podd/atom <compute_pod_atom>`
* :doc:`pod/local <compute_pod_atom>`
* :doc:`pod/global <compute_pod_atom>`
* :doc:`pressure <compute_pressure>`
* :doc:`pressure/alchemy <compute_pressure_alchemy>`
* :doc:`pressure/uef <compute_pressure_uef>`

16
doc/src/Commands_pair.rst
View File

@ -25,16 +25,16 @@ OPT.
* :doc:`none <pair_none>`
* :doc:`zero <pair_zero>`
* :doc:`hybrid (k) <pair_hybrid>`
* :doc:`hybrid/overlay (k) <pair_hybrid>`
* :doc:`hybrid/scaled <pair_hybrid>`
* :doc:`hybrid (ko) <pair_hybrid>`
* :doc:`hybrid/molecular (o) <pair_hybrid>`
* :doc:`hybrid/overlay (ko) <pair_hybrid>`
* :doc:`hybrid/scaled (o) <pair_hybrid>`
* :doc:`kim <pair_kim>`
* :doc:`list <pair_list>`
* :doc:`tracker <pair_tracker>`
*
*
*
*
* :doc:`adp (ko) <pair_adp>`
* :doc:`agni (o) <pair_agni>`
* :doc:`aip/water/2dm (t) <pair_aip_water_2dm>`
@ -94,9 +94,10 @@ OPT.
* :doc:`coul/wolf (ko) <pair_coul>`
* :doc:`coul/wolf/cs <pair_cs>`
* :doc:`dpd (giko) <pair_dpd>`
* :doc:`dpd/fdt <pair_dpd_fdt>`
* :doc:`dpd/coul/slater/long (g) <pair_dpd_coul_slater_long>`
* :doc:`dpd/ext (ko) <pair_dpd_ext>`
* :doc:`dpd/ext/tstat (ko) <pair_dpd_ext>`
* :doc:`dpd/fdt <pair_dpd_fdt>`
* :doc:`dpd/fdt/energy (k) <pair_dpd_fdt>`
* :doc:`dpd/tstat (gko) <pair_dpd>`
* :doc:`dsmc <pair_dsmc>`
@ -246,7 +247,7 @@ OPT.
* :doc:`pace (k) <pair_pace>`
* :doc:`pace/extrapolation (k) <pair_pace>`
* :doc:`pedone (o) <pair_pedone>`
* :doc:`pod <pair_pod>`
* :doc:`pod (k) <pair_pod>`
* :doc:`peri/eps <pair_peri>`
* :doc:`peri/lps (o) <pair_peri>`
* :doc:`peri/pmb (o) <pair_peri>`
@ -269,7 +270,7 @@ OPT.
* :doc:`smd/ulsph <pair_smd_ulsph>`
* :doc:`smtbq <pair_smtbq>`
* :doc:`snap (ik) <pair_snap>`
* :doc:`soft (go) <pair_soft>`
* :doc:`soft (gko) <pair_soft>`
* :doc:`sph/heatconduction (g) <pair_sph_heatconduction>`
* :doc:`sph/idealgas <pair_sph_idealgas>`
* :doc:`sph/lj (g) <pair_sph_lj>`
@ -303,6 +304,7 @@ OPT.
* :doc:`tip4p/long/soft (o) <pair_fep_soft>`
* :doc:`tri/lj <pair_tri_lj>`
* :doc:`ufm (got) <pair_ufm>`
* :doc:`uf3 (k) <pair_uf3>`
* :doc:`vashishta (gko) <pair_vashishta>`
* :doc:`vashishta/table (o) <pair_vashishta>`
* :doc:`wf/cut <pair_wf_cut>`

4
doc/src/Commands_removed.rst
View File

@ -151,10 +151,10 @@ and allow running LAMMPS with GPU acceleration.
i-PI tool
---------
.. versionchanged:: TBD
.. versionchanged:: 27June2024
The i-PI tool has been removed from the LAMMPS distribution. Instead,
instructions to install i-PI from PyPi via pip are provided.
instructions to install i-PI from PyPI via pip are provided.
restart2data tool
-----------------

3
doc/src/Developer_utils.rst
View File

@ -211,6 +211,9 @@ Argument processing
.. doxygenfunction:: bounds
:project: progguide
.. doxygenfunction:: bounds_typelabel
:project: progguide
.. doxygenfunction:: expand_args
:project: progguide

198
doc/src/Howto_bioFF.rst
View File

@ -1,6 +1,10 @@
CHARMM, AMBER, COMPASS, and DREIDING force fields
=================================================
A compact summary of the concepts, definitions, and properties of
force fields with explicit bonded interactions (like the ones discussed
in this HowTo) is given in :ref:`(Gissinger) <Typelabel2>`.
A force field has 2 parts: the formulas that define it and the
coefficients used for a particular system. Here we only discuss
formulas implemented in LAMMPS that correspond to formulas commonly used
@ -11,12 +15,42 @@ commands like :doc:`pair_coeff <pair_coeff>` or :doc:`bond_coeff
<bond_coeff>` and so on. See the :doc:`Tools <Tools>` doc page for
additional tools that can use CHARMM, AMBER, or Materials Studio
generated files to assign force field coefficients and convert their
output into LAMMPS input.
output into LAMMPS input. LAMMPS input scripts can also be generated by
`charmm-gui.org <https://charmm-gui.org/>`_.
See :ref:`(MacKerell) <howto-MacKerell>` for a description of the CHARMM
force field. See :ref:`(Cornell) <howto-Cornell>` for a description of
the AMBER force field. See :ref:`(Sun) <howto-Sun>` for a description
of the COMPASS force field.
CHARMM and AMBER
----------------
The `CHARMM force field
<https://mackerell.umaryland.edu/charmm_ff.shtml>`_ :ref:`(MacKerell)
<howto-MacKerell>` and `AMBER force field
<https://ambermd.org/AmberModels.php>`_ :ref:`(Cornell) <howto-Cornell>`
have potential energy function of the form
.. math::
V & = \sum_{bonds} E_b + \sum_{angles} \!E_a + \!\overbrace{\sum_{dihedral} \!\!E_d}^{\substack{
\text{charmm} \\
\text{charmmfsw}
}} +\!\!\! \sum_{impropers} \!\!\!E_i \\[.6em]
& \quad + \!\!\!\!\!\!\!\!\!\!\underbrace{~\sum_{pairs} \left(E_{LJ}+E_{coul}\right)}_{\substack{
\text{lj/charmm/coul/charmm} \\
\text{lj/charmm/coul/charmm/implicit} \\
\text{lj/charmm/coul/long} \\
\text{lj/charmm/coul/msm} \\
\text{lj/charmmfsw/coul/charmmfsh} \\
\text{lj/charmmfsw/coul/long}
}} \!\!\!\!\!\!\!\!+ \!\!\sum_{special}\! E_s + \!\!\!\!\sum_{residues} \!\!\!{\scriptstyle\mathrm{CMAP}(\phi,\psi)}
The terms are computed by bond styles (relationship between 2 atoms),
angle styles (between 3 atoms) , dihedral/improper styles (between 4
atoms), pair styles (non-covalently bonded pair interactions) and
special bonds. The CMAP term (see :doc:`fix cmap <fix_cmap>` command for
details) corrects for pairs of dihedral angles ("Correction MAP") to
significantly improve the structural and dynamic properties of proteins
in crystalline and solution environments :ref:`(Brooks)
<howto-Brooks>`. The AMBER force field does not include the CMAP term.
The interaction styles listed below compute force field formulas that
are consistent with common options in CHARMM or AMBER. See each
@ -31,10 +65,81 @@ command's documentation for the formula it computes.
* :doc:`pair_style <pair_charmm>` lj/charmm/coul/charmm
* :doc:`pair_style <pair_charmm>` lj/charmm/coul/charmm/implicit
* :doc:`pair_style <pair_charmm>` lj/charmm/coul/long
* :doc:`special_bonds <special_bonds>` charmm
* :doc:`special_bonds <special_bonds>` amber
The pair styles compute Lennard Jones (LJ) and Coulombic interactions
with additional switching or shifting functions that ramp the energy
and/or force smoothly to zero between an inner :math:`(a)` and outer
:math:`(b)` cutoff. The older styles with *charmm* (not *charmmfsw* or
*charmmfsh*\ ) in their name compute the LJ and Coulombic interactions
with an energy switching function (esw) S(r) which ramps the energy
smoothly to zero between the inner and outer cutoff. This can cause
irregularities in pairwise forces (due to the discontinuous second
derivative of energy at the boundaries of the switching region), which
in some cases can result in complications in energy minimization and
detectable artifacts in MD simulations.
.. grid:: 1 1 2 2
.. grid-item::
.. math::
LJ(r) &= 4 \epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} -
\left(\frac{\sigma}{r}\right)^6 \right]\\[.6em]
C(r) &= \frac{C q_i q_j}{ \epsilon r}\\[.6em]
S(r) &= \frac{ \left(b^2 - r^2\right)^2 \left(b^2 + 2r^2 - 3{a^2}\right)}
{ \left(b^2 - a^2\right)^3 }\\[.6em]
E_{LJ}(r) &= \begin{cases}
LJ(r), & r \leq a \\
LJ(r) S(r), & a < r \leq b \\
0, &r > b
\end{cases} \\[.6em]
E_{coul}(r) &= \begin{cases}
C(r), & r \leq a \\
C(r) S(r), & a < r \leq b \\
0, & r > b
\end{cases}
.. grid-item::
.. image:: img/howto_charmm_ELJ.png
:align: center
The newer styles with *charmmfsw* or *charmmfsh* in their name replace
energy switching with force switching (fsw) for LJ interactions and
force shifting (fsh) functions for Coulombic interactions
:ref:`(Steinbach) <howto-Steinbach>`
.. grid:: 1 1 2 2
.. grid-item::
.. math::
E_{LJ}(r) = & \begin{cases}
4 \epsilon \sigma^6 \left(\frac{\displaystyle\sigma
^6-r^6}{\displaystyle r^{12}}-\frac{\displaystyle\sigma ^6}{\displaystyle a^6
b^6}+\frac{\displaystyle 1}{\displaystyle a^3 b^3}\right) & r\leq a \\
\frac{\displaystyle 4 \epsilon \sigma^6 \left(\sigma ^6
\left(b^6-r^6\right)^2-b^3 r^6 \left(a^3+b^3\right)
\left(b^3-r^3\right)^2\right)}{\displaystyle b^6 r^{12}
\left(b^6-a^6\right)} & a<r \leq b\\
0, & r>b
\end{cases}\\[.6em]
E_{coul}(r) & = \begin{cases}
C(r) \frac{\displaystyle (b-r)^2}{\displaystyle r b^2}, & r \leq b \\
0, & r > b
\end{cases}
.. grid-item::
.. image:: img/howto_charmmfsw_ELJ.png
:align: center
These styles are used by LAMMPS input scripts generated by
https://charmm-gui.org/ :ref:`(Brooks) <howto-Brooks>`.
.. note::
For CHARMM, newer *charmmfsw* or *charmmfsh* styles were released in
@ -43,17 +148,33 @@ command's documentation for the formula it computes.
<pair_charmm>` and :doc:`dihedral charmm <dihedral_charmm>` doc
pages.
.. note::
The TIP3P water model is strongly recommended for use with the CHARMM
force field. In fact, `"using the SPC model with CHARMM parameters is
a bad idea"
<https://matsci.org/t/using-spc-water-with-charmm-ff/24715>`_ and `"to
enable TIP4P style water in CHARMM, you would have to write a new pair
style"
<https://matsci.org/t/hybrid-pair-styles-for-charmm-and-tip4p-ew/32609>`_
. LAMMPS input scripts generated by Solution Builder on https://charmm-gui.org
use TIP3P molecules for solvation. Any other water model can and
probably will lead to false conclusions.
COMPASS
-------
COMPASS is a general force field for atomistic simulation of common
organic molecules, inorganic small molecules, and polymers which was
developed using ab initio and empirical parameterization techniques.
See the :doc:`Tools <Tools>` page for the msi2lmp tool for creating
LAMMPS template input and data files from BIOVIA's Materials Studio
files. Please note that the msi2lmp tool is very old and largely
unmaintained, so it does not support all features of Materials Studio
provided force field files, especially additions during the last decade.
You should watch the output carefully and compare results, where
possible. See :ref:`(Sun) <howto-Sun>` for a description of the COMPASS force
field.
developed using ab initio and empirical parameterization techniques
:ref:`(Sun) <howto-Sun>`. See the :doc:`Tools <Tools>` page for the
msi2lmp tool for creating LAMMPS template input and data files from
BIOVIA's Materials Studio files. Please note that the msi2lmp tool is
very old and largely unmaintained, so it does not support all features
of Materials Studio provided force field files, especially additions
during the last decade. You should watch the output carefully and
compare results, where possible. See :ref:`(Sun) <howto-Sun>` for a
description of the COMPASS force field.
These interaction styles listed below compute force field formulas that
are consistent with the COMPASS force field. See each command's
@ -70,14 +191,21 @@ documentation for the formula it computes.
* :doc:`special_bonds <special_bonds>` lj/coul 0 0 1
DREIDING is a generic force field developed by the `Goddard group <http://www.wag.caltech.edu>`_ at Caltech and is useful for
predicting structures and dynamics of organic, biological and main-group
inorganic molecules. The philosophy in DREIDING is to use general force
constants and geometry parameters based on simple hybridization
considerations, rather than individual force constants and geometric
parameters that depend on the particular combinations of atoms involved
in the bond, angle, or torsion terms. DREIDING has an :doc:`explicit hydrogen bond term <pair_hbond_dreiding>` to describe interactions involving a
hydrogen atom on very electronegative atoms (N, O, F).
DREIDING
--------
DREIDING is a generic force field developed by the `Goddard group
<http://www.wag.caltech.edu>`_ at Caltech and is useful for predicting
structures and dynamics of organic, biological and main-group inorganic
molecules. The philosophy in DREIDING is to use general force constants
and geometry parameters based on simple hybridization considerations,
rather than individual force constants and geometric parameters that
depend on the particular combinations of atoms involved in the bond,
angle, or torsion terms. DREIDING has an :doc:`explicit hydrogen bond
term <pair_hbond_dreiding>` to describe interactions involving a
hydrogen atom on very electronegative atoms (N, O, F). Unlike CHARMM
or AMBER, the DREIDING force field has not been parameterized for
considering solvents (like water).
See :ref:`(Mayo) <howto-Mayo>` for a description of the DREIDING force field
@ -110,21 +238,31 @@ documentation for the formula it computes.
----------
.. _Typelabel2:
**(Gissinger)** J. R. Gissinger, I. Nikiforov, Y. Afshar, B. Waters, M. Choi, D. S. Karls, A. Stukowski, W. Im, H. Heinz, A. Kohlmeyer, and E. B. Tadmor, J Phys Chem B, 128, 3282-3297 (2024).
.. _howto-MacKerell:
**(MacKerell)** MacKerell, Bashford, Bellott, Dunbrack, Evanseck, Field,
Fischer, Gao, Guo, Ha, et al, J Phys Chem, 102, 3586 (1998).
**(MacKerell)** MacKerell, Bashford, Bellott, Dunbrack, Evanseck, Field, Fischer, Gao, Guo, Ha, et al (1998). J Phys Chem, 102, 3586 . https://doi.org/10.1021/jp973084f
.. _howto-Cornell:
**(Cornell)** Cornell, Cieplak, Bayly, Gould, Merz, Ferguson,
Spellmeyer, Fox, Caldwell, Kollman, JACS 117, 5179-5197 (1995).
**(Cornell)** Cornell, Cieplak, Bayly, Gould, Merz, Ferguson, Spellmeyer, Fox, Caldwell, Kollman (1995). JACS 117, 5179-5197. https://doi.org/10.1021/ja00124a002
.. _howto-Steinbach:
**(Steinbach)** Steinbach, Brooks (1994). J Comput Chem, 15, 667. https://doi.org/10.1002/jcc.540150702
.. _howto-Brooks:
**(Brooks)** Brooks, et al (2009). J Comput Chem, 30, 1545. https://onlinelibrary.wiley.com/doi/10.1002/jcc.21287
.. _howto-Sun:
**(Sun)** Sun, J. Phys. Chem. B, 102, 7338-7364 (1998).
**(Sun)** Sun (1998). J. Phys. Chem. B, 102, 7338-7364. https://doi.org/10.1021/jp980939v
.. _howto-Mayo:
**(Mayo)** Mayo, Olfason, Goddard III, J Phys Chem, 94, 8897-8909
(1990).
**(Mayo)** Mayo, Olfason, Goddard III (1990). J Phys Chem, 94, 8897-8909. https://doi.org/10.1021/j100389a010

52
doc/src/Packages_details.rst
View File

@ -52,6 +52,7 @@ page gives those details.
* :ref:`DRUDE <PKG-DRUDE>`
* :ref:`EFF <PKG-EFF>`
* :ref:`ELECTRODE <PKG-ELECTRODE>`
* :ref:`EXTRA-COMMAND <PKG-EXTRA-COMMAND>`
* :ref:`EXTRA-COMPUTE <PKG-EXTRA-COMPUTE>`
* :ref:`EXTRA-DUMP <PKG-EXTRA-DUMP>`
* :ref:`EXTRA-FIX <PKG-EXTRA-FIX>`
@ -84,6 +85,7 @@ page gives those details.
* :ref:`ML-QUIP <PKG-ML-QUIP>`
* :ref:`ML-RANN <PKG-ML-RANN>`
* :ref:`ML-SNAP <PKG-ML-SNAP>`
* :ref:`ML-UF3 <PKG-ML-UF3>`
* :ref:`MOFFF <PKG-MOFFF>`
* :ref:`MOLECULE <PKG-MOLECULE>`
* :ref:`MOLFILE <PKG-MOLFILE>`
@ -403,6 +405,7 @@ and :ref:`ASPHERE <PKG-ASPHERE>` packages are installed.
* :doc:`bond_style oxdna2/\* <bond_oxdna>`
* :doc:`bond_style oxrna2/\* <bond_oxdna>`
* :doc:`fix nve/dotc/langevin <fix_nve_dotc_langevin>`
* examples/PACKAGES/cgdna
----------
@ -676,7 +679,12 @@ DPD-BASIC package
Pair styles for the basic dissipative particle dynamics (DPD) method
and DPD thermostatting.
**Author:** Kurt Smith (U Pittsburgh), Martin Svoboda, Martin Lisal (ICPF and UJEP)
Pair style :doc:`dpd/coul/slater/long <pair_dpd_coul_slater_long>` also
includes smeared charges for coulomb interactions and thus requires the
:ref:`KSPACE <PKG-KSPACE>` package to be installed to handle the long-range
Coulomb part of the interactions.
**Authors:** Kurt Smith (U Pittsburgh), Martin Svoboda, Martin Lisal (ICPF and UJEP), Eddy Barraud (IFPEN)
**Supporting info:**
@ -685,6 +693,7 @@ and DPD thermostatting.
* :doc:`pair_style dpd/tstat <pair_dpd>`
* :doc:`pair_style dpd/ext <pair_dpd_ext>`
* :doc:`pair_style dpd/ext/tstat <pair_dpd_ext>`
* :doc:`pair_style dpd/coul/slater/long <pair_dpd_coul_slater_long>`
* examples/PACKAGES/dpd-basic
----------
@ -886,6 +895,22 @@ This package has :ref:`specific installation instructions <electrode>` on the
----------
.. _PKG-EXTRA-COMMAND:
EXTRA-COMMAND package
---------------------
**Contents:**
Additional command styles that are less commonly used.
**Supporting info:**
* src/EXTRA-COMMAND: filenames -> commands
* :doc:`general commands <Commands_all>`
----------
.. _PKG-EXTRA-COMPUTE:
EXTRA-COMPUTE package
@ -1925,6 +1950,31 @@ computes which analyze attributes of the potential.
----------
.. _PKG-ML-UF3:
ML-UF3 package
--------------
**Contents:**
A pair style for the ultra-fast force field potentials (UF3). UF3 is a
methodology for deriving a highly accurate classical potential which is
fast to evaluate and is fitted to a large archives of quantum mechanical
(DFT) data. The use of b-spline basis set in UF3 enables the rapid
evaluation of 2-body and 3-body interactions.
**Authors:** Ajinkya C Hire (University of Florida),
Hendrik Krass (University of Constance),
Matthias Rupp (Luxembourg Institute of Science and Technology),
Richard Hennig (University of Florida)
**Supporting info:**
* src/ML-UF3: filenames -> commands
* :doc:`pair_style uf3 <pair_uf3>`
* examples/uf3
* https://github.com/uf3/uf3
.. _PKG-MOFFF:
MOFFF package

10
doc/src/Packages_list.rst
View File

@ -158,6 +158,11 @@ whether an extra library is needed to build and use the package:
- :doc:`fix electrode/conp <fix_electrode>`
- PACKAGES/electrode
- no
* - :ref:`EXTRA-COMMAND <PKG-EXTRA-COMMAND>`
- additional command styles
- :doc:`general commands <Commands_all>`
- n/a
- no
* - :ref:`EXTRA-COMPUTE <PKG-EXTRA-COMPUTE>`
- additional compute styles
- :doc:`compute <compute>`
@ -318,6 +323,11 @@ whether an extra library is needed to build and use the package:
- :doc:`pair_style snap <pair_snap>`
- snap
- no
* - :ref:`ML-UF3 <PKG-ML-UF3>`
- quantum-fitted ultra fast potentials
- :doc:`pair_style uf3 <pair_uf3>`
- PACKAGES/uf3
- no
* - :ref:`MOFFF <PKG-MOFFF>`
- styles for `MOF-FF <MOFplus_>`_ force field
- :doc:`pair_style buck6d/coul/gauss <pair_buck6d_coul_gauss>`

39
doc/src/Tools.rst
View File

@ -379,7 +379,7 @@ See README file in the tools/fep directory.
i-PI tool
-------------------
.. versionchanged:: TBD
.. versionchanged:: 27June2024
The tools/i-pi directory used to contain a bundled version of the i-PI
software package for use with LAMMPS. This version, however, was
@ -389,7 +389,7 @@ The i-PI package was created and is maintained by Michele Ceriotti,
michele.ceriotti at gmail.com, to interface to a variety of molecular
dynamics codes.
i-PI is now available via PyPi using the pip package manager at:
i-PI is now available via PyPI using the pip package manager at:
https://pypi.org/project/ipi/
Here are the commands to set up a virtual environment and install
@ -728,8 +728,8 @@ CMake is required.
The LAMMPS GUI has been successfully compiled and tested on:
- Ubuntu Linux 20.04LTS x86_64 using GCC 9, Qt version 5.12
- Fedora Linux 38 x86\_64 using GCC 13 and Clang 16, Qt version 5.15LTS
- Fedora Linux 38 x86\_64 using GCC 13, Qt version 6.5LTS
- Fedora Linux 40 x86\_64 using GCC 14 and Clang 17, Qt version 5.15LTS
- Fedora Linux 40 x86\_64 using GCC 14, Qt version 6.5LTS
- Apple macOS 12 (Monterey) and macOS 13 (Ventura) with Xcode on arm64 and x86\_64, Qt version 5.15LTS
- Windows 10 and 11 x86_64 with Visual Studio 2022 and Visual C++ 14.36, Qt version 5.15LTS
- Windows 10 and 11 x86_64 with MinGW / GCC 10.0 cross-compiler on Fedora 38, Qt version 5.15LTS
@ -771,22 +771,23 @@ if necessary. When both, Qt5 and Qt6 are available, Qt6 will be preferred
unless ``-D LAMMPS_GUI_USE_QT5=yes`` is set.
It should be possible to build the LAMMPS GUI as a standalone
compilation (e.g. when LAMMPS has been compiled with traditional make),
then the CMake configuration needs to be told where to find the LAMMPS
compilation (e.g. when LAMMPS has been compiled with traditional make).
Then the CMake configuration needs to be told where to find the LAMMPS
headers and the LAMMPS library, via ``-D
LAMMPS_SOURCE_DIR=/path/to/lammps/src``. CMake will try to guess a
build folder with the LAMMPS library from that path, but it can also be
set with ``-D LAMMPS_LIB_DIR=/path/to/lammps/lib``.
Rather than linking to the LAMMPS library during compilation, it is also
possible to compile the GUI with a plugin loader library that will load
possible to compile the GUI with a plugin loader that will load
the LAMMPS library dynamically at runtime during the start of the GUI
from a shared library; e.g. ``liblammps.so`` or ``liblammps.dylib`` or
``liblammps.dll`` (depending on the operating system). This has the
advantage that the LAMMPS library can be updated LAMMPS without having
to recompile the GUI. The ABI of the LAMMPS C-library interface is very
stable and generally backward compatible. This feature is enabled by
setting ``-D LAMMPS_GUI_USE_PLUGIN=on`` and then ``-D
advantage that the LAMMPS library can be built from updated or modified
LAMMPS source without having to recompile the GUI. The ABI of the
LAMMPS C-library interface is very stable and generally backward
compatible. This feature is enabled by setting
``-D LAMMPS_GUI_USE_PLUGIN=on`` and then ``-D
LAMMPS_PLUGINLIB_DIR=/path/to/lammps/plugin/loader``. Typically, this
would be the ``examples/COUPLE/plugin`` folder of the LAMMPS
distribution.
@ -798,8 +799,8 @@ macOS
"""""
When building on macOS, the build procedure will try to manufacture a
drag-n-drop installer, LAMMPS-macOS-multiarch.dmg, when using the 'dmg'
target (i.e. ``cmake --build <build dir> --target dmg`` or ``make dmg``.
drag-n-drop installer, ``LAMMPS-macOS-multiarch.dmg``, when using the
'dmg' target (i.e. ``cmake --build <build dir> --target dmg`` or ``make dmg``.
To build multi-arch executables that will run on both, arm64 and x86_64
architectures natively, it is necessary to set the CMake variable ``-D
@ -838,12 +839,12 @@ and LAMMPS GUI can be launched from anywhere from the command line.
The standard CMake build procedure can be applied and the
``mingw-cross.cmake`` preset used. By using ``mingw64-cmake`` the CMake
command will automatically include a suitable CMake toolset file (the
regular cmake command can be used after that). After building the
libraries and executables, you can build the target 'zip'
(i.e. ``cmake --build <build dir> --target zip`` or ``make zip``
to stage all installed files into a LAMMPS_GUI folder and then
run a script to copy all required dependencies, some other files,
command will automatically include a suitable CMake toolchain file (the
regular cmake command can be used after that to modify the configuration
settings, if needed). After building the libraries and executables,
you can build the target 'zip' (i.e. ``cmake --build <build dir> --target zip``
or ``make zip`` to stage all installed files into a LAMMPS_GUI folder
and then run a script to copy all required dependencies, some other files,
and create a zip file from it.
Linux

122
doc/src/bond_oxdna.rst
View File

@ -27,6 +27,7 @@ Examples
.. code-block:: LAMMPS
# LJ units
bond_style oxdna/fene
bond_coeff * 2.0 0.25 0.7525
@ -36,6 +37,32 @@ Examples
bond_style oxrna2/fene
bond_coeff * 2.0 0.25 0.76107
bond_style oxdna/fene
bond_coeff * oxdna_lj.cgdna
# Real units
bond_style oxdna/fene
bond_coeff * 11.92337812042065 2.1295 6.409795
bond_style oxdna2/fene
bond_coeff * 11.92337812042065 2.1295 6.4430152
bond_style oxrna2/fene
bond_coeff * 11.92337812042065 2.1295 6.482800913
bond_style oxrna2/fene
bond_coeff * oxrna2_real.cgdna
.. note::
The coefficients in the above examples have to be kept fixed and
cannot be changed without reparameterizing the entire model. They are
provided in forms compatible with both *units lj* and *units real*
(see documentation of :doc:`units <units>`). These can also be read
from a potential file with correct unit style by specifying the name
of the file. Several potential files for each unit style are included
in the ``potentials`` directory of the LAMMPS distribution.
Description
"""""""""""
@ -46,15 +73,14 @@ The *oxdna/fene*, *oxdna2/fene*, and *oxrna2/fene* bond styles use the potential
E = - \frac{\epsilon}{2} \ln \left[ 1 - \left(\frac{r-r_0}{\Delta}\right)^2\right]
to define a modified finite extensible nonlinear elastic (FENE)
potential :ref:`(Ouldridge) <Ouldridge0>` to model the connectivity of the
phosphate backbone in the oxDNA/oxRNA force field for coarse-grained
potential :ref:`(Ouldridge) <Ouldridge0>` to model the connectivity of
the phosphate backbone in the oxDNA/oxRNA force field for coarse-grained
modelling of DNA/RNA.
The following coefficients must be defined for the bond type via the
:doc:`bond_coeff <bond_coeff>` command as given in the above example, or
in the data file or restart files read by the
:doc:`read_data <read_data>` or :doc:`read_restart <read_restart>`
commands:
in the data file or restart files read by the :doc:`read_data
<read_data>` or :doc:`read_restart <read_restart>` commands:
* :math:`\epsilon` (energy)
* :math:`\Delta` (distance)
@ -62,41 +88,40 @@ commands:
.. note::
The oxDNA bond style has to be used together with the
corresponding oxDNA pair styles for excluded volume interaction
*oxdna/excv* , stacking *oxdna/stk* , cross-stacking *oxdna/xstk* and
coaxial stacking interaction *oxdna/coaxstk* as well as
hydrogen-bonding interaction *oxdna/hbond* (see also documentation of
:doc:`pair_style oxdna/excv <pair_oxdna>`). For the oxDNA2
:ref:`(Snodin) <Snodin0>` bond style the analogous pair styles
*oxdna2/excv* , *oxdna2/stk* , *oxdna2/xstk* , *oxdna2/coaxstk* ,
*oxdna2/hbond* and an additional Debye-Hueckel pair style
*oxdna2/dh* have to be defined. The same applies to the oxRNA2
:ref:`(Sulc1) <Sulc01>` styles.
The coefficients in the above example have to be kept fixed and cannot
be changed without reparameterizing the entire model.
The oxDNA bond style has to be used together with the corresponding
oxDNA pair styles for excluded volume interaction *oxdna/excv* ,
stacking *oxdna/stk* , cross-stacking *oxdna/xstk* and coaxial
stacking interaction *oxdna/coaxstk* as well as hydrogen-bonding
interaction *oxdna/hbond* (see also documentation of :doc:`pair_style
oxdna/excv <pair_oxdna>`). For the oxDNA2 :ref:`(Snodin) <Snodin0>`
bond style the analogous pair styles *oxdna2/excv* , *oxdna2/stk* ,
*oxdna2/xstk* , *oxdna2/coaxstk* , *oxdna2/hbond* and an additional
Debye-Hueckel pair style *oxdna2/dh* have to be defined. The same
applies to the oxRNA2 :ref:`(Sulc1) <Sulc01>` styles.
.. note::
This bond style has to be used with the *atom_style hybrid bond ellipsoid oxdna*
(see documentation of :doc:`atom_style <atom_style>`). The *atom_style oxdna*
stores the 3'-to-5' polarity of the nucleotide strand, which is set through
the bond topology in the data file. The first (second) atom in a bond definition
is understood to point towards the 3'-end (5'-end) of the strand.
This bond style has to be used with the *atom_style hybrid bond
ellipsoid oxdna* (see documentation of :doc:`atom_style
<atom_style>`). The *atom_style oxdna* stores the 3'-to-5' polarity
of the nucleotide strand, which is set through the bond topology in
the data file. The first (second) atom in a bond definition is
understood to point towards the 3'-end (5'-end) of the strand.
.. warning::
If data files are produced with :doc:`write_data <write_data>`, then the
:doc:`newton <newton>` command should be set to *newton on* or *newton off on*.
Otherwise the data files will not have the same 3'-to-5' polarity as the
initial data file. This limitation does not apply to binary restart files
produced with :doc:`write_restart <write_restart>`.
If data files are produced with :doc:`write_data <write_data>`, then
the :doc:`newton <newton>` command should be set to *newton on* or
*newton off on*. Otherwise the data files will not have the same
3'-to-5' polarity as the initial data file. This limitation does not
apply to binary restart files produced with :doc:`write_restart
<write_restart>`.
Example input and data files for DNA and RNA duplexes can be found in
examples/PACKAGES/cgdna/examples/oxDNA/ , /oxDNA2/ and /oxRNA2/. A simple python
setup tool which creates single straight or helical DNA strands, DNA/RNA
duplexes or arrays of DNA/RNA duplexes can be found in
examples/PACKAGES/cgdna/util/.
``examples/PACKAGES/cgdna/examples/oxDNA/`, `.../oxDNA2/`` and
``.../oxRNA2/``. A simple python setup tool which creates single
straight or helical DNA strands, DNA/RNA duplexes or arrays of DNA/RNA
duplexes can be found in ``examples/PACKAGES/cgdna/util/``.
Please cite :ref:`(Henrich) <Henrich0>` in any publication that uses
this implementation. An updated documentation that contains general information
@ -113,6 +138,39 @@ and for sequence-specific hydrogen-bonding and stacking interactions
----------
Potential file reading
""""""""""""""""""""""
For each style oxdna, oxdna2 and oxrna2, the first parameter argument
can be a filename, and if it is, no further arguments should be
supplied. Therefore the following command:
.. code-block:: LAMMPS
bond_style oxdna/fene
bond_coeff * oxdna_lj.cgdna
will be interpreted as a request to read the (FENE) potential
:ref:`(Ouldridge) <Ouldridge0>` parameters from the file with the given
name. The file can define multiple potential parameters for both bonded
and pair interactions, but for the above bonded interactions there must
exist in the file a line of the form:
.. code-block:: LAMMPS
* fene epsilon delta r0
There are sample potential files for each unit style in the
``potentials`` directory of the LAMMPS distribution. The potential file
unit system must align with the units defined via the :doc:`units
<units>` command. For conversion between different *LJ* and *real* unit
systems for oxDNA, the python tool *lj2real.py* located in the
``examples/PACKAGES/cgdna/util/`` directory can be used. This tool
assumes similar file structure to the examples found in
``examples/PACKAGES/cgdna/examples/``.
----------
Restrictions
""""""""""""

4
doc/src/compute.rst
View File

@ -272,6 +272,10 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` pag
* :doc:`pe/mol/tally <compute_tally>` - potential energy between two groups of atoms separated into intermolecular and intramolecular components via the tally callback mechanism
* :doc:`pe/tally <compute_tally>` - potential energy between two groups of atoms via the tally callback mechanism
* :doc:`plasticity/atom <compute_plasticity_atom>` - Peridynamic plasticity for each atom
* :doc:`pod/atom <compute_pod_atom>` - POD descriptors for each atom
* :doc:`podd/atom <compute_pod_atom>` - derivative of POD descriptors for each atom
* :doc:`pod/local <compute_pod_atom>` - local POD descriptors and their derivatives
* :doc:`pod/global <compute_pod_atom>` - global POD descriptors and their derivatives
* :doc:`pressure <compute_pressure>` - total pressure and pressure tensor
* :doc:`pressure/alchemy <compute_pressure_alchemy>` - mixed system total pressure and pressure tensor for :doc:`fix alchemy <fix_alchemy>` runs
* :doc:`pressure/uef <compute_pressure_uef>` - pressure tensor in the reference frame of an applied flow field

145
doc/src/compute_pod_atom.rst Normal file
View File

@ -0,0 +1,145 @@
.. index:: compute pod/atom
.. index:: compute podd/atom
.. index:: compute pod/local
.. index:: compute pod/global
compute pod/atom command
========================
compute podd/atom command
=========================
compute pod/local command
=========================
compute pod/global command
==========================
Syntax
""""""
.. code-block:: LAMMPS
compute ID group-ID pod/atom param.pod coefficients.pod
compute ID group-ID podd/atom param.pod coefficients.pod
compute ID group-ID pod/local param.pod coefficients.pod
compute ID group-ID pod/global param.pod coefficients.pod
* ID, group-ID are documented in :doc:`compute <compute>` command
* pod/atom = style name of this compute command
* param.pod = the parameter file specifies parameters of the POD descriptors
* coefficients.pod = the coefficient file specifies coefficients of the POD potential
Examples
""""""""
.. code-block:: LAMMPS
compute d all pod/atom Ta_param.pod
compute dd all podd/atom Ta_param.pod
compute ldd all pod/local Ta_param.pod
compute gdd all podd/global Ta_param.pod
compute d all pod/atom Ta_param.pod Ta_coefficients.pod
compute dd all podd/atom Ta_param.pod Ta_coefficients.pod
compute ldd all pod/local Ta_param.pod Ta_coefficients.pod
compute gdd all podd/global Ta_param.pod Ta_coefficients.pod
Description
"""""""""""
.. versionadded:: 27June2024
Define a computation that calculates a set of quantities related to the
POD descriptors of the atoms in a group. These computes are used
primarily for calculating the dependence of energy and force components
on the linear coefficients in the :doc:`pod pair_style <pair_pod>`,
which is useful when training a POD potential to match target data. POD
descriptors of an atom are characterized by the radial and angular
distribution of neighbor atoms. The detailed mathematical definition is
given in the papers by :ref:`(Nguyen and Rohskopf) <Nguyen20222c>`,
:ref:`(Nguyen2023) <Nguyen20232c>`, :ref:`(Nguyen2024) <Nguyen20242c>`,
and :ref:`(Nguyen and Sema) <Nguyen20243c>`.
Compute *pod/atom* calculates the per-atom POD descriptors.
Compute *podd/atom* calculates derivatives of the per-atom POD
descriptors with respect to atom positions.
Compute *pod/local* calculates the per-atom POD descriptors and their
derivatives with respect to atom positions.
Compute *pod/global* calculates the global POD descriptors and their
derivatives with respect to atom positions.
Examples how to use Compute POD commands are found in the directory
``examples/PACKAGES/pod``.
.. warning::
All of these compute styles produce *very* large per-atom output
arrays that scale with the total number of atoms in the system.
This will result in *very* large memory consumption for systems
with a large number of atoms.
----------
Output info
"""""""""""
Compute *pod/atom* produces an 2D array of size :math:`N \times M`,
where :math:`N` is the number of atoms and :math:`M` is the number of
descriptors. Each column corresponds to a particular POD descriptor.
Compute *podd/atom* produces an 2D array of size :math:`N \times (M * 3
N)`. Each column corresponds to a particular derivative of a POD
descriptor.
Compute *pod/local* produces an 2D array of size :math:`(1 + 3N) \times
(M * N)`. The first row contains the per-atom descriptors, and the last
3N rows contain the derivatives of the per-atom descriptors with respect
to atom positions.
Compute *pod/global* produces an 2D array of size :math:`(1 + 3N) \times
(M)`. The first row contains the global descriptors, and the last 3N
rows contain the derivatives of the global descriptors with respect to
atom positions.
Restrictions
""""""""""""
These computes are part of the ML-POD package. They are only enabled
if LAMMPS was built with that package. See the :doc:`Build package
<Build_package>` page for more info.
Related commands
""""""""""""""""
:doc:`fitpod <fitpod_command>`,
:doc:`pair_style pod <pair_pod>`
Default
"""""""
none
----------
.. _Nguyen20222c:
**(Nguyen and Rohskopf)** Nguyen and Rohskopf, Journal of Computational Physics, 480, 112030, (2023).
.. _Nguyen20232c:
**(Nguyen2023)** Nguyen, Physical Review B, 107(14), 144103, (2023).
.. _Nguyen20242c:
**(Nguyen2024)** Nguyen, Journal of Computational Physics, 113102, (2024).
.. _Nguyen20243c:
**(Nguyen and Sema)** Nguyen and Sema, https://arxiv.org/abs/2405.00306, (2024).

22
doc/src/compute_rdf.rst
View File

@ -13,8 +13,8 @@ Syntax
* ID, group-ID are documented in :doc:`compute <compute>` command
* rdf = style name of this compute command
* Nbin = number of RDF bins
* itypeN = central atom type for Nth RDF histogram (see asterisk form below)
* jtypeN = distribution atom type for Nth RDF histogram (see asterisk form below)
* itypeN = central atom type for Nth RDF histogram (integer, type label, or asterisk form)
* jtypeN = distribution atom type for Nth RDF histogram (integer, type label, or asterisk form)
* zero or more keyword/value pairs may be appended
* keyword = *cutoff*
@ -96,14 +96,16 @@ is computed for :math:`g(r)` between all atom types. If one or more pairs are
listed, then a separate histogram is generated for each
*itype*,\ *jtype* pair.
The *itypeN* and *jtypeN* settings can be specified in one of two
ways. An explicit numeric value can be used, as in the fourth example
above. Or a wild-card asterisk can be used to specify a range of atom
types. This takes the form "\*" or "\*n" or "m\*" or "m\*n". If
:math:`N` is the number of atom types, then an asterisk with no numeric values
means all types from 1 to :math:`N`. A leading asterisk means all types from 1
to n (inclusive). A trailing asterisk means all types from m to :math:`N`
(inclusive). A middle asterisk means all types from m to n (inclusive).
The *itypeN* and *jtypeN* settings can be specified in one of three
ways. One or both of the types in the I,J pair can be a
:doc:`type label <Howto_type_labels>`. Or an explicit numeric value can be
used, as in the fourth example above. Or a wild-card asterisk can be used
to specify a range of atom types. This takes the form "\*" or "\*n" or
"m\*" or "m\*n". If :math:`N` is the number of atom types, then an asterisk
with no numeric values means all types from 1 to :math:`N`. A leading
asterisk means all types from 1 to n (inclusive). A trailing asterisk
means all types from m to :math:`N` (inclusive). A middle asterisk means
all types from m to n (inclusive).
If both *itypeN* and *jtypeN* are single values, as in the fourth example
above, this means that a :math:`g(r)` is computed where atoms of type *itypeN*

11
doc/src/create_atoms.rst
View File

@ -10,7 +10,7 @@ Syntax
create_atoms type style args keyword values ...
* type = atom type (1-Ntypes) of atoms to create (offset for molecule creation)
* type = atom type (1-Ntypes or type label) of atoms to create (offset for molecule creation)
* style = *box* or *region* or *single* or *mesh* or *random*
.. parsed-literal::
@ -37,7 +37,7 @@ Syntax
seed = random # seed (positive integer)
*basis* values = M itype
M = which basis atom
itype = atom type (1-N) to assign to this basis atom
itype = atom type (1-Ntypes or type label) to assign to this basis atom
*ratio* values = frac seed
frac = fraction of lattice sites (0 to 1) to populate randomly
seed = random # seed (positive integer)
@ -74,6 +74,13 @@ Examples
.. code-block:: LAMMPS
create_atoms 1 box
labelmap atom 1 Pt
create_atoms Pt box
labelmap atom 1 C 2 Si
create_atoms C region regsphere basis Si C
create_atoms 3 region regsphere basis 2 3
create_atoms 3 region regsphere basis 2 3 ratio 0.5 74637
create_atoms 3 single 0 0 5

24
doc/src/delete_bonds.rst
View File

@ -43,6 +43,9 @@ Examples
delete_bonds all bond 0*3 special
delete_bonds all stats
labelmap atom 4 hc
delete_bonds all atom hc special
Description
"""""""""""
@ -59,19 +62,20 @@ For all styles, by default, an interaction is only turned off (or on)
if all the atoms involved are in the specified group. See the *any*
keyword to change the behavior.
Several of the styles (\ *atom*, *bond*, *angle*, *dihedral*,
*improper*\ ) take a *type* as an argument. The specified *type* should
be an integer from 0 to :math:`N`, where :math:`N` is the number of relevant
Several of the styles (\ *atom*, *bond*, *angle*, *dihedral*, *improper*\ )
take a *type* as an argument. The specified *type* can be a
:doc:`type label <Howto_type_labels>`. Otherwise, the type should be an
integer from 0 to :math:`N`, where :math:`N` is the number of relevant
types (atom types, bond types, etc.). A value of 0 is only relevant for
style *bond*\ ; see details below. In all cases, a wildcard asterisk
style *bond*\ ; see details below. For numeric types, a wildcard asterisk
can be used in place of or in conjunction with the *type* argument to
specify a range of types. This takes the form "\*" or "\*n" or "m\*" or
"m\*n". If :math:`N` is the number of types, then an asterisk with no numeric
values means all types from 0 to :math:`N`. A leading asterisk means all
types from 0 to n (inclusive). A trailing asterisk means all types
from m to N (inclusive). A middle asterisk means all types from m to
n (inclusive). Note that it is fine to include a type of 0 for
non-bond styles; it will simply be ignored.
"m\*n". If :math:`N` is the number of types, then an asterisk with no
numeric values means all types from 0 to :math:`N`. A leading asterisk
means all types from 0 to n (inclusive). A trailing asterisk means all
types from m to N (inclusive). A middle asterisk means all types from m to
n (inclusive). Note that it is fine to include a type of 0 for non-bond
styles; it will simply be ignored.
For style *multi* all bond, angle, dihedral, and improper interactions
of any type, involving atoms in the group, are turned off.

28
doc/src/dump.rst
View File

@ -114,6 +114,7 @@ Syntax
proc = ID of processor that owns atom
procp1 = ID+1 of processor that owns atom
type = atom type
typelabel = atom :doc:`type label <Howto_type_labels>`
element = name of atom element, as defined by :doc:`dump_modify <dump_modify>` command
mass = atom mass
x,y,z = unscaled atom coordinates
@ -470,8 +471,9 @@ followed by one line per atom with the atom type and the :math:`x`-,
:math:`y`-, and :math:`z`-coordinate of that atom. You can use the
:doc:`dump_modify element <dump_modify>` option to change the output
from using the (numerical) atom type to an element name (or some other
label). This will help many visualization programs to guess bonds and
colors.
label). This option will help many visualization programs to guess bonds
and colors. You can use the :doc:`dump_modify types labels <dump_modify>`
option to replace numeric atom types with :doc:`type labels <Howto_type_labels>`.
.. versionadded:: 22Dec2022
@ -774,21 +776,21 @@ command creates a per-atom array with six columns:
Per-atom attributes used as arguments to the *custom* and *cfg* styles:
The *id*, *mol*, *proc*, *procp1*, *type*, *element*, *mass*, *vx*,
*vy*, *vz*, *fx*, *fy*, *fz*, *q* attributes are self-explanatory.
The *id*, *mol*, *proc*, *procp1*, *type*, *typelabel*, *element*, *mass*,
*vx*, *vy*, *vz*, *fx*, *fy*, *fz*, *q* attributes are self-explanatory.
*Id* is the atom ID. *Mol* is the molecule ID, included in the data
file for molecular systems. *Proc* is the ID of the processor (0 to
*Id* is the atom ID. *Mol* is the molecule ID, included in the data file
for molecular systems. *Proc* is the ID of the processor (0 to
:math:`N_\text{procs}-1`) that currently owns the atom. *Procp1* is the
proc ID+1, which can be convenient in place of a *type* attribute (1 to
:math:`N_\text{types}`) for coloring atoms in a visualization program.
*Type* is the atom type (1 to :math:`N_\text{types}`). *Element* is
typically the chemical name of an element, which you must assign to each
type via the :doc:`dump_modify element <dump_modify>` command. More
generally, it can be any string you wish to associated with an atom
type. *Mass* is the atom mass. The quantities *vx*, *vy*, *vz*, *fx*,
*fy*, *fz*, and *q* are components of atom velocity and force and atomic
charge.
*Type* is the atom type (1 to :math:`N_\text{types}`). *Typelabel* is the
atom :doc:`type label <Howto_type_labels>`. *Element* is typically the
chemical name of an element, which you must assign to each type via the
:doc:`dump_modify element <dump_modify>` command. More generally, it can
be any string you wish to associated with an atom type. *Mass* is the atom
mass. The quantities *vx*, *vy*, *vz*, *fx*, *fy*, *fz*, and *q* are
components of atom velocity and force and atomic charge.
There are several options for outputting atom coordinates. The *x*,
*y*, and *z* attributes write atom coordinates "unscaled", in the

14
doc/src/dump_modify.rst
View File

@ -17,7 +17,7 @@ Syntax
* one or more keyword/value pairs may be appended
* these keywords apply to various dump styles
* keyword = *append* or *at* or *balance* or *buffer* or *colname* or *delay* or *element* or *every* or *every/time* or *fileper* or *first* or *flush* or *format* or *header* or *image* or *label* or *maxfiles* or *nfile* or *pad* or *pbc* or *precision* or *region* or *refresh* or *scale* or *sfactor* or *skip* or *sort* or *tfactor* or *thermo* or *thresh* or *time* or *triclinic/general* or *units* or *unwrap*
* keyword = *append* or *at* or *balance* or *buffer* or *colname* or *delay* or *element* or *every* or *every/time* or *fileper* or *first* or *flush* or *format* or *header* or *image* or *label* or *maxfiles* or *nfile* or *pad* or *pbc* or *precision* or *region* or *refresh* or *scale* or *sfactor* or *skip* or *sort* or *tfactor* or *thermo* or *thresh* or *time* or *triclinic/general* or *types* or *units* or *unwrap*
.. parsed-literal::
@ -81,6 +81,7 @@ Syntax
these 3 args can be replaced by the word "none" to turn off thresholding
*time* arg = *yes* or *no*
*triclinic/general* arg = *yes* or *no*
*types* value = *numeric* or *labels*
*units* arg = *yes* or *no*
*unwrap* arg = *yes* or *no*
@ -849,6 +850,13 @@ The default setting is *no*\ .
----------
The *types* keyword applies only to the dump xyz style. If this keyword is
used with a value of *numeric*, then numeric atom types are printed in the
xyz file (default). If the value *labels* is specified, then
:doc:`type labels <Howto_type_labels>` are printed for atom types.
----------
The *triclinic/general* keyword only applies to the dump *atom* and
*custom* styles. It can only be used with a value of *yes* if the
simulation box was created as a general triclinic box. See the
@ -960,11 +968,11 @@ The option defaults are
* sort = id for dump styles *dcd*, *xtc*, and *xyz*
* thresh = none
* time = no
* triclinic/general no
* triclinic/general = no
* types = numeric
* units = no
* unwrap = no
* compression_level = 9 (gz variants)
* compression_level = 0 (zstd variants)
* checksum = yes (zstd variants)

741
doc/src/fitpod_command.rst
View File

@ -1,18 +1,19 @@
.. index:: fitpod
fitpod command
======================
==============
Syntax
""""""
.. code-block:: LAMMPS
fitpod Ta_param.pod Ta_data.pod
fitpod Ta_param.pod Ta_data.pod Ta_coefficients.pod
* fitpod = style name of this command
* Ta_param.pod = an input file that describes proper orthogonal descriptors (PODs)
* Ta_data.pod = an input file that specifies DFT data used to fit a POD potential
* Ta_coefficients.pod (optional) = an input file that specifies trainable coefficients of a POD potential
Examples
""""""""
@ -20,20 +21,26 @@ Examples
.. code-block:: LAMMPS
fitpod Ta_param.pod Ta_data.pod
fitpod Ta_param.pod Ta_data.pod Ta_coefficients.pod
Description
"""""""""""
.. versionadded:: 22Dec2022
Fit a machine-learning interatomic potential (ML-IAP) based on proper
orthogonal descriptors (POD). Two input files are required for this
command. The first input file describes a POD potential parameter
settings, while the second input file specifies the DFT data used for
the fitting procedure.
orthogonal descriptors (POD); please see :ref:`(Nguyen and Rohskopf)
<Nguyen20222a>`, :ref:`(Nguyen2023) <Nguyen20232a>`, :ref:`(Nguyen2024)
<Nguyen20242a>`, and :ref:`(Nguyen and Sema) <Nguyen20243a>` for details.
The fitted POD potential can be used to run MD simulations via
:doc:`pair_style pod <pair_pod>`.
The table below has one-line descriptions of all the keywords that can
be used in the first input file (i.e. ``Ta_param.pod`` in the example
above):
Two input files are required for this command. The first input file
describes a POD potential parameter settings, while the second input
file specifies the DFT data used for the fitting procedure. All keywords
except *species* have default values. If a keyword is not set in the
input file, its default value is used. The table below has one-line
descriptions of all the keywords that can be used in the first input
file (i.e. ``Ta_param.pod``)
.. list-table::
:header-rows: 1
@ -52,7 +59,7 @@ above):
- INT
- three integer constants specify boundary conditions
* - rin
- 1.0
- 0.5
- REAL
- a real number specifies the inner cut-off radius
* - rcut
@ -60,46 +67,75 @@ above):
- REAL
- a real number specifies the outer cut-off radius
* - bessel_polynomial_degree
- 3
- 4
- INT
- the maximum degree of Bessel polynomials
* - inverse_polynomial_degree
- 6
- 8
- INT
- the maximum degree of inverse radial basis functions
* - number_of_environment_clusters
- 1
- INT
- the number of clusters for environment-adaptive potentials
* - number_of_principal_components
- 2
- INT
- the number of principal components for dimensionality reduction
* - onebody
- 1
- BOOL
- turns on/off one-body potential
* - twobody_number_radial_basis_functions
- 6
- 8
- INT
- number of radial basis functions for two-body potential
* - threebody_number_radial_basis_functions
- 5
- 6
- INT
- number of radial basis functions for three-body potential
* - threebody_number_angular_basis_functions
* - threebody_angular_degree
- 5
- INT
- number of angular basis functions for three-body potential
* - fourbody_snap_twojmax
- angular degree for three-body potential
* - fourbody_number_radial_basis_functions
- 4
- INT
- number of radial basis functions for four-body potential
* - fourbody_angular_degree
- 3
- INT
- angular degree for four-body potential
* - fivebody_number_radial_basis_functions
- 0
- INT
- band limit for SNAP bispectrum components (0,2,4,6,8... allowed)
* - fourbody_snap_chemflag
- number of radial basis functions for five-body potential
* - fivebody_angular_degree
- 0
- BOOL
- turns on/off the explicit multi-element variant of the SNAP bispectrum components
* - quadratic_pod_potential
- INT
- angular degree for five-body potential
* - sixbody_number_radial_basis_functions
- 0
- BOOL
- turns on/off quadratic POD potential
- INT
- number of radial basis functions for six-body potential
* - sixbody_angular_degree
- 0
- INT
- angular degree for six-body potential
* - sevenbody_number_radial_basis_functions
- 0
- INT
- number of radial basis functions for seven-body potential
* - sevenbody_angular_degree
- 0
- INT
- angular degree for seven-body potential
Note that both the number of radial basis functions and angular degree
must decrease as the body order increases. The next table describes all
keywords that can be used in the second input file (i.e. ``Ta_data.pod``
in the example above):
All keywords except *species* have default values. If a keyword is not
set in the input file, its default value is used. The next table
describes all keywords that can be used in the second input file
(i.e. ``Ta_data.pod`` in the example above):
.. list-table::
:header-rows: 1
@ -125,6 +161,10 @@ describes all keywords that can be used in the second input file
- ""
- STRING
- specifies the path to test data files in double quotes
* - path_to_environment_configuration_set
- ""
- STRING
- specifies the path to environment configuration files in double quotes
* - fraction_training_data_set
- 1.0
- REAL
@ -133,6 +173,14 @@ describes all keywords that can be used in the second input file
- 0
- BOOL
- turns on/off randomization of the training set
* - fraction_test_data_set
- 1.0
- REAL
- a real number (<= 1.0) specifies the fraction of the test set used to validate POD
* - randomize_test_data_set
- 0
- BOOL
- turns on/off randomization of the test set
* - fitting_weight_energy
- 100.0
- REAL
@ -161,6 +209,10 @@ describes all keywords that can be used in the second input file
- 8
- INT
- number of digits after the decimal points for numbers in the coefficient file
* - group_weights
- global
- STRING
- ``table`` uses group weights defined for each group named by filename
All keywords except *path_to_training_data_set* have default values. If
a keyword is not set in the input file, its default value is used. After
@ -172,14 +224,44 @@ successful training, a number of output files are produced, if enabled:
* ``<basename>_test_analysis.pod`` reports detailed errors for all test configurations
* ``<basename>_coefficients.pod`` contains the coefficients of the POD potential
After training the POD potential, ``Ta_param.pod`` and ``<basename>_coefficients.pod``
are the two files needed to use the POD potential in LAMMPS. See
:doc:`pair_style pod <pair_pod>` for using the POD potential. Examples
about training and using POD potentials are found in the directory
lammps/examples/PACKAGES/pod.
After training the POD potential, ``Ta_param.pod`` and
``<basename>_coefficients.pod`` are the two files needed to use the POD
potential in LAMMPS. See :doc:`pair_style pod <pair_pod>` for using the
POD potential. Examples about training and using POD potentials are
found in the directory lammps/examples/PACKAGES/pod and the Github repo
https://github.com/cesmix-mit/pod-examples.
Parameterized Potential Energy Surface
""""""""""""""""""""""""""""""""""""""
Loss Function Group Weights
^^^^^^^^^^^^^^^^^^^^^^^^^^^
The *group_weights* keyword in the ``data.pod`` file is responsible for
weighting certain groups of configurations in the loss function. For
example:
.. code-block:: LAMMPS
group_weights table
Displaced_A15 100.0 1.0
Displaced_BCC 100.0 1.0
Displaced_FCC 100.0 1.0
Elastic_BCC 100.0 1.0
Elastic_FCC 100.0 1.0
GSF_110 100.0 1.0
GSF_112 100.0 1.0
Liquid 100.0 1.0
Surface 100.0 1.0
Volume_A15 100.0 1.0
Volume_BCC 100.0 1.0
Volume_FCC 100.0 1.0
This will apply an energy weight of ``100.0`` and a force weight of
``1.0`` for all groups in the ``Ta`` example. The groups are named by
their respective filename. If certain groups are left out of this table,
then the globally defined weights from the ``fitting_weight_energy`` and
``fitting_weight_force`` keywords will be used.
POD Potential
"""""""""""""
We consider a multi-element system of *N* atoms with :math:`N_{\rm e}`
unique elements. We denote by :math:`\boldsymbol r_n` and :math:`Z_n`
@ -187,535 +269,82 @@ position vector and type of an atom *n* in the system,
respectively. Note that we have :math:`Z_n \in \{1, \ldots, N_{\rm e}
\}`, :math:`\boldsymbol R = (\boldsymbol r_1, \boldsymbol r_2, \ldots,
\boldsymbol r_N) \in \mathbb{R}^{3N}`, and :math:`\boldsymbol Z = (Z_1,
Z_2, \ldots, Z_N) \in \mathbb{N}^{N}`. The potential energy surface
(PES) of the system can be expressed as a many-body expansion of the
form
Z_2, \ldots, Z_N) \in \mathbb{N}^{N}`. The total energy of the
POD potential is expressed as :math:`E(\boldsymbol R, \boldsymbol Z) =
\sum_{i=1}^N E_i(\boldsymbol R_i, \boldsymbol Z_i)`, where
.. math::
E(\boldsymbol R, \boldsymbol Z, \boldsymbol{\eta}, \boldsymbol{\mu}) \ = \ & \sum_{i} V^{(1)}(\boldsymbol r_i, Z_i, \boldsymbol \mu^{(1)} ) + \frac12 \sum_{i,j} V^{(2)}(\boldsymbol r_i, \boldsymbol r_j, Z_i, Z_j, \boldsymbol \eta, \boldsymbol \mu^{(2)}) \\
& + \frac16 \sum_{i,j,k} V^{(3)}(\boldsymbol r_i, \boldsymbol r_j, \boldsymbol r_k, Z_i, Z_j, Z_k, \boldsymbol \eta, \boldsymbol \mu^{(3)}) + \ldots
where :math:`V^{(1)}` is the one-body potential often used for
representing external field or energy of isolated elements, and the
higher-body potentials :math:`V^{(2)}, V^{(3)}, \ldots` are symmetric,
uniquely defined, and zero if two or more indices take identical values.
The superscript on each potential denotes its body order. Each *q*-body
potential :math:`V^{(q)}` depends on :math:`\boldsymbol \mu^{(q)}` which
are sets of parameters to fit the PES. Note that :math:`\boldsymbol \mu`
is a collection of all potential parameters :math:`\boldsymbol
\mu^{(1)}`, :math:`\boldsymbol \mu^{(2)}`, :math:`\boldsymbol
\mu^{(3)}`, etc, and that :math:`\boldsymbol \eta` is a set of
hyper-parameters such as inner cut-off radius :math:`r_{\rm in}` and
outer cut-off radius :math:`r_{\rm cut}`.
Interatomic potentials rely on parameters to learn relationship between
atomic environments and interactions. Since interatomic potentials are
approximations by nature, their parameters need to be set to some
reference values or fitted against data by necessity. Typically,
potential fitting finds optimal parameters, :math:`\boldsymbol \mu^*`,
to minimize a certain loss function of the predicted quantities and
data. Since the fitted potential depends on the data set used to fit it,
different data sets will yield different optimal parameters and thus
different fitted potentials. When fitting the same functional form on
*Q* different data sets, we would obtain *Q* different optimized
potentials, :math:`E(\boldsymbol R,\boldsymbol Z, \boldsymbol \eta,
\boldsymbol \mu_q^*), 1 \le q \le Q`. Consequently, there exist many
different sets of optimized parameters for empirical interatomic
potentials.
Instead of optimizing the potential parameters, inspired by the reduced
basis method :ref:`(Grepl) <Grepl20072>` for parameterized partial
differential equations, we view the parameterized PES as a parametric
manifold of potential energies
.. math::
\mathcal{M} = \{E(\boldsymbol R, \boldsymbol Z, \boldsymbol \eta, \boldsymbol \mu) \ | \ \boldsymbol \mu \in \Omega^{\boldsymbol \mu} \}
where :math:`\Omega^{\boldsymbol \mu}` is a parameter domain in which
:math:`\boldsymbol \mu` resides. The parametric manifold
:math:`\mathcal{M}` contains potential energy surfaces for all values of
:math:`\boldsymbol \mu \in \Omega^{\boldsymbol \mu}`. Therefore, the
parametric manifold yields a much richer and more transferable atomic
representation than any particular individual PES :math:`E(\boldsymbol
R, \boldsymbol Z, \boldsymbol \eta, \boldsymbol \mu^*)`.
We propose specific forms of the parameterized potentials for one-body,
two-body, and three-body interactions. We apply the Karhunen-Loeve
expansion to snapshots of the parameterized potentials to obtain sets of
orthogonal basis functions. These basis functions are aggregated
according to the chemical elements of atoms, thus leading to
multi-element proper orthogonal descriptors.
Proper Orthogonal Descriptors
"""""""""""""""""""""""""""""
Proper orthogonal descriptors are finger prints characterizing the
radial and angular distribution of a system of atoms. The detailed
mathematical definition is given in the paper by Nguyen and Rohskopf
:ref:`(Nguyen) <Nguyen20222>`.
The descriptors for the one-body interaction are used to capture energy
of isolated elements and defined as follows
.. math::
D_{ip}^{(1)} = \left\{
\begin{array}{ll}
1, & \mbox{if } Z_i = p \\
0, & \mbox{if } Z_i \neq p
\end{array}
\right.
for :math:`1 \le i \le N, 1 \le p \le N_{\rm e}`. The number of one-body
descriptors per atom is equal to the number of elements. The one-body
descriptors are independent of atom positions, but dependent on atom
types. The one-body descriptors are active only when the keyword
*onebody* is set to 1.
We adopt the usual assumption that the direct interaction between two
atoms vanishes smoothly when their distance is greater than the outer
cutoff distance :math:`r_{\rm cut}`. Furthermore, we assume that two
atoms can not get closer than the inner cutoff distance :math:`r_{\rm
in}` due to the Pauli repulsion principle. Let :math:`r \in (r_{\rm in},
r_{\rm cut})`, we introduce the following parameterized radial functions
.. math::
\phi(r, r_{\rm in}, r_{\rm cut}, \alpha, \beta) = \frac{\sin (\alpha \pi x) }{r - r_{\rm in}}, \qquad \varphi(r, \gamma) = \frac{1}{r^\gamma} ,
where the scaled distance function :math:`x` is defined below to enrich the two-body manifold
.. math::
x(r, r_{\rm in}, r_{\rm cut}, \beta) = \frac{e^{-\beta(r - r_{\rm in})/(r_{\rm cut} - r_{\rm in})} - 1}{e^{-\beta} - 1} .
We introduce the following function as a convex combination of the two functions
.. math::
\psi(r, r_{\rm in}, r_{\rm cut}, \alpha, \beta, \gamma, \kappa) = \kappa \phi(r, r_{\rm in}, r_{\rm cut}, \alpha, \beta) + (1- \kappa) \varphi(r, \gamma) .
We see that :math:`\psi` is a function of distance :math:`r`, cut-off
distances :math:`r_{\rm in}` and :math:`r_{\rm cut}`, and parameters
:math:`\alpha, \beta, \gamma, \kappa`. Together these parameters allow
the function :math:`\psi` to characterize a diverse spectrum of two-body
interactions within the cut-off interval :math:`(r_{\rm in}, r_{\rm
cut})`.
Next, we introduce the following parameterized potential
.. math::
W^{(2)}(r_{ij}, \boldsymbol \eta, \boldsymbol \mu^{(2)}) = f_{\rm c}(r_{ij}, \boldsymbol \eta) \psi(r_{ij}, \boldsymbol \eta, \boldsymbol \mu^{(2)})
where :math:`\eta_1 = r_{\rm in}, \eta_2 = r_{\rm cut}, \mu_1^{(2)} =
\alpha, \mu_2^{(2)} = \beta, \mu_3^{(2)} = \gamma`, and
:math:`\mu_4^{(2)} = \kappa`. Here the cut-off function :math:`f_{\rm
c}(r_{ij}, \boldsymbol \eta)` proposed in [refs] is used to ensure the
smooth vanishing of the potential and its derivative for :math:`r_{ij}
\ge r_{\rm cut}`:
.. math::
f_{\rm c}(r_{ij}, r_{\rm in}, r_{\rm cut}) = \exp \left(1 -\frac{1}{\sqrt{\left(1 - \frac{(r-r_{\rm in})^3}{(r_{\rm cut} - r_{\rm in})^3} \right)^2 + 10^{-6}}} \right)
Based on the parameterized potential, we form a set of snapshots as
follows. We assume that we are given :math:`N_{\rm s}` parameter tuples
:math:`\boldsymbol \mu^{(2)}_\ell, 1 \le \ell \le N_{\rm s}`. We
introduce the following set of snapshots on :math:`(r_{\rm in}, r_{\rm
cut})`:
.. math::
\xi_\ell(r_{ij}, \boldsymbol \eta) = W^{(2)}(r_{ij}, \boldsymbol \eta, \boldsymbol \mu^{(2)}_\ell), \quad \ell = 1, \ldots, N_{\rm s} .
To ensure adequate sampling of the PES for different parameters, we
choose :math:`N_{\rm s}` parameter points :math:`\boldsymbol
\mu^{(2)}_\ell = (\alpha_\ell, \beta_\ell, \gamma_\ell, \kappa_\ell), 1
\le \ell \le N_{\rm s}` as follows. The parameters :math:`\alpha \in [1,
N_\alpha]` and :math:`\gamma \in [1, N_\gamma]` are integers, where
:math:`N_\alpha` and :math:`N_\gamma` are the highest degrees for
:math:`\alpha` and :math:`\gamma`, respectively. We next choose
:math:`N_\beta` different values of :math:`\beta` in the interval
:math:`[\beta_{\min}, \beta_{\max}]`, where :math:`\beta_{\min} = 0` and
:math:`\beta_{\max} = 4`. The parameter :math:`\kappa` can be set either
0 or 1. Hence, the total number of parameter points is :math:`N_{\rm s}
= N_\alpha N_\beta + N_\gamma`. Although :math:`N_\alpha, N_\beta,
N_\gamma` can be chosen conservatively large, we find that
:math:`N_\alpha = 6, N_\beta = 3, N_\gamma = 8` are adequate for most
problems. Note that :math:`N_\alpha` and :math:`N_\gamma` correspond to
*bessel_polynomial_degree* and *inverse_polynomial_degree*,
respectively.
We employ the Karhunen-Loeve (KL) expansion to generate an orthogonal
basis set which is known to be optimal for representation of the
snapshot family :math:`\{\xi_\ell\}_{\ell=1}^{N_{\rm s}}`. The two-body
orthogonal basis functions are computed as follows
.. math::
U^{(2)}_m(r_{ij}, \boldsymbol \eta) = \sum_{\ell = 1}^{N_{\rm s}} A_{\ell m}(\boldsymbol \eta) \, \xi_\ell(r_{ij}, \boldsymbol \eta), \qquad m = 1, \ldots, N_{\rm 2b} ,
where the matrix :math:`\boldsymbol A \in \mathbb{R}^{N_{\rm s} \times
N_{\rm s}}` consists of eigenvectors of the eigenvalue problem
.. math::
\boldsymbol C \boldsymbol a = \lambda \boldsymbol a
with the entries of :math:`\boldsymbol C \in \mathbb{R}^{N_{\rm s} \times N_{\rm s}}` being given by
.. math::
C_{ij} = \frac{1}{N_{\rm s}} \int_{r_{\rm in}}^{r_{\rm cut}} \xi_i(x, \boldsymbol \eta) \xi_j(x, \boldsymbol \eta) dx, \quad 1 \le i, j \le N_{\rm s}
Note that the eigenvalues :math:`\lambda_\ell, 1 \le \ell \le N_{\rm
s}`, are ordered such that :math:`\lambda_1 \ge \lambda_2 \ge \ldots \ge
\lambda_{N_{\rm s}}`, and that the matrix :math:`\boldsymbol A` is
pe-computed and stored for any given :math:`\boldsymbol \eta`. Owing to
the rapid convergence of the KL expansion, only a small number of
orthogonal basis functions is needed to obtain accurate
approximation. The value of :math:`N_{\rm 2b}` corresponds to
*twobody_number_radial_basis_functions*.
The two-body proper orthogonal descriptors at each atom *i* are computed
by summing the orthogonal basis functions over the neighbors of atom *i*
and numerating on the atom types as follows
.. math::
D^{(2)}_{im l(p, q) }(\boldsymbol \eta) = \left\{
\begin{array}{ll}
\displaystyle \sum_{\{j | Z_j = q\}} U^{(2)}_m(r_{ij}, \boldsymbol \eta), & \mbox{if } Z_i = p \\
0, & \mbox{if } Z_i \neq p
\end{array}
\right.
for :math:`1 \le i \le N, 1 \le m \le N_{\rm 2b}, 1 \le q, p \le N_{\rm
e}`. Here :math:`l(p,q)` is a symmetric index mapping such that
.. math::
l(p,q) = \left\{
\begin{array}{ll}
q + (p-1) N_{\rm e} - p(p-1)/2, & \mbox{if } q \ge p \\
p + (q-1) N_{\rm e} - q(q-1)/2, & \mbox{if } q < p .
\end{array}
\right.
The number of two-body descriptors per atom is thus :math:`N_{\rm 2b}
N_{\rm e}(N_{\rm e}+1)/2`.
It is important to note that the orthogonal basis functions do not
depend on the atomic numbers :math:`Z_i` and :math:`Z_j`. Therefore, the
cost of evaluating the basis functions and their derivatives with
respect to :math:`r_{ij}` is independent of the number of elements
:math:`N_{\rm e}`. Consequently, even though the two-body proper
orthogonal descriptors depend on :math:`\boldsymbol Z`, their
computational complexity is independent of :math:`N_{\rm e}`.
In order to provide proper orthogonal descriptors for three-body
interactions, we need to introduce a three-body parameterized
potential. In particular, the three-body potential is defined as a
product of radial and angular functions as follows
.. math::
W^{(3)}(r_{ij}, r_{ik}, \theta_{ijk}, \boldsymbol \eta, \boldsymbol \mu^{(3)}) = \psi(r_{ij}, r_{\rm min}, r_{\rm max}, \alpha, \beta, \gamma, \kappa) f_{\rm c}(r_{ij}, r_{\rm min}, r_{\rm max}) \\
\psi(r_{ik}, r_{\rm min}, r_{\rm max}, \alpha, \beta, \gamma, \kappa) f_{\rm c}(r_{ik}, r_{\rm min}, r_{\rm max}) \\
\cos (\sigma \theta_{ijk} + \zeta)
where :math:`\sigma` is the periodic multiplicity, :math:`\zeta` is the
equilibrium angle, :math:`\boldsymbol \mu^{(3)} = (\alpha, \beta,
\gamma, \kappa, \sigma, \zeta)`. The three-body potential provides an
angular fingerprint of the atomic environment through the bond angles
:math:`\theta_{ijk}` formed with each pair of neighbors :math:`j` and
:math:`k`. Compared to the two-body potential, the three-body potential
has two extra parameters :math:`(\sigma, \zeta)` associated with the
angular component.
Let :math:`\boldsymbol \varrho = (\alpha, \beta, \gamma, \kappa)`. We
assume that we are given :math:`L_{\rm r}` parameter tuples
:math:`\boldsymbol \varrho_\ell, 1 \le \ell \le L_{\rm r}`. We
introduce the following set of snapshots on :math:`(r_{\min},
r_{\max})`:
.. math::
\zeta_\ell(r_{ij}, r_{\rm min}, r_{\rm max} ) = \psi(r_{ij}, r_{\rm min}, r_{\rm max}, \boldsymbol \varrho_\ell) f_{\rm c}(r_{ij}, r_{\rm min}, r_{\rm max}), \quad 1 \le \ell \le L_{\rm r} .
We apply the Karhunen-Loeve (KL) expansion to this set of snapshots to
obtain orthogonal basis functions as follows
.. math::
U^{r}_m(r_{ij}, r_{\rm min}, r_{\rm max} ) = \sum_{\ell = 1}^{L_{\rm r}} A_{\ell m} \, \zeta_\ell(r_{ij}, r_{\rm min}, r_{\rm max} ), \qquad m = 1, \ldots, N_{\rm r} ,
where the matrix :math:`\boldsymbol A \in \mathbb{R}^{L_{\rm r} \times L_{\rm r}}` consists
of eigenvectors of the eigenvalue problem. For the parameterized angular function,
we consider angular basis functions
.. math::
U^{a}_n(\theta_{ijk}) = \cos ((n-1) \theta_{ijk}), \qquad n = 1,\ldots, N_{\rm a},
where :math:`N_{\rm a}` is the number of angular basis functions. The orthogonal
basis functions for the parameterized potential are computed as follows
.. math::
U^{(3)}_{mn}(r_{ij}, r_{ik}, \theta_{ijk}, \boldsymbol \eta) = U^{r}_m(r_{ij}, \boldsymbol \eta) U^{r}_m(r_{ik}, \boldsymbol \eta) U^{a}_n(\theta_{ijk}),
for :math:`1 \le m \le N_{\rm r}, 1 \le n \le N_{\rm a}`. The number of three-body
orthogonal basis functions is equal to :math:`N_{\rm 3b} = N_{\rm r} N_{\rm a}` and
independent of the number of elements. The value of :math:`N_{\rm r}` corresponds to
*threebody_number_radial_basis_functions*, while that of :math:`N_{\rm a}` to
*threebody_number_angular_basis_functions*.
The three-body proper orthogonal descriptors at each atom *i*
are obtained by summing over the neighbors *j* and *k* of atom *i* as
.. math::
D^{(3)}_{imn \ell(p, q, s)}(\boldsymbol \eta) = \left\{
\begin{array}{ll}
\displaystyle \sum_{\{j | Z_j = q\}} \sum_{\{k | Z_k = s\}} U^{(3)}_{mn}(r_{ij}, r_{ik}, \theta_{ijk}, \boldsymbol \eta), & \mbox{if } Z_i = p \\
0, & \mbox{if } Z_i \neq p
\end{array}
\right.
for :math:`1 \le i \le N, 1 \le m \le N_{\rm r}, 1 \le n \le N_{\rm a}, 1 \le q, p, s \le N_{\rm e}`,
where
.. math::
\ell(p,q,s) = \left\{
\begin{array}{ll}
s + (q-1) N_{\rm e} - q(q-1)/2 + (p-1)N_{\rm e}(1+N_{\rm e})/2 , & \mbox{if } s \ge q \\
q + (s-1) N_{\rm e} - s(s-1)/2 + (p-1)N_{\rm e}(1+N_{\rm e})/2, & \mbox{if } s < q .
\end{array}
\right.
The number of three-body descriptors per atom is thus :math:`N_{\rm 3b} N_{\rm e}^2(N_{\rm e}+1)/2`.
While the number of three-body PODs is cubic function of the number of elements,
the computational complexity of the three-body PODs is independent of the number of elements.
Four-Body SNAP Descriptors
""""""""""""""""""""""""""
In addition to the proper orthogonal descriptors described above, we also employ
the spectral neighbor analysis potential (SNAP) descriptors. SNAP uses bispectrum components
to characterize the local neighborhood of each atom in a very general way. The mathematical definition
of the bispectrum calculation and its derivatives w.r.t. atom positions is described in
:doc:`compute snap <compute_sna_atom>`. In SNAP, the
total energy is decomposed into a sum over atom energies. The energy of
atom *i* is expressed as a weighted sum over bispectrum components.
.. math::
E_i^{\rm SNAP} = \sum_{k=1}^{N_{\rm 4b}} \sum_{p=1}^{N_{\rm e}} c_{kp}^{(4)} D_{ikp}^{(4)}
where the SNAP descriptors are related to the bispectrum components by
.. math::
D^{(4)}_{ikp} = \left\{
\begin{array}{ll}
\displaystyle B_{ik}, & \mbox{if } Z_i = p \\
0, & \mbox{if } Z_i \neq p
\end{array}
\right.
Here :math:`B_{ik}` is the *k*\ -th bispectrum component of atom *i*. The number of
bispectrum components :math:`N_{\rm 4b}` depends on the value of *fourbody_snap_twojmax* :math:`= 2 J_{\rm max}`
and *fourbody_snap_chemflag*. If *fourbody_snap_chemflag* = 0
then :math:`N_{\rm 4b} = (J_{\rm max}+1)(J_{\rm max}+2)(J_{\rm max}+1.5)/3`.
If *fourbody_snap_chemflag* = 1 then :math:`N_{\rm 4b} = N_{\rm e}^3 (J_{\rm max}+1)(J_{\rm max}+2)(J_{\rm max}+1.5)/3`.
The bispectrum calculation is described in more detail in :doc:`compute sna/atom <compute_sna_atom>`.
Linear Proper Orthogonal Descriptor Potentials
""""""""""""""""""""""""""""""""""""""""""""""
The proper orthogonal descriptors and SNAP descriptors are used to define the atomic energies
in the following expansion
.. math::
E_{i}(\boldsymbol \eta) = \sum_{p=1}^{N_{\rm e}} c^{(1)}_p D^{(1)}_{ip} + \sum_{m=1}^{N_{\rm 2b}} \sum_{l=1}^{N_{\rm e}(N_{\rm e}+1)/2} c^{(2)}_{ml} D^{(2)}_{iml}(\boldsymbol \eta) + \sum_{m=1}^{N_{\rm r}} \sum_{n=1}^{N_{\rm a}} \sum_{\ell=1}^{N_{\rm e}^2(N_{\rm e}+1)/2} c^{(3)}_{mn\ell} D^{(3)}_{imn\ell}(\boldsymbol \eta) + \sum_{k=1}^{N_{\rm 4b}} \sum_{p=1}^{N_{\rm e}} c_{kp}^{(4)} D_{ikp}^{(4)}(\boldsymbol \eta),
where :math:`D^{(1)}_{ip}, D^{(2)}_{iml}, D^{(3)}_{imn\ell}, D^{(4)}_{ikp}` are the one-body, two-body, three-body, four-body descriptors,
respectively, and :math:`c^{(1)}_p, c^{(2)}_{ml}, c^{(3)}_{mn\ell}, c^{(4)}_{kp}` are their respective expansion
coefficients. In a more compact notation that implies summation over descriptor indices
the atomic energies can be written as
.. math::
E_i(\boldsymbol \eta) = \sum_{m=1}^{N_{\rm e}} c^{(1)}_m D^{(1)}_{im} + \sum_{m=1}^{N_{\rm d}^{(2)}} c^{(2)}_k D^{(2)}_{im} + \sum_{m=1}^{N_{\rm d}^{(3)}} c^{(3)}_m D^{(3)}_{im} + \sum_{m=1}^{N_{\rm d}^{(4)}} c^{(4)}_m D^{(4)}_{im}
where :math:`N_{\rm d}^{(2)} = N_{\rm 2b} N_{\rm e} (N_{\rm e}+1)/2`,
:math:`N_{\rm d}^{(3)} = N_{\rm 3b} N_{\rm e}^2 (N_{\rm e}+1)/2`, and
:math:`N_{\rm d}^{(4)} = N_{\rm 4b} N_{\rm e}` are
the number of two-body, three-body, and four-body descriptors, respectively.
The potential energy is then obtained by summing local atomic energies :math:`E_i`
for all atoms :math:`i` in the system
.. math::
E(\boldsymbol \eta) = \sum_{i}^N E_{i}(\boldsymbol \eta)
Because the descriptors are one-body, two-body, and three-body terms,
the resulting POD potential is a three-body PES. We can express the potential
energy as a linear combination of the global descriptors as follows
.. math::
E(\boldsymbol \eta) = \sum_{m=1}^{N_{\rm e}} c^{(1)}_m d^{(1)}_{m} + \sum_{m=1}^{N_{\rm d}^{(2)}} c^{(2)}_m d^{(2)}_{m} + \sum_{m=1}^{N_{\rm d}^{(3)}} c^{(3)}_m d^{(3)}_{m} + \sum_{m=1}^{N_{\rm d}^{(4)}} c^{(4)}_m d^{(4)}_{m}
where the global descriptors are given by
.. math::
d_{m}^{(1)}(\boldsymbol \eta) = \sum_{i=1}^N D_{im}^{(1)}(\boldsymbol \eta), \quad d_{m}^{(2)}(\boldsymbol \eta) = \sum_{i=1}^N D_{im}^{(2)}(\boldsymbol \eta), \quad d_{m}^{(3)}(\boldsymbol \eta) = \sum_{i=1}^N D_{im}^{(3)}(\boldsymbol \eta), \quad d_{m}^{(4)}(\boldsymbol \eta) = \sum_{i=1}^N D_{im}^{(4)}(\boldsymbol \eta)
Hence, we obtain the atomic forces as
.. math::
\boldsymbol F = -\nabla E(\boldsymbol \eta) = - \sum_{m=1}^{N_{\rm d}^{(2)}} c^{(2)}_m \nabla d_m^{(2)} - \sum_{m=1}^{N_{\rm d}^{(3)}} c^{(3)}_m \nabla d_m^{(3)} - \sum_{m=1}^{N_{\rm d}^{(4)}} c^{(4)}_m \nabla d_m^{(4)}
where :math:`\nabla d_m^{(2)}`, :math:`\nabla d_m^{(3)}` and :math:`\nabla d_m^{(4)}` are derivatives of the two-body
three-body, and four-body global descriptors with respect to atom positions, respectively.
Note that since the first-body global descriptors are constant, their derivatives are zero.
Quadratic Proper Orthogonal Descriptor Potentials
"""""""""""""""""""""""""""""""""""""""""""""""""
We recall two-body PODs :math:`D^{(2)}_{ik}, 1 \le k \le N_{\rm d}^{(2)}`,
and three-body PODs :math:`D^{(3)}_{im}, 1 \le m \le N_{\rm d}^{(3)}`,
with :math:`N_{\rm d}^{(2)} = N_{\rm 2b} N_{\rm e} (N_{\rm e}+1)/2` and
:math:`N_{\rm d}^{(3)} = N_{\rm 3b} N_{\rm e}^2 (N_{\rm e}+1)/2` being
the number of descriptors per atom for the two-body PODs and three-body PODs,
respectively. We employ them to define a new set of atomic descriptors as follows
.. math::
D^{(2*3)}_{ikm} = \frac{1}{2N}\left( D^{(2)}_{ik} \sum_{j=1}^N D^{(3)}_{jm} + D^{(3)}_{im} \sum_{j=1}^N D^{(2)}_{jk} \right)
for :math:`1 \le i \le N, 1 \le k \le N_{\rm d}^{(2)}, 1 \le m \le N_{\rm d}^{(3)}`.
The new descriptors are four-body because they involve central atom :math:`i` together
with three neighbors :math:`j, k` and :math:`l`. The total number of new descriptors per atom is equal to
.. math::
N_{\rm d}^{(2*3)} = N_{\rm d}^{(2)} * N_{\rm d}^{(3)} = N_{\rm 2b} N_{\rm 3b} N_{\rm e}^3 (N_{\rm e}+1)^2/4 .
The new global descriptors are calculated as
.. math::
d^{(2*3)}_{km} = \sum_{i=1}^N D^{(2*3)}_{ikm} = \left( \sum_{i=1}^N D^{(2)}_{ik} \right) \left( \sum_{i=1}^N D^{(3)}_{im} \right) = d^{(2)}_{k} d^{(3)}_m,
for :math:`1 \le k \le N_{\rm d}^{(2)}, 1 \le m \le N_{\rm d}^{(3)}`. Hence, the gradient
of the new global descriptors with respect to atom positions is calculated as
.. math::
\nabla d^{(2*3)}_{km} = d^{(3)}_m \nabla d^{(2)}_{k} + d^{(2)}_{k} \nabla d^{(3)}_m, \quad 1 \le k \le N_{\rm d}^{(2)}, 1 \le m \le N_{\rm d}^{(3)} .
The quadratic POD potential is defined as a linear combination of the
original and new global descriptors as follows
.. math::
E^{(2*3)} = \sum_{k=1}^{N_{\rm 2d}^{(2*3)}} \sum_{m=1}^{N_{\rm 3d}^{(2*3)}} c^{(2*3)}_{km} d^{(2*3)}_{km} .
It thus follows that
.. math::
E^{(2*3)} = 0.5 \sum_{k=1}^{N_{\rm 2d}^{(2*3)}} \left( \sum_{m=1}^{N_{\rm 3d}^{(2*3)}} c^{(2*3)}_{km} d_m^{(3)} \right) d_k^{(2)} + 0.5 \sum_{m=1}^{N_{\rm 3d}^{(2*3)}} \left( \sum_{k=1}^{N_{\rm 2d}^{(2*3)}} c^{(2*3)}_{km} d_k^{(2)} \right) d_m^{(3)} ,
which is simplified to
.. math::
E^{(2*3)} = 0.5 \sum_{k=1}^{N_{\rm 2d}^{(2*3)}} b_k^{(2)} d_k^{(2)} + 0.5 \sum_{m=1}^{N_{\rm 3d}^{(2*3)}} b_m^{(3)} d_m^{(3)}
where
.. math::
b_k^{(2)} & = \sum_{m=1}^{N_{\rm 3d}^{(2*3)}} c^{(2*3)}_{km} d_m^{(3)}, \quad k = 1,\ldots, N_{\rm 2d}^{(2*3)}, \\
b_m^{(3)} & = \sum_{k=1}^{N_{\rm 2d}^{(2*3)}} c^{(2*3)}_{km} d_k^{(2)}, \quad m = 1,\ldots, N_{\rm 3d}^{(2*3)} .
The quadratic POD potential results in the following atomic forces
.. math::
\boldsymbol F^{(2*3)} = - \sum_{k=1}^{N_{\rm 2d}^{(2*3)}} \sum_{m=1}^{N_{\rm 3d}^{(2*3)}} c^{(2*3)}_{km} \nabla d^{(2*3)}_{km} .
It can be shown that
.. math::
\boldsymbol F^{(2*3)} = - \sum_{k=1}^{N_{\rm 2d}^{(2*3)}} b^{(2)}_k \nabla d_k^{(2)} - \sum_{m=1}^{N_{\rm 3d}^{(2*3)}} b^{(3)}_m \nabla d_m^{(3)} .
The calculation of the atomic forces for the quadratic POD potential
only requires the extra calculation of :math:`b_k^{(2)}` and :math:`b_m^{(3)}` which can be negligible.
As a result, the quadratic POD potential does not increase the computational complexity.
E_i(\boldsymbol R_i, \boldsymbol Z_i) \ = \ \sum_{m=1}^M c_m \mathcal{D}_{im}(\boldsymbol R_i, \boldsymbol Z_i)
Here :math:`c_m` are trainable coefficients and
:math:`\mathcal{D}_{im}(\boldsymbol R_i, \boldsymbol Z_i)` are per-atom
POD descriptors. Summing the per-atom descriptors over :math:`i` yields
the global descriptors :math:`d_m(\boldsymbol R, \boldsymbol Z) =
\sum_{i=1}^N \mathcal{D}_{im}(\boldsymbol R_i, \boldsymbol Z_i)`. It
thus follows that :math:`E(\boldsymbol R, \boldsymbol Z) = \sum_{m=1}^M
c_m d_m(\boldsymbol R, \boldsymbol Z)`.
The per-atom POD descriptors include one, two, three, four, five, six,
and seven-body descriptors, which can be specified in the first input
file. Furthermore, the per-atom POD descriptors also depend on the
number of environment clusters specified in the first input file.
Please see :ref:`(Nguyen2024) <Nguyen20242a>` and :ref:`(Nguyen and Sema)
<Nguyen20243a>` for the detailed description of the per-atom POD
descriptors.
Training
""""""""
POD potentials are trained using the least-squares regression against
A POD potential is trained using the least-squares regression against
density functional theory (DFT) data. Let :math:`J` be the number of
training configurations, with :math:`N_j` being the number of atoms in
the j-th configuration. Let :math:`\{E^{\star}_j\}_{j=1}^{J}` and
:math:`\{\boldsymbol F^{\star}_j\}_{j=1}^{J}` be the DFT energies and
forces for :math:`J` configurations. Next, we calculate the global
descriptors and their derivatives for all training configurations. Let
:math:`d_{jm}, 1 \le m \le M`, be the global descriptors associated with
the j-th configuration, where :math:`M` is the number of global
descriptors. We then form a matrix :math:`\boldsymbol A \in
\mathbb{R}^{J \times M}` with entries :math:`A_{jm} = d_{jm}/ N_j` for
:math:`j=1,\ldots,J` and :math:`m=1,\ldots,M`. Moreover, we form a
matrix :math:`\boldsymbol B \in \mathbb{R}^{\mathcal{N} \times M}` by
stacking the derivatives of the global descriptors for all training
configurations from top to bottom, where :math:`\mathcal{N} =
3\sum_{j=1}^{J} N_j`.
the j-th configuration. The training configurations are extracted from
the extended XYZ files located in a directory (i.e.,
path_to_training_data_set in the second input file). Let
:math:`\{E^{\star}_j\}_{j=1}^{J}` and :math:`\{\boldsymbol
F^{\star}_j\}_{j=1}^{J}` be the DFT energies and forces for :math:`J`
configurations. Next, we calculate the global descriptors and their
derivatives for all training configurations. Let :math:`d_{jm}, 1 \le m
\le M`, be the global descriptors associated with the j-th
configuration, where :math:`M` is the number of global descriptors. We
then form a matrix :math:`\boldsymbol A \in \mathbb{R}^{J \times M}`
with entries :math:`A_{jm} = d_{jm}/ N_j` for :math:`j=1,\ldots,J` and
:math:`m=1,\ldots,M`. Moreover, we form a matrix :math:`\boldsymbol B
\in \mathbb{R}^{\mathcal{N} \times M}` by stacking the derivatives of
the global descriptors for all training configurations from top to
bottom, where :math:`\mathcal{N} = 3\sum_{j=1}^{J} N_j`.
The coefficient vector :math:`\boldsymbol c` of the POD potential is
found by solving the following least-squares problem
.. math::
{\min}_{\boldsymbol c \in \mathbb{R}^{M}} \ w_E \|\boldsymbol A(\boldsymbol \eta) \boldsymbol c - \bar{\boldsymbol E}^{\star} \|^2 + w_F \|\boldsymbol B(\boldsymbol \eta) \boldsymbol c + \boldsymbol F^{\star} \|^2 + w_R \|\boldsymbol c \|^2,
{\min}_{\boldsymbol c \in \mathbb{R}^{M}} \ w_E \|\boldsymbol A \boldsymbol c - \bar{\boldsymbol E}^{\star} \|^2 + w_F \|\boldsymbol B \boldsymbol c + \boldsymbol F^{\star} \|^2 + w_R \|\boldsymbol c \|^2,
where :math:`w_E` and :math:`w_F` are weights for the energy
(*fitting_weight_energy*) and force (*fitting_weight_force*),
respectively; and :math:`w_R` is the regularization parameter (*fitting_regularization_parameter*). Here :math:`\bar{\boldsymbol E}^{\star} \in
\mathbb{R}^{J}` is a vector of with entries :math:`\bar{E}^{\star}_j =
E^{\star}_j/N_j` and :math:`\boldsymbol F^{\star}` is a vector of
:math:`\mathcal{N}` entries obtained by stacking :math:`\{\boldsymbol
F^{\star}_j\}_{j=1}^{J}` from top to bottom.
respectively; and :math:`w_R` is the regularization parameter
(*fitting_regularization_parameter*). Here :math:`\bar{\boldsymbol
E}^{\star} \in \mathbb{R}^{J}` is a vector of with entries
:math:`\bar{E}^{\star}_j = E^{\star}_j/N_j` and :math:`\boldsymbol
F^{\star}` is a vector of :math:`\mathcal{N}` entries obtained by
stacking :math:`\{\boldsymbol F^{\star}_j\}_{j=1}^{J}` from top to
bottom.
The training procedure is the same for both the linear and quadratic POD
potentials. However, since the quadratic POD potential has a
significantly larger number of the global descriptors, it is more
expensive to train the linear POD potential. This is because the
training of the quadratic POD potential still requires us to calculate
and store the quadratic global descriptors and their
gradient. Furthermore, the quadratic POD potential may require more
training data in order to prevent over-fitting. In order to reduce the
computational cost of fitting the quadratic POD potential and avoid
over-fitting, we can use subsets of two-body and three-body PODs for
constructing the new descriptors.
Validation
""""""""""
POD potential can be validated on a test dataset in a directory
specified by setting path_to_test_data_set in the second input file. It
is possible to validate the POD potential after the training is
complete. This is done by providing the coefficient file as an input to
:doc:`fitpod <fitpod_command>`, for example,
.. code-block:: LAMMPS
fitpod Ta_param.pod Ta_data.pod Ta_coefficients.pod
Restrictions
""""""""""""
@ -727,7 +356,11 @@ LAMMPS was built with that package. See the :doc:`Build package
Related commands
""""""""""""""""
:doc:`pair_style pod <pair_pod>`
:doc:`pair_style pod <pair_pod>`,
:doc:`compute pod/atom <compute_pod_atom>`,
:doc:`compute podd/atom <compute_pod_atom>`,
:doc:`compute pod/local <compute_pod_atom>`,
:doc:`compute pod/global <compute_pod_atom>`
Default
"""""""
@ -736,10 +369,20 @@ The keyword defaults are also given in the description of the input files.
----------
.. _Grepl20072:
.. _Nguyen20222a:
**(Grepl)** Grepl, Maday, Nguyen, and Patera, ESAIM: Mathematical Modelling and Numerical Analysis 41(3), 575-605, (2007).
**(Nguyen and Rohskopf)** Nguyen and Rohskopf, Journal of Computational Physics, 480, 112030, (2023).
.. _Nguyen20232a:
**(Nguyen2023)** Nguyen, Physical Review B, 107(14), 144103, (2023).
.. _Nguyen20242a:
**(Nguyen2024)** Nguyen, Journal of Computational Physics, 113102, (2024).
.. _Nguyen20243a:
**(Nguyen and Sema)** Nguyen and Sema, https://arxiv.org/abs/2405.00306, (2024).
.. _Nguyen20222:
**(Nguyen)** Nguyen and Rohskopf, arXiv preprint arXiv:2209.02362 (2022).

24
doc/src/fix_adapt.rst
View File

@ -21,17 +21,17 @@ Syntax
*pair* args = pstyle pparam I J v_name
pstyle = pair style name (e.g., lj/cut)
pparam = parameter to adapt over time
I,J = type pair(s) to set parameter for
I,J = type pair(s) to set parameter for (integer or type label)
v_name = variable with name that calculates value of pparam
*bond* args = bstyle bparam I v_name
bstyle = bond style name (e.g., harmonic)
bparam = parameter to adapt over time
I = type bond to set parameter for
I = type bond to set parameter for (integer or type label)
v_name = variable with name that calculates value of bparam
*angle* args = astyle aparam I v_name
astyle = angle style name (e.g., harmonic)
aparam = parameter to adapt over time
I = type angle to set parameter for
I = type angle to set parameter for (integer or type label)
v_name = variable with name that calculates value of aparam
*kspace* arg = v_name
v_name = variable with name that calculates scale factor on :math:`k`-space terms
@ -67,6 +67,9 @@ Examples
variable ramp_up equal "ramp(0.01,0.5)"
fix stretch all adapt 1 bond harmonic r0 1 v_ramp_up
labelmap atom 1 c1
fix 1 all adapt 1 pair soft a c1 c1 v_prefactor
Description
"""""""""""
@ -254,10 +257,12 @@ should be specified to indicate which type pairs to apply it to. If a global
parameter is specified, the :math:`I` and :math:`J` settings still need to be
specified, but are ignored.
Similar to the :doc:`pair_coeff command <pair_coeff>`, :math:`I` and :math:`J`
can be specified in one of two ways. Explicit numeric values can be used for
each, as in the first example above. :math:`I \le J` is required. LAMMPS sets
the coefficients for the symmetric :math:`J,I` interaction to the same values.
Similar to the :doc:`pair_coeff command <pair_coeff>`, :math:`I` and
:math:`J` can be specified in one of several ways. Explicit numeric values
can be used for each, as in the first example above. Or, one or both of
the types in the I,J pair can be a :doc:`type label <Howto_type_labels>`.
LAMMPS sets the coefficients for the symmetric :math:`J,I` interaction to
the same values.
A wild-card asterisk can be used in place of or in conjunction with
the :math:`I,J` arguments to set the coefficients for multiple pairs of atom
@ -266,8 +271,9 @@ is the number of atom types, then an asterisk with no numeric values
means all types from 1 to :math:`N`. A leading asterisk means all types from
1 to n (inclusive). A trailing asterisk means all types from m to :math:`N`
(inclusive). A middle asterisk means all types from m to n
(inclusive). Note that only type pairs with :math:`I \le J` are considered; if
asterisks imply type pairs where :math:`J < I`, they are ignored.
(inclusive). For the asterisk syntax, note that only type pairs with
:math:`I \le J` are considered; if asterisks imply type pairs where
:math:`J < I`, they are ignored.
IMPORTANT NOTE: If :doc:`pair_style hybrid or hybrid/overlay
<pair_hybrid>` is being used, then the *pstyle* will be a sub-style

22
doc/src/fix_adapt_fep.rst
View File

@ -21,13 +21,13 @@ Syntax
*pair* args = pstyle pparam I J v_name
pstyle = pair style name (e.g., lj/cut)
pparam = parameter to adapt over time
I,J = type pair(s) to set parameter for
I,J = type pair(s) to set parameter for (integer or type label)
v_name = variable with name that calculates value of pparam
*kspace* arg = v_name
v_name = variable with name that calculates scale factor on K-space terms
*atom* args = aparam v_name
aparam = parameter to adapt over time
I = type(s) to set parameter for
I = type(s) to set parameter for (integer or type label)
v_name = variable with name that calculates value of aparam
* zero or more keyword/value pairs may be appended
@ -56,6 +56,9 @@ Examples
fix 1 all adapt/fep 1 pair lj/cut epsilon * * v_scale1 coul/cut scale 3 3 v_scale2 scale yes reset yes
fix 1 all adapt/fep 10 atom diameter 1 v_size
labelmap atom 1 c1
fix 1 all adapt/fep 1 pair soft a c1 c1 v_prefactor
Example input scripts available: examples/PACKAGES/fep
@ -218,10 +221,12 @@ be specified to indicate which type pairs to apply it to. If a global
parameter is specified, the *I* and *J* settings still need to be
specified, but are ignored.
Similar to the :doc:`pair_coeff command <pair_coeff>`, I and J can be
specified in one of two ways. Explicit numeric values can be used for
each, as in the first example above. :math:`I \le J` is required. LAMMPS sets
the coefficients for the symmetric J,I interaction to the same values.
Similar to the :doc:`pair_coeff command <pair_coeff>`, :math:`I` and
:math:`J` can be specified in one of several ways. Explicit numeric values
can be used for each, as in the first example above. Or, one or both of
the types in the I,J pair can be a :doc:`type label <Howto_type_labels>`.
LAMMPS sets the coefficients for the symmetric :math:`J,I` interaction to
the same values.
A wild-card asterisk can be used in place of or in conjunction with
the :math:`I,J` arguments to set the coefficients for multiple pairs of atom
@ -230,8 +235,9 @@ the number of atom types, then an asterisk with no numeric values means
all types from 1 to :math:`N`. A leading asterisk means all types from 1 to n
(inclusive). A trailing asterisk means all types from m to :math:`N`
(inclusive). A middle asterisk means all types from m to n
(inclusive). Note that only type pairs with :math:`I \le J` are considered; if
asterisks imply type pairs where :math:`J < I`, they are ignored.
(inclusive). For the asterisk syntax, note that only type pairs with
:math:`I \le J` are considered; if asterisks imply type pairs where
:math:`J < I`, they are ignored.
IMPROTANT NOTE: If :doc:`pair_style hybrid or hybrid/overlay <pair_hybrid>` is
being used, then the *pstyle* will be a sub-style name. You must specify

23
doc/src/fix_amoeba_bitorsion.rst
View File

@ -35,7 +35,11 @@ the implementation of AMOEBA and HIPPO in LAMMPS.
Bitorsion interactions add additional potential energy contributions
to pairs of overlapping phi-psi dihedrals of amino-acids, which are
important to properly represent their conformational behavior.
important to properly represent their conformational behavior. Each
bitorsion interaction is thus defined for a 5-tuple of atoms
:math:`IJKLM` with bonds between successive atoms in the list,
i.e. two overlapping dihedral interactions for atoms :math:`IJKL` and
:math:`JKLM`.
The examples/amoeba directory has a sample input script and data file
for ubiquitin, which illustrates use of the fix amoeba/bitorsion
@ -68,14 +72,15 @@ lines:
[...]
N 3 314 315 317 318 330
The first column is an index from 1 to :math:`N` to enumerate the bitorsion
5-atom tuples; it is ignored by LAMMPS. The second column is the
*type* of the interaction; it is an index into the bitorsion force
field file. The remaining 5 columns are the atom IDs of the atoms in
the two 4-atom dihedrals that overlap to create the bitorsion 5-body
interaction. Note that the *bitorsions* and *BiTorsions* keywords for
the header and body sections match those specified in the
:doc:`read_data <read_data>` command following the data file name.
The first column is an index from 1 to :math:`N` to enumerate the
bitorsion 5-atom tuples; it is ignored by LAMMPS. The second column
is the *type* of the interaction; it is an index into the bitorsion
force field file. The remaining 5 columns are the atom IDs of the
atoms (in order) for the 5-tuple :math:`IJKLM`, as described above.
Note that the *bitorsions* and *BiTorsions* keywords for the header
and body sections match those specified in the :doc:`read_data
<read_data>` command following the data file name.
The data file should be generated by using the
tools/tinker/tinker2lmp.py conversion script which creates a LAMMPS

23
doc/src/fix_amoeba_pitorsion.rst
View File

@ -57,7 +57,7 @@ should have two lines like these in its header section:
M pitorsion types
N pitorsions
where :math:`N` is the number of pitorsion 5-body interactions and :math:`M` is
where :math:`N` is the number of pitorsion 6-body interactions and :math:`M` is
the number of pitorsion types. It should also have two sections in the body
of the data file like these with :math:`M` and :math:`N` lines each:
@ -74,21 +74,20 @@ of the data file like these with :math:`M` and :math:`N` lines each:
PiTorsions
1 1 8 10 12 18 20
2 5 18 20 22 25 27
1 1 2 4 3 20 21 24
2 5 21 23 22 37 38 41
[...]
N 3 314 315 317 318 330
N 7 27 29 28 30 35 36
For PiTorsion Coeffs, the first column is an index from 1 to :math:`M` to
enumerate the pitorsion types. The second column is the single
For PiTorsion Coeffs, the first column is an index from 1 to :math:`M`
to enumerate the pitorsion types. The second column is the single
prefactor coefficient needed for each type.
For PiTorsions, the first column is an index from 1 to :math:`N` to enumerate
the pitorsion 5-atom tuples; it is ignored by LAMMPS. The second
column is the "type" of the interaction; it is an index into the
PiTorsion Coeffs. The remaining 5 columns are the atom IDs of the
atoms in the two 4-atom dihedrals that overlap to create the pitorsion
5-body interaction.
For PiTorsions, the first column is an index from 1 to :math:`N` to
enumerate the pitorsion 6-atom tuples; it is ignored by LAMMPS. The
second column is the "type" of the interaction; it is an index into
the PiTorsion Coeffs. The remaining 6 columns are the atom IDs of the
atoms (in order) for the 6-tuple :math:`IJKLMN`, as described above.
Note that the *pitorsion types* and *pitorsions* and *PiTorsion
Coeffs* and *PiTorsions* keywords for the header and body sections of

2
doc/src/fix_atom_swap.rst
View File

@ -21,7 +21,7 @@ Syntax
.. parsed-literal::
*types* values = two or more atom types
*types* values = two or more atom types (1-Ntypes or type label)
*mu* values = chemical potential of swap types (energy units)
*ke* value = *no* or *yes*
*no* = no conservation of kinetic energy after atom swaps

2
doc/src/fix_bond_break.rst
View File

@ -13,7 +13,7 @@ Syntax
* ID, group-ID are documented in :doc:`fix <fix>` command
* bond/break = style name of this fix command
* Nevery = attempt bond breaking every this many steps
* bondtype = type of bonds to break
* bondtype = type of bonds to break (integer or type label)
* Rmax = bond longer than Rmax can break (distance units)
* zero or more keyword/value pairs may be appended
* keyword = *prob*

18
doc/src/fix_bond_create.rst
View File

@ -17,9 +17,9 @@ Syntax
* ID, group-ID are documented in :doc:`fix <fix>` command
* bond/create = style name of this fix command
* Nevery = attempt bond creation every this many steps
* itype,jtype = atoms of itype can bond to atoms of jtype
* itype,jtype = atoms of itype can bond to atoms of jtype (1-Ntypes or type label)
* Rmin = 2 atoms separated by less than Rmin can bond (distance units)
* bondtype = type of created bonds
* bondtype = type of created bonds (integer or type label)
* zero or more keyword/value pairs may be appended to args
* keyword = *iparam* or *jparam* or *prob* or *atype* or *dtype* or *itype* or *aconstrain*
@ -27,19 +27,19 @@ Syntax
*iparam* values = maxbond, newtype
maxbond = max # of bonds of bondtype the itype atom can have
newtype = change the itype atom to this type when maxbonds exist
newtype = change the itype atom to this type when maxbonds exist (1-Ntypes or type label)
*jparam* values = maxbond, newtype
maxbond = max # of bonds of bondtype the jtype atom can have
newtype = change the jtype atom to this type when maxbonds exist
newtype = change the jtype atom to this type when maxbonds exist (1-Ntypes or type label)
*prob* values = fraction seed
fraction = create a bond with this probability if otherwise eligible
seed = random number seed (positive integer)
*atype* value = angletype
angletype = type of created angles
angletype = type of created angles (integer or type label)
*dtype* value = dihedraltype
dihedraltype = type of created dihedrals
dihedraltype = type of created dihedrals (integer or type label)
*itype* value = impropertype
impropertype = type of created impropers
impropertype = type of created impropers (integer or type label)
*aconstrain* value = amin amax
amin = minimal angle at which new bonds can be created
amax = maximal angle at which new bonds can be created
@ -54,6 +54,10 @@ Examples
fix 5 all bond/create 1 3 3 0.8 1 prob 0.5 85784 iparam 2 3 atype 1 dtype 2
fix 5 all bond/create/angle 10 1 2 1.122 1 aconstrain 120 180 prob 1 4928459 iparam 2 1 jparam 2 2
labelmap atom 1 c1 2 n2
labelmap bond 1 c1-n2
fix 5 all bond/create 10 c1 n2 0.8 c1-n2
Description
"""""""""""

11
doc/src/fix_charge_regulation.rst
View File

@ -13,8 +13,8 @@ Syntax
* ID, group-ID are documented in fix command
* charge/regulation = style name of this fix command
* cation_type = atom type of free cations
* anion_type = atom type of free anions
* cation_type = atom type of free cations (integer or type label)
* anion_type = atom type of free anions (integer or type label)
* zero or more keyword/value pairs may be appended
@ -27,8 +27,8 @@ Syntax
*pIp* value = activity (effective concentration) of free cations (in the -log10 representation)
*pIm* value = activity (effective concentration) of free anions (in the -log10 representation)
*pKs* value = solvent self-dissociation constant (in the -log10 representation)
*acid_type* = atom type of acid groups
*base_type* = atom type of base groups
*acid_type* = atom type of acid groups (integer or type label)
*base_type* = atom type of base groups (integer or type label)
*lunit_nm* value = unit length used by LAMMPS (# in the units of nanometers)
*temp* value = temperature
*tempfixid* value = fix ID of temperature thermostat
@ -51,6 +51,9 @@ Examples
fix chareg all charge/regulation 1 2 acid_type 3 base_type 4 pKa 5.0 pKb 6.0 pH 7.0 pIp 3.0 pIm 3.0 nevery 200 nmc 200 seed 123 tempfixid fT
fix chareg all charge/regulation 1 2 pIp 3 pIm 3 onlysalt yes 2 -1 seed 123 tag yes temp 1.0
labelmap atom 1 H+ 2 OH-
fix chareg all charge/regulation H+ OH- pIp 3 pIm 3 onlysalt yes 2 -1 seed 123 tag yes temp 1.0
Description
"""""""""""

2
doc/src/fix_deform.rst
View File

@ -64,6 +64,8 @@ Syntax
effectively an engineering shear strain rate
*erate* value = R
R = engineering shear strain rate (1/time units)
*erate/rescale* value = R (ONLY available in :doc:`fix deform/pressure <fix_deform_pressure>` command)
R = engineering shear strain rate (1/time units)
*trate* value = R
R = true shear strain rate (1/time units)
*wiggle* values = A Tp

5
doc/src/fix_deposit.rst
View File

@ -13,7 +13,7 @@ Syntax
* ID, group-ID are documented in :doc:`fix <fix>` command
* deposit = style name of this fix command
* N = # of atoms or molecules to insert
* type = atom type to assign to inserted atoms (offset for molecule insertion)
* type = atom type (1-Ntypes or type label) to assign to inserted atoms (offset for molecule insertion)
* M = insert a single atom or molecule every M steps
* seed = random # seed (positive integer)
* one or more keyword/value pairs may be appended to args
@ -76,6 +76,9 @@ Examples
fix 4 sputter deposit 1000 2 500 12235 region sphere vz -1.0 -1.0 target 5.0 5.0 0.0 units lattice
fix 5 insert deposit 200 2 100 777 region disk gaussian 5.0 5.0 9.0 1.0 units box
labelmap atom 1 Au
fix 4 sputter deposit 1000 Au 500 12235 region sphere vz -1.0 -1.0 target 5.0 5.0 0.0 units lattice
Description
"""""""""""

71
doc/src/fix_gcmc.rst
View File

@ -15,7 +15,7 @@ Syntax
* N = invoke this fix every N steps
* X = average number of GCMC exchanges to attempt every N steps
* M = average number of MC moves to attempt every N steps
* type = atom type for inserted atoms (must be 0 if mol keyword used)
* type = atom type (1-Ntypes or type label) for inserted atoms (must be 0 if mol keyword used)
* seed = random # seed (positive integer)
* T = temperature of the ideal gas reservoir (temperature units)
* mu = chemical potential of the ideal gas reservoir (energy units)
@ -45,7 +45,7 @@ Syntax
*group* value = group-ID
group-ID = group-ID for inserted atoms (string)
*grouptype* values = type group-ID
type = atom type (int)
type = atom type (1-Ntypes or type label)
group-ID = group-ID for inserted atoms (string)
*intra_energy* value = intramolecular energy (energy units)
*tfac_insert* value = scale up/down temperature of inserted atoms (unitless)
@ -62,52 +62,47 @@ Examples
fix 3 water gcmc 10 100 100 0 3456543 3.0 -2.5 0.1 mol my_one_water maxangle 180 full_energy
fix 4 my_gas gcmc 1 10 10 1 123456543 300.0 -12.5 1.0 region disk
labelmap atom 1 Li
fix 2 ion gcmc 10 1000 1000 Li 29494 298.0 -0.5 0.01
Description
"""""""""""
This fix performs grand canonical Monte Carlo (GCMC) exchanges of
atoms or molecules with an imaginary ideal gas
reservoir at the specified T and chemical potential (mu) as discussed
in :ref:`(Frenkel) <Frenkel2>`. It also
attempts Monte Carlo (MC) moves (translations and molecule
rotations) within the simulation cell or
region. If used with the :doc:`fix nvt <fix_nh>`
This fix performs grand canonical Monte Carlo (GCMC) exchanges of atoms or
molecules with an imaginary ideal gas reservoir at the specified T and
chemical potential (mu) as discussed in :ref:`(Frenkel) <Frenkel2>`. It
also attempts Monte Carlo (MC) moves (translations and molecule rotations)
within the simulation cell or region. If used with the :doc:`fix nvt <fix_nh>`
command, simulations in the grand canonical ensemble (muVT, constant
chemical potential, constant volume, and constant temperature) can be
performed. Specific uses include computing isotherms in microporous
materials, or computing vapor-liquid coexistence curves.
Every N timesteps the fix attempts both GCMC exchanges
(insertions or deletions) and MC moves of gas atoms or molecules.
On those timesteps, the average number of attempted GCMC exchanges is X,
while the average number of attempted MC moves is M.
For GCMC exchanges of either molecular or atomic gasses,
these exchanges can be either deletions or insertions,
with equal probability.
Every N timesteps the fix attempts both GCMC exchanges (insertions or
deletions) and MC moves of gas atoms or molecules. On those timesteps, the
average number of attempted GCMC exchanges is X, while the average number
of attempted MC moves is M. For GCMC exchanges of either molecular or
atomic gasses, these exchanges can be either deletions or insertions, with
equal probability.
The possible choices for MC moves are translation of an atom,
translation of a molecule, and rotation of a molecule.
The relative amounts of each are determined by the optional
*mcmoves* keyword (see below).
The default behavior is as follows.
If the *mol* keyword is used, only molecule translations
and molecule rotations are performed with equal probability.
Conversely, if the *mol* keyword is not used, only atom
translations are performed.
M should typically be
chosen to be approximately equal to the expected number of gas atoms
or molecules of the given type within the simulation cell or region,
which will result in roughly one MC move per atom or molecule
per MC cycle.
The possible choices for MC moves are translation of an atom, translation
of a molecule, and rotation of a molecule. The relative amounts of each are
determined by the optional *mcmoves* keyword (see below). The default
behavior is as follows. If the *mol* keyword is used, only molecule
translations and molecule rotations are performed with equal probability.
Conversely, if the *mol* keyword is not used, only atom translations are
performed. M should typically be chosen to be approximately equal to the
expected number of gas atoms or molecules of the given type within the
simulation cell or region, which will result in roughly one MC move per
atom or molecule per MC cycle.
All inserted particles are always added to two groups: the default
group "all" and the fix group specified in the fix command.
In addition, particles are also added to any groups
specified by the *group* and *grouptype* keywords. If inserted
particles are individual atoms, they are assigned the atom type given
by the type argument. If they are molecules, the type argument has no
effect and must be set to zero. Instead, the type of each atom in the
inserted molecule is specified in the file read by the
All inserted particles are always added to two groups: the default group
"all" and the fix group specified in the fix command. In addition,
particles are also added to any groups specified by the *group* and
*grouptype* keywords. If inserted particles are individual atoms, they are
assigned the atom type given by the type argument. If they are molecules,
the type argument has no effect and must be set to zero. Instead, the type
of each atom in the inserted molecule is specified in the file read by the
:doc:`molecule <molecule>` command.
.. note::

2
doc/src/fix_ipi.rst
View File

@ -80,7 +80,7 @@ Obtaining i-PI
""""""""""""""
Here are the commands to set up a virtual environment and install
i-PI into it with all its dependencies via the PyPi repository and
i-PI into it with all its dependencies via the PyPI repository and
the pip package manager.
.. code-block:: sh

5
doc/src/fix_mol_swap.rst
View File

@ -14,7 +14,7 @@ Syntax
* atom/swap = style name of this fix command
* N = invoke this fix every N steps
* X = number of swaps to attempt every N steps
* itype,jtype = two atom types to swap with each other
* itype,jtype = two atom types (1-Ntypes or type label) to swap with each other
* seed = random # seed (positive integer)
* T = scaling temperature of the MC swaps (temperature units)
* zero or more keyword/value pairs may be appended to args
@ -34,6 +34,9 @@ Examples
fix 2 all mol/swap 100 1 2 3 29494 300.0 ke no
fix mySwap fluid mol/swap 500 10 1 2 482798 1.0
labelmap atom 1 A 2 B
fix mySwap fluid mol/swap 500 10 A B 482798 1.0
Description
"""""""""""

31
doc/src/fix_sgcmc.rst
View File

@ -15,7 +15,7 @@ Syntax
* every_nsteps = number of MD steps between MC cycles
* swap_fraction = fraction of a full MC cycle carried out at each call (a value of 1.0 will perform as many trial moves as there are atoms)
* temperature = temperature that enters Boltzmann factor in Metropolis criterion (usually the same as MD temperature)
* deltamu = chemical potential difference(s) (`N-1` values must be provided, with `N` being the number of elements)
* deltamu = `N-1` chemical potential differences :math:`\mu_1-\mu_2, \ldots, \mu_1-\mu_N` (`N` is the number of atom types)
* Zero or more keyword/value pairs may be appended to fix definition line:
.. parsed-literal::
@ -23,7 +23,7 @@ Syntax
keyword = *variance* or *randseed* or *window_moves* or *window_size*
*variance* kappa conc1 [conc2] ... [concN]
kappa = variance constraint parameter
conc1,conc2,... = target concentration(s) in the range 0.0-1.0 (*N-1* values must be provided, with *N* being the number of elements)
`c_2`, `c_3`,..., `c_N` = `N-1` target concentration fractions
*randseed* N
N = seed for pseudo random number generator
*window_moves* N
@ -90,11 +90,10 @@ the simulation, e.g., to speed up equilibration at low temperatures.
------------
The parameter *deltamu* is used to set the chemical potential difference
in the SGC MC algorithm (see Eq. 16 in :ref:`Sadigh1 <Sadigh1>`). By
convention it is the difference of the chemical potentials of elements
`B`, `C` ..., with respect to element A. When the simulation includes
`N` elements, `N-1` values must be specified.
The parameter *deltamu* is used to set the chemical potential differences
in the SGC MC algorithm (see Eq. 16 in :ref:`Sadigh1 <Sadigh1>`).
The `N-1` differences are defined as :math:`\mu_1-\mu_2, \ldots, \mu_1-\mu_N`,
where `N` is the number of atom types.
------------
@ -104,12 +103,12 @@ the effective average constraint in the parallel VC-SGC MC algorithm
(parameter :math:`\delta\mu_0` in Eq. (20) of :ref:`Sadigh1
<Sadigh1>`). The parameter *kappa* specifies the variance constraint
(see Eqs. (20-21) in :ref:`Sadigh1 <Sadigh1>`).
The parameter *conc* sets the target concentration (parameter
:math:`c_0` in Eqs. (20-21) of :ref:`Sadigh1 <Sadigh1>`). The atomic
concentrations refer to components `B`, `C` ..., with `A` being set
automatically. When the simulation includes `N` elements, `N-1`
concentration values must be specified.
The parameter *conc* sets the `N-1` target atomic concentration
fractions (parameter :math:`c_0` in Eqs. (20-21) of :ref:`Sadigh1 <Sadigh1>`)
:math:`0 \le c_2, \ldots, c_N \le 1`, with
:math:`c_1 = 1 - \Sigma_{i=2}^N c_i`.
When the simulation includes `N` atom types (elements),
`N-1` concentration values must be specified.
------------
@ -143,10 +142,10 @@ components of the vector represent the following quantities:
* 1 = The absolute number of accepted trial swaps during the last MC step
* 2 = The absolute number of rejected trial swaps during the last MC step
* 3 = The current global concentration of species *A* (= number of atoms of type 1 / total number of atoms)
* 4 = The current global concentration of species *B* (= number of atoms of type 2 / total number of atoms)
* 3 = Current global concentration `c_1` (= number of atoms of type 1 / total number of atoms)
* 4 = Current global concentration `c_2` (= number of atoms of type 2 / total number of atoms)
* ...
* N+2: The current global concentration of species *X* (= number of atoms of type *N* / total number of atoms)
* N+2: Current global concentration `c_N` (= number of atoms of type *N* / total number of atoms)
The vector values calculated by this fix are "intensive".

12
doc/src/fix_wall_gran.rst
View File

@ -115,6 +115,18 @@ friction and twisting friction supported by the :doc:`pair_style granular <pair_
supported for walls. These are discussed in greater detail on the doc
page for :doc:`pair_style granular <pair_granular>`.
.. note::
When *fstyle* *granular* is specified, the associated *fstyle_params* are taken as
those for a wall/particle interaction. For example, for the *hertz/material* normal
contact model with :math:`E = 960` and :math:`\nu = 0.2`, the effective Young's
modulus for a wall/particle interaction is computed as
:math:`E_{eff} = \frac{960}{2(1-0.2^2)} = 500`. Any pair coefficients defined by
:doc:`pair_style granular <pair_granular>` are not taken into consideration. To
model different wall/particle interactions for particles of different material
types, the user may define multiple fix wall/gran commands operating on separate
groups (e.g. based on particle type) each with a different wall/particle effective
Young's modulus.
Note that you can choose a different force styles and/or different
values for the wall/particle coefficients than for particle/particle
interactions. E.g. if you wish to model the wall as a different

5
doc/src/fix_widom.rst
View File

@ -14,7 +14,7 @@ Syntax
* widom = style name of this fix command
* N = invoke this fix every N steps
* M = number of Widom insertions to attempt every N steps
* type = atom type for inserted atoms (must be 0 if mol keyword used)
* type = atom type (1-Ntypes or type label) for inserted atoms (must be 0 if mol keyword used)
* seed = random # seed (positive integer)
* T = temperature of the system (temperature units)
* zero or more keyword/value pairs may be appended to args
@ -38,6 +38,9 @@ Examples
fix 2 gas widom 1 50000 1 19494 2.0
fix 3 water widom 1000 100 0 29494 300.0 mol h2omol full_energy
labelmap atom 1 Li
fix 2 ion widom 1 50000 Li 19494 2.0
Description
"""""""""""

12
doc/src/group.rst
View File

@ -20,13 +20,13 @@ Syntax
*empty* = no args
*region* args = region-ID
*type* or *id* or *molecule*
args = list of one or more atom types, atom IDs, or molecule IDs
any entry in list can be a sequence formatted as A:B or A:B:C where
args = list of one or more atom types (1-Ntypes or type label), atom IDs, or molecule IDs
any numeric entry in list can be a sequence formatted as A:B or A:B:C where
A = starting index, B = ending index,
C = increment between indices, 1 if not specified
args = logical value
logical = "<" or "<=" or ">" or ">=" or "==" or "!="
value = an atom type or atom ID or molecule ID (depending on *style*\ )
value = an atom type (1-Ntypes or type label) or atom ID or molecule ID (depending on *style*\ )
args = logical value1 value2
logical = "<>"
value1,value2 = atom types or atom IDs or molecule IDs (depending on *style*\ )
@ -52,6 +52,7 @@ Examples
group edge region regstrip
group water type 3 4
group water type OW HT
group sub id 10 25 50
group sub id 10 25 50 500:1000
group sub id 100:10000:10
@ -119,7 +120,7 @@ three styles can use arguments specified in one of two formats.
The first format is a list of values (types or IDs). For example, the
second command in the examples above puts all atoms of type 3 or 4 into
the group named *water*\ . Each entry in the list can be a
the group named *water*\ . Each numeric entry in the list can be a
colon-separated sequence ``A:B`` or ``A:B:C``, as in two of the examples
above. A "sequence" generates a sequence of values (types or IDs),
with an optional increment. The first example with ``500:1000`` has the
@ -135,7 +136,8 @@ except ``<>`` take a single argument. The third example above adds all
atoms with IDs from 1 to 150 to the group named *sub*\ . The logical ``<>``
means "between" and takes 2 arguments. The fourth example above adds all
atoms belonging to molecules with IDs from 50 to 250 (inclusive) to
the group named polyA.
the group named polyA. For the *type* style, type labels are converted into
numeric types before being evaluated.
The *variable* style evaluates a variable to determine which atoms to
add to the group. It must be an :doc:`atom-style variable <variable>`

64
doc/src/group2ndx.rst
View File

@ -34,32 +34,66 @@ Description
Write or read a Gromacs style index file in text format that associates
atom IDs with the corresponding group definitions. This index file can be
used with in combination with Gromacs analysis tools or to import group
definitions into the :doc:`fix colvars <fix_colvars>` input file. It can
also be used to save and restore group definitions for static groups.
definitions into the :doc:`fix colvars <fix_colvars>` input file.
It can also be used to save and restore group definitions for static groups
using the individual atom IDs. This may be important if the original
group definition depends on a region or otherwise on the geometry and thus
cannot be easily recreated.
Another application would be to import atom groups defined for Gromacs
simulation into LAMMPS. When translating Gromacs topology and geometry
data to LAMMPS.
The *group2ndx* command will write group definitions to an index file.
Without specifying any group IDs, all groups will be written to the index
file. When specifying group IDs, only those groups will be written to the
index file. In order to follow the Gromacs conventions, the group *all*
will be renamed to *System* in the index file.
Without specifying any group IDs, all groups will be written to the
index file. When specifying group IDs, only those groups will be
written to the index file. In order to follow the Gromacs conventions,
the group *all* will be renamed to *System* in the index file.
The *ndx2group* command will create of update group definitions from those
stored in an index file. Without specifying any group IDs, all groups except
*System* will be read from the index file and the corresponding groups
recreated. If a group of the same name already exists, it will be completely
reset. When specifying group IDs, those groups, if present, will be read
from the index file and restored.
The *ndx2group* command will create of update group definitions from
those stored in an index file. Without specifying any group IDs, all
groups except *System* will be read from the index file and the
corresponding groups recreated. If a group of the same name already
exists, it will be completely reset. When specifying group IDs, those
groups, if present, will be read from the index file and restored.
File Format
"""""""""""
The file format is equivalent and compatible with what is produced by
the `Gromacs make_ndx command <https://manual.gromacs.org/current/onlinehelp/gmx-make_ndx.html>`_.
and follows the `Gromacs definition of an ndx file <https://manual.gromacs.org/current/reference-manual/file-formats.html#ndx>`_
Each group definition begins with the group name in square brackets with
blanks, e.g. \[ water \] and is then followed by the list of atom
indices, which may be spread over multiple lines. Here is a small
example file:
.. code-block:: ini
[ Oxygen ]
1 4 7
[ Hydrogen ]
2 3 5 6
8 9
[ Water ]
1 2 3 4 5 6 7 8 9
The index file defines 3 groups: Oxygen, Hydrogen, and Water and the
latter happens to be the union of the first two.
----------
Restrictions
""""""""""""
This command requires that atoms have atom IDs, since this is the
These commands require that atoms have atom IDs, since this is the
information that is written to the index file.
These commands are part of the COLVARS package. They are only
enabled if LAMMPS was built with that package. See the :doc:`Build package <Build_package>` page for more info.
These commands are part of the EXTRA-COMMAND package. They are only
enabled if LAMMPS was built with that package. See the
:doc:`Build package <Build_package>` page for more info.
Related commands
""""""""""""""""

BIN
doc/src/img/howto_charmm_ELJ.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 50 KiB

BIN
doc/src/img/howto_charmmfsw_ELJ.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 47 KiB

7
doc/src/labelmap.rst
View File

@ -24,6 +24,7 @@ Examples
.. code-block:: LAMMPS
labelmap atom 1 c1 2 hc 3 cp 4 nt
labelmap atom 3 carbon 4 'c3"' 5 "c1'" 6 "c#"
labelmap atom $(label2type(atom,carbon)) C # change type label from 'carbon' to 'C'
labelmap clear
@ -43,8 +44,8 @@ The label map can also be defined by the :doc:`read_data <read_data>`
command when it reads these sections in a data file: Atom Type Labels,
Bond Type Labels, etc. See the :doc:`Howto type labels
<Howto_type_labels>` doc page for a general discussion of how type
labels can be used. See :ref:`(Gissinger) <Typelabel>` for a discussion
of the type label implementation in LAMMPS and its uses.
labels can be used. See :ref:`(Gissinger) <Typelabel1>` for a
discussion of the type label implementation in LAMMPS and its uses.
Valid type labels can contain any alphanumeric character, but must not
start with a number, a '#', or a '*' character. They can contain other
@ -102,6 +103,6 @@ none
-----------
.. _Typelabel:
.. _Typelabel1:
**(Gissinger)** J. R. Gissinger, I. Nikiforov, Y. Afshar, B. Waters, M. Choi, D. S. Karls, A. Stukowski, W. Im, H. Heinz, A. Kohlmeyer, and E. B. Tadmor, J Phys Chem B, 128, 3282-3297 (2024).

42
doc/src/pair_charmm.rst
View File

@ -112,26 +112,22 @@ Description
These pair styles compute Lennard Jones (LJ) and Coulombic
interactions with additional switching or shifting functions that ramp
the energy and/or force smoothly to zero between an inner and outer
cutoff. They are implementations of the widely used CHARMM force
field used in the `CHARMM <https://www.charmm.org>`_ MD code (and
others). See :ref:`(MacKerell) <pair-MacKerell>` for a description of the
CHARMM force field.
cutoff. They implement the widely used CHARMM force field, see
:doc:`Howto discussion on biomolecular force fields <Howto_bioFF>` for
details.
The styles with *charmm* (not *charmmfsw* or *charmmfsh*\ ) in their
name are the older, original LAMMPS implementations. They compute the
LJ and Coulombic interactions with an energy switching function (esw,
shown in the formula below as S(r)), which ramps the energy smoothly
to zero between the inner and outer cutoff. This can cause
irregularities in pairwise forces (due to the discontinuous second
derivative of energy at the boundaries of the switching region), which
in some cases can result in detectable artifacts in an MD simulation.
LJ and Coulombic interactions with an energy switching function which
ramps the energy smoothly to zero between the inner and outer cutoff.
This can cause irregularities in pairwise forces (due to the discontinuous
second derivative of energy at the boundaries of the switching region),
which in some cases can result in detectable artifacts in an MD simulation.
The newer styles with *charmmfsw* or *charmmfsh* in their name replace
the energy switching with force switching (fsw) and force shifting
(fsh) functions, for LJ and Coulombic interactions respectively.
These follow the formulas and description given in
:ref:`(Steinbach) <Steinbach>` and :ref:`(Brooks) <Brooks1>` to minimize these
artifacts.
.. note::
@ -152,26 +148,6 @@ artifacts.
the CHARMM force field energies and forces, when using one of these
two CHARMM pair styles.
.. math::
E = & LJ(r) \qquad \qquad \qquad r < r_{\rm in} \\
= & S(r) * LJ(r) \qquad \qquad r_{\rm in} < r < r_{\rm out} \\
= & 0 \qquad \qquad \qquad \qquad r > r_{\rm out} \\
E = & C(r) \qquad \qquad \qquad r < r_{\rm in} \\
= & S(r) * C(r) \qquad \qquad r_{\rm in} < r < r_{\rm out} \\
= & 0 \qquad \qquad \qquad \qquad r > r_{\rm out} \\
LJ(r) = & 4 \epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} -
\left(\frac{\sigma}{r}\right)^6 \right] \\
C(r) = & \frac{C q_i q_j}{ \epsilon r} \\
S(r) = & \frac{ \left[r_{\rm out}^2 - r^2\right]^2
\left[r_{\rm out}^2 + 2r^2 - 3{r_{\rm in}^2}\right]}
{ \left[r_{\rm out}^2 - {r_{\rm in}}^2\right]^3 }
where S(r) is the energy switching function mentioned above for the
*charmm* styles. See the :ref:`(Steinbach) <Steinbach>` paper for the
functional forms of the force switching and force shifting functions
used in the *charmmfsw* and *charmmfsh* styles.
When using the *lj/charmm/coul/charmm styles*, both the LJ and
Coulombic terms require an inner and outer cutoff. They can be the
same for both formulas or different depending on whether 2 or 4

196
doc/src/pair_dpd_coul_slater_long.rst Normal file
View File

@ -0,0 +1,196 @@
.. index:: pair_style dpd/coul/slater/long
.. index:: pair_style dpd/coul/slater/long/gpu
pair_style dpd/coul/slater/long command
=======================================
Accelerator Variants: *dpd/coul/slater/long/gpu*
Syntax
""""""
.. code-block:: LAMMPS
pair_style dpd/coul/slater/long T cutoff_DPD seed lambda cutoff_coul
* T = temperature (temperature units)
* cutoff_DPD = global cutoff for DPD interactions (distance units)
* seed = random # seed (positive integer)
* lambda = decay length of the charge (distance units)
* cutoff_coul = global cutoff for Coulombic interactions (distance units)
Examples
""""""""
.. code-block:: LAMMPS
pair_style dpd/coul/slater/long 1.0 2.5 34387 0.25 3.0
pair_coeff 1 1 78.0 4.5 # not charged by default
pair_coeff 2 2 78.0 4.5 yes
Description
"""""""""""
.. versionadded:: 27June2024
Style *dpd/coul/slater/long* computes a force field for dissipative
particle dynamics (DPD) following the exposition in :ref:`(Groot)
<Groot5>`. It also allows for the use of charged particles in the
model by adding a long-range Coulombic term to the DPD interactions.
The short-range portion of the Coulombics is calculated by this pair
style. The long-range Coulombics are computed by use of the
:doc:`kspace_style <kspace_style>` command, e.g. using the Ewald or
PPPM styles.
Coulombic forces in mesoscopic models such as DPD employ potentials
without explicit excluded-volume interactions. The goal is to prevent
artificial ionic pair formation by including a charge distribution in
the Coulomb potential, following the formulation in :ref:`(Melchor1)
<Melchor1>`.
.. note::
This pair style is effectively the combination of the
:doc:`pair_style dpd <pair_dpd>` and :doc:`pair_style
coul/slater/long <pair_coul_slater>` commands, but should be more
efficient (especially on GPUs) than using :doc:`pair_style
hybrid/overlay dpd coul/slater/long <pair_hybrid>`. That is
particularly true for the GPU package version of the pair style since
this version is compatible with computing neighbor lists on the GPU
instead of the CPU as is required for hybrid styles.
In the charged DPD model, the force on bead I due to bead J is given
as a sum of 4 terms:
.. math::
\vec{f} = & (F^C + F^D + F^R + F^E) \hat{r_{ij}} \\
F^C = & A w(r) \qquad \qquad \qquad \qquad \qquad r < r_{DPD} \\
F^D = & - \gamma w^2(r) (\hat{r_{ij}} \bullet \vec{v}_{ij}) \qquad \qquad r < r_{DPD} \\
F^R = & \sigma w(r) \alpha (\Delta t)^{-1/2} \qquad \qquad \qquad r < r_{DPD} \\
w(r) = & 1 - \frac{r}{r_{DPD}} \\
F^E = & \frac{C q_iq_j}{\epsilon r^2} \left( 1- exp\left( \frac{2r_{ij}}{\lambda} \right) \left( 1 + \frac{2r_{ij}}{\lambda} \left( 1 + \frac{r_{ij}}{\lambda} \right)\right) \right)
where :math:`F^C` is a conservative force, :math:`F^D` is a
dissipative force, :math:`F^R` is a random force, and :math:`F^E` is
an electrostatic force. :math:`\hat{r_{ij}}` is a unit vector in the
direction :math:`r_i - r_j`, :math:`\vec{v}_{ij}` is the vector
difference in velocities of the two atoms :math:`\vec{v}_i -
\vec{v}_j`, :math:`\alpha` is a Gaussian random number with zero mean
and unit variance, *dt* is the timestep size, and :math:`w(r)` is a
weighting factor that varies between 0 and 1.
:math:`\sigma` is set equal to :math:`\sqrt{2 k_B T \gamma}`, where
:math:`k_B` is the Boltzmann constant and *T* is the temperature
parameter in the pair_style command.
:math:`r_{DPD}` is the pairwise cutoff for the first 3 DPD terms in
the formula as specified by *cutoff_DPD*. For the :math:`F^E` term,
pairwise interactions within the specified *cutoff_coul* distance are
computed directly; interactions beyond that distance are computed in
reciprocal space. *C* is the same Coulomb conversion factor used in
the Coulombic formulas described on the :doc:`pair_coul <pair_coul>`
doc page.
The following parameters must be defined for each pair of atoms types
via the :doc:`pair_coeff <pair_coeff>` command as in the examples
above, or in the data file or restart files read by the
:doc:`read_data <read_data>` or :doc:`read_restart <read_restart>`
commands:
* A (force units)
* :math:`\gamma` (force/velocity units)
* is_charged (optional boolean, default = no)
The *is_charged* parameter is optional and can be specified as *yes* or
*no*. *Yes* should be used for interactions between two types of
charged particles. *No* is the default and should be used for
interactions between two types of particles when one or both are
uncharged.
----------
.. include:: accel_styles.rst
----------
Mixing, shift, table, tail correction, restart, rRESPA info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
This pair style does not support mixing. Thus, coefficients for all
I,J pairs must be specified explicitly.
This pair style does not support the :doc:`pair_modify <pair_modify>`
shift option for the energy of the pair interaction.
The :doc:`pair_modify <pair_modify>` table option is not relevant
for this pair style.
This pair style does not support the :doc:`pair_modify <pair_modify>`
tail option for adding long-range tail corrections to energy and
pressure.
This pair style writes its information to :doc:`binary restart files
<restart>`, so pair_style and pair_coeff commands do not need to be
specified in an input script that reads a restart file. Note that the
user-specified random number seed is stored in the restart file, so
when a simulation is restarted, each processor will re-initialize its
random number generator the same way it did initially. This means the
random forces will be random, but will not be the same as they would
have been if the original simulation had continued past the restart
time.
This pair style can only be used via the *pair* keyword of the
:doc:`run_style respa <run_style>` command. It does not support the
*inner*, *middle*, *outer* keywords.
----------
Restrictions
""""""""""""
This style is part of the DPD-BASIC package. It is only enabled if
LAMMPS was built with that package. See the :doc:`Build package
<Build_package>` page for more info.
The default frequency for rebuilding neighbor lists is every 10 steps
(see the :doc:`neigh_modify <neigh_modify>` command). This may be too
infrequent since particles move rapidly and can overlap by large
amounts. If this setting yields a non-zero number of "dangerous"
reneighborings (printed at the end of a simulation), you should
experiment with forcing reneighboring more often and see if system
energies/trajectories change.
This pair style requires use of the :doc:`comm_modify vel yes
<comm_modify>` command so that velocities are stored by ghost atoms.
This pair style also requires use of a long-range solvers from the
KSPACE package.
This pair style will not restart exactly when using the
:doc:`read_restart <read_restart>` command, though they should provide
statistically similar results. This is because the forces they compute
depend on atom velocities. See the :doc:`read_restart <read_restart>`
command for more details.
Related commands
""""""""""""""""
:doc:`pair_style dpd <pair_dpd>`, :doc:`pair_style coul/slater/long <pair_coul_slater>`,
Default
"""""""
For the pair_coeff command, the default is is_charged = no.
----------
.. _Groot5:
**(Groot)** Groot and Warren, J Chem Phys, 107, 4423-35 (1997).
.. _Melchor1:
**(Melchor)** Gonzalez-Melchor, Mayoral, Velazquez, and Alejandre, J Chem Phys, 125, 224107 (2006).

30
doc/src/pair_granular.rst
View File

@ -111,7 +111,7 @@ For the *hertz* model, the normal component of force is given by:
\mathbf{F}_{ne, Hertz} = k_n R_{eff}^{1/2}\delta_{ij}^{3/2} \mathbf{n}
Here, :math:`R_{eff} = \frac{R_i R_j}{R_i + R_j}` is the effective
Here, :math:`R_{eff} = R = \frac{R_i R_j}{R_i + R_j}` is the effective
radius, denoted for simplicity as *R* from here on. For *hertz*, the
units of the spring constant :math:`k_n` are *force*\ /\ *length*\ \^2, or
equivalently *pressure*\ .
@ -120,13 +120,14 @@ For the *hertz/material* model, the force is given by:
.. math::
\mathbf{F}_{ne, Hertz/material} = \frac{4}{3} E_{eff} R_{eff}^{1/2}\delta_{ij}^{3/2} \mathbf{n}
\mathbf{F}_{ne, Hertz/material} = \frac{4}{3} E_{eff} R^{1/2}\delta_{ij}^{3/2} \mathbf{n}
Here, :math:`E_{eff} = E = \left(\frac{1-\nu_i^2}{E_i} + \frac{1-\nu_j^2}{E_j}\right)^{-1}` is the effective Young's
modulus, with :math:`\nu_i, \nu_j` the Poisson ratios of the particles of
types *i* and *j*\ . Note that if the elastic modulus and the shear
modulus of the two particles are the same, the *hertz/material* model
is equivalent to the *hertz* model with :math:`k_n = 4/3 E_{eff}`
Here, :math:`E_{eff} = E = \left(\frac{1-\nu_i^2}{E_i} + \frac{1-\nu_j^2}{E_j}\right)^{-1}`
is the effective Young's modulus, with :math:`\nu_i, \nu_j` the Poisson ratios
of the particles of types *i* and *j*. :math:`E_{eff}` is denoted as *E* from here on.
Note that if the elastic modulus and the shear modulus of the two particles are the
same, the *hertz/material* model is equivalent to the *hertz* model with
:math:`k_n = 4/3 E`
The *dmt* model corresponds to the
:ref:`(Derjaguin-Muller-Toporov) <DMT1975>` cohesive model, where the force
@ -270,7 +271,8 @@ where :math:`k_n = \frac{4}{3} E_{eff}` for the *hertz/material* model. Since
*coeff_restitution* accounts for the effective mass, effective radius, and
pairwise overlaps (except when used with the *hooke* normal model) when calculating
the damping coefficient, it accurately reproduces the specified coefficient of
restitution for both monodisperse and polydisperse particle pairs.
restitution for both monodisperse and polydisperse particle pairs. This damping
model is not compatible with cohesive normal models such as *JKR* or *DMT*.
The total normal force is computed as the sum of the elastic and
damping components:
@ -441,11 +443,11 @@ discussion above. To match the Mindlin solution, one should set
G_{eff} = \left(\frac{2-\nu_i}{G_i} + \frac{2-\nu_j}{G_j}\right)^{-1}
where :math:`G` is the shear modulus, related to Young's modulus :math:`E`
and Poisson's ratio :math:`\nu` by :math:`G = E/(2(1+\nu))`. This can also be
achieved by specifying *NULL* for :math:`k_t`, in which case a
normal contact model that specifies material parameters :math:`E` and
:math:`\nu` is required (e.g. *hertz/material*, *dmt* or *jkr*\ ). In this
where :math:`G_i` is the shear modulus of a particle of type :math:`i`, related to Young's
modulus :math:`E_i` and Poisson's ratio :math:`\nu_i` by :math:`G_i = E_i/(2(1+\nu_i))`.
This can also be achieved by specifying *NULL* for :math:`k_t`, in which case a
normal contact model that specifies material parameters :math:`E_i` and
:math:`\nu_i` is required (e.g. *hertz/material*, *dmt* or *jkr*\ ). In this
case, mixing of the shear modulus for different particle types *i* and
*j* is done according to the formula above.
@ -575,7 +577,7 @@ opposite torque on each particle, according to:
.. math::
\tau_{roll,i} = R_{eff} \mathbf{n} \times \mathbf{F}_{roll}
\tau_{roll,i} = R \mathbf{n} \times \mathbf{F}_{roll}
.. math::

54
doc/src/pair_hybrid.rst
View File

@ -1,28 +1,41 @@
.. index:: pair_style hybrid
.. index:: pair_style hybrid/kk
.. index:: pair_style hybrid/omp
.. index:: pair_style hybrid/molecular
.. index:: pair_style hybrid/molecular/omp
.. index:: pair_style hybrid/overlay
.. index:: pair_style hybrid/overlay/omp
.. index:: pair_style hybrid/overlay/kk
.. index:: pair_style hybrid/scaled
.. index:: pair_style hybrid/scaled/omp
pair_style hybrid command
=========================
Accelerator Variants: *hybrid/kk*
Accelerator Variants: *hybrid/kk*, *hybrid/omp*
pair_style hybrid/molecular command
===================================
Accelerator Variant: *hybrid/molecular/omp*
pair_style hybrid/overlay command
=================================
Accelerator Variants: *hybrid/overlay/kk*
Accelerator Variants: *hybrid/overlay/kk*, *hybrid/overlay/omp*
pair_style hybrid/scaled command
==================================
Accelerator Variant: *hybrid/scaled/omp*
Syntax
""""""
.. code-block:: LAMMPS
pair_style hybrid style1 args style2 args ...
pair_style hybrid/molecular factor1 style1 args factor2 style 2 args
pair_style hybrid/overlay style1 args style2 args ...
pair_style hybrid/scaled factor1 style1 args factor2 style 2 args ...
@ -47,6 +60,10 @@ Examples
pair_coeff * * tersoff Si.tersoff Si
pair_coeff * * sw Si.sw Si
pair_style hybrid/molecular lj/cut 2.5 lj/cut 2.5
pair_coeff * * lj/cut 1 1.0 1.0
pair_coeff * * lj/cut 2 1.5 1.0
variable one equal ramp(1.0,0.0)
variable two equal 1.0-v_one
pair_style hybrid/scaled v_one lj/cut 2.5 v_two morse 2.5
@ -56,17 +73,26 @@ Examples
Description
"""""""""""
The *hybrid*, *hybrid/overlay*, and *hybrid/scaled* styles enable the
use of multiple pair styles in one simulation. With the *hybrid* style,
exactly one pair style is assigned to each pair of atom types. With the
*hybrid/overlay* and *hybrid/scaled* styles, one or more pair styles can
be assigned to each pair of atom types. The assignment of pair styles
to type pairs is made via the :doc:`pair_coeff <pair_coeff>` command.
The major difference between the *hybrid/overlay* and *hybrid/scaled*
styles is that the *hybrid/scaled* adds a scale factor for each
sub-style contribution to forces, energies and stresses. Because of the
added complexity, the *hybrid/scaled* style has more overhead and thus
may be slower than *hybrid/overlay*.
The *hybrid*, *hybrid/overlay*, *hybrid/molecular*, and *hybrid/scaled*
styles enable the use of multiple pair styles in one simulation. With
the *hybrid* style, exactly one pair style is assigned to each pair of
atom types. With the *hybrid/overlay* and *hybrid/scaled* styles, one
or more pair styles can be assigned to each pair of atom types. With
the hybrid/molecular style, pair styles are assigned to either intra-
or inter-molecular interactions.
The assignment of pair styles to type pairs is made via the
:doc:`pair_coeff <pair_coeff>` command. The major difference between
the *hybrid/overlay* and *hybrid/scaled* styles is that the
*hybrid/scaled* adds a scale factor for each sub-style contribution to
forces, energies and stresses. Because of the added complexity, the
*hybrid/scaled* style has more overhead and thus may be slower than
*hybrid/overlay*.
The *hybrid/molecular* pair style accepts *only* two sub-styles: the
first is assigned to intra-molecular interactions (i.e. both atoms
have the same molecule ID), the second to inter-molecular interactions
(i.e. interacting atoms have different molecule IDs).
Here are two examples of hybrid simulations. The *hybrid* style could
be used for a simulation of a metal droplet on a LJ surface. The metal
@ -476,6 +502,8 @@ the same or else LAMMPS will generate an error.
Pair style *hybrid/scaled* currently only works for non-accelerated
pair styles and pair styles from the OPT package.
Pair style *hybrid/molecular* is not compatible with manybody potentials.
When using pair styles from the GPU package they must not be listed
multiple times. LAMMPS will detect this and abort.

182
doc/src/pair_oxdna.rst
View File

@ -37,18 +37,19 @@ Syntax
*oxdna/stk* args = seq T xi kappa 6.0 0.4 0.9 0.32 0.75 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65
seq = seqav (for average sequence stacking strength) or seqdep (for sequence-dependent stacking strength)
T = temperature (oxDNA units, 0.1 = 300 K)
xi = 1.3448 (temperature-independent coefficient in stacking strength)
kappa = 2.6568 (coefficient of linear temperature dependence in stacking strength)
T = temperature (LJ units: 0.1 = 300 K, real units: 300 = 300 K)
xi = 1.3448 (LJ units) or 8.01727944817084 (real units), temperature-independent coefficient in stacking strength
kappa = 2.6568 (LJ units) or 0.005279604 (real units), coefficient of linear temperature dependence in stacking strength
*oxdna/hbond* args = seq eps 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
seq = seqav (for average sequence base-pairing strength) or seqdep (for sequence-dependent base-pairing strength)
eps = 1.077 (between base pairs A-T and C-G) or 0 (all other pairs)
eps = 1.077 (LJ units) or 6.42073911784652 (real units), average hydrogen bonding strength between A-T and C-G Watson-Crick base pairs, 0 between all other pairs
Examples
""""""""
.. code-block:: LAMMPS
# LJ units
pair_style hybrid/overlay oxdna/excv oxdna/stk oxdna/hbond oxdna/xstk oxdna/coaxstk
pair_coeff * * oxdna/excv 2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32
pair_coeff * * oxdna/stk seqdep 0.1 1.3448 2.6568 6.0 0.4 0.9 0.32 0.75 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65
@ -58,55 +59,105 @@ Examples
pair_coeff * * oxdna/xstk 47.5 0.575 0.675 0.495 0.655 2.25 0.791592653589793 0.58 1.7 1.0 0.68 1.7 1.0 0.68 1.5 0 0.65 1.7 0.875 0.68 1.7 0.875 0.68
pair_coeff * * oxdna/coaxstk 46.0 0.4 0.6 0.22 0.58 2.0 2.541592653589793 0.65 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 -0.65 2.0 -0.65
pair_style hybrid/overlay oxdna/excv oxdna/stk oxdna/hbond oxdna/xstk oxdna/coaxstk
pair_coeff * * oxdna/excv oxdna_lj.cgdna
pair_coeff * * oxdna/stk seqav 0.1 1.3448 2.6568 oxdna_lj.cgdna
pair_coeff * * oxdna/hbond seqav oxdna_lj.cgdna
pair_coeff 1 4 oxdna/hbond seqav oxdna_lj.cgdna
pair_coeff 2 3 oxdna/hbond seqav oxdna_lj.cgdna
pair_coeff * * oxdna/xstk oxdna_lj.cgdna
pair_coeff * * oxdna/coaxstk oxdna_lj.cgdna
# Real units
pair_style hybrid/overlay oxdna/excv oxdna/stk oxdna/hbond oxdna/xstk oxdna/coaxstk
pair_coeff * * oxdna/excv 11.92337812042065 5.9626 5.74965 11.92337812042065 4.38677 4.259 11.92337812042065 2.81094 2.72576
pair_coeff * * oxdna/stk seqdep 300.0 8.01727944817084 0.005279604 0.70439070204273 3.4072 7.6662 2.72576 6.3885 1.3 0.0 0.8 0.9 0.0 0.95 0.9 0.0 0.95 2.0 0.65 2.0 0.65
pair_coeff * * oxdna/hbond seqdep 0.0 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592654 0.7 4.0 1.570796327 0.45 4.0 1.570796327 0.45
pair_coeff 1 4 oxdna/hbond seqdep 6.42073911784652 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592654 0.7 4.0 1.570796327 0.45 4.0 1.570796327 0.45
pair_coeff 2 3 oxdna/hbond seqdep 6.42073911784652 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592654 0.7 4.0 1.570796327 0.45 4.0 1.570796327 0.45
pair_coeff * * oxdna/xstk 3.9029021145006 4.89785 5.74965 4.21641 5.57929 2.25 0.791592654 0.58 1.7 1.0 0.68 1.7 1.0 0.68 1.5 0.0 0.65 1.7 0.875 0.68 1.7 0.875 0.68
pair_coeff * * oxdna/coaxstk 3.77965257404268 3.4072 5.1108 1.87396 4.94044 2.0 2.541592654 0.65 1.3 0.0 0.8 0.9 0.0 0.95 0.9 0.0 0.95 2.0 -0.65 2.0 -0.65
pair_style hybrid/overlay oxdna/excv oxdna/stk oxdna/hbond oxdna/xstk oxdna/coaxstk
pair_coeff * * oxdna/excv oxdna_real.cgdna
pair_coeff * * oxdna/stk seqav 300.0 8.01727944817084 0.005279604 oxdna_real.cgdna
pair_coeff * * oxdna/hbond seqav oxdna_real.cgdna
pair_coeff 1 4 oxdna/hbond seqav oxdna_real.cgdna
pair_coeff 2 3 oxdna/hbond seqav oxdna_real.cgdna
pair_coeff * * oxdna/xstk oxdna_real.cgdna
pair_coeff * * oxdna/coaxstk oxdna_real.cgdna
.. note::
The coefficients in the above examples are provided in forms
compatible with both *units lj* and *units real* (see documentation
of :doc:`units <units>`). These can also be read from a potential
file with correct unit style by specifying the name of the
file. Several potential files for each unit style are included in the
``potentials`` directory of the LAMMPS distribution.
Description
"""""""""""
The *oxdna* pair styles compute the pairwise-additive parts of the oxDNA force field
for coarse-grained modelling of DNA. The effective interaction between the nucleotides consists of potentials for the
excluded volume interaction *oxdna/excv*, the stacking *oxdna/stk*, cross-stacking *oxdna/xstk*
and coaxial stacking interaction *oxdna/coaxstk* as well
as the hydrogen-bonding interaction *oxdna/hbond* between complementary pairs of nucleotides on
opposite strands. Average sequence or sequence-dependent stacking and base-pairing strengths
are supported :ref:`(Sulc) <Sulc1>`. Quasi-unique base-pairing between nucleotides can be achieved by using
more complementary pairs of atom types like 5-8 and 6-7, 9-12 and 10-11, 13-16 and 14-15, etc.
This prevents the hybridization of in principle complementary bases within Ntypes/4 bases
up and down along the backbone.
The *oxdna* pair styles compute the pairwise-additive parts of the oxDNA
force field for coarse-grained modelling of DNA. The effective
interaction between the nucleotides consists of potentials for the
excluded volume interaction *oxdna/excv*, the stacking *oxdna/stk*,
cross-stacking *oxdna/xstk* and coaxial stacking interaction
*oxdna/coaxstk* as well as the hydrogen-bonding interaction
*oxdna/hbond* between complementary pairs of nucleotides on opposite
strands. Average sequence or sequence-dependent stacking and
base-pairing strengths are supported :ref:`(Sulc) <Sulc1>`. Quasi-unique
base-pairing between nucleotides can be achieved by using more
complementary pairs of atom types like 5-8 and 6-7, 9-12 and 10-11,
13-16 and 14-15, etc. This prevents the hybridization of in principle
complementary bases within Ntypes/4 bases up and down along the
backbone.
The exact functional form of the pair styles is rather complex.
The individual potentials consist of products of modulation factors,
which themselves are constructed from a number of more basic potentials
(Morse, Lennard-Jones, harmonic angle and distance) as well as quadratic smoothing and modulation terms.
We refer to :ref:`(Ouldridge-DPhil) <Ouldridge-DPhil1>` and :ref:`(Ouldridge) <Ouldridge1>`
for a detailed description of the oxDNA force field.
The exact functional form of the pair styles is rather complex. The
individual potentials consist of products of modulation factors, which
themselves are constructed from a number of more basic potentials
(Morse, Lennard-Jones, harmonic angle and distance) as well as quadratic
smoothing and modulation terms. We refer to :ref:`(Ouldridge-DPhil)
<Ouldridge-DPhil1>` and :ref:`(Ouldridge) <Ouldridge1>` for a detailed
description of the oxDNA force field.
.. note::
These pair styles have to be used together with the related oxDNA bond style
*oxdna/fene* for the connectivity of the phosphate backbone (see also documentation of
:doc:`bond_style oxdna/fene <bond_oxdna>`). Most of the coefficients
in the above example have to be kept fixed and cannot be changed without reparameterizing the entire model.
Exceptions are the first four coefficients after *oxdna/stk* (seq=seqdep, T=0.1, xi=1.3448 and kappa=2.6568 in the above example)
and the first coefficient after *oxdna/hbond* (seq=seqdep in the above example).
When using a Langevin thermostat, e.g. through :doc:`fix langevin <fix_langevin>`
or :doc:`fix nve/dotc/langevin <fix_nve_dotc_langevin>`
the temperature coefficients have to be matched to the one used in the fix.
These pair styles have to be used together with the related oxDNA
bond style *oxdna/fene* for the connectivity of the phosphate
backbone (see also documentation of :doc:`bond_style oxdna/fene
<bond_oxdna>`). Most of the coefficients in the above example have to
be kept fixed and cannot be changed without reparameterizing the
entire model. Exceptions are the first four coefficients after
*oxdna/stk* (seq=seqdep, T=0.1, xi=1.3448 and kappa=2.6568 and
corresponding *real unit* equivalents in the above examples) and the
first coefficient after *oxdna/hbond* (seq=seqdep in the above
example). When using a Langevin thermostat, e.g. through :doc:`fix
langevin <fix_langevin>` or :doc:`fix nve/dotc/langevin
<fix_nve_dotc_langevin>` the temperature coefficients have to be
matched to the one used in the fix.
.. note::
These pair styles have to be used with the *atom_style hybrid bond ellipsoid oxdna*
(see documentation of :doc:`atom_style <atom_style>`). The *atom_style oxdna*
stores the 3'-to-5' polarity of the nucleotide strand, which is set through
the bond topology in the data file. The first (second) atom in a bond definition
is understood to point towards the 3'-end (5'-end) of the strand.
These pair styles have to be used with the *atom_style hybrid bond
ellipsoid oxdna* (see documentation of :doc:`atom_style
<atom_style>`). The *atom_style oxdna* stores the 3'-to-5' polarity
of the nucleotide strand, which is set through the bond topology in
the data file. The first (second) atom in a bond definition is
understood to point towards the 3'-end (5'-end) of the strand.
Example input and data files for DNA duplexes can be found in examples/PACKAGES/cgdna/examples/oxDNA/ and /oxDNA2/.
A simple python setup tool which creates single straight or helical DNA strands,
DNA duplexes or arrays of DNA duplexes can be found in examples/PACKAGES/cgdna/util/.
Example input and data files for DNA duplexes can be found in
``examples/PACKAGES/cgdna/examples/oxDNA/`` and ``.../oxDNA2/``. A
simple python setup tool which creates single straight or helical DNA
strands, DNA duplexes or arrays of DNA duplexes can be found in
``examples/PACKAGES/cgdna/util/``.
Please cite :ref:`(Henrich) <Henrich1>` in any publication that uses
this implementation. An updated documentation that contains general information
on the model, its implementation and performance as well as the structure of
the data and input file can be found `here <PDF/CG-DNA.pdf>`_.
this implementation. An updated documentation that contains general
information on the model, its implementation and performance as well as
the structure of the data and input file can be found `here
<PDF/CG-DNA.pdf>`_.
Please cite also the relevant oxDNA publications
:ref:`(Ouldridge) <Ouldridge1>`,
@ -115,6 +166,57 @@ and :ref:`(Sulc) <Sulc1>`.
----------
Potential file reading
""""""""""""""""""""""
For each pair style above the first non-modifiable argument can be a
filename, and if it is, no further arguments should be
supplied. Therefore the following command:
.. code-block:: LAMMPS
pair_coeff 1 4 oxdna/hbond seqav oxdna_lj.cgdna
will be interpreted as a request to read the corresponding hydrogen
bonding potential parameters from the file with the given name. The file
can define multiple potential parameters for both bonded and pair
interactions, but for the example pair interaction above there must
exist in the file a line of the form:
.. code-block:: LAMMPS
1 4 hbond <coefficients>
If potential customization is required, the potential file reading can
be mixed with the manual specification of the potential parameters. For
example, the following command:
.. code-block:: LAMMPS
pair_style hybrid/overlay oxdna/excv oxdna/stk oxdna/hbond oxdna/xstk oxdna/coaxstk
pair_coeff * * oxdna/excv oxdna_lj.cgdna
pair_coeff * * oxdna/stk seqav 0.1 1.3448 2.6568 6.0 0.4 0.9 0.32 0.75 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65
pair_coeff * * oxdna/hbond seqav oxdna_lj.cgdna
pair_coeff 1 4 oxdna/hbond seqav oxdna_lj.cgdna
pair_coeff 2 3 oxdna/hbond seqav oxdna_lj.cgdna
pair_coeff * * oxdna/xstk oxdna_lj.cgdna
pair_coeff * * oxdna/coaxstk 46.0 0.4 0.6 0.22 0.58 2.0 2.541592653589793 0.65 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 -0.65 2.0 -0.65
will read the stacking and coaxial stacking potential parameters from
the manual specification and all others from the potential file
*oxdna_lj.cgdna*.
There are sample potential files for each unit style in the
``potentials`` directory of the LAMMPS distribution. The potential file
unit system must align with the units defined via the :doc:`units
<units>` command. For conversion between different *LJ* and *real* unit
systems for oxDNA, the python tool *lj2real.py* located in the
``examples/PACKAGES/cgdna/util/`` directory can be used. This tool
assumes similar file structure to the examples found in
``examples/PACKAGES/cgdna/examples/``.
----------
Restrictions
""""""""""""

189
doc/src/pair_oxdna2.rst
View File

@ -41,14 +41,14 @@ Syntax
*oxdna2/stk* args = seq T xi kappa 6.0 0.4 0.9 0.32 0.75 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65
seq = seqav (for average sequence stacking strength) or seqdep (for sequence-dependent stacking strength)
T = temperature (oxDNA units, 0.1 = 300 K)
xi = 1.3523 (temperature-independent coefficient in stacking strength)
kappa = 2.6717 (coefficient of linear temperature dependence in stacking strength)
T = temperature (LJ units: 0.1 = 300 K, real units: 300 = 300 K)
xi = 1.3523 (LJ units) or 8.06199211612242 (real units), temperature-independent coefficient in stacking strength
kappa = 2.6717 (LJ units) or 0.005309213 (real units), coefficient of linear temperature dependence in stacking strength
*oxdna2/hbond* args = seq eps 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
seq = seqav (for average sequence base-pairing strength) or seqdep (for sequence-dependent base-pairing strength)
eps = 1.0678 (between base pairs A-T and C-G) or 0 (all other pairs)
eps = 1.0678 (LJ units) or 6.36589157849259 (real units), average hydrogen bonding strength between A-T and C-G Watson-Crick base pairs, 0 between all other pairs
*oxdna2/dh* args = T rhos qeff
T = temperature (oxDNA units, 0.1 = 300 K)
T = temperature (LJ units: 0.1 = 300 K, real units: 300 = 300 K)
rhos = salt concentration (mole per litre)
qeff = 0.815 (effective charge in elementary charges)
@ -57,6 +57,7 @@ Examples
.. code-block:: LAMMPS
# LJ units
pair_style hybrid/overlay oxdna2/excv oxdna2/stk oxdna2/hbond oxdna2/xstk oxdna2/coaxstk oxdna2/dh
pair_coeff * * oxdna2/excv 2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32
pair_coeff * * oxdna2/stk seqdep 0.1 1.3523 2.6717 6.0 0.4 0.9 0.32 0.75 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65
@ -67,61 +68,169 @@ Examples
pair_coeff * * oxdna2/coaxstk 58.5 0.4 0.6 0.22 0.58 2.0 2.891592653589793 0.65 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 40.0 3.116592653589793
pair_coeff * * oxdna2/dh 0.1 0.5 0.815
pair_style hybrid/overlay oxdna2/excv oxdna2/stk oxdna2/hbond oxdna2/xstk oxdna2/coaxstk oxdna2/dh
pair_coeff * * oxdna2/excv oxdna2_lj.cgdna
pair_coeff * * oxdna2/stk seqdep 0.1 1.3523 2.6717 oxdna2_lj.cgdna
pair_coeff * * oxdna2/hbond seqdep oxdna2_lj.cgdna
pair_coeff 1 4 oxdna2/hbond seqdep oxdna2_lj.cgdna
pair_coeff 2 3 oxdna2/hbond seqdep oxdna2_lj.cgdna
pair_coeff * * oxdna2/xstk oxdna2_lj.cgdna
pair_coeff * * oxdna2/coaxstk oxdna2_lj.cgdna
pair_coeff * * oxdna2/dh 0.1 0.5 oxdna2_lj.cgdna
# Real units
pair_style hybrid/overlay oxdna2/excv oxdna2/stk oxdna2/hbond oxdna2/xstk oxdna2/coaxstk oxdna2/dh
pair_coeff * * oxdna2/excv 11.92337812042065 5.9626 5.74965 11.92337812042065 4.38677 4.259 11.92337812042065 2.81094 2.72576
pair_coeff * * oxdna2/stk seqdep 300.0 8.06199211612242 0.005309213 0.70439070204273 3.4072 7.6662 2.72576 6.3885 1.3 0.0 0.8 0.9 0.0 0.95 0.9 0.0 0.95 2.0 0.65 2.0 0.65
pair_coeff * * oxdna2/hbond seqdep 0.0 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592654 0.7 4.0 1.570796327 0.45 4.0 1.570796327 0.45
pair_coeff 1 4 oxdna2/hbond seqdep 6.36589157849259 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592654 0.7 4.0 1.570796327 0.45 4.0 1.570796327 0.45
pair_coeff 2 3 oxdna2/hbond seqdep 6.36589157849259 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592654 0.7 4.0 1.570796327 0.45 4.0 1.570796327 0.45
pair_coeff * * oxdna2/xstk 3.9029021145006 4.89785 5.74965 4.21641 5.57929 2.25 0.791592654 0.58 1.7 1.0 0.68 1.7 1.0 0.68 1.5 0.0 0.65 1.7 0.875 0.68 1.7 0.875 0.68
pair_coeff * * oxdna2/coaxstk 4.80673207785863 3.4072 5.1108 1.87396 4.94044 2.0 2.891592653589793 0.65 1.3 0.0 0.8 0.9 0.0 0.95 0.9 0.0 0.95 40.0 3.116592653589793
pair_coeff * * oxdna2/dh 300.0 0.5 0.815
pair_style hybrid/overlay oxdna2/excv oxdna2/stk oxdna2/hbond oxdna2/xstk oxdna2/coaxstk oxdna2/dh
pair_coeff * * oxdna2/excv oxdna2_real.cgdna
pair_coeff * * oxdna2/stk seqdep 300.0 8.06199211612242 0.005309213 oxdna2_real.cgdna
pair_coeff * * oxdna2/hbond seqdep oxdna2_real.cgdna
pair_coeff 1 4 oxdna2/hbond seqdep oxdna2_real.cgdna
pair_coeff 2 3 oxdna2/hbond seqdep oxdna2_real.cgdna
pair_coeff * * oxdna2/xstk oxdna2_real.cgdna
pair_coeff * * oxdna2/coaxstk oxdna2_real.cgdna
pair_coeff * * oxdna2/dh 300.0 0.5 oxdna2_real.cgdna
.. note::
The coefficients in the above examples are provided in forms
compatible with both *units lj* and *units real* (see documentation
of :doc:`units <units>`). These can also be read from a potential
file with correct unit style by specifying the name of the
file. Several potential files for each unit style are included in the
``potentials`` directory of the LAMMPS distribution.
Description
"""""""""""
The *oxdna2* pair styles compute the pairwise-additive parts of the oxDNA force field
for coarse-grained modelling of DNA. The effective interaction between the nucleotides consists of potentials for the
excluded volume interaction *oxdna2/excv*, the stacking *oxdna2/stk*, cross-stacking *oxdna2/xstk*
and coaxial stacking interaction *oxdna2/coaxstk*, electrostatic Debye-Hueckel interaction *oxdna2/dh*
as well as the hydrogen-bonding interaction *oxdna2/hbond* between complementary pairs of nucleotides on
opposite strands. Average sequence or sequence-dependent stacking and base-pairing strengths
are supported :ref:`(Sulc) <Sulc2>`. Quasi-unique base-pairing between nucleotides can be achieved by using
more complementary pairs of atom types like 5-8 and 6-7, 9-12 and 10-11, 13-16 and 14-15, etc.
This prevents the hybridization of in principle complementary bases within Ntypes/4 bases
The *oxdna2* pair styles compute the pairwise-additive parts of the
oxDNA force field for coarse-grained modelling of DNA. The effective
interaction between the nucleotides consists of potentials for the
excluded volume interaction *oxdna2/excv*, the stacking *oxdna2/stk*,
cross-stacking *oxdna2/xstk* and coaxial stacking interaction
*oxdna2/coaxstk*, electrostatic Debye-Hueckel interaction *oxdna2/dh* as
well as the hydrogen-bonding interaction *oxdna2/hbond* between
complementary pairs of nucleotides on opposite strands. Average sequence
or sequence-dependent stacking and base-pairing strengths are supported
:ref:`(Sulc) <Sulc2>`. Quasi-unique base-pairing between nucleotides can
be achieved by using more complementary pairs of atom types like 5-8 and
6-7, 9-12 and 10-11, 13-16 and 14-15, etc. This prevents the
hybridization of in principle complementary bases within Ntypes/4 bases
up and down along the backbone.
The exact functional form of the pair styles is rather complex.
The individual potentials consist of products of modulation factors,
which themselves are constructed from a number of more basic potentials
(Morse, Lennard-Jones, harmonic angle and distance) as well as quadratic smoothing and modulation terms.
We refer to :ref:`(Snodin) <Snodin2>` and the original oxDNA publications :ref:`(Ouldridge-DPhil) <Ouldridge-DPhil2>`
and :ref:`(Ouldridge) <Ouldridge2>` for a detailed description of the oxDNA2 force field.
The exact functional form of the pair styles is rather complex. The
individual potentials consist of products of modulation factors, which
themselves are constructed from a number of more basic potentials
(Morse, Lennard-Jones, harmonic angle and distance) as well as quadratic
smoothing and modulation terms. We refer to :ref:`(Snodin) <Snodin2>`
and the original oxDNA publications :ref:`(Ouldridge-DPhil)
<Ouldridge-DPhil2>` and :ref:`(Ouldridge) <Ouldridge2>` for a detailed
description of the oxDNA2 force field.
.. note::
These pair styles have to be used together with the related oxDNA2 bond style
*oxdna2/fene* for the connectivity of the phosphate backbone (see also documentation of
:doc:`bond_style oxdna2/fene <bond_oxdna>`). Most of the coefficients
in the above example have to be kept fixed and cannot be changed without reparameterizing the entire model.
Exceptions are the first four coefficients after *oxdna2/stk* (seq=seqdep, T=0.1, xi=1.3523 and kappa=2.6717 in the above example),
the first coefficient after *oxdna2/hbond* (seq=seqdep in the above example) and the three coefficients
after *oxdna2/dh* (T=0.1, rhos=0.5, qeff=0.815 in the above example). When using a Langevin thermostat
e.g. through :doc:`fix langevin <fix_langevin>` or :doc:`fix nve/dotc/langevin <fix_nve_dotc_langevin>`
the temperature coefficients have to be matched to the one used in the fix.
These pair styles have to be used together with the related oxDNA2
bond style *oxdna2/fene* for the connectivity of the phosphate
backbone (see also documentation of :doc:`bond_style oxdna2/fene
<bond_oxdna>`). Most of the coefficients in the above example have to
be kept fixed and cannot be changed without reparameterizing the
entire model. Exceptions are the first four coefficients after
*oxdna2/stk* (seq=seqdep, T=0.1, xi=1.3523 and kappa=2.6717 and
corresponding *real unit* equivalents in the above examples). the
first coefficient after *oxdna2/hbond* (seq=seqdep in the above
example) and the three coefficients after *oxdna2/dh* (T=0.1,
rhos=0.5, qeff=0.815 in the above example). When using a Langevin
thermostat e.g. through :doc:`fix langevin <fix_langevin>` or
:doc:`fix nve/dotc/langevin <fix_nve_dotc_langevin>` the temperature
coefficients have to be matched to the one used in the fix.
.. note::
These pair styles have to be used with the *atom_style hybrid bond ellipsoid oxdna*
(see documentation of :doc:`atom_style <atom_style>`). The *atom_style oxdna*
stores the 3'-to-5' polarity of the nucleotide strand, which is set through
the bond topology in the data file. The first (second) atom in a bond definition
is understood to point towards the 3'-end (5'-end) of the strand.
These pair styles have to be used with the *atom_style hybrid bond
ellipsoid oxdna* (see documentation of :doc:`atom_style
<atom_style>`). The *atom_style oxdna* stores the 3'-to-5' polarity
of the nucleotide strand, which is set through the bond topology in
the data file. The first (second) atom in a bond definition is
understood to point towards the 3'-end (5'-end) of the strand.
Example input and data files for DNA duplexes can be found in examples/PACKAGES/cgdna/examples/oxDNA/ and /oxDNA2/.
A simple python setup tool which creates single straight or helical DNA strands,
DNA duplexes or arrays of DNA duplexes can be found in examples/PACKAGES/cgdna/util/.
Example input and data files for DNA duplexes can be found in
``examples/PACKAGES/cgdna/examples/oxDNA/`` and ``.../oxDNA2/``. A
simple python setup tool which creates single straight or helical DNA
strands, DNA duplexes or arrays of DNA duplexes can be found in
``examples/PACKAGES/cgdna/util/``.
Please cite :ref:`(Henrich) <Henrich2>` in any publication that uses
this implementation. An updated documentation that contains general information
on the model, its implementation and performance as well as the structure of
the data and input file can be found `here <PDF/CG-DNA.pdf>`_.
this implementation. An updated documentation that contains general
information on the model, its implementation and performance as well as
the structure of the data and input file can be found `here
<PDF/CG-DNA.pdf>`_.
Please cite also the relevant oxDNA2 publications
:ref:`(Snodin) <Snodin2>` and :ref:`(Sulc) <Sulc2>`.
----------
Potential file reading
""""""""""""""""""""""
For each pair style above the first non-modifiable argument can be a
filename (with exception of Debye-Hueckel, for which the effective
charge argument can be a filename), and if it is, no further arguments
should be supplied. Therefore the following command:
.. code-block:: LAMMPS
pair_coeff 1 4 oxdna2/hbond seqdep oxdna_real.cgdna
will be interpreted as a request to read the corresponding hydrogen
bonding potential parameters from the file with the given name. The
file can define multiple potential parameters for both bonded and pair
interactions, but for the example pair interaction above there must
exist in the file a line of the form:
.. code-block:: LAMMPS
1 4 hbond <coefficients>
If potential customization is required, the potential file reading can
be mixed with the manual specification of the potential parameters. For
example, the following command:
.. code-block:: LAMMPS
pair_style hybrid/overlay oxdna2/excv oxdna2/stk oxdna2/hbond oxdna2/xstk oxdna2/coaxstk oxdna2/dh
pair_coeff * * oxdna2/excv 2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32
pair_coeff * * oxdna2/stk seqdep 0.1 1.3523 2.6717 oxdna2_lj.cgdna
pair_coeff * * oxdna2/hbond seqdep oxdna2_lj.cgdna
pair_coeff 1 4 oxdna2/hbond seqdep oxdna2_lj.cgdna
pair_coeff 2 3 oxdna2/hbond seqdep oxdna2_lj.cgdna
pair_coeff * * oxdna2/xstk oxdna2_lj.cgdna
pair_coeff * * oxdna2/coaxstk oxdna2_lj.cgdna
pair_coeff * * oxdna2/dh 0.1 0.5 0.815
will read the excluded volume and Debye-Hueckel effective charge *qeff*
parameters from the manual specification and all others from the
potential file *oxdna2_lj.cgdna*.
There are sample potential files for each unit style in the ``potentials``
directory of the LAMMPS distribution. The potential file unit system
must align with the units defined via the :doc:`units <units>`
command. For conversion between different *LJ* and *real* unit systems
for oxDNA, the python tool *lj2real.py* located in the
``examples/PACKAGES/cgdna/util/`` directory can be used. This tool assumes
similar file structure to the examples found in
``examples/PACKAGES/cgdna/examples/``.
----------
Restrictions
""""""""""""

197
doc/src/pair_oxrna2.rst
View File

@ -41,14 +41,14 @@ Syntax
*oxrna2/stk* args = seq T xi kappa 6.0 0.43 0.93 0.35 0.78 0.9 0 0.95 0.9 0 0.95 1.3 0 0.8 1.3 0 0.8 2.0 0.65 2.0 0.65
seq = seqav (for average sequence stacking strength) or seqdep (for sequence-dependent stacking strength)
T = temperature (oxDNA units, 0.1 = 300 K)
xi = 1.40206 (temperature-independent coefficient in stacking strength)
kappa = 2.77 (coefficient of linear temperature dependence in stacking strength)
T = temperature (LJ units: 0.1 = 300 K, real units: 300 = 300 K)
xi = 1.40206 (LJ units) or 8.35864576375849 (real units), temperature-independent coefficient in stacking strength
kappa = 2.77 (LJ units) or 0.005504556 (real units), coefficient of linear temperature dependence in stacking strength
*oxrna2/hbond* args = seq eps 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
seq = seqav (for average sequence base-pairing strength) or seqdep (for sequence-dependent base-pairing strength)
eps = 0.870439 (between base pairs A-T, C-G and G-T) or 0 (all other pairs)
eps = 0.870439 (LJ units) or 5.18928666388042 (real units), average hydrogen bonding strength between A-U and C-G Watson-Crick and G-U wobble base pairs, 0 between all other pairs
*oxrna2/dh* args = T rhos qeff
T = temperature (oxDNA units, 0.1 = 300 K)
T = temperature (LJ units: 0.1 = 300 K, real units: 300 = 300 K)
rhos = salt concentration (mole per litre)
qeff = 1.02455 (effective charge in elementary charges)
@ -57,6 +57,7 @@ Examples
.. code-block:: LAMMPS
# LJ units
pair_style hybrid/overlay oxrna2/excv oxrna2/stk oxrna2/hbond oxrna2/xstk oxrna2/coaxstk oxrna2/dh
pair_coeff * * oxrna2/excv 2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32
pair_coeff * * oxrna2/stk seqdep 0.1 1.40206 2.77 6.0 0.43 0.93 0.35 0.78 0.9 0 0.95 0.9 0 0.95 1.3 0 0.8 1.3 0 0.8 2.0 0.65 2.0 0.65
@ -68,58 +69,168 @@ Examples
pair_coeff * * oxrna2/coaxstk 80 0.5 0.6 0.42 0.58 2.0 2.592 0.65 1.3 0.151 0.8 0.9 0.685 0.95 0.9 0.685 0.95 2.0 -0.65 2.0 -0.65
pair_coeff * * oxrna2/dh 0.1 0.5 1.02455
pair_style hybrid/overlay oxrna2/excv oxrna2/stk oxrna2/hbond oxrna2/xstk oxrna2/coaxstk oxrna2/dh
pair_coeff * * oxrna2/excv oxrna2_lj.cgdna
pair_coeff * * oxrna2/stk seqdep 0.1 1.40206 2.77 oxrna2_lj.cgdna
pair_coeff * * oxrna2/hbond seqdep oxrna2_lj.cgdna
pair_coeff 1 4 oxrna2/hbond seqdep oxrna2_lj.cgdna
pair_coeff 2 3 oxrna2/hbond seqdep oxrna2_lj.cgdna
pair_coeff 3 4 oxrna2/hbond seqdep oxrna2_lj.cgdna
pair_coeff * * oxrna2/xstk oxrna2_lj.cgdna
pair_coeff * * oxrna2/coaxstk oxrna2_lj.cgdna
pair_coeff * * oxrna2/dh 0.1 0.5 oxrna2_lj.cgdna
# Real units
pair_style hybrid/overlay oxrna2/excv oxrna2/stk oxrna2/hbond oxrna2/xstk oxrna2/coaxstk oxrna2/dh
pair_coeff * * oxrna2/excv 11.92337812042065 5.9626 5.74965 11.92337812042065 4.38677 4.259 11.92337812042065 2.81094 2.72576
pair_coeff * * oxrna2/stk seqdep 300.0 8.35864576375849 0.005504556 0.70439070204273 3.66274 7.92174 2.9813 6.64404 0.9 0.0 0.95 0.9 0.0 0.95 1.3 0.0 0.8 1.3 0.0 0.8 2.0 0.65 2.0 0.65
pair_coeff * * oxrna2/hbond seqdep 0.0 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 1 4 oxrna2/hbond seqdep 5.18928666388042 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 2 3 oxrna2/hbond seqdep 5.18928666388042 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 3 4 oxrna2/hbond seqdep 5.18928666388042 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff * * oxrna2/xstk 4.92690859644113 4.259 5.1108 3.57756 4.94044 2.25 0.505 0.58 1.7 1.266 0.68 1.7 1.266 0.68 1.7 0.309 0.68 1.7 0.309 0.68
pair_coeff * * oxrna2/coaxstk 6.57330882442206 4.259 5.1108 3.57756 4.94044 2.0 2.592 0.65 1.3 0.151 0.8 0.9 0.685 0.95 0.9 0.685 0.95 2.0 -0.65 2.0 -0.65
pair_coeff * * oxrna2/dh 300.0 0.5 1.02455
pair_style hybrid/overlay oxrna2/excv oxrna2/stk oxrna2/hbond oxrna2/xstk oxrna2/coaxstk oxrna2/dh
pair_coeff * * oxrna2/excv oxrna2_real.cgdna
pair_coeff * * oxrna2/stk seqdep 300.0 8.35864576375849 0.005504556 oxrna2_real.cgdna
pair_coeff * * oxrna2/hbond seqdep oxrna2_real.cgdna
pair_coeff 1 4 oxrna2/hbond seqdep oxrna2_real.cgdna
pair_coeff 2 3 oxrna2/hbond seqdep oxrna2_real.cgdna
pair_coeff 3 4 oxrna2/hbond seqdep oxrna2_real.cgdna
pair_coeff * * oxrna2/xstk oxrna2_real.cgdna
pair_coeff * * oxrna2/coaxstk oxrna2_real.cgdna
pair_coeff * * oxrna2/dh 300.0 0.5 oxrna2_real.cgdna
.. note::
The coefficients in the above examples are provided in forms
compatible with both *units lj* and *units real* (see documentation
of :doc:`units <units>`). These can also be read from a potential
file with correct unit style by specifying the name of the
file. Several potential files for each unit style are included in the
``potentials`` directory of the LAMMPS distribution.
Description
"""""""""""
The *oxrna2* pair styles compute the pairwise-additive parts of the oxDNA force field
for coarse-grained modelling of DNA. The effective interaction between the nucleotides consists of potentials for the
excluded volume interaction *oxrna2/excv*, the stacking *oxrna2/stk*, cross-stacking *oxrna2/xstk*
and coaxial stacking interaction *oxrna2/coaxstk*, electrostatic Debye-Hueckel interaction *oxrna2/dh*
as well as the hydrogen-bonding interaction *oxrna2/hbond* between complementary pairs of nucleotides on
opposite strands. Average sequence or sequence-dependent stacking and base-pairing strengths
are supported :ref:`(Sulc2) <Sulc32>`. Quasi-unique base-pairing between nucleotides can be achieved by using
more complementary pairs of atom types like 5-8 and 6-7, 9-12 and 10-11, 13-16 and 14-15, etc.
This prevents the hybridization of in principle complementary bases within Ntypes/4 bases
The *oxrna2* pair styles compute the pairwise-additive parts of the
oxDNA force field for coarse-grained modelling of RNA. The effective
interaction between the nucleotides consists of potentials for the
excluded volume interaction *oxrna2/excv*, the stacking *oxrna2/stk*,
cross-stacking *oxrna2/xstk* and coaxial stacking interaction
*oxrna2/coaxstk*, electrostatic Debye-Hueckel interaction *oxrna2/dh* as
well as the hydrogen-bonding interaction *oxrna2/hbond* between
complementary pairs of nucleotides on opposite strands. Average sequence
or sequence-dependent stacking and base-pairing strengths are supported
:ref:`(Sulc2) <Sulc32>`. Quasi-unique base-pairing between nucleotides
can be achieved by using more complementary pairs of atom types like 5-8
and 6-7, 9-12 and 10-11, 13-16 and 14-15, etc. This prevents the
hybridization of in principle complementary bases within Ntypes/4 bases
up and down along the backbone.
The exact functional form of the pair styles is rather complex.
The individual potentials consist of products of modulation factors,
which themselves are constructed from a number of more basic potentials
(Morse, Lennard-Jones, harmonic angle and distance) as well as quadratic smoothing and modulation terms.
We refer to :ref:`(Sulc1) <Sulc31>` and the original oxDNA publications :ref:`(Ouldridge-DPhil) <Ouldridge-DPhil3>`
and :ref:`(Ouldridge) <Ouldridge3>` for a detailed description of the oxRNA2 force field.
The exact functional form of the pair styles is rather complex. The
individual potentials consist of products of modulation factors, which
themselves are constructed from a number of more basic potentials
(Morse, Lennard-Jones, harmonic angle and distance) as well as quadratic
smoothing and modulation terms. We refer to :ref:`(Sulc1) <Sulc31>` and
the original oxDNA publications :ref:`(Ouldridge-DPhil)
<Ouldridge-DPhil3>` and :ref:`(Ouldridge) <Ouldridge3>` for a detailed
description of the oxRNA2 force field.
.. note::
These pair styles have to be used together with the related oxDNA2 bond style
*oxrna2/fene* for the connectivity of the phosphate backbone (see also documentation of
:doc:`bond_style oxrna2/fene <bond_oxdna>`). Most of the coefficients
in the above example have to be kept fixed and cannot be changed without reparameterizing the entire model.
Exceptions are the first four coefficients after *oxrna2/stk* (seq=seqdep, T=0.1, xi=1.40206 and kappa=2.77 in the above example),
the first coefficient after *oxrna2/hbond* (seq=seqdep in the above example) and the three coefficients
after *oxrna2/dh* (T=0.1, rhos=0.5, qeff=1.02455 in the above example). When using a Langevin thermostat
e.g. through :doc:`fix langevin <fix_langevin>` or :doc:`fix nve/dotc/langevin <fix_nve_dotc_langevin>`
the temperature coefficients have to be matched to the one used in the fix.
These pair styles have to be used together with the related oxDNA2
bond style *oxrna2/fene* for the connectivity of the phosphate
backbone (see also documentation of :doc:`bond_style oxrna2/fene
<bond_oxdna>`). Most of the coefficients in the above example have to
be kept fixed and cannot be changed without reparameterizing the
entire model. Exceptions are the first four coefficients after
*oxrna2/stk* (seq=seqdep, T=0.1, xi=1.40206 and kappa=2.77 and
corresponding *real unit* equivalents in the above examples), the
first coefficient after *oxrna2/hbond* (seq=seqdep in the above
example) and the three coefficients after *oxrna2/dh* (T=0.1,
rhos=0.5, qeff=1.02455 in the above example). When using a Langevin
thermostat e.g. through :doc:`fix langevin <fix_langevin>` or
:doc:`fix nve/dotc/langevin <fix_nve_dotc_langevin>` the temperature
coefficients have to be matched to the one used in the fix.
.. note::
These pair styles have to be used with the *atom_style hybrid bond ellipsoid oxdna*
(see documentation of :doc:`atom_style <atom_style>`). The *atom_style oxdna*
stores the 3'-to-5' polarity of the nucleotide strand, which is set through
the bond topology in the data file. The first (second) atom in a bond definition
is understood to point towards the 3'-end (5'-end) of the strand.
These pair styles have to be used with the *atom_style hybrid bond
ellipsoid oxdna* (see documentation of :doc:`atom_style
<atom_style>`). The *atom_style oxdna* stores the 3'-to-5' polarity
of the nucleotide strand, which is set through the bond topology in
the data file. The first (second) atom in a bond definition is
understood to point towards the 3'-end (5'-end) of the strand.
Example input and data files for DNA duplexes can be found in examples/PACKAGES/cgdna/examples/oxDNA/ and /oxDNA2/.
A simple python setup tool which creates single straight or helical DNA strands,
DNA duplexes or arrays of DNA duplexes can be found in examples/PACKAGES/cgdna/util/.
Example input and data files for DNA duplexes can be found in
``examples/PACKAGES/cgdna/examples/oxDNA/`` and ``.../oxDNA2/``. A simple python
setup tool which creates single straight or helical DNA strands, DNA
duplexes or arrays of DNA duplexes can be found in
``examples/PACKAGES/cgdna/util/``.
Please cite :ref:`(Henrich) <Henrich3>` in any publication that uses
this implementation. The article contains general information
on the model, its implementation and performance as well as the structure of
the data and input file. The preprint version of the article can be found
`here <PDF/CG-DNA.pdf>`_.
Please cite also the relevant oxRNA2 publications
:ref:`(Sulc1) <Sulc31>` and :ref:`(Sulc2) <Sulc32>`.
this implementation. The article contains general information on the
model, its implementation and performance as well as the structure of
the data and input file. The preprint version of the article can be
found `here <PDF/CG-DNA.pdf>`_. Please cite also the relevant oxRNA2
publications :ref:`(Sulc1) <Sulc31>` and :ref:`(Sulc2) <Sulc32>`.
----------
Potential file reading
""""""""""""""""""""""
For each pair style above the first non-modifiable argument can be a
filename (with exception of Debye-Hueckel, for which the effective
charge argument can be a filename), and if it is, no further arguments
should be supplied. Therefore the following command:
.. code-block:: LAMMPS
pair_coeff 3 4 oxrna2/hbond seqdep oxrna2_lj.cgdna
will be interpreted as a request to read the corresponding hydrogen
bonding potential parameters from the file with the given name. The
file can define multiple potential parameters for both bonded and pair
interactions, but for the example pair interaction above there must
exist in the file a line of the form:
.. code-block:: LAMMPS
3 4 hbond <coefficients>
If potential customization is required, the potential file reading can
be mixed with the manual specification of the potential parameters. For
example, the following command:
.. code-block:: LAMMPS
pair_style hybrid/overlay oxrna2/excv oxrna2/stk oxrna2/hbond oxrna2/xstk oxrna2/coaxstk oxrna2/dh
pair_coeff * * oxrna2/excv 2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32
pair_coeff * * oxrna2/stk seqdep 0.1 1.40206 2.77 oxrna2_lj.cgdna
pair_coeff * * oxrna2/hbond seqdep oxrna2_lj.cgdna
pair_coeff 1 4 oxrna2/hbond seqdep oxrna2_lj.cgdna
pair_coeff 2 3 oxrna2/hbond seqdep oxrna2_lj.cgdna
pair_coeff 3 4 oxrna2/hbond seqdep oxrna2_lj.cgdna
pair_coeff * * oxrna2/xstk oxrna2_lj.cgdna
pair_coeff * * oxrna2/coaxstk oxrna2_lj.cgdna
pair_coeff * * oxrna2/dh 0.1 0.5 1.02455
will read the excluded volume and Debye-Hueckel effective charge *qeff*
parameters from the manual specification and all others from the
potential file *oxrna2_lj.cgdna*.
There are sample potential files for each unit style in the
``potentials`` directory of the LAMMPS distribution. The potential file
unit system must align with the units defined via the :doc:`units
<units>` command. For conversion between different *LJ* and *real* unit
systems for oxDNA, the python tool *lj2real.py* located in the
``examples/PACKAGES/cgdna/util/`` directory can be used. This tool
assumes similar file structure to the examples found in
``examples/PACKAGES/cgdna/examples/``.
----------

89
doc/src/pair_pod.rst
View File

@ -1,8 +1,11 @@
.. index:: pair_style pod
.. index:: pair_style pod/kk
pair_style pod command
========================
Accelerator Variants: *pod/kk*
Syntax
""""""
@ -24,29 +27,33 @@ Description
.. versionadded:: 22Dec2022
Pair style *pod* defines the proper orthogonal descriptor (POD)
potential :ref:`(Nguyen) <Nguyen20221>`. The mathematical definition of
the POD potential is described from :doc:`fitpod <fitpod_command>`, which is
used to fit the POD potential to *ab initio* energy and force data.
potential :ref:`(Nguyen and Rohskopf) <Nguyen20222b>`,
:ref:`(Nguyen2023) <Nguyen20232b>`, :ref:`(Nguyen2024) <Nguyen20242b>`,
and :ref:`(Nguyen and Sema) <Nguyen20243b>`. The :doc:`fitpod
<fitpod_command>` is used to fit the POD potential.
Only a single pair_coeff command is used with the *pod* style which
specifies a POD parameter file followed by a coefficient file.
specifies a POD parameter file followed by a coefficient file, a
projection matrix file, and a centroid file.
The coefficient file (``Ta_coefficients.pod``) contains coefficients for the
POD potential. The top of the coefficient file can contain any number of
blank and comment lines (start with #), but follows a strict format
after that. The first non-blank non-comment line must contain:
The POD parameter file (``Ta_param.pod``) can contain blank and comment
lines (start with #) anywhere. Each non-blank non-comment line must
contain one keyword/value pair. See :doc:`fitpod <fitpod_command>` for
the description of all the keywords that can be assigned in the
parameter file.
* POD_coefficients: *ncoeff*
The coefficient file (``Ta_coefficients.pod``) contains coefficients for
the POD potential. The top of the coefficient file can contain any
number of blank and comment lines (start with #), but follows a strict
format after that. The first non-blank non-comment line must contain:
This is followed by *ncoeff* coefficients, one per line. The coefficient
* model_coefficients: *ncoeff* *nproj* *ncentroid*
This is followed by *ncoeff* coefficients, *nproj* projection matrix entries,
and *ncentroid* centroid coordinates, one per line. The coefficient
file is generated after training the POD potential using :doc:`fitpod
<fitpod_command>`.
The POD parameter file (``Ta_param.pod``) can contain blank and comment lines
(start with #) anywhere. Each non-blank non-comment line must contain
one keyword/value pair. See :doc:`fitpod <fitpod_command>` for the description
of all the keywords that can be assigned in the parameter file.
As an example, if a LAMMPS indium phosphide simulation has 4 atoms
types, with the first two being indium and the third and fourth being
phophorous, the pair_coeff command would look like this:
@ -67,7 +74,33 @@ the *hybrid* pair style. The NULL values are placeholders for atom
types that will be used with other potentials.
Examples about training and using POD potentials are found in the
directory lammps/examples/PACKAGES/pod.
directory lammps/examples/PACKAGES/pod and the Github repo https://github.com/cesmix-mit/pod-examples.
----------
Mixing, shift, table, tail correction, restart, rRESPA info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
For atom type pairs I,J and I != J, where types I and J correspond to
two different element types, mixing is performed by LAMMPS with
user-specifiable parameters as described above. You never need to
specify a pair_coeff command with I != J arguments for this style.
This pair style does not support the :doc:`pair_modify <pair_modify>`
shift, table, and tail options.
This pair style does not write its information to :doc:`binary restart
files <restart>`, since it is stored in potential files. Thus, you need
to re-specify the pair_style and pair_coeff commands in an input script
that reads a restart file.
This pair style can only be used via the *pair* keyword of the
:doc:`run_style respa <run_style>` command. It does not support the
*inner*, *middle*, *outer* keywords.
----------
.. include:: accel_styles.rst
----------
@ -78,12 +111,14 @@ This style is part of the ML-POD package. It is only enabled if LAMMPS
was built with that package. See the :doc:`Build package
<Build_package>` page for more info.
This pair style does not compute per-atom energies and per-atom stresses.
Related commands
""""""""""""""""
:doc:`fitpod <fitpod_command>`,
:doc:`compute pod/atom <compute_pod_atom>`,
:doc:`compute podd/atom <compute_pod_atom>`,
:doc:`compute pod/local <compute_pod_atom>`,
:doc:`compute pod/global <compute_pod_atom>`
Default
"""""""
@ -92,6 +127,20 @@ none
----------
.. _Nguyen20221:
.. _Nguyen20222b:
**(Nguyen and Rohskopf)** Nguyen and Rohskopf, Journal of Computational Physics, 480, 112030, (2023).
.. _Nguyen20232b:
**(Nguyen2023)** Nguyen, Physical Review B, 107(14), 144103, (2023).
.. _Nguyen20242b:
**(Nguyen2024)** Nguyen, Journal of Computational Physics, 113102, (2024).
.. _Nguyen20243b:
**(Nguyen and Sema)** Nguyen and Sema, https://arxiv.org/abs/2405.00306, (2024).
**(Nguyen)** Nguyen and Rohskopf, arXiv preprint arXiv:2209.02362 (2022).

3
doc/src/pair_soft.rst
View File

@ -1,11 +1,12 @@
.. index:: pair_style soft
.. index:: pair_style soft/gpu
.. index:: pair_style soft/kk
.. index:: pair_style soft/omp
pair_style soft command
=======================
Accelerator Variants: *soft/gpu*, *soft/omp*
Accelerator Variants: *soft/gpu*, *soft/kk*, *soft/omp*
Syntax
""""""

3
doc/src/pair_style.rst
View File

@ -108,6 +108,7 @@ accelerated styles exist.
* :doc:`none <pair_none>` - turn off pairwise interactions
* :doc:`hybrid <pair_hybrid>` - multiple styles of pairwise interactions
* :doc:`hybrid/molecular <pair_hybrid>` - different pair styles for intra- and inter-molecular interactions
* :doc:`hybrid/overlay <pair_hybrid>` - multiple styles of superposed pairwise interactions
* :doc:`hybrid/scaled <pair_hybrid>` - multiple styles of scaled superposed pairwise interactions
* :doc:`zero <pair_zero>` - neighbor list but no interactions
@ -171,6 +172,7 @@ accelerated styles exist.
* :doc:`coul/wolf <pair_coul>` - Coulomb via Wolf potential
* :doc:`coul/wolf/cs <pair_cs>` - Coulomb via Wolf potential with core/shell adjustments
* :doc:`dpd <pair_dpd>` - dissipative particle dynamics (DPD)
* :doc:`dpd/coul/slater/long <pair_dpd_coul_slater_long>` - dissipative particle dynamics (DPD) with electrostatic interactions
* :doc:`dpd/ext <pair_dpd_ext>` - generalized force field for DPD
* :doc:`dpd/ext/tstat <pair_dpd_ext>` - pairwise DPD thermostatting with generalized force field
* :doc:`dpd/fdt <pair_dpd_fdt>` - DPD for constant temperature and pressure
@ -382,6 +384,7 @@ accelerated styles exist.
* :doc:`tracker <pair_tracker>` - monitor information about pairwise interactions
* :doc:`tri/lj <pair_tri_lj>` - LJ potential between triangles
* :doc:`ufm <pair_ufm>` -
* :doc:`uf3 <pair_uf3>` - UF3 machine-learning potential
* :doc:`vashishta <pair_vashishta>` - Vashishta 2-body and 3-body potential
* :doc:`vashishta/table <pair_vashishta>` -
* :doc:`wf/cut <pair_wf_cut>` - Wang-Frenkel Potential for short-ranged interactions

213
doc/src/pair_uf3.rst Normal file
View File

@ -0,0 +1,213 @@
.. index:: pair_style uf3
.. index:: pair_style uf3/kk
pair_style uf3 command
======================
Accelerator Variants: *uf3/kk*
Syntax
""""""
.. code-block:: LAMMPS
pair_style style BodyFlag
* style = *uf3* or *uf3/kk*
.. parsed-literal::
BodyFlag = Indicates whether to calculate only 2-body or 2 and 3-body interactions. Possible values: 2 or 3
Examples
""""""""
.. code-block:: LAMMPS
pair_style uf3 3
pair_coeff * * Nb.uf3 Nb
pair_style uf3 2
pair_coeff * * NbSn.uf3 Nb Sn
pair_style uf3 3
pair_coeff * * NbSn.uf3 Nb Sn
Description
"""""""""""
.. versionadded:: 27June2024
The *uf3* style computes the :ref:`Ultra-Fast Force Fields (UF3)
<Xie23>` potential, a machine-learning interatomic potential. In UF3,
the total energy of the system is defined via two- and three-body
interactions:
.. math::
E & = \sum_{i,j} V_2(r_{ij}) + \sum_{i,j,k} V_3 (r_{ij},r_{ik},r_{jk}) \\
V_2(r_{ij}) & = \sum_{n=0}^N c_n B_n(r_{ij}) \\
V_3 (r_{ij},r_{ik},r_{jk}) & = \sum_{l=0}^{N_l} \sum_{m=0}^{N_m} \sum_{n=0}^{N_n} c_{l,m,n} B_l(r_{ij}) B_m(r_{ik}) B_n(r_{jk})
where :math:`V_2(r_{ij})` and :math:`V_3 (r_{ij},r_{ik},r_{jk})` are the
two- and three-body interactions, respectively. For the two-body the
summation is over all neighbors J and for the three-body the summation
is over all neighbors J and K of atom I within a cutoff distance
determined from the potential files. :math:`B_n(r_{ij})` are the cubic
b-spline basis, :math:`c_n` and :math:`c_{l,m,n}` are the machine-learned
interaction parameters and :math:`N`, :math:`N_l`, :math:`N_m`, and
:math:`N_n` denote the number of basis functions per spline or tensor
spline dimension.
With *uf3* style only a single pair_coeff command is used to indicate the
UF3 LAMMPS potential file containing all the two- and three-body interactions
followed by N additional arguments specifying the mapping of UF3 elements to
LAMMPS atom types, where N is the number of LAMMPS atom types:
* UF3 LAMMPS potential file
* N elements names = mapping of UF3 elements to atom types
As an example, if a LAMMPS simulation contains 2 atom types (elements
'A' and 'B'), the pair_coeff command will be:
.. code-block:: LAMMPS
pair_style uf3 3
pair_coeff * * AB.uf3 A B
The AB.uf3 file should contain all two-body (A-A, A-B, B-B) and three-body
(A-A-A, A-A-B, A-B-B, B-A-A, B-A-B, B-B-B).
If a value of "2" is specified in the :code:`pair_style uf3` command,
only the two-body potentials are needed. For 3-body interaction the
first atom type is the central atom. We recommend using the
:code:`generate_uf3_lammps_pots.py` script (found `here
<https://github.com/uf3/uf3/tree/develop/lammps_plugin/scripts>`_) for
generating the UF3 LAMMPS potential file from the UF3 JSON potentials.
----------
UF3 LAMMPS potential file in the *potentials* directory of the LAMMPS
distribution have a ".uf3" suffix. The interaction block in UF3 LAMMPS potential
file should start with :code:`#UF3 POT` and end with :code:`#` characters.
Following shows the format of a generic 2-body and 3-body potential block in
UF3 LAMMPS potential file-
.. code-block:: LAMMPS
#UF3 POT UNITS: units DATE: POT_GEN_DATE AUTHOR: AUTHOR_NAME CITATION: CITE
2B ELEMENT1 ELEMENT2 LEADING_TRIM TRAILING_TRIM
Rij_CUTOFF NUM_OF_KNOTS
BSPLINE_KNOTS
NUM_OF_COEFF
COEFF
#
#UF3 POT UNITS: units DATE: POT_GEN_DATE AUTHOR: AUTHOR_NAME CITATION: CITE
3B ELEMENT1 ELEMENT2 ELEMENT3 LEADING_TRIM TRAILING_TRIM
Rjk_CUTOFF Rik_CUTOFF Rij_CUTOFF NUM_OF_KNOTS_JK NUM_OF_KNOTS_IK NUM_OF_KNOTS_IJ
BSPLINE_KNOTS_FOR_JK
BSPLINE_KNOTS_FOR_IK
BSPLINE_KNOTS_FOR_IJ
SHAPE_OF_COEFF_MATRIX[I][J][K]
COEFF_MATRIX[0][0][K]
COEFF_MATRIX[0][1][K]
COEFF_MATRIX[0][2][K]
.
.
.
COEFF_MATRIX[1][0][K]
COEFF_MATRIX[1][1][K]
COEFF_MATRIX[1][2][K]
.
.
.
#
The second line indicates whether the block contains data for 2-body
(:code:`2B`) or 3-body (:code:`3B`) interaction. This is followed by element
combination interaction, :code:`LEADING_TRIM` and :code:`TRAILING_TRIM`
number on the same line. The current implementation is only tested for
:code:`LEADING_TRIM=0` and :code:`TRAILING_TRIM=3`.
If other values are used LAMMPS is terminated after issuing an error message.
The :code:`Rij_CUTOFF` sets the 2-body cutoff for the interaction described
by the potential block. :code:`NUM_OF_KNOTS` is the number of knots
(or the length of the knot vector) present on the very next line. The
:code:`BSPLINE_KNOTS` line should contain all the knots in ascending order.
:code:`NUM_OF_COEFF` is the number of coefficients in the :code:`COEFF` line.
All the numbers in the BSPLINE_KNOTS and COEFF line should be space-separated.
Similar to the 2-body potential block, the third line sets the cutoffs and
length of the knots. The cutoff distance between atom-type I and J is
:code:`Rij_CUTOFF`, atom-type I and K is :code:`Rik_CUTOFF` and between
J and K is :code:`Rjk_CUTOFF`.
.. note::
The current implementation only works for UF3 potentials with cutoff
distances for 3-body interactions that follows
:code:`2Rij_CUTOFF=2Rik_CUTOFF=Rjk_CUTOFF` relation.
The :code:`BSPLINE_KNOTS_FOR_JK`, :code:`BSPLINE_KNOTS_FOR_IK`, and
:code:`BSPLINE_KNOTS_FOR_IJ` lines (note the order) contain the knots in
increasing order for atoms J and K, I and K, and atoms I and J
respectively. The number of knots is defined by the
:code:`NUM_OF_KNOTS_*` characters in the previous line. The shape of
the coefficient matrix is defined on the
:code:`SHAPE_OF_COEFF_MATRIX[I][J][K]` line followed by the columns of
the coefficient matrix, one per line, as shown above. For example, if
the coefficient matrix has the shape of 8x8x13, then
:code:`SHAPE_OF_COEFF_MATRIX[I][J][K]` will be :code:`8 8 13` followed
by 64 (8x8) lines each containing 13 coefficients separated by space.
----------
.. include:: accel_styles.rst
----------
Mixing, shift, table, tail correction, restart, rRESPA info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
For atom type pairs I,J and I != J, where types I and J correspond to
two different element types, mixing is performed by LAMMPS as described
above from values in the potential file.
This pair style does not support the :doc:`pair_modify <pair_modify>`
shift, table, and tail options.
This pair style does not write its information to :doc:`binary restart
files <restart>`, since it is stored in potential file.
This pair style can only be used via the *pair* keyword of the
:doc:`run_style respa <run_style>` command. It does not support the
*inner*, *middle*, *outer* keywords.
Restrictions
""""""""""""
The 'uf3' pair style is part of the ML-UF3 package. It is only enabled
if LAMMPS was built with that package. See the :doc:`Build package
<Build_package>` page for more info.
This pair style requires the :doc:`newton <newton>` setting to be "on".
The UF3 LAMMPS potential file provided with LAMMPS (see the potentials
directory) are parameterized for metal :doc:`units <units>`.
The single() function of 'uf3' pair style only return the 2-body
interaction energy.
Related commands
""""""""""""""""
:doc:`pair_coeff <pair_coeff>`
Default
"""""""
none
----------
.. _Xie23:
**(Xie23)** Xie, S.R., Rupp, M. & Hennig, R.G. Ultra-fast interpretable machine-learning potentials. npj Comput Mater 9, 162 (2023). https://doi.org/10.1038/s41524-023-01092-7

8
doc/src/run_style.rst
View File

@ -327,10 +327,12 @@ Restrictions
""""""""""""
The *verlet/split* style can only be used if LAMMPS was built with the
REPLICA package. Correspondingly the *respa/omp* style is available
REPLICA package. Correspondingly the *respa/omp* style is available
only if the OPENMP package was included. See the :doc:`Build package
<Build_package>` page for more info. It is not compatible with
kspace styles from the INTEL package.
<Build_package>` page for more info.
Run style *verlet/split* is not compatible with kspace styles from
the INTEL package and it is not compatible with any TIP4P styles.
Whenever using rRESPA, the user should experiment with trade-offs in
speed and accuracy for their system, and verify that they are

2
doc/src/variable.rst
View File

@ -957,7 +957,7 @@ of points, equally spaced by 1 in their x coordinate: (1,V1), (2,V2),
length N. The returned value is the slope of the line. If the line
has a single point or is vertical, it returns 1.0e20.
.. versionadded:: TBD
.. versionadded:: 27June2024
The sort(x) and rsort(x) functions sort the data of the input vector by
their numeric value: sort(x) sorts in ascending order, rsort(x) sorts

3
doc/utils/requirements.txt
View File

@ -1,6 +1,7 @@
Sphinx >= 5.3.0, <7.3
Sphinx >= 5.3.0, <8.0
sphinxcontrib-spelling
sphinxcontrib-jquery
sphinx-design
git+https://github.com/akohlmey/sphinx-fortran@parallel-read
sphinx-tabs>=3.4.1
breathe

1
doc/utils/sphinx-config/conf.py.in
View File

@ -57,6 +57,7 @@ extensions = [
'table_from_list',
'tab_or_note',
'breathe',
'sphinx_design'
]
images_config = {

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

@ -1173,6 +1173,7 @@ finitecutflag
Finnis
Fiorin
fitpod
fivebody
fixID
fj
Fji
@ -1434,6 +1435,7 @@ Hendrik
Henin
Henkelman
Henkes
Hennig
henrich
Henrich
Hermitian
@ -1595,6 +1597,7 @@ interlayer
intermolecular
interoperable
Interparticle
interpretable
interstitials
intertube
Intr
@ -1817,6 +1820,7 @@ Koziol
Kp
kradius
Kraker
Krass
Kraus
Kremer
Kress
@ -2894,6 +2898,7 @@ Pmolrotate
Pmoltrans
pN
png
podd
Podhorszki
Poiseuille
poisson
@ -2959,6 +2964,7 @@ Priya
proc
Proc
procs
procgrid
progguide
Prony
ps
@ -3271,6 +3277,7 @@ Rudranarayan
Rudzinski
Runge
runtime
Rupp
Rutuparna
rx
rxd
@ -3365,6 +3372,7 @@ setmask
Setmask
setpoint
setvel
sevenbody
sfftw
sfree
Sg
@ -3415,6 +3423,7 @@ SiO
Siochi
Sirk
Sival
sixbody
sizeI
sizeJ
sizex
@ -3730,7 +3739,6 @@ tokyo
tol
tomic
toolchain
toolset
topologies
Toporov
Torder
@ -3808,6 +3816,7 @@ typeJ
typelabel
typeN
typesafe
typestr
Tz
Tzou
ub
@ -3820,6 +3829,8 @@ uChem
uCond
uef
UEF
uf
uf3
ufm
Uhlenbeck
Ui

0
examples/PACKAGES/cgdna/examples/oxDNA/duplex1/data.duplex1 → examples/PACKAGES/cgdna/examples/lj_units/oxDNA/duplex1/data.duplex1
View File

0
examples/PACKAGES/cgdna/examples/oxDNA/duplex1/in.duplex1 → examples/PACKAGES/cgdna/examples/lj_units/oxDNA/duplex1/in.duplex1
View File

0
examples/PACKAGES/cgdna/examples/oxDNA/duplex1/log.2Jul21.duplex1.g++.1 → examples/PACKAGES/cgdna/examples/lj_units/oxDNA/duplex1/log.2Jul21.duplex1.g++.1
View File

0
examples/PACKAGES/cgdna/examples/oxDNA/duplex1/log.2Jul21.duplex1.g++.4 → examples/PACKAGES/cgdna/examples/lj_units/oxDNA/duplex1/log.2Jul21.duplex1.g++.4
View File

0
examples/PACKAGES/cgdna/examples/oxDNA/duplex2/data.duplex2 → examples/PACKAGES/cgdna/examples/lj_units/oxDNA/duplex2/data.duplex2
View File

0
examples/PACKAGES/cgdna/examples/oxDNA/duplex2/in.duplex2 → examples/PACKAGES/cgdna/examples/lj_units/oxDNA/duplex2/in.duplex2
View File

0
examples/PACKAGES/cgdna/examples/oxDNA/duplex2/log.2Jul21.duplex2.g++.1 → examples/PACKAGES/cgdna/examples/lj_units/oxDNA/duplex2/log.2Jul21.duplex2.g++.1
View File

0
examples/PACKAGES/cgdna/examples/oxDNA/duplex2/log.2Jul21.duplex2.g++.4 → examples/PACKAGES/cgdna/examples/lj_units/oxDNA/duplex2/log.2Jul21.duplex2.g++.4
View File

0
examples/PACKAGES/cgdna/examples/oxDNA2/duplex1/data.duplex1 → examples/PACKAGES/cgdna/examples/lj_units/oxDNA/potential_file/data.duplex1
View File

70
examples/PACKAGES/cgdna/examples/lj_units/oxDNA/potential_file/in.duplex1 Normal file
View File

@ -0,0 +1,70 @@
variable number equal 1
variable ofreq equal 1000
variable efreq equal 1000
variable T equal 0.1
units lj
dimension 3
newton on
boundary p p p
atom_style hybrid bond ellipsoid oxdna
atom_modify sort 0 1.0
# Pair interactions require lists of neighbours to be calculated
neighbor 2.0 bin
neigh_modify every 1 delay 0 check yes
read_data data.duplex1
set atom * mass 3.1575
group all type 1 4
# oxDNA bond interactions - FENE backbone
bond_style oxdna/fene
bond_coeff * oxdna_lj.cgdna
special_bonds lj 0 1 1
# oxDNA pair interactions
pair_style hybrid/overlay oxdna/excv oxdna/stk oxdna/hbond oxdna/xstk oxdna/coaxstk
pair_coeff * * oxdna/excv oxdna_lj.cgdna
pair_coeff * * oxdna/stk seqav 0.1 1.3448 2.6568 oxdna_lj.cgdna
pair_coeff * * oxdna/hbond seqav oxdna_lj.cgdna
pair_coeff 1 4 oxdna/hbond seqav oxdna_lj.cgdna
pair_coeff 2 3 oxdna/hbond seqav oxdna_lj.cgdna
pair_coeff * * oxdna/xstk oxdna_lj.cgdna
pair_coeff * * oxdna/coaxstk oxdna_lj.cgdna
# NVE ensemble
fix 1 all nve/asphere
#fix 2 all langevin ${T} ${T} 2.5 457145 angmom 10
timestep 1e-5
#comm_style tiled
fix 3 all balance 1000 1.03 shift xyz 10 1.03
comm_modify cutoff 3.8
compute quat all property/atom quatw quati quatj quatk
compute erot all erotate/asphere
compute ekin all ke
compute epot all pe
variable erot equal c_erot
variable ekin equal c_ekin
variable epot equal c_epot
variable etot equal c_erot+c_ekin+c_epot
fix 5 all print ${efreq} "$(step) ekin = ${ekin} | erot = ${erot} | epot = ${epot} | etot = ${etot}" screen yes
dump out all custom ${ofreq} out.${number}.lammpstrj id mol type x y z ix iy iz vx vy vz c_quat[1] c_quat[2] c_quat[3] c_quat[4] angmomx angmomy angmomz
dump_modify out sort id
dump_modify out format line "%d %d %d %22.15le %22.15le %22.15le %d %d %d %22.15le %22.15le %22.15le %22.15le %22.15le %22.15le %22.15le %22.15le %22.15le %22.15le"
run 1000000
write_data last_config.${number}.* nocoeff
#write_restart last_config.${number}.*

1113
examples/PACKAGES/cgdna/examples/lj_units/oxDNA/potential_file/log.19May24.duplex1.g++.1 Normal file
View File

File diff suppressed because it is too large Load Diff

1113
examples/PACKAGES/cgdna/examples/lj_units/oxDNA/potential_file/log.19May24.duplex1.g++.4 Normal file
View File

File diff suppressed because it is too large Load Diff

1
examples/PACKAGES/cgdna/examples/lj_units/oxDNA/potential_file/oxdna_lj.cgdna Symbolic link
View File

@ -0,0 +1 @@
../../../../../../../potentials/oxdna_lj.cgdna

0
examples/PACKAGES/cgdna/examples/oxDNA2/dsring/data.dsring → examples/PACKAGES/cgdna/examples/lj_units/oxDNA2/dsring/data.dsring
View File

0
examples/PACKAGES/cgdna/examples/oxDNA2/dsring/in.dsring → examples/PACKAGES/cgdna/examples/lj_units/oxDNA2/dsring/in.dsring
View File

0
examples/PACKAGES/cgdna/examples/oxDNA2/dsring/log.2Jul21.dsring.g++.1 → examples/PACKAGES/cgdna/examples/lj_units/oxDNA2/dsring/log.2Jul21.dsring.g++.1
View File

0
examples/PACKAGES/cgdna/examples/oxDNA2/dsring/log.2Jul21.dsring.g++.4 → examples/PACKAGES/cgdna/examples/lj_units/oxDNA2/dsring/log.2Jul21.dsring.g++.4
View File

68
examples/PACKAGES/cgdna/examples/lj_units/oxDNA2/duplex1/data.duplex1 Normal file
View File

@ -0,0 +1,68 @@
LAMMPS data file via write_data, version 27 May 2021
10 atoms
4 atom types
8 bonds
1 bond types
10 ellipsoids
-20 20 xlo xhi
-20 20 ylo yhi
-20 20 zlo zhi
Masses
1 3.1575
2 3.1575
3 3.1575
4 3.1575
Atoms # hybrid
1 1 -0.33741452300167507 -0.43708835412476305 0.6450685042019271 1 1 3.7269849963023267 0 0 0
2 2 -0.32142606102826937 -0.7137743037592722 1.1817366147004618 1 1 3.7269849963023267 0 0 0
3 3 -0.130363628207774 -0.9147144801536078 1.62581312195109 1 1 3.7269849963023267 0 0 0
4 4 0.16795127962282844 -0.9808507459807022 2.0894908590909003 1 1 3.7269849963023267 0 0 0
5 1 0.46370423490634166 -0.7803347954883079 2.4251986815515827 1 1 3.7269849963023267 0 0 0
6 4 -0.4462950185476711 0.09062163051035639 2.4668941268777607 2 1 3.7269849963023267 0 0 0
7 1 -0.03377054097560965 0.20979847489755046 2.078208732038921 2 1 3.7269849963023267 0 0 0
8 2 0.3297325391466579 0.17657587120899895 1.7206328374934152 2 1 3.7269849963023267 0 0 0
9 3 0.6063699309305985 0.04682595158675571 1.2335049647817748 2 1 3.7269849963023267 0 0 0
10 4 0.8003979559814726 -0.364393011459011 0.9884025318908612 2 1 3.7269849963023267 0 0 0
Velocities
1 0.320321385294804 -0.13632815939410442 -0.029398292452023418 0.3064009492028237 -0.15808560233691588 0.35878007201886397
2 0.16868594667473025 -0.04950805613064363 0.15811007290373785 -0.07666583909321756 -0.0008074676325318194 -0.21475821019816385
3 -0.22924557018300165 0.08381876748892438 -0.0919832851533896 0.4039387481683193 0.6123610642545824 -0.11063432848545783
4 -0.22186204313310393 0.04952817499985707 -0.0693642101605732 -0.1358248430264938 0.4118493572385653 -0.056529305922687775
5 0.08156431270087049 -0.2564594759800144 0.1724544416027875 0.05439894663158808 0.09338481510384318 0.3205408219238883
6 0.03598698404367743 -0.04868642973674152 0.02860105256592922 0.04007709957283992 -0.317943400069374 0.36438025397586354
7 -0.00822868972307372 0.047514932936351305 -0.027726409420297023 0.18356392696891796 -0.49877294396308003 0.06993313839189567
8 -0.07177147672242379 0.1718272727853115 0.39056151182616994 -0.16728362538690794 -0.47839708820957955 -0.17897249005947627
9 -0.1748638855727651 -0.0781638161351808 0.0560181565271157 -0.28062568580131014 0.2405396522734162 -0.4311598357169048
10 0.18870318272756448 -0.1066780134639517 0.12610657946741227 -0.05740397100183697 0.36748833227892685 0.1498775724372025
Bonds
1 1 1 2
2 1 2 3
3 1 3 4
4 1 4 5
5 1 6 7
6 1 7 8
7 1 8 9
8 1 9 10
Ellipsoids
1 1.173984503142341 1.173984503142341 1.173984503142341 0.9890278201757743 0.01779228232037064 -0.14337734159225404 0.030827642240801516
2 1.173984503142341 1.173984503142341 1.173984503142341 0.939687458852748 0.04174166924055095 -0.023337773785056866 0.338674565089608
3 1.173984503142341 1.173984503142341 1.173984503142341 0.8210113150655425 0.03012140921736572 0.017666019956944813 0.5698429897612057
4 1.173984503142341 1.173984503142341 1.173984503142341 0.6623662858285051 -0.028186343967346823 0.022942552517501488 0.7482981175276918
5 1.173984503142341 1.173984503142341 1.173984503142341 0.3601488726765216 0.0513614985821682 0.0724224158335286 0.9286602067807472
6 1.173984503142341 1.173984503142341 1.173984503142341 0.11941234710084649 0.9244660117493703 -0.35317942248051865 -0.07979711784524246
7 1.173984503142341 1.173984503142341 1.173984503142341 -0.17949125421205164 0.7412884899431119 -0.6379094464220707 0.1065166771202199
8 1.173984503142341 1.173984503142341 1.173984503142341 -0.10483691088405202 0.5508895999584645 -0.8250090480220789 0.06992811634525403
9 1.173984503142341 1.173984503142341 1.173984503142341 0.07777239911646 -0.3724087549185288 0.9103052384821374 -0.1631181963720798
10 1.173984503142341 1.173984503142341 1.173984503142341 0.16279109707978262 0.027148630125149613 0.9849325709665359 -0.0516705065113425

0
examples/PACKAGES/cgdna/examples/oxDNA2/duplex1/in.duplex1 → examples/PACKAGES/cgdna/examples/lj_units/oxDNA2/duplex1/in.duplex1
View File

0
examples/PACKAGES/cgdna/examples/oxDNA2/duplex1/log.2Jul21.duplex1.g++.1 → examples/PACKAGES/cgdna/examples/lj_units/oxDNA2/duplex1/log.2Jul21.duplex1.g++.1
View File

0
examples/PACKAGES/cgdna/examples/oxDNA2/duplex1/log.2Jul21.duplex1.g++.4 → examples/PACKAGES/cgdna/examples/lj_units/oxDNA2/duplex1/log.2Jul21.duplex1.g++.4
View File

0
examples/PACKAGES/cgdna/examples/oxDNA2/duplex2/data.duplex2 → examples/PACKAGES/cgdna/examples/lj_units/oxDNA2/duplex2/data.duplex2
View File

0
examples/PACKAGES/cgdna/examples/oxDNA2/duplex2/in.duplex2 → examples/PACKAGES/cgdna/examples/lj_units/oxDNA2/duplex2/in.duplex2
View File

0
examples/PACKAGES/cgdna/examples/oxDNA2/duplex2/log.2Jul21.duplex2.g++.1 → examples/PACKAGES/cgdna/examples/lj_units/oxDNA2/duplex2/log.2Jul21.duplex2.g++.1
View File

0
examples/PACKAGES/cgdna/examples/oxDNA2/duplex2/log.2Jul21.duplex2.g++.4 → examples/PACKAGES/cgdna/examples/lj_units/oxDNA2/duplex2/log.2Jul21.duplex2.g++.4
View File

0
examples/PACKAGES/cgdna/examples/oxDNA2/duplex3/data.duplex3 → examples/PACKAGES/cgdna/examples/lj_units/oxDNA2/duplex3/data.duplex3
View File

0
examples/PACKAGES/cgdna/examples/oxDNA2/duplex3/in.duplex3 → examples/PACKAGES/cgdna/examples/lj_units/oxDNA2/duplex3/in.duplex3
View File

0
examples/PACKAGES/cgdna/examples/oxDNA2/duplex3/log.2Jul21.duplex3.g++.1 → examples/PACKAGES/cgdna/examples/lj_units/oxDNA2/duplex3/log.2Jul21.duplex3.g++.1
View File

0
examples/PACKAGES/cgdna/examples/oxDNA2/duplex3/log.2Jul21.duplex3.g++.4 → examples/PACKAGES/cgdna/examples/lj_units/oxDNA2/duplex3/log.2Jul21.duplex3.g++.4
View File

68
examples/PACKAGES/cgdna/examples/lj_units/oxDNA2/potential_file/data.duplex1 Normal file
View File

@ -0,0 +1,68 @@
LAMMPS data file via write_data, version 27 May 2021
10 atoms
4 atom types
8 bonds
1 bond types
10 ellipsoids
-20 20 xlo xhi
-20 20 ylo yhi
-20 20 zlo zhi
Masses
1 3.1575
2 3.1575
3 3.1575
4 3.1575
Atoms # hybrid
1 1 -0.33741452300167507 -0.43708835412476305 0.6450685042019271 1 1 3.7269849963023267 0 0 0
2 2 -0.32142606102826937 -0.7137743037592722 1.1817366147004618 1 1 3.7269849963023267 0 0 0
3 3 -0.130363628207774 -0.9147144801536078 1.62581312195109 1 1 3.7269849963023267 0 0 0
4 4 0.16795127962282844 -0.9808507459807022 2.0894908590909003 1 1 3.7269849963023267 0 0 0
5 1 0.46370423490634166 -0.7803347954883079 2.4251986815515827 1 1 3.7269849963023267 0 0 0
6 4 -0.4462950185476711 0.09062163051035639 2.4668941268777607 2 1 3.7269849963023267 0 0 0
7 1 -0.03377054097560965 0.20979847489755046 2.078208732038921 2 1 3.7269849963023267 0 0 0
8 2 0.3297325391466579 0.17657587120899895 1.7206328374934152 2 1 3.7269849963023267 0 0 0
9 3 0.6063699309305985 0.04682595158675571 1.2335049647817748 2 1 3.7269849963023267 0 0 0
10 4 0.8003979559814726 -0.364393011459011 0.9884025318908612 2 1 3.7269849963023267 0 0 0
Velocities
1 0.320321385294804 -0.13632815939410442 -0.029398292452023418 0.3064009492028237 -0.15808560233691588 0.35878007201886397
2 0.16868594667473025 -0.04950805613064363 0.15811007290373785 -0.07666583909321756 -0.0008074676325318194 -0.21475821019816385
3 -0.22924557018300165 0.08381876748892438 -0.0919832851533896 0.4039387481683193 0.6123610642545824 -0.11063432848545783
4 -0.22186204313310393 0.04952817499985707 -0.0693642101605732 -0.1358248430264938 0.4118493572385653 -0.056529305922687775
5 0.08156431270087049 -0.2564594759800144 0.1724544416027875 0.05439894663158808 0.09338481510384318 0.3205408219238883
6 0.03598698404367743 -0.04868642973674152 0.02860105256592922 0.04007709957283992 -0.317943400069374 0.36438025397586354
7 -0.00822868972307372 0.047514932936351305 -0.027726409420297023 0.18356392696891796 -0.49877294396308003 0.06993313839189567
8 -0.07177147672242379 0.1718272727853115 0.39056151182616994 -0.16728362538690794 -0.47839708820957955 -0.17897249005947627
9 -0.1748638855727651 -0.0781638161351808 0.0560181565271157 -0.28062568580131014 0.2405396522734162 -0.4311598357169048
10 0.18870318272756448 -0.1066780134639517 0.12610657946741227 -0.05740397100183697 0.36748833227892685 0.1498775724372025
Bonds
1 1 1 2
2 1 2 3
3 1 3 4
4 1 4 5
5 1 6 7
6 1 7 8
7 1 8 9
8 1 9 10
Ellipsoids
1 1.173984503142341 1.173984503142341 1.173984503142341 0.9890278201757743 0.01779228232037064 -0.14337734159225404 0.030827642240801516
2 1.173984503142341 1.173984503142341 1.173984503142341 0.939687458852748 0.04174166924055095 -0.023337773785056866 0.338674565089608
3 1.173984503142341 1.173984503142341 1.173984503142341 0.8210113150655425 0.03012140921736572 0.017666019956944813 0.5698429897612057
4 1.173984503142341 1.173984503142341 1.173984503142341 0.6623662858285051 -0.028186343967346823 0.022942552517501488 0.7482981175276918
5 1.173984503142341 1.173984503142341 1.173984503142341 0.3601488726765216 0.0513614985821682 0.0724224158335286 0.9286602067807472
6 1.173984503142341 1.173984503142341 1.173984503142341 0.11941234710084649 0.9244660117493703 -0.35317942248051865 -0.07979711784524246
7 1.173984503142341 1.173984503142341 1.173984503142341 -0.17949125421205164 0.7412884899431119 -0.6379094464220707 0.1065166771202199
8 1.173984503142341 1.173984503142341 1.173984503142341 -0.10483691088405202 0.5508895999584645 -0.8250090480220789 0.06992811634525403
9 1.173984503142341 1.173984503142341 1.173984503142341 0.07777239911646 -0.3724087549185288 0.9103052384821374 -0.1631181963720798
10 1.173984503142341 1.173984503142341 1.173984503142341 0.16279109707978262 0.027148630125149613 0.9849325709665359 -0.0516705065113425

72
examples/PACKAGES/cgdna/examples/lj_units/oxDNA2/potential_file/in.duplex1 Normal file
View File

@ -0,0 +1,72 @@
variable number equal 1
variable ofreq equal 1000
variable efreq equal 1000
variable T equal 0.1
variable rhos equal 0.2
units lj
dimension 3
newton on
boundary p p p
atom_style hybrid bond ellipsoid oxdna
atom_modify sort 0 1.0
# Pair interactions require lists of neighbours to be calculated
neighbor 2.0 bin
neigh_modify every 1 delay 0 check yes
read_data data.duplex1
set atom * mass 3.1575
group all type 1 4
# oxDNA2 bond interactions - FENE backbone
bond_style oxdna2/fene
bond_coeff * oxdna2_lj.cgdna
special_bonds lj 0 1 1
# oxDNA2 pair interactions
pair_style hybrid/overlay oxdna2/excv oxdna2/stk oxdna2/hbond oxdna2/xstk oxdna2/coaxstk oxdna2/dh
pair_coeff * * oxdna2/excv oxdna2_lj.cgdna
pair_coeff * * oxdna2/stk seqdep 0.1 1.3523 2.6717 oxdna2_lj.cgdna
pair_coeff * * oxdna2/hbond seqdep oxdna2_lj.cgdna
pair_coeff 1 4 oxdna2/hbond seqdep oxdna2_lj.cgdna
pair_coeff 2 3 oxdna2/hbond seqdep oxdna2_lj.cgdna
pair_coeff * * oxdna2/xstk oxdna2_lj.cgdna
pair_coeff * * oxdna2/coaxstk oxdna2_lj.cgdna
pair_coeff * * oxdna2/dh 0.1 0.5 oxdna2_lj.cgdna
# NVE ensemble
fix 1 all nve/asphere
#fix 2 all langevin ${T} ${T} 2.5 457145 angmom 10
timestep 1e-5
#comm_style tiled
fix 3 all balance 1000 1.03 shift xyz 10 1.03
comm_modify cutoff 3.8
compute quat all property/atom quatw quati quatj quatk
compute erot all erotate/asphere
compute ekin all ke
compute epot all pe
variable erot equal c_erot
variable ekin equal c_ekin
variable epot equal c_epot
variable etot equal c_erot+c_ekin+c_epot
fix 5 all print ${efreq} "$(step) ekin = ${ekin} | erot = ${erot} | epot = ${epot} | etot = ${etot}" screen yes
dump out all custom ${ofreq} out.${number}.lammpstrj id mol type x y z ix iy iz vx vy vz c_quat[1] c_quat[2] c_quat[3] c_quat[4] angmomx angmomy angmomz
dump_modify out sort id
dump_modify out format line "%d %d %d %22.15le %22.15le %22.15le %d %d %d %22.15le %22.15le %22.15le %22.15le %22.15le %22.15le %22.15le %22.15le %22.15le %22.15le"
run 1000000
write_data last_config.${number}.* nocoeff
#write_restart last_config.${number}.*

Some files were not shown because too many files have changed in this diff Show More