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

Merge pull request #2508 from athomps/mliappy3

Browse Source 
Add support for Python based models to the MLIAP package
This commit is contained in:
Axel Kohlmeyer
2020-12-29 09:59:43 -05:00
committed by GitHub
parent ab5a1f229e 5dc868ec30
commit 427d037747
47 changed files with 1630 additions and 87 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
cmake/CMakeLists.txt
View File

@ -373,7 +373,7 @@ else()
set(CUDA_REQUEST_PIC) set(CUDA_REQUEST_PIC)
endif() endif()
foreach(PKG_WITH_INCL KSPACE PYTHON VORONOI USER-COLVARS USER-MOLFILE USER-NETCDF USER-PLUMED USER-QMMM foreach(PKG_WITH_INCL KSPACE PYTHON MLIAP VORONOI USER-COLVARS USER-MOLFILE USER-NETCDF USER-PLUMED USER-QMMM
USER-QUIP USER-SCAFACOS USER-SMD USER-VTK KIM LATTE MESSAGE MSCG COMPRESS) USER-QUIP USER-SCAFACOS USER-SMD USER-VTK KIM LATTE MESSAGE MSCG COMPRESS)
if(PKG_${PKG_WITH_INCL}) if(PKG_${PKG_WITH_INCL})
include(Packages/${PKG_WITH_INCL}) include(Packages/${PKG_WITH_INCL})

30
cmake/Modules/FindCythonize.cmake Normal file
View File

@ -0,0 +1,30 @@
# Find the Cythonize tool.
#
# This code sets the following variables:
#
# Cythonize_EXECUTABLE
#
# adapted from https://github.com/cmarshall108/cython-cmake-example/blob/master/cmake/FindCython.cmake
#=============================================================================
if(CMAKE_VERSION VERSION_LESS 3.12)
find_package(PythonInterp 3.6 QUIET) # Deprecated since version 3.12
if(PYTHONINTERP_FOUND)
set(Python3_EXECUTABLE ${PYTHON_EXECUTABLE})
endif()
else()
find_package(Python3 3.6 COMPONENTS Interpreter QUIET)
endif()
# Use the Cython executable that lives next to the Python executable
# if it is a local installation.
if(Python3_EXECUTABLE)
get_filename_component(_python_path ${Python3_EXECUTABLE} PATH)
find_program(Cythonize_EXECUTABLE
NAMES cythonize3 cythonize cythonize.bat
HINTS ${_python_path})
endif()
include(FindPackageHandleStandardArgs)
FIND_PACKAGE_HANDLE_STANDARD_ARGS(Cythonize REQUIRED_VARS Cythonize_EXECUTABLE)
mark_as_advanced(Cythonize_EXECUTABLE)

1
cmake/Modules/LAMMPSUtils.cmake
View File

@ -50,6 +50,7 @@ function(check_for_autogen_files source_dir)
file(GLOB SRC_AUTOGEN_FILES ${source_dir}/style_*.h) file(GLOB SRC_AUTOGEN_FILES ${source_dir}/style_*.h)
file(GLOB SRC_AUTOGEN_PACKAGES ${source_dir}/packages_*.h) file(GLOB SRC_AUTOGEN_PACKAGES ${source_dir}/packages_*.h)
list(APPEND SRC_AUTOGEN_FILES ${SRC_AUTOGEN_PACKAGES} ${source_dir}/lmpinstalledpkgs.h ${source_dir}/lmpgitversion.h) list(APPEND SRC_AUTOGEN_FILES ${SRC_AUTOGEN_PACKAGES} ${source_dir}/lmpinstalledpkgs.h ${source_dir}/lmpgitversion.h)
list(APPEND SRC_AUTOGEN_FILES ${SRC_AUTOGEN_PACKAGES} ${source_dir}/mliap_model_python_couple.h ${source_dir}/mliap_model_python_couple.cpp)
foreach(_SRC ${SRC_AUTOGEN_FILES}) foreach(_SRC ${SRC_AUTOGEN_FILES})
get_filename_component(FILENAME "${_SRC}" NAME) get_filename_component(FILENAME "${_SRC}" NAME)
if(EXISTS ${source_dir}/${FILENAME}) if(EXISTS ${source_dir}/${FILENAME})

31
cmake/Modules/Packages/MLIAP.cmake Normal file
View File

@ -0,0 +1,31 @@
# if PYTHON package is included we may also include Python support in MLIAP
set(MLIAP_ENABLE_PYTHON_DEFAULT OFF)
if(PKG_PYTHON)
find_package(Cythonize)
if(Cythonize_FOUND)
set(MLIAP_ENABLE_PYTHON_DEFAULT ON)
endif()
endif()
option(MLIAP_ENABLE_PYTHON "Build MLIAP package with Python support" ${MLIAP_ENABLE_PYTHON_DEFAULT})
if(MLIAP_ENABLE_PYTHON)
find_package(Cythonize REQUIRED)
if(NOT PKG_PYTHON)
message(FATAL_ERROR "Must enable PYTHON package for including Python support in MLIAP")
endif()
set(MLIAP_BINARY_DIR ${CMAKE_BINARY_DIR}/cython)
set(MLIAP_CYTHON_SRC ${LAMMPS_SOURCE_DIR}/MLIAP/mliap_model_python_couple.pyx)
get_filename_component(MLIAP_CYTHON_BASE ${MLIAP_CYTHON_SRC} NAME_WE)
file(MAKE_DIRECTORY ${MLIAP_BINARY_DIR})
add_custom_command(OUTPUT ${MLIAP_BINARY_DIR}/${MLIAP_CYTHON_BASE}.cpp ${MLIAP_BINARY_DIR}/${MLIAP_CYTHON_BASE}.h
COMMAND ${CMAKE_COMMAND} -E copy_if_different ${MLIAP_CYTHON_SRC} ${MLIAP_BINARY_DIR}/${MLIAP_CYTHON_BASE}.pyx
COMMAND ${Cythonize_EXECUTABLE} -3 ${MLIAP_BINARY_DIR}/${MLIAP_CYTHON_BASE}.pyx
WORKING_DIRECTORY ${MLIAP_BINARY_DIR}
MAIN_DEPENDENCY ${MLIAP_CYTHON_SRC}
COMMENT "Generating C++ sources with cythonize...")
target_compile_definitions(lammps PRIVATE -DMLIAP_PYTHON)
target_sources(lammps PRIVATE ${MLIAP_BINARY_DIR}/${MLIAP_CYTHON_BASE}.cpp)
target_include_directories(lammps PRIVATE ${MLIAP_BINARY_DIR})
endif()

49
doc/src/Build_extras.rst
View File

@ -37,6 +37,7 @@ This is the list of packages that may require additional steps.
* :ref:`KOKKOS <kokkos>` * :ref:`KOKKOS <kokkos>`
* :ref:`LATTE <latte>` * :ref:`LATTE <latte>`
* :ref:`MESSAGE <message>` * :ref:`MESSAGE <message>`
* :ref:`MLIAP <mliap>`
* :ref:`MSCG <mscg>` * :ref:`MSCG <mscg>`
* :ref:`OPT <opt>` * :ref:`OPT <opt>`
* :ref:`POEMS <poems>` * :ref:`POEMS <poems>`
@ -770,6 +771,54 @@ be installed on your system.
---------- ----------
.. _mliap:
MLIAP package
---------------------------
Building the MLIAP package requires including the :ref:`SNAP <PKG-SNAP>`
package. There will be an error message if this requirement is not satisfied.
Using the *mliappy* model also requires enabling Python support, which
in turn requires the :ref:`PYTHON <PKG-PYTHON>`
package **and** requires you have the `cython <https://cython.org>`_ software
installed and with it a working ``cythonize`` command. This feature requires
compiling LAMMPS with Python version 3.6 or later.
.. tabs::
.. tab:: CMake build
.. code-block:: bash
-D MLIAP_ENABLE_PYTHON=value # enable mliappy model (default is autodetect)
Without this setting, CMake will check whether it can find a
suitable Python version and the ``cythonize`` command and choose
the default accordingly. During the build procedure the provided
.pyx file(s) will be automatically translated to C++ code and compiled.
Please do **not** run ``cythonize`` manually in the ``src/MLIAP`` folder,
as that can lead to compilation errors if Python support is not enabled.
If you did by accident, please remove the generated .cpp and .h files.
.. tab:: Traditional make
The build uses the ``lib/python/Makefile.mliap_python`` file in the
compile/link process to add a rule to update the files generated by
the ``cythonize`` command in case the corresponding .pyx file(s) were
modified. You may need to modify ``lib/python/Makefile.lammps``
if the LAMMPS build fails.
To manually enforce building MLIAP with Python support enabled,
you can add
``-DMLIAP_PYTHON`` to the ``LMP_INC`` variable in your machine makefile.
You may have to manually run the ``cythonize`` command on .pyx file(s)
in the ``src`` folder, if this is not automatically done during
installing the MLIAP package. Please do **not** run ``cythonize``
in the ``src/MLIAP`` folder, as that can lead to compilation errors
if Python support is not enabled.
If you did by accident, please remove the generated .cpp and .h files.
----------
.. _mscg: .. _mscg:
MSCG package MSCG package

20
doc/src/Packages_details.rst
View File

@ -662,19 +662,31 @@ MLIAP package
**Contents:** **Contents:**
A general interface for machine-learning interatomic potentials. A general interface for machine-learning interatomic potentials, including PyTorch.
**Install:** **Install:**
To use this package, also the :ref:`SNAP package <PKG-SNAP>` needs to be installed. To use this package, also the :ref:`SNAP package <PKG-SNAP>` package needs
to be installed. To make the *mliappy* model available, also the
:ref:`PYTHON package <PKG-PYTHON>` package needs to be installed, the version of
Python must be 3.6 or later, and the `cython <https://cython.org/>`_ software
must be installed.
**Author:** Aidan Thompson (Sandia). **Author:** Aidan Thompson (Sandia), Nicholas Lubbers (LANL).
**Supporting info:** **Supporting info:**
* src/MLIAP: filenames -> commands * src/MLIAP: filenames -> commands
* src/MLIAP/README
* :doc:`pair_style mliap <pair_mliap>` * :doc:`pair_style mliap <pair_mliap>`
* examples/mliap * :doc:`compute_style mliap <compute_mliap>`
* examples/mliap (see README)
When built with the *mliappy* model this package includes an extension for
coupling with Python models, including PyTorch. In this case, the Python
interpreter linked to LAMMPS will need the ``cython`` and ``numpy`` modules
installed. The provided examples build models with PyTorch, which would
therefore also needs to be installed to run those examples.
---------- ----------

19
doc/src/compute_mliap.rst
View File

