diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt index 696a1b8cdd..12b713f15c 100644 --- a/cmake/CMakeLists.txt +++ b/cmake/CMakeLists.txt @@ -373,7 +373,7 @@ else() set(CUDA_REQUEST_PIC) 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) if(PKG_${PKG_WITH_INCL}) include(Packages/${PKG_WITH_INCL}) diff --git a/cmake/Modules/FindCythonize.cmake b/cmake/Modules/FindCythonize.cmake new file mode 100644 index 0000000000..9a37904ea7 --- /dev/null +++ b/cmake/Modules/FindCythonize.cmake @@ -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) diff --git a/cmake/Modules/LAMMPSUtils.cmake b/cmake/Modules/LAMMPSUtils.cmake index 339ba867bd..37275843fa 100644 --- a/cmake/Modules/LAMMPSUtils.cmake +++ b/cmake/Modules/LAMMPSUtils.cmake @@ -50,6 +50,7 @@ function(check_for_autogen_files source_dir) file(GLOB SRC_AUTOGEN_FILES ${source_dir}/style_*.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}/mliap_model_python_couple.h ${source_dir}/mliap_model_python_couple.cpp) foreach(_SRC ${SRC_AUTOGEN_FILES}) get_filename_component(FILENAME "${_SRC}" NAME) if(EXISTS ${source_dir}/${FILENAME}) diff --git a/cmake/Modules/Packages/MLIAP.cmake b/cmake/Modules/Packages/MLIAP.cmake new file mode 100644 index 0000000000..d3f601a1e1 --- /dev/null +++ b/cmake/Modules/Packages/MLIAP.cmake @@ -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() diff --git a/doc/src/Build_extras.rst b/doc/src/Build_extras.rst index 4dd17d71f8..78165f0150 100644 --- a/doc/src/Build_extras.rst +++ b/doc/src/Build_extras.rst @@ -37,6 +37,7 @@ This is the list of packages that may require additional steps. * :ref:`KOKKOS ` * :ref:`LATTE ` * :ref:`MESSAGE ` + * :ref:`MLIAP ` * :ref:`MSCG ` * :ref:`OPT ` * :ref:`POEMS ` @@ -770,6 +771,54 @@ be installed on your system. ---------- +.. _mliap: + +MLIAP package +--------------------------- + +Building the MLIAP package requires including the :ref:`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 ` +package **and** requires you have the `cython `_ 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 package diff --git a/doc/src/Packages_details.rst b/doc/src/Packages_details.rst index 2b2cb90aff..e044adfcb3 100644 --- a/doc/src/Packages_details.rst +++ b/doc/src/Packages_details.rst @@ -662,19 +662,31 @@ MLIAP package **Contents:** -A general interface for machine-learning interatomic potentials. +A general interface for machine-learning interatomic potentials, including PyTorch. **Install:** -To use this package, also the :ref:`SNAP package ` needs to be installed. +To use this package, also the :ref:`SNAP package ` package needs +to be installed. To make the *mliappy* model available, also the +:ref:`PYTHON package ` package needs to be installed, the version of +Python must be 3.6 or later, and the `cython `_ software +must be installed. -**Author:** Aidan Thompson (Sandia). +**Author:** Aidan Thompson (Sandia), Nicholas Lubbers (LANL). **Supporting info:** * src/MLIAP: filenames -> commands +* src/MLIAP/README * :doc:`pair_style mliap ` -* examples/mliap +* :doc:`compute_style 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. ---------- diff --git a/doc/src/compute_mliap.rst b/doc/src/compute_mliap.rst index fcc336e682..579086ad4a 100644 --- a/doc/src/compute_mliap.rst +++ b/doc/src/compute_mliap.rst @@ -18,7 +18,7 @@ Syntax .. parsed-literal:: *model* values = style - style = *linear* or *quadratic* + style = *linear* or *quadratic* or *mliappy* *descriptor* values = style filename style = *sna* 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 *model* and *descriptor* in either order. -The *model* keyword is followed by a model style, currently limited to -either *linear* or *quadratic*. +The *model* keyword is followed by the model style (*linear*, *quadratic* or *mliappy*). +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. -Currently the only descriptor style is *sna*, indicating the bispectrum component -descriptors used by the Spectral Neighbor Analysis Potential (SNAP) potentials of -:doc:`pair_style snap `. +The compute currently supports just one descriptor style, but it is +is straightforward to add new descriptor styles. +The SNAP descriptor style *sna* is the same as that used by :doc:`pair_style snap `, +including the linear, quadratic, and chem variants. A single additional argument specifies the descriptor filename containing the parameters and setting used by the SNAP descriptor. The descriptor filename usually ends in the *.mliap.descriptor* extension. @@ -162,9 +164,10 @@ potentials, see the examples in `FitSNAP `_. Restrictions """""""""""" -This compute 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 +This compute 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 requires building LAMMPS with the SNAP package. +The *mliappy* model requires building LAMMPS with the PYTHON package. See the :doc:`Build package ` doc page for more info. Related commands diff --git a/doc/src/pair_mliap.rst b/doc/src/pair_mliap.rst index 09295cd8be..af0cc9e855 100644 --- a/doc/src/pair_mliap.rst +++ b/doc/src/pair_mliap.rst @@ -16,7 +16,7 @@ Syntax .. parsed-literal:: *model* values = style filename - style = *linear* or *quadratic* + style = *linear* or *quadratic* or *mliappy* filename = name of file containing model definitions *descriptor* values = style filename style = *sna* @@ -40,12 +40,15 @@ definitions of the interatomic potential functional form (*model*) and the geometric quantities that characterize the atomic positions (*descriptor*). By defining *model* and *descriptor* separately, it is possible to use many different models with a given descriptor, -or many different descriptors with a given model. Currently, the pair_style -supports just two models, *linear* and *quadratic*, -and one descriptor, *sna*, the SNAP descriptor used by :doc:`pair_style snap `, including the linear, quadratic, -and chem variants. Work is currently underway to extend -the interface to handle neural network energy models, -and it is also straightforward to add new descriptor styles. +or many different descriptors with a given model. The +pair style currently supports just one descriptor style, but it is +is straightforward to add new descriptor styles. +The SNAP descriptor style *sna* is the same as that used by :doc:`pair_style snap `, +including the linear, quadratic, and chem variants. +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 of energy, force, and stress w.r.t. model parameters. This information can be accessed using the related :doc:`compute mliap ` command. @@ -59,9 +62,8 @@ that specify the mapping of MLIAP element names to LAMMPS atom types, where N is the number of LAMMPS atom types. -The *model* keyword is followed by a model style, currently limited to -either *linear* or *quadratic*. In both cases, -this is followed by a single argument specifying the model filename containing the +The *model* keyword is followed by the model style. This is followed +by a single argument specifying the model filename containing the parameters for a set of elements. The model filename usually ends in the *.mliap.model* extension. It may contain parameters for many elements. The only requirement is that it @@ -82,6 +84,16 @@ for the :doc:`pair_style snap ` coefficient file. Specifically, the line containing the element weight and radius is omitted, 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. Currently the only descriptor style is *sna*, indicating the bispectrum component 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 """""""""""" -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 requires building LAMMPS with the SNAP package. +The *mliappy* model requires building LAMMPS with the PYTHON package. See the :doc:`Build package ` doc page for more info. + Related commands """""""""""""""" diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index de2065c1e3..0dff97459a 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -558,6 +558,7 @@ Cygwin cylindrically Cyrot cyrstals +cython Daivis Dammak dampflag diff --git a/examples/mliap/README b/examples/mliap/README new file mode 100644 index 0000000000..af848e42ba --- /dev/null +++ b/examples/mliap/README @@ -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. + diff --git a/examples/mliap/Ta06A.mliap.pytorch b/examples/mliap/Ta06A.mliap.pytorch new file mode 100644 index 0000000000..5489688a82 --- /dev/null +++ b/examples/mliap/Ta06A.mliap.pytorch @@ -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 + diff --git a/examples/mliap/convert_mliap_Ta06A.py b/examples/mliap/convert_mliap_Ta06A.py new file mode 100644 index 0000000000..8a7466743f --- /dev/null +++ b/examples/mliap/convert_mliap_Ta06A.py @@ -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") diff --git a/examples/mliap/in.mliap.pytorch.Ta06A b/examples/mliap/in.mliap.pytorch.Ta06A new file mode 100644 index 0000000000..59cdcc8c81 --- /dev/null +++ b/examples/mliap/in.mliap.pytorch.Ta06A @@ -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} + diff --git a/examples/mliap/in.mliap.pytorch.relu1hidden b/examples/mliap/in.mliap.pytorch.relu1hidden new file mode 100644 index 0000000000..f19ffb608d --- /dev/null +++ b/examples/mliap/in.mliap.pytorch.relu1hidden @@ -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} + diff --git a/examples/mliap/in.mliap.snap.Ta06A b/examples/mliap/in.mliap.snap.Ta06A index 3d94d5c9fc..7b9b5b47d1 100644 --- a/examples/mliap/in.mliap.snap.Ta06A +++ b/examples/mliap/in.mliap.snap.Ta06A @@ -1,4 +1,4 @@ -# Demonstrate MLIAP interface to kinear SNAP potential +# Demonstrate MLIAP interface to linear SNAP potential # Initialize simulation diff --git a/examples/mliap/log.04Dec20.mliap.pytorch.Ta06A.g++.1 b/examples/mliap/log.04Dec20.mliap.pytorch.Ta06A.g++.1 new file mode 100644 index 0000000000..f56a79650e --- /dev/null +++ b/examples/mliap/log.04Dec20.mliap.pytorch.Ta06A.g++.1 @@ -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 diff --git a/examples/mliap/log.04Dec20.mliap.pytorch.Ta06A.g++.4 b/examples/mliap/log.04Dec20.mliap.pytorch.Ta06A.g++.4 new file mode 100644 index 0000000000..7f1bc818d1 --- /dev/null +++ b/examples/mliap/log.04Dec20.mliap.pytorch.Ta06A.g++.4 @@ -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 diff --git a/examples/mliap/mliap_pytorch_Ta06A.py b/examples/mliap/mliap_pytorch_Ta06A.py new file mode 100644 index 0000000000..c4387e21a3 --- /dev/null +++ b/examples/mliap/mliap_pytorch_Ta06A.py @@ -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) diff --git a/examples/mliap/relu1hidden.mliap.pytorch b/examples/mliap/relu1hidden.mliap.pytorch new file mode 100644 index 0000000000..32ac4c85a3 --- /dev/null +++ b/examples/mliap/relu1hidden.mliap.pytorch @@ -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 + diff --git a/examples/mliap/relu1hidden.mliap.pytorch.model.pt b/examples/mliap/relu1hidden.mliap.pytorch.model.pt new file mode 100644 index 0000000000..ecc04341a7 Binary files /dev/null and b/examples/mliap/relu1hidden.mliap.pytorch.model.pt differ diff --git a/lib/python/Makefile.mliap_python b/lib/python/Makefile.mliap_python new file mode 100644 index 0000000000..00483a4a6a --- /dev/null +++ b/lib/python/Makefile.mliap_python @@ -0,0 +1,3 @@ + +../mliap_model_python_couple.cpp: ../mliap_model_python_couple.pyx + cythonize -3 ../mliap_model_python_couple.cpp diff --git a/python/install.py b/python/install.py index a6b69c1ee6..6765c33925 100644 --- a/python/install.py +++ b/python/install.py @@ -98,19 +98,23 @@ os.chdir(os.path.dirname(args.package)) from distutils.core import setup from distutils.sysconfig import get_python_lib 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: sys.argv = ["setup.py","install"] # as if had run "python setup.py install" - setup(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'], - data_files = [(os.path.join(get_python_lib(), 'lammps'), [args.lib])]) + setup_kwargs['data_files']=[(os.path.join(get_python_lib(), 'lammps'), [args.lib])] + setup(**setup_kwargs) except: tryuser=True print ("Installation into global site-packages folder failed.\nTrying user folder %s now." % site.USER_SITE) @@ -118,14 +122,7 @@ except: if tryuser: try: sys.argv = ["setup.py","install","--user"] # as if had run "python setup.py install --user" - setup(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'], - data_files = [(os.path.join(site.USER_SITE, 'lammps'), [args.lib])]) + setup_kwargs['data_files']=[(os.path.join(site.USER_SITE, 'lammps'), [args.lib])] + setup(**setup_kwargs) except: print("Installation into user site package folder failed.") diff --git a/python/lammps/mliap/__init__.py b/python/lammps/mliap/__init__.py new file mode 100644 index 0000000000..0d63cc810b --- /dev/null +++ b/python/lammps/mliap/__init__.py @@ -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 diff --git a/python/lammps/mliap/loader.py b/python/lammps/mliap/loader.py new file mode 100644 index 0000000000..b875d8d834 --- /dev/null +++ b/python/lammps/mliap/loader.py @@ -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) + diff --git a/python/lammps/mliap/pytorch.py b/python/lammps/mliap/pytorch.py new file mode 100644 index 0000000000..aa2cf1c97c --- /dev/null +++ b/python/lammps/mliap/pytorch.py @@ -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) diff --git a/python/setup.py b/python/setup.py index 9be04138d5..aff0b14671 100644 --- a/python/setup.py +++ b/python/setup.py @@ -22,5 +22,5 @@ setup( url = "https://lammps.sandia.gov", description = "LAMMPS Molecular Dynamics Python package", license = "GPL", - packages=["lammps"] + packages=["lammps","lammps.mliap"], ) diff --git a/src/Depend.sh b/src/Depend.sh index e629ca9197..f77d435fc5 100755 --- a/src/Depend.sh +++ b/src/Depend.sh @@ -106,6 +106,10 @@ if (test $1 = "PERI") then depend USER-OMP fi +if (test $1 = "PYTHON") then + depend MLIAP +fi + if (test $1 = "RIGID") then depend KOKKOS depend USER-OMP @@ -114,6 +118,7 @@ fi if (test $1 = "SNAP") then depend KOKKOS + depend MLIAP fi if (test $1 = "USER-CGSDK") then diff --git a/src/MLIAP/Install.sh b/src/MLIAP/Install.sh new file mode 100755 index 0000000000..e6ba4042a9 --- /dev/null +++ b/src/MLIAP/Install.sh @@ -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 diff --git a/src/MLIAP/README b/src/MLIAP/README index 18ce3297e6..59eb5e34b3 100644 --- a/src/MLIAP/README +++ b/src/MLIAP/README @@ -8,12 +8,15 @@ definitions of the interatomic potential functional form (*model*) and the geometric quantities that characterize the atomic positions (*descriptor*). By defining *model* and *descriptor* separately, it is possible to use many different models with a given descriptor, -or many different descriptors with a given model. Currently, the pair_style -supports just two models, *linear* and *quadratic*, -and one descriptor, *sna*, the SNAP descriptor, including the -linear, quadratic, and chem variants. Work is currently underway to extend -the interface to handle neural network energy models, -and it is also straightforward to add new descriptor styles. +or many different descriptors with a given model. The pair_style +supports the following models: *linear*, *quadratic*, and +*mliappy* (general Python interface to things like PyTorch, see below +for build instructions). +It currently supports only one class of descriptors, +*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, and stress tensor w.r.t. model parameters. @@ -27,4 +30,28 @@ reference potential. To see how this command 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 diff --git a/src/MLIAP/compute_mliap.cpp b/src/MLIAP/compute_mliap.cpp index 34d086462f..0842e421c0 100644 --- a/src/MLIAP/compute_mliap.cpp +++ b/src/MLIAP/compute_mliap.cpp @@ -11,12 +11,20 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing author: Aidan Thompson (SNL) +------------------------------------------------------------------------- */ + #include #include "mliap_data.h" #include "mliap_model_linear.h" #include "mliap_model_quadratic.h" #include "mliap_descriptor_snap.h" +#ifdef MLIAP_PYTHON +#include "mliap_model_python.h" +#endif + #include "compute_mliap.h" #include "atom.h" #include "update.h" @@ -65,7 +73,14 @@ ComputeMLIAP::ComputeMLIAP(LAMMPS *lmp, int narg, char **arg) : } else if (strcmp(arg[iarg+1],"quadratic") == 0) { model = new MLIAPModelQuadratic(lmp); 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; } else if (strcmp(arg[iarg],"descriptor") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute mliap command"); diff --git a/src/MLIAP/mliap_data.cpp b/src/MLIAP/mliap_data.cpp index e6502d001a..c310f0efc7 100644 --- a/src/MLIAP/mliap_data.cpp +++ b/src/MLIAP/mliap_data.cpp @@ -11,6 +11,10 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing author: Aidan Thompson (SNL) +------------------------------------------------------------------------- */ + #include "mliap_data.h" #include "atom.h" @@ -25,7 +29,8 @@ MLIAPData::MLIAPData(LAMMPS *lmp, int gradgradflag_in, int *map_in, class MLIAPModel* model_in, class MLIAPDescriptor* descriptor_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), numneighs(nullptr), iatoms(nullptr), ielems(nullptr), jatoms(nullptr), jelems(nullptr), rij(nullptr), graddesc(nullptr), model(nullptr), descriptor(nullptr), list(nullptr) @@ -64,6 +69,7 @@ MLIAPData::~MLIAPData() { memory->destroy(betas); memory->destroy(descriptors); + memory->destroy(eatoms); memory->destroy(gamma_row_index); memory->destroy(gamma_col_index); memory->destroy(gamma); @@ -133,6 +139,7 @@ void MLIAPData::generate_neighdata(NeighList* list_in, int eflag_in, int vflag_i if (natoms_max < natoms) { memory->grow(betas,natoms,ndescriptors,"MLIAPData:betas"); memory->grow(descriptors,natoms,ndescriptors,"MLIAPData:descriptors"); + memory->grow(eatoms,natoms,"MLIAPData:eatoms"); natoms_max = natoms; } @@ -267,6 +274,7 @@ double MLIAPData::memory_usage() bytes += natoms*ndescriptors*sizeof(int); // betas bytes += natoms*ndescriptors*sizeof(int); // descriptors + bytes += natoms*sizeof(double); // eatoms bytes += natomneigh_max*sizeof(int); // iatoms bytes += natomneigh_max*sizeof(int); // ielems diff --git a/src/MLIAP/mliap_data.h b/src/MLIAP/mliap_data.h index ffac9ccd4c..7fb60bd723 100644 --- a/src/MLIAP/mliap_data.h +++ b/src/MLIAP/mliap_data.h @@ -34,8 +34,10 @@ class MLIAPData : protected Pointers { int yoffset, zoffset; int ndims_force, ndims_virial; 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* eatoms; // energies for all atoms in list + double energy; // energy int ndescriptors; // number of descriptors int nparams; // number of model parameters per element int nelements; // number of elements diff --git a/src/MLIAP/mliap_descriptor.cpp b/src/MLIAP/mliap_descriptor.cpp index b3e0cebd2e..ac168dc94b 100644 --- a/src/MLIAP/mliap_descriptor.cpp +++ b/src/MLIAP/mliap_descriptor.cpp @@ -11,6 +11,10 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing author: Aidan Thompson (SNL) +------------------------------------------------------------------------- */ + #include "mliap_descriptor.h" using namespace LAMMPS_NS; diff --git a/src/MLIAP/mliap_descriptor_snap.cpp b/src/MLIAP/mliap_descriptor_snap.cpp index cfbc51cb5c..2cb88de6b4 100644 --- a/src/MLIAP/mliap_descriptor_snap.cpp +++ b/src/MLIAP/mliap_descriptor_snap.cpp @@ -11,6 +11,10 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing author: Aidan Thompson (SNL) +------------------------------------------------------------------------- */ + #include "mliap_descriptor_snap.h" #include "atom.h" diff --git a/src/MLIAP/mliap_model.cpp b/src/MLIAP/mliap_model.cpp index 04b2ac663a..aa198e6085 100644 --- a/src/MLIAP/mliap_model.cpp +++ b/src/MLIAP/mliap_model.cpp @@ -11,6 +11,10 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing author: Aidan Thompson (SNL) +------------------------------------------------------------------------- */ + #include "mliap_model.h" #include "comm.h" @@ -28,13 +32,9 @@ using namespace LAMMPS_NS; MLIAPModel::MLIAPModel(LAMMPS* lmp, char* coefffilename) : Pointers(lmp) { - coeffelem = nullptr; - if (coefffilename) read_coeffs(coefffilename); - else { - nparams = 0; - nelements = 0; - ndescriptors = 0; - } + nparams = 0; + nelements = 0; + ndescriptors = 0; nonlinearflag = 0; } @@ -42,9 +42,9 @@ MLIAPModel::MLIAPModel(LAMMPS* lmp, char* coefffilename) : Pointers(lmp) MLIAPModel::~MLIAPModel() { - memory->destroy(coeffelem); } + /* ---------------------------------------------------------------------- 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 @@ -165,7 +180,7 @@ void MLIAPModel::read_coeffs(char *coefffilename) memory usage ------------------------------------------------------------------------- */ -double MLIAPModel::memory_usage() +double MLIAPModelSimple::memory_usage() { double bytes = 0; diff --git a/src/MLIAP/mliap_model.h b/src/MLIAP/mliap_model.h index b2f7312897..96cfff0a3d 100644 --- a/src/MLIAP/mliap_model.h +++ b/src/MLIAP/mliap_model.h @@ -30,15 +30,26 @@ public: virtual void compute_gradgrads(class MLIAPData*)=0; virtual void compute_force_gradients(class MLIAPData*)=0; virtual void init(); - virtual double memory_usage(); + virtual double memory_usage()=0; 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 nparams; // number of parameters per element 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 + virtual void read_coeffs(char *); + }; } diff --git a/src/MLIAP/mliap_model_linear.cpp b/src/MLIAP/mliap_model_linear.cpp index b3224284e5..3c8f9aa39a 100644 --- a/src/MLIAP/mliap_model_linear.cpp +++ b/src/MLIAP/mliap_model_linear.cpp @@ -11,6 +11,10 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing author: Aidan Thompson (SNL) +------------------------------------------------------------------------- */ + #include "mliap_model_linear.h" #include "pair_mliap.h" #include "mliap_data.h" @@ -21,7 +25,7 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ MLIAPModelLinear::MLIAPModelLinear(LAMMPS* lmp, char* coefffilename) : - MLIAPModel(lmp, coefffilename) + MLIAPModelSimple(lmp, coefffilename) { if (nparams > 0) ndescriptors = nparams - 1; } @@ -51,8 +55,9 @@ int MLIAPModelLinear::get_nparams() void MLIAPModelLinear::compute_gradients(MLIAPData* data) { + data->energy = 0.0; + for (int ii = 0; ii < data->natoms; ii++) { - const int i = data->iatoms[ii]; const int ielem = data->ielems[ii]; double* coeffi = coeffelem[ielem]; @@ -74,7 +79,8 @@ void MLIAPModelLinear::compute_gradients(MLIAPData* data) for (int icoeff = 0; icoeff < data->ndescriptors; icoeff++) etmp += coeffi[icoeff+1]*data->descriptors[ii][icoeff]; - data->pairmliap->e_tally(i,etmp); + data->energy += etmp; + data->eatoms[ii] = etmp; } } } diff --git a/src/MLIAP/mliap_model_linear.h b/src/MLIAP/mliap_model_linear.h index 326407faa8..89c69d3537 100644 --- a/src/MLIAP/mliap_model_linear.h +++ b/src/MLIAP/mliap_model_linear.h @@ -18,7 +18,7 @@ namespace LAMMPS_NS { -class MLIAPModelLinear : public MLIAPModel { +class MLIAPModelLinear : public MLIAPModelSimple { public: MLIAPModelLinear(LAMMPS*, char* = nullptr); ~MLIAPModelLinear(); diff --git a/src/MLIAP/mliap_model_python.cpp b/src/MLIAP/mliap_model_python.cpp new file mode 100644 index 0000000000..e47356390f --- /dev/null +++ b/src/MLIAP/mliap_model_python.cpp @@ -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 +#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 diff --git a/src/MLIAP/mliap_model_python.h b/src/MLIAP/mliap_model_python.h new file mode 100644 index 0000000000..57110f36e2 --- /dev/null +++ b/src/MLIAP/mliap_model_python.h @@ -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 + diff --git a/src/MLIAP/mliap_model_python_couple.pyx b/src/MLIAP/mliap_model_python_couple.pyx new file mode 100644 index 0000000000..1f04964ef6 --- /dev/null +++ b/src/MLIAP/mliap_model_python_couple.pyx @@ -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( 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 = 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( &data.betas[0][0]) + desc_np = np.asarray( &data.descriptors[0][0]) + elem_np = np.asarray( &data.ielems[0]) + en_np = np.asarray( &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 = energy + return diff --git a/src/MLIAP/mliap_model_quadratic.cpp b/src/MLIAP/mliap_model_quadratic.cpp index 481417078a..915ec30435 100644 --- a/src/MLIAP/mliap_model_quadratic.cpp +++ b/src/MLIAP/mliap_model_quadratic.cpp @@ -11,6 +11,10 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing author: Aidan Thompson (SNL) +------------------------------------------------------------------------- */ + #include "mliap_model_quadratic.h" #include "pair_mliap.h" #include "mliap_data.h" @@ -22,8 +26,9 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ 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; nonlinearflag = 1; } @@ -52,8 +57,9 @@ int MLIAPModelQuadratic::get_nparams() void MLIAPModelQuadratic::compute_gradients(MLIAPData* data) { + data->energy = 0.0; + for (int ii = 0; ii < data->natoms; ii++) { - const int i = data->iatoms[ii]; const int ielem = data->ielems[ii]; double* coeffi = coeffelem[ielem]; @@ -99,7 +105,8 @@ void MLIAPModelQuadratic::compute_gradients(MLIAPData* data) etmp += coeffi[k++]*bveci*bvecj; } } - data->pairmliap->e_tally(i,etmp); + data->energy += etmp; + data->eatoms[ii] = etmp; } } } diff --git a/src/MLIAP/mliap_model_quadratic.h b/src/MLIAP/mliap_model_quadratic.h index b41df18a07..28f1732f9b 100644 --- a/src/MLIAP/mliap_model_quadratic.h +++ b/src/MLIAP/mliap_model_quadratic.h @@ -18,7 +18,7 @@ namespace LAMMPS_NS { -class MLIAPModelQuadratic : public MLIAPModel { +class MLIAPModelQuadratic : public MLIAPModelSimple { public: MLIAPModelQuadratic(LAMMPS*, char* = nullptr); ~MLIAPModelQuadratic(); diff --git a/src/MLIAP/pair_mliap.cpp b/src/MLIAP/pair_mliap.cpp index a422271f19..0c6f3903f6 100644 --- a/src/MLIAP/pair_mliap.cpp +++ b/src/MLIAP/pair_mliap.cpp @@ -11,12 +11,19 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing author: Aidan Thompson (SNL) +------------------------------------------------------------------------- */ + #include "pair_mliap.h" #include "mliap_data.h" #include "mliap_model_linear.h" #include "mliap_model_quadratic.h" #include "mliap_descriptor_snap.h" +#ifdef MLIAP_PYTHON +#include "mliap_model_python.h" +#endif #include "atom.h" #include "error.h" @@ -27,6 +34,7 @@ #include #include +#include "error.h" using namespace LAMMPS_NS; @@ -65,6 +73,17 @@ PairMLIAP::~PairMLIAP() 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); data->generate_neighdata(list, eflag, vflag); @@ -78,6 +97,8 @@ void PairMLIAP::compute(int eflag, int vflag) model->compute_gradients(data); + e_tally(data); + // calculate force contributions beta_i*dB_i/dR_j descriptor->compute_forces(data); @@ -107,6 +128,7 @@ void PairMLIAP::allocate() void PairMLIAP::settings(int narg, char ** arg) { + if (narg < 4) 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"); model = new MLIAPModelQuadratic(lmp,arg[iarg+2]); 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"); modelflag = 1; } 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->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_atom) eatom[i] += ei; + if (eflag_global) eng_vdwl += data->energy; + if (eflag_atom) + for (int ii = 0; ii < data->natoms; ii++) { + const int i = data->iatoms[ii]; + eatom[i] += data->eatoms[ii]; + } } /* ---------------------------------------------------------------------- diff --git a/src/MLIAP/pair_mliap.h b/src/MLIAP/pair_mliap.h index 61acbbb278..c31634e923 100644 --- a/src/MLIAP/pair_mliap.h +++ b/src/MLIAP/pair_mliap.h @@ -31,7 +31,7 @@ public: virtual void compute(int, int); void settings(int, char **); virtual void coeff(int, char **); - void e_tally(int, double); + void e_tally(class MLIAPData*); void v_tally(int, int, double*, double*); virtual void init_style(); virtual double init_one(int, int); diff --git a/src/PYTHON/python_impl.cpp b/src/PYTHON/python_impl.cpp index 1b3fabfa62..42b1135ceb 100644 --- a/src/PYTHON/python_impl.cpp +++ b/src/PYTHON/python_impl.cpp @@ -26,6 +26,14 @@ #include #include // 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; enum{NONE,INT,DOUBLE,STRING,PTR}; @@ -47,11 +55,17 @@ PythonImpl::PythonImpl(LAMMPS *lmp) : Pointers(lmp) nfunc = 0; pfuncs = nullptr; - // one-time initialization of Python interpreter // pyMain stores pointer to main module 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(); // only needed for Python 2.x and Python 3 < 3.7 diff --git a/src/Purge.list b/src/Purge.list index 0251f923be..2853ce2a7c 100644 --- a/src/Purge.list +++ b/src/Purge.list @@ -49,6 +49,8 @@ packages_ntopo.h # other auto-generated files lmpinstalledpkgs.h lmpgitversion.h +mliap_model_python_couple.cpp +mliap_model_python_couple.h # removed on 9 Sep 2020 mergesort.h # renamed on 8 May 2020