@ -18,7 +18,7 @@ Syntax
.. parsed-literal:: .. parsed-literal::
*model* values = style *model* values = style
style = *linear* or *quadratic* style = *linear* or *quadratic* or *mliappy*
*descriptor* values = style filename *descriptor* values = style filename
style = *sna* style = *sna*
filename = name of file containing descriptor definitions filename = name of file containing descriptor definitions
@ -56,13 +56,15 @@ and it is also straightforward to add new descriptor styles.
The compute *mliap* command must be followed by two keywords The compute *mliap* command must be followed by two keywords
*model* and *descriptor* in either order. *model* and *descriptor* in either order.
The *model* keyword is followed by a model style, currently limited to The *model* keyword is followed by the model style (*linear*, *quadratic* or *mliappy*).
either *linear* or *quadratic*. The *mliappy* model is only available
if lammps is built with MLIAPPY package.
The *descriptor* keyword is followed by a descriptor style, and additional arguments. The *descriptor* keyword is followed by a descriptor style, and additional arguments.
Currently the only descriptor style is *sna*, indicating the bispectrum component The compute currently supports just one descriptor style, but it is
descriptors used by the Spectral Neighbor Analysis Potential (SNAP) potentials of is straightforward to add new descriptor styles.
:doc:`pair_style snap <pair_snap>`. The SNAP descriptor style *sna* is the same as that used by :doc:`pair_style snap <pair_snap>`,
including the linear, quadratic, and chem variants.
A single additional argument specifies the descriptor filename A single additional argument specifies the descriptor filename
containing the parameters and setting used by the SNAP descriptor. containing the parameters and setting used by the SNAP descriptor.
The descriptor filename usually ends in the *.mliap.descriptor* extension. The descriptor filename usually ends in the *.mliap.descriptor* extension.
@ -162,9 +164,10 @@ potentials, see the examples in `FitSNAP <https://github.com/FitSNAP/FitSNAP>`_.
Restrictions Restrictions
"""""""""""" """"""""""""
This compute is part of the MLIAP package. It is only enabled if This compute is part of the MLIAP package. It is only enabled if LAMMPS
LAMMPS was built with that package. In addition, building LAMMPS with the MLIAP package was built with that package. In addition, building LAMMPS with the MLIAP package
requires building LAMMPS with the SNAP package. requires building LAMMPS with the SNAP package.
The *mliappy* model requires building LAMMPS with the PYTHON package.
See the :doc:`Build package <Build_package>` doc page for more info. See the :doc:`Build package <Build_package>` doc page for more info.
Related commands Related commands

36
doc/src/pair_mliap.rst
View File

@ -16,7 +16,7 @@ Syntax
.. parsed-literal:: .. parsed-literal::
*model* values = style filename *model* values = style filename
style = *linear* or *quadratic* style = *linear* or *quadratic* or *mliappy*
filename = name of file containing model definitions filename = name of file containing model definitions
*descriptor* values = style filename *descriptor* values = style filename
style = *sna* style = *sna*
@ -40,12 +40,15 @@ definitions of the interatomic potential functional form (*model*)
and the geometric quantities that characterize the atomic positions and the geometric quantities that characterize the atomic positions
(*descriptor*). By defining *model* and *descriptor* separately, (*descriptor*). By defining *model* and *descriptor* separately,
it is possible to use many different models with a given descriptor, it is possible to use many different models with a given descriptor,
or many different descriptors with a given model. Currently, the pair_style or many different descriptors with a given model. The
supports just two models, *linear* and *quadratic*, pair style currently supports just one descriptor style, but it is
and one descriptor, *sna*, the SNAP descriptor used by :doc:`pair_style snap <pair_snap>`, including the linear, quadratic, is straightforward to add new descriptor styles.
and chem variants. Work is currently underway to extend The SNAP descriptor style *sna* is the same as that used by :doc:`pair_style snap <pair_snap>`,
the interface to handle neural network energy models, including the linear, quadratic, and chem variants.
and it is also straightforward to add new descriptor styles. The available models are *linear*, *quadratic*, and *mliappy*.
The *mliappy* style can be used to couple python models,
e.g. PyTorch neural network energy models, and requires building
LAMMPS with the PYTHON package (see below).
In order to train a model, it is useful to know the gradient or derivative In order to train a model, it is useful to know the gradient or derivative
of energy, force, and stress w.r.t. model parameters. This information of energy, force, and stress w.r.t. model parameters. This information
can be accessed using the related :doc:`compute mliap <compute_mliap>` command. can be accessed using the related :doc:`compute mliap <compute_mliap>` command.
@ -59,9 +62,8 @@ that specify the mapping of MLIAP
element names to LAMMPS atom types, element names to LAMMPS atom types,
where N is the number of LAMMPS atom types. where N is the number of LAMMPS atom types.
The *model* keyword is followed by a model style, currently limited to The *model* keyword is followed by the model style. This is followed
either *linear* or *quadratic*. In both cases, by a single argument specifying the model filename containing the
this is followed by a single argument specifying the model filename containing the
parameters for a set of elements. parameters for a set of elements.
The model filename usually ends in the *.mliap.model* extension. The model filename usually ends in the *.mliap.model* extension.
It may contain parameters for many elements. The only requirement is that it It may contain parameters for many elements. The only requirement is that it
@ -82,6 +84,16 @@ for the :doc:`pair_style snap <pair_snap>` coefficient file.
Specifically, the line containing the element weight and radius is omitted, Specifically, the line containing the element weight and radius is omitted,
since these are handled by the *descriptor*. since these are handled by the *descriptor*.
Notes on mliappy models:
When the *model* keyword is *mliappy*, the filename should end in '.pt',
'.pth' for pytorch models, or be a pickle file. To load a model from
memory (i.e. an existing python object), specify the filename as
"LATER", and then call `lammps.mliap.load_model(model)` from python
before using the pair style. When using lammps via the library mode, you will need to call
`lammps.mliappy.activate_mliappy(lmp)` on the active lammps object
before the pair style is defined. This call locates and loads the mliap-specific
python module that is built into lammps.
The *descriptor* keyword is followed by a descriptor style, and additional arguments. The *descriptor* keyword is followed by a descriptor style, and additional arguments.
Currently the only descriptor style is *sna*, indicating the bispectrum component Currently the only descriptor style is *sna*, indicating the bispectrum component
descriptors used by the Spectral Neighbor Analysis Potential (SNAP) potentials of descriptors used by the Spectral Neighbor Analysis Potential (SNAP) potentials of
@ -138,11 +150,13 @@ This pair style can only be used via the *pair* keyword of the
Restrictions Restrictions
"""""""""""" """"""""""""
This style is part of the MLIAP package. It is only enabled if LAMMPS This pair style is part of the MLIAP package. It is only enabled if LAMMPS
was built with that package. In addition, building LAMMPS with the MLIAP package was built with that package. In addition, building LAMMPS with the MLIAP package
requires building LAMMPS with the SNAP package. requires building LAMMPS with the SNAP package.
The *mliappy* model requires building LAMMPS with the PYTHON package.
See the :doc:`Build package <Build_package>` doc page for more info. See the :doc:`Build package <Build_package>` doc page for more info.
Related commands Related commands
"""""""""""""""" """"""""""""""""

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

@ -558,6 +558,7 @@ Cygwin
cylindrically cylindrically
Cyrot Cyrot
cyrstals cyrstals
cython
Daivis Daivis
Dammak Dammak
dampflag dampflag

103
examples/mliap/README Normal file
View File

@ -0,0 +1,103 @@
This directory contains multiple examples of
machine-learning potentials defined using the
MLIAP package in LAMMPS. The input files
are described below.
in.mliap.snap.Ta06A
-------------------
Run linear SNAP, equivalent to examples/snap/in.snap.Ta06A
in.mliap.snap.WBe.PRB2019
-------------------------
Run linear SNAP, equivalent to examples/snap/in.snap.WBe.PRB2019
in.mliap.snap.quadratic
-----------------------
Run quadratic SNAP
in.mliap.snap.chem
------------------
Run EME-SNAP, equivalent to examples/snap/in.snap.InP.JCPA2020
in.mliap.snap.compute
---------------------
Generate the A matrix, the gradients (w.r.t. coefficients)
of total potential energy, forces, and stress tensor for
linear SNAP, equivalent to in.snap.compute
in.mliap.quadratic.compute
--------------------------
Generate the A matrix, the gradients (w.r.t. coefficients)
of total potential energy, forces, and stress tensor for
for quadratic SNAP, equivalent to in.snap.compute.quadratic
in.mliap.pytorch.Ta06A
-----------------------
This reproduces the output of in.mliap.snap.Ta06A above,
but using the Python coupling to PyTorch.
This example can be run in two different ways:
1: Running a LAMMPS executable: in.mliap.pytorch.Ta06A
First run ``python convert_mliap_Ta06A.py``. It creates
a PyTorch energy model that replicates the
SNAP Ta06A potential and saves it in the file
"Ta06A.mliap.pytorch.model.pt".
You can then run the example as follows
`lmp -in in.mliap.pytorch.Ta06A -echo both`
The resultant log.lammps output should be identical to that generated
by in.mliap.snap.Ta06A.
If this fails, see the instructions for building the MLIAP package
with Python support enabled. Also, confirm that the
LAMMPS Python embedded Python interpreter is
working by running ../examples/in.python.
2: Running a Python script: mliap_pytorch_Ta06A.py
Before testing this, ensure that the previous method
(running a LAMMPS executable) works.
You can run the example in serial:
`python mliap_pytorch_Ta06A.py`
or in parallel:
`mpirun -np 4 python mliap_pytorch_Ta06A.py`
The resultant log.lammps output should be identical to that generated
by in.mliap.snap.Ta06A and in.mliap.pytorch.Ta06A.
Not all Python installations support this mode of operation.
It requires that the Python interpreter be initialized. If not,
the script will exit with an error message.
in.mliap.pytorch.relu1hidden
----------------------------
This example demonstrates a simple neural network potential
using PyTorch and SNAP descriptors.
`lmp -in in.mliap.pytorch.relu1hidden -echo both`
It was trained on just the energy component (no forces) of
the data used in the original SNAP Ta06A potential for
tantalum (Thompson, Swiler, Trott, Foiles, Tucker,
J Comp Phys, 285, 316 (2015).). Because of the very small amount
of energy training data, it uses just 1 hidden layer with
a ReLU activation function. It is not expected to be
very accurate for forces.
NOTE: Unlike the previous example, this example uses
a pre-built PyTorch file `Ta06A.mliap.pytorch.model.pt`.
It is read using `torch.load`,
which implicitly uses the Python `pickle` module.
This is known to be insecure. It is possible to construct malicious
pickle data that will execute arbitrary code during unpickling. Never
load data that could have come from an untrusted source, or that
could have been tampered with. Only load data you trust.

18
examples/mliap/Ta06A.mliap.pytorch Normal file
View File

@ -0,0 +1,18 @@
# DATE: 2014-09-05 UNITS: metal CONTRIBUTOR: Aidan Thompson athomps@sandia.gov CITATION: Thompson, Swiler, Trott, Foiles and Tucker, arxiv.org, 1409.3880 (2014)
# Definition of SNAP potential Ta_Cand06A
# Assumes 1 LAMMPS atom type
variable zblcutinner equal 4
variable zblcutouter equal 4.8
variable zblz equal 73
# Specify hybrid with SNAP, ZBL
pair_style hybrid/overlay &
zbl ${zblcutinner} ${zblcutouter} &
mliap model mliappy Ta06A.mliap.pytorch.model.pt &
descriptor sna Ta06A.mliap.descriptor
pair_coeff 1 1 zbl ${zblz} ${zblz}
pair_coeff * * mliap Ta

26
examples/mliap/convert_mliap_Ta06A.py Normal file
View File

@ -0,0 +1,26 @@
import sys
import numpy as np
import torch
# torch.nn.modules useful for defining a MLIAPPY model.
from lammps.mliap.pytorch import TorchWrapper, IgnoreElems
# Read coefficients
coeffs = np.genfromtxt("Ta06A.mliap.model",skip_header=6)
# Write coefficients to a pytorch linear model
bias = coeffs[0]
weights = coeffs[1:]
lin = torch.nn.Linear(weights.shape[0],1)
lin.to(torch.float64)
with torch.autograd.no_grad():
lin.weight.set_(torch.from_numpy(weights).unsqueeze(0))
lin.bias.set_(torch.as_tensor(bias,dtype=torch.float64).unsqueeze(0))
# Wrap the pytorch model for usage with mliappy coupling.
model = IgnoreElems(lin) # The linear module does not use the types.
n_descriptors = lin.weight.shape[1]
n_elements = 1
linked_model = TorchWrapper(model,n_descriptors=n_descriptors,n_elements=n_elements)
torch.save(linked_model,"Ta06A.mliap.pytorch.model.pt")

53
examples/mliap/in.mliap.pytorch.Ta06A Normal file
View File

@ -0,0 +1,53 @@
# Demonstrate MLIAP/PyTorch interface to linear SNAP potential
# Initialize simulation
variable nsteps index 100
variable nrep equal 4
variable a equal 3.316
units metal
# generate the box and atom positions using a BCC lattice
variable nx equal ${nrep}
variable ny equal ${nrep}
variable nz equal ${nrep}
boundary p p p
lattice bcc $a
region box block 0 ${nx} 0 ${ny} 0 ${nz}
create_box 1 box
create_atoms 1 box
mass 1 180.88
# choose potential
include Ta06A.mliap.pytorch
# Setup output
compute eatom all pe/atom
compute energy all reduce sum c_eatom
compute satom all stress/atom NULL
compute str all reduce sum c_satom[1] c_satom[2] c_satom[3]
variable press equal (c_str[1]+c_str[2]+c_str[3])/(3*vol)
thermo_style custom step temp epair c_energy etotal press v_press
thermo 10
thermo_modify norm yes
# Set up NVE run
timestep 0.5e-3
neighbor 1.0 bin
neigh_modify once no every 1 delay 0 check yes
# Run MD
velocity all create 300.0 4928459 loop geom
fix 1 all nve
run ${nsteps}

53
examples/mliap/in.mliap.pytorch.relu1hidden Normal file
View File

@ -0,0 +1,53 @@
# Demonstrate MLIAP interface to linear SNAP potential
# Initialize simulation
variable nsteps index 100
variable nrep equal 4
variable a equal 3.316
units metal
# generate the box and atom positions using a BCC lattice
variable nx equal ${nrep}
variable ny equal ${nrep}
variable nz equal ${nrep}
boundary p p p
lattice bcc $a
region box block 0 ${nx} 0 ${ny} 0 ${nz}
create_box 1 box
create_atoms 1 box
mass 1 180.88
# choose potential
include relu1hidden.mliap.pytorch
# Setup output
compute eatom all pe/atom
compute energy all reduce sum c_eatom
compute satom all stress/atom NULL
compute str all reduce sum c_satom[1] c_satom[2] c_satom[3]
variable press equal (c_str[1]+c_str[2]+c_str[3])/(3*vol)
thermo_style custom step temp epair c_energy etotal press v_press
thermo 10
thermo_modify norm yes
# Set up NVE run
timestep 0.5e-3
neighbor 1.0 bin
neigh_modify once no every 1 delay 0 check yes
# Run MD
velocity all create 300.0 4928459 loop geom
fix 1 all nve
run ${nsteps}

2
examples/mliap/in.mliap.snap.Ta06A
View File

@ -1,4 +1,4 @@
# Demonstrate MLIAP interface to kinear SNAP potential # Demonstrate MLIAP interface to linear SNAP potential
# Initialize simulation # Initialize simulation

157
examples/mliap/log.04Dec20.mliap.pytorch.Ta06A.g++.1 Normal file
View File

@ -0,0 +1,157 @@
LAMMPS (30 Nov 2020)
using 48 OpenMP thread(s) per MPI task
# Demonstrate MLIAP/PyTorch interface to linear SNAP potential
# Initialize simulation
variable nsteps index 100
variable nrep equal 4
variable a equal 3.316
units metal
# generate the box and atom positions using a BCC lattice
variable nx equal ${nrep}
variable nx equal 4
variable ny equal ${nrep}
variable ny equal 4
variable nz equal ${nrep}
variable nz equal 4
boundary p p p
lattice bcc $a
lattice bcc 3.316
Lattice spacing in x,y,z = 3.3160000 3.3160000 3.3160000
region box block 0 ${nx} 0 ${ny} 0 ${nz}
region box block 0 4 0 ${ny} 0 ${nz}
region box block 0 4 0 4 0 ${nz}
region box block 0 4 0 4 0 4
create_box 1 box
Created orthogonal box = (0.0000000 0.0000000 0.0000000) to (13.264000 13.264000 13.264000)
1 by 1 by 1 MPI processor grid
create_atoms 1 box
Created 128 atoms
create_atoms CPU = 0.002 seconds
mass 1 180.88
# choose potential
include Ta06A.mliap.pytorch
# DATE: 2014-09-05 UNITS: metal CONTRIBUTOR: Aidan Thompson athomps@sandia.gov CITATION: Thompson, Swiler, Trott, Foiles and Tucker, arxiv.org, 1409.3880 (2014)
# Definition of SNAP potential Ta_Cand06A
# Assumes 1 LAMMPS atom type
variable zblcutinner equal 4
variable zblcutouter equal 4.8
variable zblz equal 73
# Specify hybrid with SNAP, ZBL
pair_style hybrid/overlay zbl ${zblcutinner} ${zblcutouter} mliap model mliappy Ta06A.mliap.pytorch.model.pkl descriptor sna Ta06A.mliap.descriptor
pair_style hybrid/overlay zbl 4 ${zblcutouter} mliap model mliappy Ta06A.mliap.pytorch.model.pkl descriptor sna Ta06A.mliap.descriptor
pair_style hybrid/overlay zbl 4 4.8 mliap model mliappy Ta06A.mliap.pytorch.model.pkl descriptor sna Ta06A.mliap.descriptor
Loading python model complete.
Reading potential file Ta06A.mliap.descriptor with DATE: 2014-09-05
SNAP keyword rcutfac 4.67637
SNAP keyword twojmax 6
SNAP keyword nelems 1
SNAP keyword elems Ta
SNAP keyword radelems 0.5
SNAP keyword welems 1
SNAP keyword rfac0 0.99363
SNAP keyword rmin0 0
SNAP keyword bzeroflag 0
pair_coeff 1 1 zbl ${zblz} ${zblz}
pair_coeff 1 1 zbl 73 ${zblz}
pair_coeff 1 1 zbl 73 73
pair_coeff * * mliap Ta
# Setup output
compute eatom all pe/atom
compute energy all reduce sum c_eatom
compute satom all stress/atom NULL
compute str all reduce sum c_satom[1] c_satom[2] c_satom[3]
variable press equal (c_str[1]+c_str[2]+c_str[3])/(3*vol)
thermo_style custom step temp epair c_energy etotal press v_press
thermo 10
thermo_modify norm yes
# Set up NVE run
timestep 0.5e-3
neighbor 1.0 bin
neigh_modify once no every 1 delay 0 check yes
# Run MD
velocity all create 300.0 4928459 loop geom
fix 1 all nve
run ${nsteps}
run 100
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 5.8
ghost atom cutoff = 5.8
binsize = 2.9, bins = 5 5 5
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair zbl, perpetual, half/full from (2)
attributes: half, newton on
pair build: halffull/newton
stencil: none
bin: none
(2) pair mliap, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 159.8 | 159.8 | 159.8 Mbytes
Step Temp E_pair c_energy TotEng Press v_press
0 300 -11.85157 -11.85157 -11.813095 2717.1661 -2717.1661
10 296.01467 -11.851059 -11.851059 -11.813095 2697.4796 -2697.4796
20 284.53666 -11.849587 -11.849587 -11.813095 2289.1527 -2289.1527
30 266.51577 -11.847275 -11.847275 -11.813095 1851.7131 -1851.7131
40 243.05007 -11.844266 -11.844266 -11.813095 1570.684 -1570.684
50 215.51032 -11.840734 -11.840734 -11.813094 1468.1899 -1468.1899
60 185.48331 -11.836883 -11.836883 -11.813094 1524.8757 -1524.8757
70 154.6736 -11.832931 -11.832931 -11.813094 1698.3351 -1698.3351
80 124.79303 -11.829099 -11.829099 -11.813094 1947.0715 -1947.0715
90 97.448054 -11.825592 -11.825592 -11.813094 2231.9563 -2231.9563
100 74.035418 -11.822589 -11.822589 -11.813094 2515.8526 -2515.8526
Loop time of 2.00236 on 48 procs for 100 steps with 128 atoms
Performance: 2.157 ns/day, 11.124 hours/ns, 49.941 timesteps/s
288.8% CPU use with 1 MPI tasks x 48 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 1.9998 | 1.9998 | 1.9998 | 0.0 | 99.87
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.0011814 | 0.0011814 | 0.0011814 | 0.0 | 0.06
Output | 0.00059724 | 0.00059724 | 0.00059724 | 0.0 | 0.03
Modify | 0.00047352 | 0.00047352 | 0.00047352 | 0.0 | 0.02
Other | | 0.0003468 | | | 0.02
Nlocal: 128.000 ave 128 max 128 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 727.000 ave 727 max 727 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 3712.00 ave 3712 max 3712 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 7424.00 ave 7424 max 7424 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 7424
Ave neighs/atom = 58.000000
Neighbor list builds = 0
Dangerous builds = 0
Total wall time: 0:00:03

157
examples/mliap/log.04Dec20.mliap.pytorch.Ta06A.g++.4 Normal file
View File

@ -0,0 +1,157 @@
LAMMPS (30 Nov 2020)
using 48 OpenMP thread(s) per MPI task
# Demonstrate MLIAP/PyTorch interface to linear SNAP potential
# Initialize simulation
variable nsteps index 100
variable nrep equal 4
variable a equal 3.316
units metal
# generate the box and atom positions using a BCC lattice
variable nx equal ${nrep}
variable nx equal 4
variable ny equal ${nrep}
variable ny equal 4
variable nz equal ${nrep}
variable nz equal 4
boundary p p p
lattice bcc $a
lattice bcc 3.316
Lattice spacing in x,y,z = 3.3160000 3.3160000 3.3160000
region box block 0 ${nx} 0 ${ny} 0 ${nz}
region box block 0 4 0 ${ny} 0 ${nz}
region box block 0 4 0 4 0 ${nz}
region box block 0 4 0 4 0 4
create_box 1 box
Created orthogonal box = (0.0000000 0.0000000 0.0000000) to (13.264000 13.264000 13.264000)
1 by 2 by 2 MPI processor grid
create_atoms 1 box
Created 128 atoms
create_atoms CPU = 0.002 seconds
mass 1 180.88
# choose potential
include Ta06A.mliap.pytorch
# DATE: 2014-09-05 UNITS: metal CONTRIBUTOR: Aidan Thompson athomps@sandia.gov CITATION: Thompson, Swiler, Trott, Foiles and Tucker, arxiv.org, 1409.3880 (2014)
# Definition of SNAP potential Ta_Cand06A
# Assumes 1 LAMMPS atom type
variable zblcutinner equal 4
variable zblcutouter equal 4.8
variable zblz equal 73
# Specify hybrid with SNAP, ZBL
pair_style hybrid/overlay zbl ${zblcutinner} ${zblcutouter} mliap model mliappy Ta06A.mliap.pytorch.model.pkl descriptor sna Ta06A.mliap.descriptor
pair_style hybrid/overlay zbl 4 ${zblcutouter} mliap model mliappy Ta06A.mliap.pytorch.model.pkl descriptor sna Ta06A.mliap.descriptor
pair_style hybrid/overlay zbl 4 4.8 mliap model mliappy Ta06A.mliap.pytorch.model.pkl descriptor sna Ta06A.mliap.descriptor
Loading python model complete.
Reading potential file Ta06A.mliap.descriptor with DATE: 2014-09-05
SNAP keyword rcutfac 4.67637
SNAP keyword twojmax 6
SNAP keyword nelems 1
SNAP keyword elems Ta
SNAP keyword radelems 0.5
SNAP keyword welems 1
SNAP keyword rfac0 0.99363
SNAP keyword rmin0 0
SNAP keyword bzeroflag 0
pair_coeff 1 1 zbl ${zblz} ${zblz}
pair_coeff 1 1 zbl 73 ${zblz}
pair_coeff 1 1 zbl 73 73
pair_coeff * * mliap Ta
# Setup output
compute eatom all pe/atom
compute energy all reduce sum c_eatom
compute satom all stress/atom NULL
compute str all reduce sum c_satom[1] c_satom[2] c_satom[3]
variable press equal (c_str[1]+c_str[2]+c_str[3])/(3*vol)
thermo_style custom step temp epair c_energy etotal press v_press
thermo 10
thermo_modify norm yes
# Set up NVE run
timestep 0.5e-3
neighbor 1.0 bin
neigh_modify once no every 1 delay 0 check yes
# Run MD
velocity all create 300.0 4928459 loop geom
fix 1 all nve
run ${nsteps}
run 100
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 5.8
ghost atom cutoff = 5.8
binsize = 2.9, bins = 5 5 5
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair zbl, perpetual, half/full from (2)
attributes: half, newton on
pair build: halffull/newton
stencil: none
bin: none
(2) pair mliap, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 159.7 | 159.7 | 159.7 Mbytes
Step Temp E_pair c_energy TotEng Press v_press
0 300 -11.85157 -11.85157 -11.813095 2717.1661 -2717.1661
10 296.01467 -11.851059 -11.851059 -11.813095 2697.4796 -2697.4796
20 284.53666 -11.849587 -11.849587 -11.813095 2289.1527 -2289.1527
30 266.51577 -11.847275 -11.847275 -11.813095 1851.7131 -1851.7131
40 243.05007 -11.844266 -11.844266 -11.813095 1570.684 -1570.684
50 215.51032 -11.840734 -11.840734 -11.813094 1468.1899 -1468.1899
60 185.48331 -11.836883 -11.836883 -11.813094 1524.8757 -1524.8757
70 154.6736 -11.832931 -11.832931 -11.813094 1698.3351 -1698.3351
80 124.79303 -11.829099 -11.829099 -11.813094 1947.0715 -1947.0715
90 97.448054 -11.825592 -11.825592 -11.813094 2231.9563 -2231.9563
100 74.035418 -11.822589 -11.822589 -11.813094 2515.8526 -2515.8526
Loop time of 0.562802 on 192 procs for 100 steps with 128 atoms
Performance: 7.676 ns/day, 3.127 hours/ns, 177.682 timesteps/s
99.7% CPU use with 4 MPI tasks x 48 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.53583 | 0.54622 | 0.55401 | 0.9 | 97.05
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.0071442 | 0.01491 | 0.025289 | 5.4 | 2.65
Output | 0.00092525 | 0.00095771 | 0.0010166 | 0.0 | 0.17
Modify | 0.00014479 | 0.00015043 | 0.00015893 | 0.0 | 0.03
Other | | 0.0005624 | | | 0.10
Nlocal: 32.0000 ave 32 max 32 min
Histogram: 4 0 0 0 0 0 0 0 0 0
Nghost: 431.000 ave 431 max 431 min
Histogram: 4 0 0 0 0 0 0 0 0 0
Neighs: 928.000 ave 928 max 928 min
Histogram: 4 0 0 0 0 0 0 0 0 0
FullNghs: 1856.00 ave 1856 max 1856 min
Histogram: 4 0 0 0 0 0 0 0 0 0
Total # of neighbors = 7424
Ave neighs/atom = 58.000000
Neighbor list builds = 0
Dangerous builds = 0
Total wall time: 0:00:02

104
examples/mliap/mliap_pytorch_Ta06A.py Normal file
View File

@ -0,0 +1,104 @@
# Demonstrate how to load a model from the python side.
# This is essentially the same as in.mliap.pytorch.Ta06A
# except that python is the driving program, and lammps
# is in library mode.
before_loading =\
"""# Demonstrate MLIAP/PyTorch interface to linear SNAP potential
# Initialize simulation
variable nsteps index 100
variable nrep equal 4
variable a equal 3.316
units metal
# generate the box and atom positions using a BCC lattice
variable nx equal ${nrep}
variable ny equal ${nrep}
variable nz equal ${nrep}
boundary p p p
lattice bcc $a
region box block 0 ${nx} 0 ${ny} 0 ${nz}
create_box 1 box
create_atoms 1 box
mass 1 180.88
# choose potential
# DATE: 2014-09-05 UNITS: metal CONTRIBUTOR: Aidan Thompson athomps@sandia.gov CITATION: Thompson, Swiler, Trott, Foiles and Tucker, arxiv.org, 1409.3880 (2014)
# Definition of SNAP potential Ta_Cand06A
# Assumes 1 LAMMPS atom type
variable zblcutinner equal 4
variable zblcutouter equal 4.8
variable zblz equal 73
# Specify hybrid with SNAP, ZBL
pair_style hybrid/overlay &
zbl ${zblcutinner} ${zblcutouter} &
mliap model mliappy LATER &
descriptor sna Ta06A.mliap.descriptor
pair_coeff 1 1 zbl ${zblz} ${zblz}
pair_coeff * * mliap Ta
"""
after_loading =\
"""
# Setup output
compute eatom all pe/atom
compute energy all reduce sum c_eatom
compute satom all stress/atom NULL
compute str all reduce sum c_satom[1] c_satom[2] c_satom[3]
variable press equal (c_str[1]+c_str[2]+c_str[3])/(3*vol)
thermo_style custom step temp epair c_energy etotal press v_press
thermo 10
thermo_modify norm yes
# Set up NVE run
timestep 0.5e-3
neighbor 1.0 bin
neigh_modify once no every 1 delay 0 check yes
# Run MD
velocity all create 300.0 4928459 loop geom
fix 1 all nve
run ${nsteps}
"""
import lammps
lmp = lammps.lammps(cmdargs=['-echo','both'])
# Before defining the pair style, one must do the following:
import lammps.mliap
lammps.mliap.activate_mliappy(lmp)
# Otherwise, when running lammps in library mode,
# you will get an error:
# "ERROR: Loading MLIAPPY coupling module failure."
# Setup the simulation and declare an empty model
# by specifying model filename as "LATER"
lmp.commands_string(before_loading)
# Define the model however you like. In this example
# we load it from disk:
import torch
model = torch.load('Ta06A.mliap.pytorch.model.pt')
# Connect the PyTorch model to the mliap pair style.
lammps.mliap.load_model(model)
# run the simulation with the mliap pair style
lmp.commands_string(after_loading)

18
examples/mliap/relu1hidden.mliap.pytorch Normal file
View File

@ -0,0 +1,18 @@
# DATE: 2014-09-05 UNITS: metal CONTRIBUTOR: Aidan Thompson athomps@sandia.gov CITATION: Thompson, Swiler, Trott, Foiles and Tucker, arxiv.org, 1409.3880 (2014)
# Definition of SNAP potential Ta_Cand06A
# Assumes 1 LAMMPS atom type
variable zblcutinner equal 4
variable zblcutouter equal 4.8
variable zblz equal 73
# Specify hybrid with SNAP, ZBL
pair_style hybrid/overlay &
zbl ${zblcutinner} ${zblcutouter} &
mliap model mliappy relu1hidden.mliap.pytorch.model.pt &
descriptor sna Ta06A.mliap.descriptor
pair_coeff 1 1 zbl ${zblz} ${zblz}
pair_coeff * * mliap Ta

BIN
examples/mliap/relu1hidden.mliap.pytorch.model.pt Normal file
View File

Binary file not shown.

3
lib/python/Makefile.mliap_python Normal file
View File

@ -0,0 +1,3 @@
../mliap_model_python_couple.cpp: ../mliap_model_python_couple.pyx
cythonize -3 ../mliap_model_python_couple.cpp

35
python/install.py
View File

@ -98,19 +98,23 @@ os.chdir(os.path.dirname(args.package))
from distutils.core import setup from distutils.core import setup
from distutils.sysconfig import get_python_lib from distutils.sysconfig import get_python_lib
import site import site
tryuser=False
#Arguments common to global or user install -- everything but data_files
setup_kwargs= dict(name="lammps",
version=verstr,
author="Steve Plimpton",
author_email="sjplimp@sandia.gov",
url="https://lammps.sandia.gov",
description="LAMMPS Molecular Dynamics Python package",
license="GPL",
packages=["lammps","lammps.mliap"],
)
tryuser=False
try: try:
sys.argv = ["setup.py","install"] # as if had run "python setup.py install" sys.argv = ["setup.py","install"] # as if had run "python setup.py install"
setup(name = "lammps", setup_kwargs['data_files']=[(os.path.join(get_python_lib(), 'lammps'), [args.lib])]
version = verstr, setup(**setup_kwargs)
author = "Steve Plimpton",
author_email = "sjplimp@sandia.gov",
url = "https://lammps.sandia.gov",
description = "LAMMPS Molecular Dynamics Python package",
license = "GPL",
packages=['lammps'],
data_files = [(os.path.join(get_python_lib(), 'lammps'), [args.lib])])
except: except:
tryuser=True tryuser=True
print ("Installation into global site-packages folder failed.\nTrying user folder %s now." % site.USER_SITE) print ("Installation into global site-packages folder failed.\nTrying user folder %s now." % site.USER_SITE)
@ -118,14 +122,7 @@ except:
if tryuser: if tryuser:
try: try:
sys.argv = ["setup.py","install","--user"] # as if had run "python setup.py install --user" sys.argv = ["setup.py","install","--user"] # as if had run "python setup.py install --user"
setup(name = "lammps", setup_kwargs['data_files']=[(os.path.join(site.USER_SITE, 'lammps'), [args.lib])]
version = verstr, setup(**setup_kwargs)
author = "Steve Plimpton",
author_email = "sjplimp@sandia.gov",
url = "https://lammps.sandia.gov",
description = "LAMMPS Molecular Dynamics Python package",
license = "GPL",
packages=['lammps'],
data_files = [(os.path.join(site.USER_SITE, 'lammps'), [args.lib])])
except: except:
print("Installation into user site package folder failed.") print("Installation into user site package folder failed.")

13
python/lammps/mliap/__init__.py Normal file
View File

@ -0,0 +1,13 @@
# Check compatiblity of this build with the python shared library.
# If this fails, lammps will segfault because its library will
# try to improperly start up a new interpreter.
import sysconfig
import ctypes
library = sysconfig.get_config_vars('INSTSONAME')[0]
pylib = ctypes.CDLL(library)
if not pylib.Py_IsInitialized():
raise RuntimeError("This interpreter is not compatible with python-based mliap for LAMMPS.")
del sysconfig, ctypes, library, pylib
from .loader import load_model, activate_mliappy

52
python/lammps/mliap/loader.py Normal file
View File

@ -0,0 +1,52 @@
# ----------------------------------------------------------------------
# LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
# http://lammps.sandia.gov, Sandia National Laboratories
# Steve Plimpton, sjplimp@sandia.gov
#
# Copyright (2003) Sandia Corporation. Under the terms of Contract
# DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
# certain rights in this software. This software is distributed under
# the GNU General Public License.
#
# See the README file in the top-level LAMMPS directory.
# -------------------------------------------------------------------------
# ----------------------------------------------------------------------
# Contributing author: Nicholas Lubbers (LANL)
# -------------------------------------------------------------------------
import sys
import importlib.util
import importlib.machinery
def activate_mliappy(lmp):
try:
# Begin Importlib magic to find the embedded python module
# This is needed because the filename for liblammps does not
# match the spec for normal python modules, wherein
# file names match with PyInit function names.
# Also, python normally doesn't look for extensions besides '.so'
# We fix both of these problems by providing an explict
# path to the extension module 'mliap_model_python_couple' in
path = lmp.lib._name
loader = importlib.machinery.ExtensionFileLoader('mliap_model_python_couple', path)
spec = importlib.util.spec_from_loader('mliap_model_python_couple', loader)
module = importlib.util.module_from_spec(spec)
sys.modules['mliap_model_python_couple'] = module
spec.loader.exec_module(module)
# End Importlib magic to find the embedded python module
except Exception as ee:
raise ImportError("Could not load MLIAP python coupling module.") from ee
def load_model(model):
try:
import mliap_model_python_couple
except ImportError as ie:
raise ImportError("MLIAP python module must be activated before loading\n"
"the pair style. Call lammps.mliap.activate_mliappy(lmp)."
) from ie
mliap_model_python_couple.load_from_python(model)

65
python/lammps/mliap/pytorch.py Normal file
View File

@ -0,0 +1,65 @@
# ----------------------------------------------------------------------
# LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
# http://lammps.sandia.gov, Sandia National Laboratories
# Steve Plimpton, sjplimp@sandia.gov
#
# Copyright (2003) Sandia Corporation. Under the terms of Contract
# DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
# certain rights in this software. This software is distributed under
# the GNU General Public License.
#
# See the README file in the top-level LAMMPS directory.
# -------------------------------------------------------------------------
# ----------------------------------------------------------------------
# Contributing author: Nicholas Lubbers (LANL)
# -------------------------------------------------------------------------
import numpy as np
import torch
def calc_n_params(model):
return sum(p.nelement() for p in model.parameters())
class TorchWrapper(torch.nn.Module):
def __init__(self, model,n_descriptors,n_elements,n_params=None,device=None,dtype=torch.float64):
super().__init__()
self.model = model
self.device = device
self.dtype = dtype
# Put model on device and convert to dtype
self.to(self.dtype)
self.to(self.device)
if n_params is None:
n_params = calc_n_params(model)
self.n_params = n_params
self.n_descriptors = n_descriptors
self.n_elements = n_elements
def forward(self, elems, bispectrum, beta, energy):
bispectrum = torch.from_numpy(bispectrum).to(dtype=self.dtype, device=self.device).requires_grad_(True)
elems = torch.from_numpy(elems).to(dtype=torch.long, device=self.device) - 1
with torch.autograd.enable_grad():
energy_nn = self.model(bispectrum, elems)
if energy_nn.ndim > 1:
energy_nn = energy_nn.flatten()
beta_nn = torch.autograd.grad(energy_nn.sum(), bispectrum)[0]
beta[:] = beta_nn.detach().cpu().numpy().astype(np.float64)
energy[:] = energy_nn.detach().cpu().numpy().astype(np.float64)
class IgnoreElems(torch.nn.Module):
def __init__(self,subnet):
super().__init__()
self.subnet = subnet
def forward(self,bispectrum,elems):
return self.subnet(bispectrum)

2
python/setup.py
View File

@ -22,5 +22,5 @@ setup(
url = "https://lammps.sandia.gov", url = "https://lammps.sandia.gov",
description = "LAMMPS Molecular Dynamics Python package", description = "LAMMPS Molecular Dynamics Python package",
license = "GPL", license = "GPL",
packages=["lammps"] packages=["lammps","lammps.mliap"],
) )

5
src/Depend.sh
View File

@ -106,6 +106,10 @@ if (test $1 = "PERI") then
depend USER-OMP depend USER-OMP
fi fi
if (test $1 = "PYTHON") then
depend MLIAP
fi
if (test $1 = "RIGID") then if (test $1 = "RIGID") then
depend KOKKOS depend KOKKOS
depend USER-OMP depend USER-OMP
@ -114,6 +118,7 @@ fi
if (test $1 = "SNAP") then if (test $1 = "SNAP") then
depend KOKKOS depend KOKKOS
depend MLIAP
fi fi
if (test $1 = "USER-CGSDK") then if (test $1 = "USER-CGSDK") then

67
src/MLIAP/Install.sh Executable file
View File

@ -0,0 +1,67 @@
# Install/unInstall package files in LAMMPS
# mode = 0/1/2 for uninstall/install/update
mode=$1
# arg1 = file, arg2 = file it depends on
# enforce using portable C locale
LC_ALL=C
export LC_ALL
action () {
if (test $mode = 0) then
rm -f ../$1
elif (! cmp -s $1 ../$1) then
if (test -z "$2" || test -e ../$2) then
cp $1 ..
if (test $mode = 2) then
echo " updating src/$1"
fi
fi
elif (test -n "$2") then
if (test ! -e ../$2) then
rm -f ../$1
fi
fi
}
# enforce package dependency
if (test $1 = 1 || test $1 = 2) then
if (test ! -e ../sna.h) then
echo "Must install SNAP package to use MLIAP package"
exit 1
fi
fi
# all package files with no dependencies
for file in *.cpp *.h *.pyx; do
test -f ${file} && action $file
done
# edit 2 Makefile.package files to include/exclude package info
if (test $1 = 1) then
if (test "$(type cythonize 2> /dev/null)" != "" && test -e ../python_impl.cpp) then
if (test -e ../Makefile.package) then
sed -i -e 's|^PKG_INC =[ \t]*|&-DMLIAP_PYTHON |' ../Makefile.package
fi
if (test -e ../Makefile.package.settings) then
sed -i -e '/^include.*python.*mliap_python.*$/d' ../Makefile.package.settings
# multiline form needed for BSD sed on Macs
sed -i -e '4 i \
include ..\/..\/lib\/python\/Makefile.mliap_python
' ../Makefile.package.settings
fi
cythonize -3 ../mliap_model_python_couple.pyx
fi
elif (test $1 = 0) then
if (test -e ../Makefile.package) then
sed -i -e 's/[^ \t]*-DMLIAP_PYTHON[^ \t]* //g' ../Makefile.package
fi
rm -f ../mliap_model_python_couple.cpp ../mliap_model_python_couple.h
sed -i -e '/^include.*python.*mliap_python.*$/d' ../Makefile.package.settings
fi

41
src/MLIAP/README
View File

@ -8,12 +8,15 @@ definitions of the interatomic potential functional form (*model*)
and the geometric quantities that characterize the atomic positions and the geometric quantities that characterize the atomic positions
(*descriptor*). By defining *model* and *descriptor* separately, (*descriptor*). By defining *model* and *descriptor* separately,
it is possible to use many different models with a given descriptor, it is possible to use many different models with a given descriptor,
or many different descriptors with a given model. Currently, the pair_style or many different descriptors with a given model. The pair_style
supports just two models, *linear* and *quadratic*, supports the following models: *linear*, *quadratic*, and
and one descriptor, *sna*, the SNAP descriptor, including the *mliappy* (general Python interface to things like PyTorch, see below
linear, quadratic, and chem variants. Work is currently underway to extend for build instructions).
the interface to handle neural network energy models, It currently supports only one class of descriptors,
and it is also straightforward to add new descriptor styles. *sna*, the SNAP descriptors, including the
linear, quadratic, and chem variants.
It is straightforward to add new descriptor and model
styles.
The mliap compute style provides gradients of the energy, force, The mliap compute style provides gradients of the energy, force,
and stress tensor w.r.t. model parameters. and stress tensor w.r.t. model parameters.
@ -27,4 +30,28 @@ reference potential.
To see how this command To see how this command
can be used within a Python workflow to train machine-learning interatomic can be used within a Python workflow to train machine-learning interatomic
potentials, see the examples in FitSNAP https://github.com/FitSNAP/FitSNAP>. potentials, see the examples in FitSNAP https://github.com/FitSNAP/FitSNAP.
*Additional instructions for building MLIAP with Python support enabled*
The *mliappy* energy model requires that the MLIAP package
be compiled with Python support enabled.
This extension, written by Nick Lubbers (LANL),
provides a coupling to PyTorch and other Python modules.
This should be automatically
enabled by default if the prerequisite software is installed. It can be
enforced during CMake configuration by setting the variable
MLIAP_ENABLE_PYTHON=yes or for conventional build by adding -DMLIAP_PYTHON
to the LMP_INC variable in your makefile and running the cythonize script
on the .pyx file(s) copied to the src folder.
This requires to also install the PYTHON package and have the cython
(https://cython.org) software installed. During configuration/compilation
the cythonize script will be used to convert the provided .pyx file(s)
to C++ code. Please do not run the cythonize script in the src/MLIAP folder.
If you have done so by accident, you need to delete the generated .cpp and .h
file(s) in the src/MLIAP folder or there may be problems during compilation.
More information on building LAMMPS with this package is here:
https://lammps.sandia.gov/doc/Build_extras.html#mliap

17
src/MLIAP/compute_mliap.cpp
View File

@ -11,12 +11,20 @@
See the README file in the top-level LAMMPS directory. See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Aidan Thompson (SNL)
------------------------------------------------------------------------- */
#include <cstring> #include <cstring>
#include "mliap_data.h" #include "mliap_data.h"
#include "mliap_model_linear.h" #include "mliap_model_linear.h"
#include "mliap_model_quadratic.h" #include "mliap_model_quadratic.h"
#include "mliap_descriptor_snap.h" #include "mliap_descriptor_snap.h"
#ifdef MLIAP_PYTHON
#include "mliap_model_python.h"
#endif
#include "compute_mliap.h" #include "compute_mliap.h"
#include "atom.h" #include "atom.h"
#include "update.h" #include "update.h"
@ -65,7 +73,14 @@ ComputeMLIAP::ComputeMLIAP(LAMMPS *lmp, int narg, char **arg) :
} else if (strcmp(arg[iarg+1],"quadratic") == 0) { } else if (strcmp(arg[iarg+1],"quadratic") == 0) {
model = new MLIAPModelQuadratic(lmp); model = new MLIAPModelQuadratic(lmp);
iarg += 2; iarg += 2;
} else error->all(FLERR,"Illegal compute mliap command"); }
#ifdef MLIAP_PYTHON
else if (strcmp(arg[iarg+1],"mliappy") == 0) {
model = new MLIAPModelPython(lmp);
iarg += 2;
}
#endif
else error->all(FLERR,"Illegal compute mliap command");
modelflag = 1; modelflag = 1;
} else if (strcmp(arg[iarg],"descriptor") == 0) { } else if (strcmp(arg[iarg],"descriptor") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal compute mliap command"); if (iarg+2 > narg) error->all(FLERR,"Illegal compute mliap command");

10
src/MLIAP/mliap_data.cpp
View File

@ -11,6 +11,10 @@
See the README file in the top-level LAMMPS directory. See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Aidan Thompson (SNL)
------------------------------------------------------------------------- */
#include "mliap_data.h" #include "mliap_data.h"
#include "atom.h" #include "atom.h"
@ -25,7 +29,8 @@ MLIAPData::MLIAPData(LAMMPS *lmp, int gradgradflag_in, int *map_in,
class MLIAPModel* model_in, class MLIAPModel* model_in,
class MLIAPDescriptor* descriptor_in, class MLIAPDescriptor* descriptor_in,
class PairMLIAP* pairmliap_in) : class PairMLIAP* pairmliap_in) :
Pointers(lmp), gradforce(nullptr), betas(nullptr), descriptors(nullptr), gamma(nullptr), Pointers(lmp), gradforce(nullptr), betas(nullptr),
descriptors(nullptr), eatoms(nullptr), gamma(nullptr),
gamma_row_index(nullptr), gamma_col_index(nullptr), egradient(nullptr), gamma_row_index(nullptr), gamma_col_index(nullptr), egradient(nullptr),
numneighs(nullptr), iatoms(nullptr), ielems(nullptr), jatoms(nullptr), jelems(nullptr), numneighs(nullptr), iatoms(nullptr), ielems(nullptr), jatoms(nullptr), jelems(nullptr),
rij(nullptr), graddesc(nullptr), model(nullptr), descriptor(nullptr), list(nullptr) rij(nullptr), graddesc(nullptr), model(nullptr), descriptor(nullptr), list(nullptr)
@ -64,6 +69,7 @@ MLIAPData::~MLIAPData()
{ {
memory->destroy(betas); memory->destroy(betas);
memory->destroy(descriptors); memory->destroy(descriptors);
memory->destroy(eatoms);
memory->destroy(gamma_row_index); memory->destroy(gamma_row_index);
memory->destroy(gamma_col_index); memory->destroy(gamma_col_index);
memory->destroy(gamma); memory->destroy(gamma);
@ -133,6 +139,7 @@ void MLIAPData::generate_neighdata(NeighList* list_in, int eflag_in, int vflag_i
if (natoms_max < natoms) { if (natoms_max < natoms) {
memory->grow(betas,natoms,ndescriptors,"MLIAPData:betas"); memory->grow(betas,natoms,ndescriptors,"MLIAPData:betas");
memory->grow(descriptors,natoms,ndescriptors,"MLIAPData:descriptors"); memory->grow(descriptors,natoms,ndescriptors,"MLIAPData:descriptors");
memory->grow(eatoms,natoms,"MLIAPData:eatoms");
natoms_max = natoms; natoms_max = natoms;
} }
@ -267,6 +274,7 @@ double MLIAPData::memory_usage()
bytes += natoms*ndescriptors*sizeof(int); // betas bytes += natoms*ndescriptors*sizeof(int); // betas
bytes += natoms*ndescriptors*sizeof(int); // descriptors bytes += natoms*ndescriptors*sizeof(int); // descriptors
bytes += natoms*sizeof(double); // eatoms
bytes += natomneigh_max*sizeof(int); // iatoms bytes += natomneigh_max*sizeof(int); // iatoms
bytes += natomneigh_max*sizeof(int); // ielems bytes += natomneigh_max*sizeof(int); // ielems

2
src/MLIAP/mliap_data.h
View File

@ -36,6 +36,8 @@ class MLIAPData : protected Pointers {
double **gradforce; double **gradforce;
double** betas; // betas for all atoms in list double** betas; // betas for all atoms in list
double** descriptors; // descriptors for all atoms in list double** descriptors; // descriptors for all atoms in list
double* eatoms; // energies for all atoms in list
double energy; // energy
int ndescriptors; // number of descriptors int ndescriptors; // number of descriptors
int nparams; // number of model parameters per element int nparams; // number of model parameters per element
int nelements; // number of elements int nelements; // number of elements

4
src/MLIAP/mliap_descriptor.cpp
View File

@ -11,6 +11,10 @@
See the README file in the top-level LAMMPS directory. See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Aidan Thompson (SNL)
------------------------------------------------------------------------- */
#include "mliap_descriptor.h" #include "mliap_descriptor.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;

4
src/MLIAP/mliap_descriptor_snap.cpp
View File

@ -11,6 +11,10 @@
See the README file in the top-level LAMMPS directory. See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Aidan Thompson (SNL)
------------------------------------------------------------------------- */
#include "mliap_descriptor_snap.h" #include "mliap_descriptor_snap.h"
#include "atom.h" #include "atom.h"

29
src/MLIAP/mliap_model.cpp
View File

@ -11,6 +11,10 @@
See the README file in the top-level LAMMPS directory. See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Aidan Thompson (SNL)
------------------------------------------------------------------------- */
#include "mliap_model.h" #include "mliap_model.h"
#include "comm.h" #include "comm.h"
@ -28,13 +32,9 @@ using namespace LAMMPS_NS;
MLIAPModel::MLIAPModel(LAMMPS* lmp, char* coefffilename) : Pointers(lmp) MLIAPModel::MLIAPModel(LAMMPS* lmp, char* coefffilename) : Pointers(lmp)
{ {
coeffelem = nullptr;
if (coefffilename) read_coeffs(coefffilename);
else {
nparams = 0; nparams = 0;
nelements = 0; nelements = 0;
ndescriptors = 0; ndescriptors = 0;
}
nonlinearflag = 0; nonlinearflag = 0;
} }
@ -42,9 +42,9 @@ MLIAPModel::MLIAPModel(LAMMPS* lmp, char* coefffilename) : Pointers(lmp)
MLIAPModel::~MLIAPModel() MLIAPModel::~MLIAPModel()
{ {
memory->destroy(coeffelem);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
placeholder placeholder
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
@ -73,7 +73,22 @@ void MLIAPModel::set_ndescriptors(int ndescriptors_in)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void MLIAPModel::read_coeffs(char *coefffilename)
MLIAPModelSimple::MLIAPModelSimple(LAMMPS* lmp, char* coefffilename) : MLIAPModel(lmp, coefffilename)
{
coeffelem = nullptr;
if (coefffilename) read_coeffs(coefffilename);
}
/* ---------------------------------------------------------------------- */
MLIAPModelSimple::~MLIAPModelSimple()
{
memory->destroy(coeffelem);
}
void MLIAPModelSimple::read_coeffs(char *coefffilename)
{ {
// open coefficient file on proc 0 // open coefficient file on proc 0
@ -165,7 +180,7 @@ void MLIAPModel::read_coeffs(char *coefffilename)
memory usage memory usage
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
double MLIAPModel::memory_usage() double MLIAPModelSimple::memory_usage()
{ {
double bytes = 0; double bytes = 0;

17
src/MLIAP/mliap_model.h
View File

@ -30,15 +30,26 @@ public:
virtual void compute_gradgrads(class MLIAPData*)=0; virtual void compute_gradgrads(class MLIAPData*)=0;
virtual void compute_force_gradients(class MLIAPData*)=0; virtual void compute_force_gradients(class MLIAPData*)=0;
virtual void init(); virtual void init();
virtual double memory_usage(); virtual double memory_usage()=0;
int nelements; // # of unique elements int nelements; // # of unique elements
int nonlinearflag; // 1 if gradient() requires escriptors int nonlinearflag; // 1 if gradient() requires descriptors
int ndescriptors; // number of descriptors int ndescriptors; // number of descriptors
int nparams; // number of parameters per element int nparams; // number of parameters per element
protected: protected:
void read_coeffs(char *); virtual void read_coeffs(char *)=0;
};
class MLIAPModelSimple : public MLIAPModel {
public:
MLIAPModelSimple(LAMMPS*, char*);
~MLIAPModelSimple();
virtual double memory_usage();
protected:
double **coeffelem; // element coefficients double **coeffelem; // element coefficients
virtual void read_coeffs(char *);
}; };
} }

12
src/MLIAP/mliap_model_linear.cpp
View File

@ -11,6 +11,10 @@
See the README file in the top-level LAMMPS directory. See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Aidan Thompson (SNL)
------------------------------------------------------------------------- */
#include "mliap_model_linear.h" #include "mliap_model_linear.h"
#include "pair_mliap.h" #include "pair_mliap.h"
#include "mliap_data.h" #include "mliap_data.h"
@ -21,7 +25,7 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
MLIAPModelLinear::MLIAPModelLinear(LAMMPS* lmp, char* coefffilename) : MLIAPModelLinear::MLIAPModelLinear(LAMMPS* lmp, char* coefffilename) :
MLIAPModel(lmp, coefffilename) MLIAPModelSimple(lmp, coefffilename)
{ {
if (nparams > 0) ndescriptors = nparams - 1; if (nparams > 0) ndescriptors = nparams - 1;
} }
@ -51,8 +55,9 @@ int MLIAPModelLinear::get_nparams()
void MLIAPModelLinear::compute_gradients(MLIAPData* data) void MLIAPModelLinear::compute_gradients(MLIAPData* data)
{ {
data->energy = 0.0;
for (int ii = 0; ii < data->natoms; ii++) { for (int ii = 0; ii < data->natoms; ii++) {
const int i = data->iatoms[ii];
const int ielem = data->ielems[ii]; const int ielem = data->ielems[ii];
double* coeffi = coeffelem[ielem]; double* coeffi = coeffelem[ielem];
@ -74,7 +79,8 @@ void MLIAPModelLinear::compute_gradients(MLIAPData* data)
for (int icoeff = 0; icoeff < data->ndescriptors; icoeff++) for (int icoeff = 0; icoeff < data->ndescriptors; icoeff++)
etmp += coeffi[icoeff+1]*data->descriptors[ii][icoeff]; etmp += coeffi[icoeff+1]*data->descriptors[ii][icoeff];
data->pairmliap->e_tally(i,etmp); data->energy += etmp;
data->eatoms[ii] = etmp;
} }
} }
} }

2
src/MLIAP/mliap_model_linear.h
View File

@ -18,7 +18,7 @@
namespace LAMMPS_NS { namespace LAMMPS_NS {
class MLIAPModelLinear : public MLIAPModel { class MLIAPModelLinear : public MLIAPModelSimple {
public: public:
MLIAPModelLinear(LAMMPS*, char* = nullptr); MLIAPModelLinear(LAMMPS*, char* = nullptr);
~MLIAPModelLinear(); ~MLIAPModelLinear();

202
src/MLIAP/mliap_model_python.cpp Normal file
View File

@ -0,0 +1,202 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Nicholas Lubbers (LANL)
------------------------------------------------------------------------- */
#ifdef MLIAP_PYTHON
#include <Python.h>
#include "mliap_model_python.h"
#include "mliap_model_python_couple.h"
#include "pair_mliap.h"
#include "mliap_data.h"
#include "error.h"
#include "utils.h"
#include "lmppython.h"
#include "python_compat.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
MLIAPModelPython::MLIAPModelPython(LAMMPS* lmp, char* coefffilename) :
MLIAPModel(lmp, coefffilename)
{
model_loaded = 0;
python->init();
PyGILState_STATE gstate = PyGILState_Ensure();
PyObject * pyMain = PyImport_AddModule("__main__");
if (!pyMain) {
PyGILState_Release(gstate);
error->all(FLERR,"Could not initialize embedded Python");
}
PyObject* coupling_module = PyImport_ImportModule("mliap_model_python_couple");
if (!coupling_module) {
PyErr_Print();
PyErr_Clear();
PyGILState_Release(gstate);
error->all(FLERR,"Loading MLIAPPY coupling module failure.");
}
// Recipe from lammps/src/pair_python.cpp :
// add current directory to PYTHONPATH
PyObject * py_path = PySys_GetObject((char *)"path");
PyList_Append(py_path, PY_STRING_FROM_STRING("."));
// if LAMMPS_POTENTIALS environment variable is set, add it to PYTHONPATH as well
const char * potentials_path = getenv("LAMMPS_POTENTIALS");
if (potentials_path != NULL) {
PyList_Append(py_path, PY_STRING_FROM_STRING(potentials_path));
}
PyGILState_Release(gstate);
if (coefffilename) read_coeffs(coefffilename);
nonlinearflag=1;
}
/* ---------------------------------------------------------------------- */
MLIAPModelPython::~MLIAPModelPython(){
MLIAPPY_unload_model(this);
}
/* ----------------------------------------------------------------------
get number of parameters
---------------------------------------------------------------------- */
int MLIAPModelPython::get_nparams()
{
return nparams;
}
void MLIAPModelPython::read_coeffs(char * fname)
{
PyGILState_STATE gstate = PyGILState_Ensure();
int loaded = MLIAPPY_load_model(this, fname);
if (PyErr_Occurred()) {
PyErr_Print();
PyErr_Clear();
PyGILState_Release(gstate);
error->all(FLERR,"Loading python model failure.");
}
PyGILState_Release(gstate);
if (loaded) {
this->connect_param_counts();
}
else {
utils::logmesg(lmp,"Loading python model deferred.\n");
}
}
// Finalize loading of the model.
void MLIAPModelPython::connect_param_counts()
{
PyGILState_STATE gstate = PyGILState_Ensure();
nelements = MLIAPPY_nelements(this);
nparams = MLIAPPY_nparams(this);
ndescriptors = MLIAPPY_ndescriptors(this);
if (PyErr_Occurred()) {
PyErr_Print();
PyErr_Clear();
PyGILState_Release(gstate);
error->all(FLERR,"Loading python model failure.");
}
PyGILState_Release(gstate);
model_loaded = 1;
utils::logmesg(lmp,"Loading python model complete.\n");
}
/* ----------------------------------------------------------------------
Calculate model gradients w.r.t descriptors
for each atom beta_i = dE(B_i)/dB_i
---------------------------------------------------------------------- */
void MLIAPModelPython::compute_gradients(MLIAPData* data)
{
if (not model_loaded) {
error->all(FLERR,"Model not loaded.");
}
PyGILState_STATE gstate = PyGILState_Ensure();
MLIAPPY_compute_gradients(this, data);
if (PyErr_Occurred()) {
PyErr_Print();
PyErr_Clear();
PyGILState_Release(gstate);
error->all(FLERR,"Running python model failure.");
}
PyGILState_Release(gstate);
}
/* ----------------------------------------------------------------------
Calculate model double gradients w.r.t descriptors and parameters
for each atom energy gamma_lk = d2E(B)/dB_k/dsigma_l,
where sigma_l is a parameter, B_k a descriptor,
and atom subscript i is omitted
gamma is in CSR format:
nnz = number of non-zero values
gamma_row_index[inz] = l indices, 0 <= l < nparams
gamma_col_indexiinz] = k indices, 0 <= k < ndescriptors
gamma[i][inz] = non-zero values, 0 <= inz < nnz
egradient is derivative of energy w.r.t. parameters
---------------------------------------------------------------------- */
void MLIAPModelPython::compute_gradgrads(class MLIAPData* data)
{
error->all(FLERR,"compute_gradgrads not implemented");
}
/* ----------------------------------------------------------------------
calculate gradients of forces w.r.t. parameters
egradient is derivative of energy w.r.t. parameters
---------------------------------------------------------------------- */
void MLIAPModelPython::compute_force_gradients(class MLIAPData* data)
{
error->all(FLERR,"compute_force_gradients not implemented");
}
/* ----------------------------------------------------------------------
count the number of non-zero entries in gamma matrix
---------------------------------------------------------------------- */
int MLIAPModelPython::get_gamma_nnz(class MLIAPData* data)
{
// todo: get_gamma_nnz
return 0;
}
double MLIAPModelPython::memory_usage()
{
// todo: get approximate memory usage in coupling code.
return 0;
}
#endif

47
src/MLIAP/mliap_model_python.h Normal file
View File

@ -0,0 +1,47 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifndef LMP_MLIAP_MODEL_PYTHON_H
#define LMP_MLIAP_MODEL_PYTHON_H
#include "mliap_model.h"
namespace LAMMPS_NS {
class MLIAPModelPython : public MLIAPModel {
public:
MLIAPModelPython(LAMMPS*, char* = NULL);
~MLIAPModelPython();
virtual int get_nparams();
virtual int get_gamma_nnz(class MLIAPData*);
virtual void compute_gradients(class MLIAPData*);
virtual void compute_gradgrads(class MLIAPData*);
virtual void compute_force_gradients(class MLIAPData*);
virtual double memory_usage();
void connect_param_counts(); // If possible convert this to protected/private and
// and figure out how to declare cython fn
// load_from_python as a friend.
int model_loaded;
protected:
virtual void read_coeffs(char *);
private:
};
}
#endif

120
src/MLIAP/mliap_model_python_couple.pyx Normal file
View File

@ -0,0 +1,120 @@
# cython: language_level=3
# distutils: language = c++
cimport cython
import pickle
# For converting C arrays to numpy arrays
import numpy as np
# For converting void * to integer for tracking object identity
from libc.stdint cimport uintptr_t
from libcpp.string cimport string
cdef extern from "mliap_data.h" namespace "LAMMPS_NS":
cdef cppclass MLIAPData:
# Array shapes
int natoms
int ndescriptors
# Input data
int * ielems # types for all atoms in list
double ** descriptors # descriptors for all atoms in list
# Output data to write to
double ** betas # betas for all atoms in list
double * eatoms # energy for all atoms in list
double energy
cdef extern from "mliap_model_python.h" namespace "LAMMPS_NS":
cdef cppclass MLIAPModelPython:
void connect_param_counts()
class MLIAPPYModelNotLinked(Exception): pass
LOADED_MODELS = {}
cdef object c_id(MLIAPModelPython * c_model):
"""
Use python-style id of object to keep track of identity.
Note, this is probably not a perfect general strategy but it should work fine with LAMMPS pair styles.
"""
return int(<uintptr_t> c_model)
cdef object retrieve(MLIAPModelPython * c_model) with gil:
try:
model = LOADED_MODELS[c_id(c_model)]
except KeyError as ke:
raise KeyError("Model has not been loaded.") from ke
if model is None:
raise MLIAPPYModelNotLinked("Model not linked, connect the model from the python side.")
return model
cdef public int MLIAPPY_load_model(MLIAPModelPython * c_model, char* fname) with gil:
str_fname = fname.decode('utf-8') # Python 3 only; not Python 2 not supported.
if str_fname == "LATER":
model = None
returnval = 0
else:
if str_fname.endswith(".pt") or str_fname.endswith('.pth'):
import torch
model = torch.load(str_fname)
else:
with open(str_fname,'rb') as pfile:
model = pickle.load(pfile)
returnval = 1
LOADED_MODELS[c_id(c_model)] = model
return returnval
def load_from_python(model):
unloaded_models = [k for k, v in LOADED_MODELS.items() if v is None]
num_models = len(unloaded_models)
cdef MLIAPModelPython * lmp_model
if num_models == 0:
raise ValueError("No model in the waiting area.")
elif num_models > 1:
raise ValueError("Model is amibguous, more than one model in waiting area.")
else:
c_id = unloaded_models[0]
LOADED_MODELS[c_id]=model
lmp_model = <MLIAPModelPython *> <uintptr_t> c_id
lmp_model.connect_param_counts()
cdef public void MLIAPPY_unload_model(MLIAPModelPython * c_model) with gil:
del LOADED_MODELS[c_id(c_model)]
cdef public int MLIAPPY_nparams(MLIAPModelPython * c_model) with gil:
return int(retrieve(c_model).n_params)
cdef public int MLIAPPY_nelements(MLIAPModelPython * c_model) with gil:
return int(retrieve(c_model).n_elements)
cdef public int MLIAPPY_ndescriptors(MLIAPModelPython * c_model) with gil:
return int(retrieve(c_model).n_descriptors)
cdef public void MLIAPPY_compute_gradients(MLIAPModelPython * c_model, MLIAPData * data) with gil:
model = retrieve(c_model)
n_d = data.ndescriptors
n_a = data.natoms
# Make numpy arrays from pointers
beta_np = np.asarray(<double[:n_a,:n_d] > &data.betas[0][0])
desc_np = np.asarray(<double[:n_a,:n_d]> &data.descriptors[0][0])
elem_np = np.asarray(<int[:n_a]> &data.ielems[0])
en_np = np.asarray(<double[:n_a]> &data.eatoms[0])
# Invoke python model on numpy arrays.
model(elem_np,desc_np,beta_np,en_np)
# Get the total energy from the atom energy.
energy = np.sum(en_np)
data.energy = <double> energy
return

13
src/MLIAP/mliap_model_quadratic.cpp
View File

@ -11,6 +11,10 @@
See the README file in the top-level LAMMPS directory. See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Aidan Thompson (SNL)
------------------------------------------------------------------------- */
#include "mliap_model_quadratic.h" #include "mliap_model_quadratic.h"
#include "pair_mliap.h" #include "pair_mliap.h"
#include "mliap_data.h" #include "mliap_data.h"
@ -22,8 +26,9 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
MLIAPModelQuadratic::MLIAPModelQuadratic(LAMMPS* lmp, char* coefffilename) : MLIAPModelQuadratic::MLIAPModelQuadratic(LAMMPS* lmp, char* coefffilename) :
MLIAPModel(lmp, coefffilename) MLIAPModelSimple(lmp, coefffilename)
{ {
if (coefffilename) read_coeffs(coefffilename);
if (nparams > 0) ndescriptors = sqrt(2*nparams)-1; if (nparams > 0) ndescriptors = sqrt(2*nparams)-1;
nonlinearflag = 1; nonlinearflag = 1;
} }
@ -52,8 +57,9 @@ int MLIAPModelQuadratic::get_nparams()
void MLIAPModelQuadratic::compute_gradients(MLIAPData* data) void MLIAPModelQuadratic::compute_gradients(MLIAPData* data)
{ {
data->energy = 0.0;
for (int ii = 0; ii < data->natoms; ii++) { for (int ii = 0; ii < data->natoms; ii++) {
const int i = data->iatoms[ii];
const int ielem = data->ielems[ii]; const int ielem = data->ielems[ii];
double* coeffi = coeffelem[ielem]; double* coeffi = coeffelem[ielem];
@ -99,7 +105,8 @@ void MLIAPModelQuadratic::compute_gradients(MLIAPData* data)
etmp += coeffi[k++]*bveci*bvecj; etmp += coeffi[k++]*bveci*bvecj;
} }
} }
data->pairmliap->e_tally(i,etmp); data->energy += etmp;
data->eatoms[ii] = etmp;
} }
} }
} }

2
src/MLIAP/mliap_model_quadratic.h
View File

@ -18,7 +18,7 @@
namespace LAMMPS_NS { namespace LAMMPS_NS {
class MLIAPModelQuadratic : public MLIAPModel { class MLIAPModelQuadratic : public MLIAPModelSimple {
public: public:
MLIAPModelQuadratic(LAMMPS*, char* = nullptr); MLIAPModelQuadratic(LAMMPS*, char* = nullptr);
~MLIAPModelQuadratic(); ~MLIAPModelQuadratic();

45
src/MLIAP/pair_mliap.cpp
View File

@ -11,12 +11,19 @@
See the README file in the top-level LAMMPS directory. See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Aidan Thompson (SNL)
------------------------------------------------------------------------- */
#include "pair_mliap.h" #include "pair_mliap.h"
#include "mliap_data.h" #include "mliap_data.h"
#include "mliap_model_linear.h" #include "mliap_model_linear.h"
#include "mliap_model_quadratic.h" #include "mliap_model_quadratic.h"
#include "mliap_descriptor_snap.h" #include "mliap_descriptor_snap.h"
#ifdef MLIAP_PYTHON
#include "mliap_model_python.h"
#endif
#include "atom.h" #include "atom.h"
#include "error.h" #include "error.h"
@ -27,6 +34,7 @@
#include <cmath> #include <cmath>
#include <cstring> #include <cstring>
#include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
@ -65,6 +73,17 @@ PairMLIAP::~PairMLIAP()
void PairMLIAP::compute(int eflag, int vflag) void PairMLIAP::compute(int eflag, int vflag)
{ {
// consistency checks
if (data->ndescriptors != model->ndescriptors) {
error->all(FLERR,"Incompatible model and descriptor descriptor count");
};
if (data->nelements != model->nelements) {
error->all(FLERR,"Incompatible model and descriptor element count");
};
ev_init(eflag,vflag); ev_init(eflag,vflag);
data->generate_neighdata(list, eflag, vflag); data->generate_neighdata(list, eflag, vflag);
@ -78,6 +97,8 @@ void PairMLIAP::compute(int eflag, int vflag)
model->compute_gradients(data); model->compute_gradients(data);
e_tally(data);
// calculate force contributions beta_i*dB_i/dR_j // calculate force contributions beta_i*dB_i/dR_j
descriptor->compute_forces(data); descriptor->compute_forces(data);
@ -107,6 +128,7 @@ void PairMLIAP::allocate()
void PairMLIAP::settings(int narg, char ** arg) void PairMLIAP::settings(int narg, char ** arg)
{ {
if (narg < 4) if (narg < 4)
error->all(FLERR,"Illegal pair_style command"); error->all(FLERR,"Illegal pair_style command");
@ -130,6 +152,12 @@ void PairMLIAP::settings(int narg, char ** arg)
if (iarg+3 > narg) error->all(FLERR,"Illegal pair_style mliap command"); if (iarg+3 > narg) error->all(FLERR,"Illegal pair_style mliap command");
model = new MLIAPModelQuadratic(lmp,arg[iarg+2]); model = new MLIAPModelQuadratic(lmp,arg[iarg+2]);
iarg += 3; iarg += 3;
#ifdef MLIAP_PYTHON
} else if (strcmp(arg[iarg+1],"mliappy") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal pair_style mliap command");
model = new MLIAPModelPython(lmp,arg[iarg+2]);
iarg += 3;
#endif
} else error->all(FLERR,"Illegal pair_style mliap command"); } else error->all(FLERR,"Illegal pair_style mliap command");
modelflag = 1; modelflag = 1;
} else if (strcmp(arg[iarg],"descriptor") == 0) { } else if (strcmp(arg[iarg],"descriptor") == 0) {
@ -211,22 +239,21 @@ void PairMLIAP::coeff(int narg, char **arg)
data = new MLIAPData(lmp, gradgradflag, map, model, descriptor, this); data = new MLIAPData(lmp, gradgradflag, map, model, descriptor, this);
data->init(); data->init();
// consistency checks
if (data->ndescriptors != model->ndescriptors)
error->all(FLERR,"Incompatible model and descriptor definitions");
if (data->nelements != model->nelements)
error->all(FLERR,"Incompatible model and descriptor definitions");
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
add energy of atom i to global and per-atom energy add energies to eng_vdwl and per-atom energy
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void PairMLIAP::e_tally(int i, double ei) void PairMLIAP::e_tally(MLIAPData* data)
{ {
if (eflag_global) eng_vdwl += ei; if (eflag_global) eng_vdwl += data->energy;
if (eflag_atom) eatom[i] += ei; if (eflag_atom)
for (int ii = 0; ii < data->natoms; ii++) {
const int i = data->iatoms[ii];
eatom[i] += data->eatoms[ii];
}
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------

2
src/MLIAP/pair_mliap.h
View File

@ -31,7 +31,7 @@ public:
virtual void compute(int, int); virtual void compute(int, int);
void settings(int, char **); void settings(int, char **);
virtual void coeff(int, char **); virtual void coeff(int, char **);
void e_tally(int, double); void e_tally(class MLIAPData*);
void v_tally(int, int, double*, double*); void v_tally(int, int, double*, double*);
virtual void init_style(); virtual void init_style();
virtual double init_one(int, int); virtual double init_one(int, int);

16
src/PYTHON/python_impl.cpp
View File

@ -26,6 +26,14 @@
#include <cstring> #include <cstring>
#include <Python.h> // IWYU pragma: export #include <Python.h> // IWYU pragma: export
#ifdef MLIAP_PYTHON
#include "mliap_model_python.h"
// The above should somehow really be included in the next file.
// We could get around this with cython --capi-reexport-cincludes
// However, that exposes -too many- headers.
#include "mliap_model_python_couple.h"
#endif
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
enum{NONE,INT,DOUBLE,STRING,PTR}; enum{NONE,INT,DOUBLE,STRING,PTR};
@ -47,11 +55,17 @@ PythonImpl::PythonImpl(LAMMPS *lmp) : Pointers(lmp)
nfunc = 0; nfunc = 0;
pfuncs = nullptr; pfuncs = nullptr;
// one-time initialization of Python interpreter // one-time initialization of Python interpreter
// pyMain stores pointer to main module // pyMain stores pointer to main module
external_interpreter = Py_IsInitialized(); external_interpreter = Py_IsInitialized();
#ifdef MLIAP_PYTHON
// Inform python intialization scheme of the mliappy module.
// This -must- happen before python is initialized.
int err = PyImport_AppendInittab("mliap_model_python_couple", PyInit_mliap_model_python_couple);
if (err) error->all(FLERR,"Could not register MLIAPPY embedded python module.");
#endif
Py_Initialize(); Py_Initialize();
// only needed for Python 2.x and Python 3 < 3.7 // only needed for Python 2.x and Python 3 < 3.7

2
src/Purge.list
View File

@ -49,6 +49,8 @@ packages_ntopo.h
# other auto-generated files # other auto-generated files
lmpinstalledpkgs.h lmpinstalledpkgs.h
lmpgitversion.h lmpgitversion.h
mliap_model_python_couple.cpp
mliap_model_python_couple.h
# removed on 9 Sep 2020 # removed on 9 Sep 2020
mergesort.h mergesort.h
# renamed on 8 May 2020 # renamed on 8 May 2020