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

Merge pull request #3458 from Boogie3D/mliappy_unified

Browse Source 
MLIAP Unified Interface
This commit is contained in:
Axel Kohlmeyer
2022-09-26 20:33:17 -04:00
committed by GitHub
parent 10068a0fde 1dc3c8fc26
commit be98b5a168
31 changed files with 1760 additions and 138 deletions
Download Patch File Download Diff File Expand all files Collapse all files

20
cmake/Modules/Packages/ML-IAP.cmake
View File

@ -25,16 +25,18 @@ if(MLIAP_ENABLE_PYTHON)
endif()
set(MLIAP_BINARY_DIR ${CMAKE_BINARY_DIR}/cython)
set(MLIAP_CYTHON_SRC ${LAMMPS_SOURCE_DIR}/ML-IAP/mliap_model_python_couple.pyx)
get_filename_component(MLIAP_CYTHON_BASE ${MLIAP_CYTHON_SRC} NAME_WE)
file(GLOB MLIAP_CYTHON_SRC ${LAMMPS_SOURCE_DIR}/ML-IAP/*.pyx)
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...")
foreach(MLIAP_CYTHON_FILE ${MLIAP_CYTHON_SRC})
get_filename_component(MLIAP_CYTHON_BASE ${MLIAP_CYTHON_FILE} NAME_WE)
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_FILE} ${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_FILE}
COMMENT "Generating C++ sources with cythonize...")
target_sources(lammps PRIVATE ${MLIAP_BINARY_DIR}/${MLIAP_CYTHON_BASE}.cpp)
endforeach()
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()

7
doc/src/Library_utility.rst
View File

@ -19,6 +19,7 @@ functions. They do not directly call the LAMMPS library.
- :cpp:func:`lammps_force_timeout`
- :cpp:func:`lammps_has_error`
- :cpp:func:`lammps_get_last_error_message`
- :cpp:func:`lammps_python_api_version`
The :cpp:func:`lammps_free` function is a clean-up function to free
memory that the library had allocated previously via other function
@ -100,3 +101,9 @@ where such memory buffers were allocated that require the use of
.. doxygenfunction:: lammps_get_last_error_message
:project: progguide
-----------------------
.. doxygenfunction:: lammps_python_api_version
:project: progguide

82
doc/src/pair_mliap.rst
View File

@ -10,8 +10,8 @@ Syntax
pair_style mliap ... keyword values ...
* two keyword/value pairs must be appended
* keyword = *model* or *descriptor*
* one or two keyword/value pairs must be appended
* keyword = *model* or *descriptor* or *unified*
.. parsed-literal::
@ -21,6 +21,9 @@ Syntax
*descriptor* values = style filename
style = *sna* or *so3*
filename = name of file containing descriptor definitions
*unified* values = filename ghostneigh_flag
filename = name of file containing serialized unified Python object
ghostneigh_flag = 0/1 to turn off/on inclusion of ghost neighbors in neighbors list
Examples
""""""""
@ -30,20 +33,24 @@ Examples
pair_style mliap model linear InP.mliap.model descriptor sna InP.mliap.descriptor
pair_style mliap model quadratic W.mliap.model descriptor sna W.mliap.descriptor
pair_style mliap model nn Si.nn.mliap.model descriptor so3 Si.nn.mliap.descriptor
pair_style mliap unified mliap_unified_lj_Ar.pkl 0
pair_coeff * * In P
Description
"""""""""""
Pair style *mliap* provides a general interface to families of
machine-learning interatomic potentials. It allows separate definitions
machine-learning interatomic potentials. It allows separate 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. The pair style currently supports only *sna* and *so3*
descriptor styles, but it is is straightforward to add new descriptor
styles.
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. The pair style currently supports only
*sna* and *so3* descriptor styles, but it is straightforward to add new
descriptor styles. By using the *unified* keyword, it is possible to
define a Python model that combines functionalities of both *model* and
*descriptor*.
The SNAP descriptor style *sna* is the same as that used by
:doc:`pair_style snap <pair_snap>`, including the linear, quadratic, and
@ -61,11 +68,11 @@ basis :ref:`(Bartok) <Bartok2013>` and :ref:`(Zagaceta) <Zagaceta2020>`.
The available models are *linear* and *nn*.
The pair_style *mliap* command must be followed by two keywords *model*
and *descriptor* in either order. A single *pair_coeff* command is also
required. The first 2 arguments must be \* \* so as to span all LAMMPS
atom types. This is followed by a list of N arguments that specify the
mapping of MLIAP element names to LAMMPS atom types, where N is the
number of LAMMPS atom types.
and *descriptor* in either order, or the one keyword *unified*. A
single *pair_coeff* command is also required. The first 2 arguments
must be \* \* so as to span all LAMMPS atom types. This is followed by
a list of N arguments 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 the model style. This is followed by
a single argument specifying the model filename containing the
@ -113,14 +120,15 @@ The detail of *nn* module implementation can be found at :ref:`(Yanxon) <Yanxon2
.. admonition:: 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.
When the *model* keyword is *mliappy*, if the filename ends in '.pt',
or '.pth', it will be loaded using pytorch; otherwise, it will be
loaded as 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 two descriptor styles are available: *sna* and *so3*.
@ -167,6 +175,36 @@ to specify the path for these *model* and *descriptor* files.
larger systems since their size depends on the total number of
neighbors per MPI process.
.. versionadded:: TBD
The *unified* keyword is followed by an argument specifying the
filename containing the serialized unified Python object and the "ghostneigh" toggle
(0/1) to disable/enable the construction of neighbors lists including
neighbors of ghost atoms. If the filename ends in '.pt', or '.pth', it will be loaded
using pytorch; otherwise, it will be loaded as a pickle file.
If ghostneigh is enabled, it is recommended to set :doc:`comm_modify <comm_modify>`
cutoff manually, such as in the following example.
.. code-block:: LAMMPS
variable ninteractions equal 2
variable cutdist equal 7.5
variable skin equal 1.0
variable commcut equal (${ninteractions}*${cutdist})+${skin}
neighbor ${skin} bin
comm_modify cutoff ${commcut}
.. note::
To load a model from memory
(i.e. an existing python object), call `lammps.mliap.load_unified(unified)`
from python, and then specify the filename as "EXISTS". 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.
Mixing, shift, table, tail correction, restart, rRESPA info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""

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

@ -1206,6 +1206,7 @@ gewald
Gezelter
Gflop
gfortran
ghostneigh
ghostwhite
Giacomo
gif

54
examples/mliap/README
View File

@ -1,6 +1,6 @@
This directory contains multiple examples of
machine-learning potentials defined using the
MLIAP package in LAMMPS. The input files
ML-IAP package in LAMMPS. The input files
are described below.
in.mliap.snap.Ta06A
@ -52,7 +52,7 @@ You can then run the example as follows
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
If this fails, see the instructions for building the ML-IAP package
with Python support enabled. Also, confirm that the
LAMMPS Python embedded Python interpreter is
working by running ../examples/in.python.
@ -109,4 +109,54 @@ in.mliap.nn.cu
-------------------------
Run a neural network potential for Cu, a combination of SNAP descriptors and the "nn" model style
in.mliap.unified.lj.Ar
-----------------------
This reproduces the output of examples/melt/in.melt,
but using the ``unified`` keyword and
the Python coupling to PyTorch.
This example can be run in two different ways:
1: Running a LAMMPS executable: in.mliap.unified.lj.Ar
You can run the example as follows
`lmp -in in.mliap.unified.lj.Ar -echo both`
The resultant thermo output should be identical to that generated
by in.melt.
The input first runs some python code that creates an ML-IAP unified model
that replicates the ``lj/cut`` pair style with parameters as defined in
examples/melt/in/melt and saves it in the file "mliap_unified_lj_Ar.pkl".
If this fails, see the instructions for building the ML-IAP package
with Python support enabled. Also, confirm that the LAMMPS Python
embedded Python interpreter is working by running ../examples/in.python
and that the LAMMPS python module is installed and working.
2: Running a Python script: mliap_unified_lj_Ar.py
Before testing this, ensure that the previous method
(running a LAMMPS executable) works.
You can run the example in serial:
`python mliap_unified_lj_Ar.py`
or in parallel:
`mpirun -np 4 python mliap_unified_lj_Ar.py`
The resultant log.lammps output should be identical to that generated
by in.melt and in.mliap.unified.lj.Ar.
in.mliap.so3.Ni_Mo
------------------
Example of linear model with SO3 descriptors for Ni/Mo
in.mliap.so3.nn.Si
------------------
Example of NN model with SO3 descriptors for Si

53
examples/mliap/in.mliap.unified.lj.Ar Normal file
View File

@ -0,0 +1,53 @@
# create pickle of unified model
python create_pickle here """
import lammps
import lammps.mliap
from lammps.mliap.mliap_unified_lj import MLIAPUnifiedLJ
def create_pickle():
unified = MLIAPUnifiedLJ(["Ar"])
unified.pickle('mliap_unified_lj_Ar.pkl')
"""
python create_pickle invoke
# 3d Lennard-Jones melt
units lj
atom_style atomic
lattice fcc 0.8442
region box block 0 10 0 10 0 10
create_box 1 box
create_atoms 1 box
mass 1 1.0
velocity all create 3.0 87287 loop geom
pair_style mliap unified mliap_unified_lj_Ar.pkl 0
pair_coeff * * Ar
neighbor 0.3 bin
neigh_modify every 20 delay 0 check no
fix 1 all nve
#dump id all atom 50 dump.melt
#dump 2 all image 25 image.*.jpg type type &
# axes yes 0.8 0.02 view 60 -30
#dump_modify 2 pad 3
#dump 3 all movie 1 movie.mpg type type &
# axes yes 0.8 0.02 view 60 -30
#dump_modify 3 pad 3
#dump 4 all custom 1 forces.xyz fx fy fz
thermo 50
run 250
# clean up
shell rm -f mliap_unified_lj_Ar.pkl

107
examples/mliap/log.23Sep22.mliap.unified.lj.Ar.g++.1 Normal file
View File

@ -0,0 +1,107 @@
LAMMPS (15 Sep 2022)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
using 1 OpenMP thread(s) per MPI task
# create pickle of unified model
python create_pickle here """
import lammps
import lammps.mliap
from lammps.mliap.mliap_unified_lj import MLIAPUnifiedLJ
def create_pickle():
unified = MLIAPUnifiedLJ(["Ar"])
unified.pickle('mliap_unified_lj_Ar.pkl')
"""
python create_pickle invoke
# 3d Lennard-Jones melt
units lj
atom_style atomic
lattice fcc 0.8442
Lattice spacing in x,y,z = 1.6795962 1.6795962 1.6795962
region box block 0 10 0 10 0 10
create_box 1 box
Created orthogonal box = (0 0 0) to (16.795962 16.795962 16.795962)
1 by 1 by 1 MPI processor grid
create_atoms 1 box
Created 4000 atoms
using lattice units in orthogonal box = (0 0 0) to (16.795962 16.795962 16.795962)
create_atoms CPU = 0.001 seconds
mass 1 1.0
velocity all create 3.0 87287 loop geom
pair_style mliap unified mliap_unified_lj_Ar.pkl 0
pair_coeff * * Ar
neighbor 0.3 bin
neigh_modify every 20 delay 0 check no
fix 1 all nve
#dump id all atom 50 dump.melt
#dump 2 all image 25 image.*.jpg type type # axes yes 0.8 0.02 view 60 -30
#dump_modify 2 pad 3
#dump 3 all movie 1 movie.mpg type type # axes yes 0.8 0.02 view 60 -30
#dump_modify 3 pad 3
#dump 4 all custom 1 forces.xyz fx fy fz
thermo 50
run 250
Neighbor list info ...
update: every = 20 steps, delay = 0 steps, check = no
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 2.8
ghost atom cutoff = 2.8
binsize = 1.4, bins = 12 12 12
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) 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) = 12.7 | 12.7 | 12.7 Mbytes
Step Temp E_pair E_mol TotEng Press
0 3 -6.7733681 0 -2.2744931 -3.7033504
50 1.6842865 -4.8082494 0 -2.2824513 5.5666131
100 1.6712577 -4.7875609 0 -2.281301 5.6613913
150 1.6444751 -4.7471034 0 -2.2810074 5.8614211
200 1.6471542 -4.7509053 0 -2.2807916 5.8805431
250 1.6645597 -4.7774327 0 -2.2812174 5.7526089
Loop time of 3.48563 on 1 procs for 250 steps with 4000 atoms
Performance: 30984.374 tau/day, 71.723 timesteps/s
79.5% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 3.3348 | 3.3348 | 3.3348 | 0.0 | 95.67
Neigh | 0.13055 | 0.13055 | 0.13055 | 0.0 | 3.75
Comm | 0.0089845 | 0.0089845 | 0.0089845 | 0.0 | 0.26
Output | 0.00018194 | 0.00018194 | 0.00018194 | 0.0 | 0.01
Modify | 0.0090026 | 0.0090026 | 0.0090026 | 0.0 | 0.26
Other | | 0.002151 | | | 0.06
Nlocal: 4000 ave 4000 max 4000 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 5506 ave 5506 max 5506 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 0 ave 0 max 0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 303576 ave 303576 max 303576 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 303576
Ave neighs/atom = 75.894
Neighbor list builds = 12
Dangerous builds not checked
# clean up
shell rm -f mliap_unified_lj_Ar.pkl
Total wall time: 0:00:03

107
examples/mliap/log.23Sep22.mliap.unified.lj.Ar.g++.4 Normal file
View File

@ -0,0 +1,107 @@
LAMMPS (15 Sep 2022)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
using 1 OpenMP thread(s) per MPI task
# create pickle of unified model
python create_pickle here """
import lammps
import lammps.mliap
from lammps.mliap.mliap_unified_lj import MLIAPUnifiedLJ
def create_pickle():
unified = MLIAPUnifiedLJ(["Ar"])
unified.pickle('mliap_unified_lj_Ar.pkl')
"""
python create_pickle invoke
# 3d Lennard-Jones melt
units lj
atom_style atomic
lattice fcc 0.8442
Lattice spacing in x,y,z = 1.6795962 1.6795962 1.6795962
region box block 0 10 0 10 0 10
create_box 1 box
Created orthogonal box = (0 0 0) to (16.795962 16.795962 16.795962)
1 by 2 by 2 MPI processor grid
create_atoms 1 box
Created 4000 atoms
using lattice units in orthogonal box = (0 0 0) to (16.795962 16.795962 16.795962)
create_atoms CPU = 0.000 seconds
mass 1 1.0
velocity all create 3.0 87287 loop geom
pair_style mliap unified mliap_unified_lj_Ar.pkl 0
pair_coeff * * Ar
neighbor 0.3 bin
neigh_modify every 20 delay 0 check no
fix 1 all nve
#dump id all atom 50 dump.melt
#dump 2 all image 25 image.*.jpg type type # axes yes 0.8 0.02 view 60 -30
#dump_modify 2 pad 3
#dump 3 all movie 1 movie.mpg type type # axes yes 0.8 0.02 view 60 -30
#dump_modify 3 pad 3
#dump 4 all custom 1 forces.xyz fx fy fz
thermo 50
run 250
Neighbor list info ...
update: every = 20 steps, delay = 0 steps, check = no
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 2.8
ghost atom cutoff = 2.8
binsize = 1.4, bins = 12 12 12
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) 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) = 5.774 | 5.774 | 5.774 Mbytes
Step Temp E_pair E_mol TotEng Press
0 3 -6.7733681 0 -2.2744931 -3.7033504
50 1.6842865 -4.8082494 0 -2.2824513 5.5666131
100 1.6712577 -4.7875609 0 -2.281301 5.6613913
150 1.6444751 -4.7471034 0 -2.2810074 5.8614211
200 1.6471542 -4.7509053 0 -2.2807916 5.8805431
250 1.6645597 -4.7774327 0 -2.2812174 5.7526089
Loop time of 1.18059 on 4 procs for 250 steps with 4000 atoms
Performance: 91479.359 tau/day, 211.758 timesteps/s
78.3% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 1.0795 | 1.1007 | 1.1197 | 1.7 | 93.23
Neigh | 0.035618 | 0.03627 | 0.037413 | 0.4 | 3.07
Comm | 0.019611 | 0.038646 | 0.060389 | 9.4 | 3.27
Output | 0.00013034 | 0.0001421 | 0.00017083 | 0.0 | 0.01
Modify | 0.0037233 | 0.0037641 | 0.0038121 | 0.1 | 0.32
Other | | 0.001086 | | | 0.09
Nlocal: 1000 ave 1008 max 987 min
Histogram: 1 0 0 0 0 0 1 0 1 1
Nghost: 2711.25 ave 2728 max 2693 min
Histogram: 1 0 0 0 0 2 0 0 0 1
Neighs: 0 ave 0 max 0 min
Histogram: 4 0 0 0 0 0 0 0 0 0
FullNghs: 75894 ave 77414 max 74214 min
Histogram: 1 0 0 1 0 0 0 1 0 1
Total # of neighbors = 303576
Ave neighs/atom = 75.894
Neighbor list builds = 12
Dangerous builds not checked
# clean up
shell rm -f mliap_unified_lj_Ar.pkl
Total wall time: 0:00:01

101
examples/mliap/mliap_numpy_Ta06A.py Normal file
View File

@ -0,0 +1,101 @@
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 numpy as np
class LinearModel():
def __init__(self,file):
coeffs = np.genfromtxt(file,skip_header=6)
self.bias = coeffs[0]
self.weights = coeffs[1:]
self.n_params = len(coeffs)
self.n_descriptors = len(self.weights)
self.n_elements = 1
def __call__(self,elems,bispectrum,beta,energy):
energy[:] = bispectrum @ self.weights + self.bias
beta[:] = self.weights
mymodel = LinearModel("Ta06A.mliap.model")
import lammps
lmp = lammps.lammps(cmdargs=['-echo','both'])
import lammps.mliap
lammps.mliap.activate_mliappy(lmp)
lmp.commands_string(before_loading)
lammps.mliap.load_model(mymodel)
lmp.commands_string(after_loading)

65
examples/mliap/mliap_unified_lj_Ar.py Normal file
View File

@ -0,0 +1,65 @@
# Demonstrate how to load a unified model from python.
# This is essentially the same as in.mliap.unified.lj.Ar
# except that python is the driving program, and lammps
# is in library mode.
before_loading =\
"""# 3d Lennard-Jones melt
units lj
atom_style atomic
lattice fcc 0.8442
region box block 0 10 0 10 0 10
create_box 1 box
create_atoms 1 box
mass 1 1.0
velocity all create 3.0 87287 loop geom
"""
after_loading =\
"""
pair_style mliap unified EXISTS
pair_coeff * * Ar
neighbor 0.3 bin
neigh_modify every 20 delay 0 check no
fix 1 all nve
thermo 50
run 250
"""
import lammps
lmp = lammps.lammps(cmdargs=['-echo','both'])
# Before defining the pair style, you must do the following:
import lammps.mliap
lammps.mliap.activate_mliappy(lmp)
# Otherwise, when running lammps in library mode,
# you will get an error
# Setup the simulation to just before declaring
# the mliap unified pair style
lmp.commands_string(before_loading)
# Define the model however you like. In this example
# we simply import the unified L-J example from mliap
from lammps.mliap.mliap_unified_lj import MLIAPUnifiedLJ
unified = MLIAPUnifiedLJ(["Ar"])
# You can also load the model from a pickle file.
# import pickle
# with open('mliap_unified_lj_Ar.pkl', 'rb') as pfile:
# unified = pickle.load(pfile)
# Connect the L-J model to the mliap unified pair style.
lammps.mliap.load_unified(unified)
# Run the simulation with the mliap unified pair style
# Use pre-loaded model by specifying model filename as "EXISTS"
lmp.commands_string(after_loading)

9
python/lammps/mliap/__init__.py
View File

@ -5,6 +5,7 @@
import sysconfig
import ctypes
import platform
import warnings
py_ver = sysconfig.get_config_vars('VERSION')[0]
OS_name = platform.system()
@ -25,8 +26,10 @@ except Exception as e:
raise OSError("Unable to locate python shared library") from e
if not pylib.Py_IsInitialized():
raise RuntimeError("This interpreter is not compatible with python-based mliap for LAMMPS.")
warnings.warn("This interpreter is not compatible with python-based MLIAP for LAMMPS. "
"Attempting to activate the MLIAP-python coupling from python may result "
"in undefined behavior.")
else:
from .loader import load_model, load_unified, activate_mliappy
del sysconfig, ctypes, library, pylib
from .loader import load_model, activate_mliappy

67
python/lammps/mliap/loader.py
View File

@ -17,27 +17,58 @@
import sys
import importlib
import importlib.util
import importlib.machinery
import importlib.abc
from ctypes import pythonapi, c_int, c_void_p, py_object
# This dynamic loader imports a python module embedded in a shared library.
# The default value of api_version is 1013 because it has been stable since 2006.
class DynamicLoader(importlib.abc.Loader):
def __init__(self,module_name,library,api_version=1013):
self.api_version = api_version
attr = "PyInit_"+module_name
initfunc = getattr(library,attr)
# c_void_p is standin for PyModuleDef *
initfunc.restype = c_void_p
initfunc.argtypes = ()
self.module_def = initfunc()
def create_module(self, spec):
createfunc = pythonapi.PyModule_FromDefAndSpec2
# c_void_p is standin for PyModuleDef *
createfunc.argtypes = c_void_p, py_object, c_int
createfunc.restype = py_object
module = createfunc(self.module_def, spec, self.api_version)
return module
def exec_module(self, module):
execfunc = pythonapi.PyModule_ExecDef
# c_void_p is standin for PyModuleDef *
execfunc.argtypes = py_object, c_void_p
execfunc.restype = c_int
result = execfunc(module, self.module_def)
if result<0:
raise ImportError()
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
library = lmp.lib
module_names = ["mliap_model_python_couple", "mliap_unified_couple"]
api_version = library.lammps_python_api_version()
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
for module_name in module_names:
# Make Machinery
loader = DynamicLoader(module_name,library,api_version)
spec = importlib.util.spec_from_loader(module_name,loader)
# Do the import
module = importlib.util.module_from_spec(spec)
sys.modules[module_name] = module
spec.loader.exec_module(module)
except Exception as ee:
raise ImportError("Could not load ML-IAP python coupling module.") from ee
@ -50,3 +81,11 @@ def load_model(model):
) from ie
mliap_model_python_couple.load_from_python(model)
def load_unified(model):
try:
import mliap_unified_couple
except ImportError as ie:
raise ImportError("ML-IAP python module must be activated before loading\n"
"the pair style. Call lammps.mliap.activate_mliappy(lmp)."
) from ie
mliap_unified_couple.load_from_python(model)

28
python/lammps/mliap/mliap_unified_abc.py Normal file
View File

@ -0,0 +1,28 @@
from abc import ABC, abstractmethod
import pickle
class MLIAPUnified(ABC):
"""Abstract base class for MLIAPUnified."""
def __init__(self):
self.interface = None
self.element_types = None
self.ndescriptors = None
self.nparams = None
self.rcutfac = None
@abstractmethod
def compute_gradients(self, data):
"""Compute gradients."""
@abstractmethod
def compute_descriptors(self, data):
"""Compute descriptors."""
@abstractmethod
def compute_forces(self, data):
"""Compute forces."""
def pickle(self, fname):
with open(fname, 'wb') as fp:
pickle.dump(self, fp)

44
python/lammps/mliap/mliap_unified_lj.py Normal file
View File

@ -0,0 +1,44 @@
from .mliap_unified_abc import MLIAPUnified
import numpy as np
class MLIAPUnifiedLJ(MLIAPUnified):
"""Test implementation for MLIAPUnified."""
def __init__(self, element_types, epsilon=1.0, sigma=1.0, rcutfac=1.25):
super().__init__()
self.element_types = element_types
self.ndescriptors = 1
self.nparams = 3
# Mimicking the LJ pair-style:
# pair_style lj/cut 2.5
# pair_coeff * * 1 1
self.epsilon = epsilon
self.sigma = sigma
self.rcutfac = rcutfac
def compute_gradients(self, data):
"""Test compute_gradients."""
def compute_descriptors(self, data):
"""Test compute_descriptors."""
def compute_forces(self, data):
"""Test compute_forces."""
eij, fij = self.compute_pair_ef(data)
data.update_pair_energy(eij)
data.update_pair_forces(fij)
def compute_pair_ef(self, data):
rij = data.rij
r2inv = 1.0 / np.sum(rij ** 2, axis=1)
r6inv = r2inv * r2inv * r2inv
lj1 = 4.0 * self.epsilon * self.sigma**12
lj2 = 4.0 * self.epsilon * self.sigma**6
eij = r6inv * (lj1 * r6inv - lj2)
fij = r6inv * (3.0 * lj2 - 6.0 * lj2 * r6inv) * r2inv
fij = fij[:, np.newaxis] * rij
return eij, fij

14
src/ML-IAP/Install.sh
View File

@ -42,6 +42,7 @@ done
# Install cython pyx file only if also Python is available
action mliap_model_python_couple.pyx python_impl.cpp
action mliap_unified_couple.pyx python_impl.cpp
# edit 2 Makefile.package files to include/exclude package info
@ -57,7 +58,7 @@ if (test $1 = 1) then
include ..\/..\/lib\/python\/Makefile.mliap_python
' ../Makefile.package.settings
fi
cythonize -3 ../mliap_model_python_couple.pyx
cythonize -3 ../mliap_model_python_couple.pyx ../mliap_unified_couple.pyx
fi
elif (test $1 = 0) then
@ -65,7 +66,8 @@ 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
rm -f ../mliap_model_python_couple.cpp ../mliap_model_python_couple.h \
../mliap_unified_couple.cpp ../mliap_unified_couple.h
sed -i -e '/^include.*python.*mliap_python.*$/d' ../Makefile.package.settings
elif (test $1 = 2) then
@ -73,7 +75,8 @@ elif (test $1 = 2) 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
rm -f ../mliap_model_python_couple.cpp ../mliap_model_python_couple.h \
../mliap_unified_couple.cpp ../mliap_unified_couple.h
sed -i -e '/^include.*python.*mliap_python.*$/d' ../Makefile.package.settings
if (test -e ../Makefile.package) then
sed -i -e 's|^PKG_INC =[ \t]*|&-DMLIAP_PYTHON |' ../Makefile.package
@ -85,8 +88,9 @@ elif (test $1 = 2) then
include ..\/..\/lib\/python\/Makefile.mliap_python
' ../Makefile.package.settings
fi
cythonize -3 ../mliap_model_python_couple.pyx
cythonize -3 ../mliap_model_python_couple.pyx ../mliap_unified_couple.pyx
else
rm -f ../mliap_model_python_couple.cpp ../mliap_model_python_couple.h
rm -f ../mliap_model_python_couple.cpp ../mliap_model_python_couple.h \
../mliap_unified_couple.cpp ../mliap_unified_couple.h
fi
fi

49
src/ML-IAP/mliap_data.cpp
View File

@ -31,11 +31,12 @@ 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),
Pointers(lmp), f(nullptr), 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)
numneighs(nullptr), iatoms(nullptr), pair_i(nullptr), ielems(nullptr),
jatoms(nullptr), jelems(nullptr), elems(nullptr), rij(nullptr),
graddesc(nullptr), model(nullptr), descriptor(nullptr), list(nullptr)
{
gradgradflag = gradgradflag_in;
map = map_in;
@ -85,10 +86,12 @@ MLIAPData::~MLIAPData()
memory->destroy(gradforce);
memory->destroy(iatoms);
memory->destroy(pair_i);
memory->destroy(ielems);
memory->destroy(numneighs);
memory->destroy(jatoms);
memory->destroy(jelems);
memory->destroy(elems);
memory->destroy(rij);
memory->destroy(graddesc);
}
@ -107,6 +110,7 @@ void MLIAPData::init()
void MLIAPData::generate_neighdata(NeighList* list_in, int eflag_in, int vflag_in)
{
list = list_in;
f = atom->f;
double **x = atom->x;
int *type = atom->type;
@ -114,21 +118,25 @@ void MLIAPData::generate_neighdata(NeighList* list_in, int eflag_in, int vflag_i
int *numneigh = list->numneigh;
int **firstneigh = list->firstneigh;
// grow nmax gradforce array if necessary
// grow nmax gradforce, elems arrays if necessary
if (atom->nmax > nmax) {
nmax = atom->nmax;
memory->grow(gradforce,nmax,size_gradforce,
"MLIAPData:gradforce");
memory->grow(elems,nmax,"MLIAPData:elems");
}
// clear gradforce array
// clear gradforce, elems arrays
int nall = atom->nlocal + atom->nghost;
for (int i = 0; i < nall; i++)
ntotal = nall;
for (int i = 0; i < nall; i++) {
for (int j = 0; j < size_gradforce; j++) {
gradforce[i][j] = 0.0;
}
elems[i] = 0;
}
// grow arrays if necessary
@ -153,8 +161,9 @@ void MLIAPData::generate_neighdata(NeighList* list_in, int eflag_in, int vflag_i
grow_neigharrays();
npairs = 0;
int ij = 0;
for (int ii = 0; ii < nlistatoms; ii++) {
for (int ii = 0; ii < natomneigh; ii++) {
const int i = ilist[ii];
const double xtmp = x[i][0];
@ -178,6 +187,7 @@ void MLIAPData::generate_neighdata(NeighList* list_in, int eflag_in, int vflag_i
const int jelem = map[jtype];
if (rsq < descriptor->cutsq[ielem][jelem]) {
pair_i[ij] = i;
jatoms[ij] = j;
jelems[ij] = jelem;
rij[ij][0] = delx;
@ -190,6 +200,12 @@ void MLIAPData::generate_neighdata(NeighList* list_in, int eflag_in, int vflag_i
iatoms[ii] = i;
ielems[ii] = ielem;
numneighs[ii] = ninside;
npairs += ninside;
}
for (int i = 0; i < nall; i++) {
const int itype = type[i];
elems[i] = map[itype];
}
eflag = eflag_in;
@ -204,12 +220,14 @@ void MLIAPData::grow_neigharrays()
{
// grow neighbor atom arrays if necessary
if (natomneigh_max < nlistatoms) {
memory->grow(iatoms,nlistatoms,"MLIAPData:iatoms");
memory->grow(ielems,nlistatoms,"MLIAPData:ielems");
memory->grow(numneighs,nlistatoms,"MLIAPData:numneighs");
natomneigh_max = nlistatoms;
natomneigh = list->inum;
if (list->ghost == 1)
natomneigh += list->gnum;
if (natomneigh_max < natomneigh) {
memory->grow(iatoms,natomneigh,"MLIAPData:iatoms");
memory->grow(ielems,natomneigh,"MLIAPData:ielems");
memory->grow(numneighs,natomneigh,"MLIAPData:numneighs");
natomneigh_max = natomneigh;
}
// grow neighbor arrays if necessary
@ -221,7 +239,7 @@ void MLIAPData::grow_neigharrays()
int *type = atom->type;
int nneigh = 0;
for (int ii = 0; ii < nlistatoms; ii++) {
for (int ii = 0; ii < natomneigh; ii++) {
const int i = ilist[ii];
const double xtmp = x[i][0];
@ -249,6 +267,7 @@ void MLIAPData::grow_neigharrays()
}
if (nneigh_max < nneigh) {
memory->grow(pair_i,nneigh,"MLIAPData:pair_i");
memory->grow(jatoms,nneigh,"MLIAPData:jatoms");
memory->grow(jelems,nneigh,"MLIAPData:jelems");
memory->grow(rij,nneigh,3,"MLIAPData:rij");
@ -264,6 +283,7 @@ double MLIAPData::memory_usage()
bytes += (double)nelements*nparams*sizeof(double); // egradient
bytes += (double)nmax*size_gradforce*sizeof(double); // gradforce
bytes += (double)nmax*sizeof(int); // elems
if (gradgradflag == 1) {
bytes += (double)natomgamma_max*
@ -282,6 +302,7 @@ double MLIAPData::memory_usage()
bytes += (double)natomneigh_max*sizeof(int); // ielems
bytes += (double)natomneigh_max*sizeof(int); // numneighs
bytes += (double)nneigh_max*sizeof(int); // pair_i
bytes += (double)nneigh_max*sizeof(int); // jatoms
bytes += (double)nneigh_max*sizeof(int); // jelems
bytes += (double)nneigh_max*3*sizeof(double); // rij"

8
src/ML-IAP/mliap_data.h
View File

@ -35,6 +35,7 @@ class MLIAPData : protected Pointers {
int size_gradforce;
int yoffset, zoffset;
int ndims_force, ndims_virial;
double **f;
double **gradforce;
double **betas; // betas for all atoms in list
double **descriptors; // descriptors for all atoms in list
@ -56,15 +57,20 @@ class MLIAPData : protected Pointers {
// data structures for mliap neighbor list
// only neighbors strictly inside descriptor cutoff
int nlistatoms; // current number of atoms in neighborlist
int ntotal; // total number of owned and ghost atoms on this proc
int nlistatoms; // current number of atoms in local atom lists
int nlistatoms_max; // allocated size of descriptor array
int natomneigh; // current number of atoms and ghosts in atom neighbor arrays
int natomneigh_max; // allocated size of atom neighbor arrays
int *numneighs; // neighbors count for each atom
int *iatoms; // index of each atom
int *ielems; // element of each atom
int nneigh_max; // number of ij neighbors allocated
int npairs; // number of ij neighbor pairs
int *pair_i; // index of each i atom for each ij pair
int *jatoms; // index of each neighbor
int *jelems; // element of each neighbor
int *elems; // element of each atom in or not in the neighborlist
double **rij; // distance vector of each neighbor
double ***graddesc; // descriptor gradient w.r.t. each neighbor
int eflag; // indicates if energy is needed

3
src/ML-IAP/mliap_descriptor.cpp
View File

@ -25,7 +25,7 @@ using namespace LAMMPS_NS;
MLIAPDescriptor::MLIAPDescriptor(LAMMPS *lmp) :
Pointers(lmp), ndescriptors(0), nelements(0), elements(nullptr), cutsq(nullptr),
radelem(nullptr), wjelem(nullptr)
cutghost(nullptr), radelem(nullptr), wjelem(nullptr)
{
cutmax = 0.0;
}
@ -37,6 +37,7 @@ MLIAPDescriptor::~MLIAPDescriptor()
for (int i = 0; i < nelements; i++) delete[] elements[i];
delete[] elements;
memory->destroy(cutsq);
memory->destroy(cutghost);
memory->destroy(radelem);
memory->destroy(wjelem);
}

15
src/ML-IAP/mliap_descriptor.h
View File

@ -29,13 +29,14 @@ class MLIAPDescriptor : protected Pointers {
virtual void init() = 0;
virtual double memory_usage();
int ndescriptors; // number of descriptors
int nelements; // # of unique elements
char **elements; // names of unique elements
double **cutsq; // nelem x nelem rcutsq values
double cutmax; // maximum cutoff needed
double *radelem; // element radii
double *wjelem; // elements weights
int ndescriptors; // number of descriptors
int nelements; // # of unique elements
char **elements; // names of unique elements
double **cutsq; // nelem x nelem rcutsq values
double **cutghost; // cutoff for each ghost pair
double cutmax; // maximum cutoff needed
double *radelem; // element radii
double *wjelem; // elements weights
protected:
};

16
src/ML-IAP/mliap_model_python.h
View File

@ -21,20 +21,20 @@ namespace LAMMPS_NS {
class MLIAPModelPython : public MLIAPModel {
public:
MLIAPModelPython(LAMMPS *, char * = nullptr);
~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();
~MLIAPModelPython() override;
int get_nparams() override;
int get_gamma_nnz(class MLIAPData *) override;
void compute_gradients(class MLIAPData *) override;
void compute_gradgrads(class MLIAPData *) override;
void compute_force_gradients(class MLIAPData *) override;
double memory_usage() override;
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 *);
void read_coeffs(char *) override;
private:
};

291
src/ML-IAP/mliap_unified.cpp Normal file
View File

@ -0,0 +1,291 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, 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: Steven Ray Anaya (LANL)
------------------------------------------------------------------------- */
#ifdef MLIAP_PYTHON
#include "mliap_unified.h"
#include <Python.h>
#include "error.h"
#include "lmppython.h"
#include "memory.h"
#include "mliap_data.h"
#include "mliap_unified_couple.h"
#include "pair_mliap.h"
#include "python_compat.h"
#include "utils.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
MLIAPDummyDescriptor::MLIAPDummyDescriptor(LAMMPS *_lmp) : MLIAPDescriptor(_lmp) {}
MLIAPDummyDescriptor::~MLIAPDummyDescriptor()
{
// manually decrement borrowed reference from Python
Py_DECREF(unified_interface);
}
/* ----------------------------------------------------------------------
invoke compute_descriptors from Cython interface
---------------------------------------------------------------------- */
void MLIAPDummyDescriptor::compute_descriptors(class MLIAPData *data)
{
PyGILState_STATE gstate = PyGILState_Ensure();
compute_descriptors_python(unified_interface, data);
if (PyErr_Occurred()) {
PyErr_Print();
PyErr_Clear();
PyGILState_Release(gstate);
lmp->error->all(FLERR, "Running mliappy unified compute_descriptors failure.");
}
PyGILState_Release(gstate);
}
/* ----------------------------------------------------------------------
invoke compute_forces from Cython interface
---------------------------------------------------------------------- */
void MLIAPDummyDescriptor::compute_forces(class MLIAPData *data)
{
PyGILState_STATE gstate = PyGILState_Ensure();
compute_forces_python(unified_interface, data);
if (PyErr_Occurred()) {
PyErr_Print();
PyErr_Clear();
PyGILState_Release(gstate);
lmp->error->all(FLERR, "Running mliappy unified compute_forces failure.");
}
PyGILState_Release(gstate);
}
// not implemented
void MLIAPDummyDescriptor::compute_force_gradients(class MLIAPData *)
{
error->all(FLERR, "compute_force_gradients not implemented");
}
// not implemented
void MLIAPDummyDescriptor::compute_descriptor_gradients(class MLIAPData *)
{
error->all(FLERR, "compute_descriptor_gradients not implemented");
}
void MLIAPDummyDescriptor::init()
{
memory->create(radelem, nelements, "mliap_dummy_descriptor:radelem");
for (int ielem = 0; ielem < nelements; ielem++) { radelem[ielem] = 1; }
double cut;
cutmax = 0.0;
memory->create(cutsq, nelements, nelements, "mliap/descriptor/dummy:cutsq");
memory->create(cutghost, nelements, nelements, "mliap/descriptor/dummy:cutghost");
for (int ielem = 0; ielem < nelements; ielem++) {
// rcutfac set from python, is global cutoff for all elements
cut = 2.0 * radelem[ielem] * rcutfac;
if (cut > cutmax) cutmax = cut;
cutsq[ielem][ielem] = cut * cut;
cutghost[ielem][ielem] = cut * cut;
for (int jelem = ielem + 1; jelem < nelements; jelem++) {
cut = (radelem[ielem] + radelem[jelem]) * rcutfac;
cutsq[ielem][jelem] = cutsq[jelem][ielem] = cut * cut;
cutghost[ielem][jelem] = cutghost[jelem][ielem] = cut * cut;
}
}
}
void MLIAPDummyDescriptor::set_elements(char **elems, int nelems)
{
nelements = nelems;
elements = new char *[nelems];
for (int i = 0; i < nelems; i++) { elements[i] = utils::strdup(elems[i]); }
}
/* ---------------------------------------------------------------------- */
MLIAPDummyModel::MLIAPDummyModel(LAMMPS *lmp, char *coefffilename) : MLIAPModel(lmp, coefffilename)
{
nonlinearflag = 1;
}
MLIAPDummyModel::~MLIAPDummyModel()
{
// manually decrement borrowed reference from Python
Py_DECREF(unified_interface);
}
int MLIAPDummyModel::get_nparams()
{
return nparams;
}
int MLIAPDummyModel::get_gamma_nnz(class MLIAPData *)
{
// TODO: get_gamma_nnz
return 0;
}
/* ----------------------------------------------------------------------
invoke compute_gradients from Cython interface
---------------------------------------------------------------------- */
void MLIAPDummyModel::compute_gradients(class MLIAPData *data)
{
PyGILState_STATE gstate = PyGILState_Ensure();
compute_gradients_python(unified_interface, data);
if (PyErr_Occurred()) {
PyErr_Print();
PyErr_Clear();
PyGILState_Release(gstate);
lmp->error->all(FLERR, "Running mliappy unified compute_gradients failure.");
}
PyGILState_Release(gstate);
}
// not implemented
void MLIAPDummyModel::compute_gradgrads(class MLIAPData *)
{
error->all(FLERR, "compute_gradgrads not implemented");
}
// not implemented
void MLIAPDummyModel::compute_force_gradients(class MLIAPData *)
{
error->all(FLERR, "compute_force_gradients not implemented");
}
/* ----------------------------------------------------------------------
memory usage unclear due to Cython/Python implementation
---------------------------------------------------------------------- */
double MLIAPDummyModel::memory_usage()
{
// TODO: implement memory usage in Cython(?)
return 0;
}
// not implemented
void MLIAPDummyModel::read_coeffs(char *)
{
error->all(FLERR, "read_coeffs not implemented");
}
/* ----------------------------------------------------------------------
build the unified interface object, connect to dummy model and descriptor
---------------------------------------------------------------------- */
MLIAPBuildUnified_t LAMMPS_NS::build_unified(char *unified_fname, MLIAPData *data, LAMMPS *lmp,
char *coefffilename)
{
lmp->python->init();
PyGILState_STATE gstate = PyGILState_Ensure();
PyObject *pyMain = PyImport_AddModule("__main__");
if (!pyMain) {
PyGILState_Release(gstate);
lmp->error->all(FLERR, "Could not initialize embedded Python");
}
PyObject *unified_module = PyImport_ImportModule("mliap_unified_couple");
if (PyErr_Occurred()) {
PyErr_Print();
PyErr_Clear();
PyGILState_Release(gstate);
lmp->error->all(FLERR, "Loading mliappy unified module failure.");
}
// Connect dummy model, dummy descriptor, data to Python unified
MLIAPDummyModel *model = new MLIAPDummyModel(lmp, coefffilename);
MLIAPDummyDescriptor *descriptor = new MLIAPDummyDescriptor(lmp);
PyObject *unified_interface = mliap_unified_connect(unified_fname, model, descriptor);
if (PyErr_Occurred()) {
PyErr_Print();
PyErr_Clear();
PyGILState_Release(gstate);
lmp->error->all(FLERR, "Running mliappy unified module failure.");
}
// Borrowed references must be manually incremented
model->unified_interface = unified_interface;
Py_INCREF(unified_interface);
descriptor->unified_interface = unified_interface;
Py_INCREF(unified_interface);
PyGILState_Release(gstate);
MLIAPBuildUnified_t build = {data, descriptor, model};
return build;
}
/* ---------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
set energy for ij atom pairs
---------------------------------------------------------------------- */
void LAMMPS_NS::update_pair_energy(MLIAPData *data, double *eij)
{
double e_total = 0.0;
const auto nlistatoms = data->nlistatoms;
for (int ii = 0; ii < nlistatoms; ii++) data->eatoms[ii] = 0;
for (int ii = 0; ii < data->npairs; ii++) {
int i = data->pair_i[ii];
double e = 0.5 * eij[ii];
// must not count any contribution where i is not a local atom
if (i < nlistatoms) {
data->eatoms[i] += e;
e_total += e;
}
}
data->energy = e_total;
}
/* ----------------------------------------------------------------------
set forces for ij atom pairs
---------------------------------------------------------------------- */
void LAMMPS_NS::update_pair_forces(MLIAPData *data, double *fij)
{
const auto nlistatoms = data->nlistatoms;
double **f = data->f;
for (int ii = 0; ii < data->npairs; ii++) {
int ii3 = ii * 3;
int i = data->pair_i[ii];
int j = data->jatoms[ii];
// must not count any contribution where i is not a local atom
if (i < nlistatoms) {
f[i][0] += fij[ii3];
f[i][1] += fij[ii3 + 1];
f[i][2] += fij[ii3 + 2];
f[j][0] -= fij[ii3];
f[j][1] -= fij[ii3 + 1];
f[j][2] -= fij[ii3 + 2];
if (data->vflag) data->pairmliap->v_tally(i, j, &fij[ii3], data->rij[ii]);
}
}
}
#endif

69
src/ML-IAP/mliap_unified.h Normal file
View File

@ -0,0 +1,69 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, 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_UNIFIED_H
#define LMP_MLIAP_UNIFIED_H
#include "mliap_data.h"
#include "mliap_descriptor.h"
#include "mliap_model.h"
#include <Python.h>
namespace LAMMPS_NS {
class MLIAPDummyDescriptor : public MLIAPDescriptor {
public:
MLIAPDummyDescriptor(LAMMPS *);
~MLIAPDummyDescriptor() override;
void compute_descriptors(class MLIAPData *) override;
void compute_forces(class MLIAPData *) override;
void compute_force_gradients(class MLIAPData *) override;
void compute_descriptor_gradients(class MLIAPData *) override;
void init() override;
void set_elements(char **, int);
PyObject *unified_interface; // MLIAPUnifiedInterface
double rcutfac;
};
class MLIAPDummyModel : public MLIAPModel {
public:
MLIAPDummyModel(LAMMPS *, char * = nullptr);
~MLIAPDummyModel() override;
int get_nparams() override;
int get_gamma_nnz(class MLIAPData *) override;
void compute_gradients(class MLIAPData *) override;
void compute_gradgrads(class MLIAPData *) override;
void compute_force_gradients(class MLIAPData *) override;
double memory_usage() override;
PyObject *unified_interface;
protected:
void read_coeffs(char *) override;
};
struct MLIAPBuildUnified_t {
MLIAPData *data;
MLIAPDummyDescriptor *descriptor;
MLIAPDummyModel *model;
};
MLIAPBuildUnified_t build_unified(char *, MLIAPData *, LAMMPS *, char * = NULL);
void update_pair_energy(MLIAPData *, double *);
void update_pair_forces(MLIAPData *, double *);
} // namespace LAMMPS_NS
#endif

393
src/ML-IAP/mliap_unified_couple.pyx Normal file
View File

@ -0,0 +1,393 @@
# cython: language_level=3
# distutils: language = c++
import pickle
import numpy as np
import lammps.mliap
cimport cython
from cpython.ref cimport PyObject
from libc.stdlib cimport malloc, free
cdef extern from "lammps.h" namespace "LAMMPS_NS":
cdef cppclass LAMMPS:
pass
cdef extern from "mliap_data.h" namespace "LAMMPS_NS":
cdef cppclass MLIAPData:
# ----- may not need -----
int size_array_rows
int size_array_cols
int natoms
int yoffset
int zoffset
int ndims_force
int ndims_virial
# -END- may not need -END-
int size_gradforce
# ----- write only -----
double ** f
double ** gradforce
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
# -END- write only -END-
int ndescriptors # number of descriptors
int nparams # number of model parameters per element
int nelements # number of elements
# data structures for grad-grad list (gamma)
# ----- ignore for now -----
int gamma_nnz # number of non-zero entries in gamma
double ** gamma # gamma element
int ** gamma_row_index # row (parameter) index
int ** gamma_col_index # column (descriptor) index
double * egradient # energy gradient w.r.t. parameters
# -END- ignore for now -END-
# data structures for mliap neighbor list
# only neighbors strictly inside descriptor cutoff
int ntotal # total number of owned and ghost atoms on this proc
int nlistatoms # current number of atoms in local atom lists
int natomneigh # current number of atoms and ghosts in atom neighbor arrays
int * numneighs # neighbors count for each atom
int * iatoms # index of each atom
int * pair_i # index of each i atom for each ij pair
int * ielems # element of each atom
int nneigh_max # number of ij neighbors allocated
int npairs # number of ij neighbor pairs
int * jatoms # index of each neighbor
int * jelems # element of each neighbor
int * elems # element of each atom in or not in the neighborlist
double ** rij # distance vector of each neighbor
# ----- write only -----
double *** graddesc # descriptor gradient w.r.t. each neighbor
# -END- write only -END-
int eflag # indicates if energy is needed
int vflag # indicates if virial is needed
cdef extern from "mliap_unified.h" namespace "LAMMPS_NS":
cdef cppclass MLIAPDummyDescriptor:
MLIAPDummyDescriptor(PyObject *, LAMMPS *) except +
int ndescriptors # number of descriptors
int nelements # # of unique elements
char **elements # names of unique elements
double cutmax # maximum cutoff needed
double rcutfac
double *radelem # element radii
void compute_descriptors(MLIAPData *)
void compute_forces(MLIAPData *)
void set_elements(char **, int)
cdef cppclass MLIAPDummyModel:
MLIAPDummyModel(PyObject *, LAMMPS *, char * = NULL) except +
int ndescriptors # number of descriptors
int nparams # number of parameters per element
int nelements; # # of unique elements
void compute_gradients(MLIAPData *)
cdef void update_pair_energy(MLIAPData *, double *) except +
cdef void update_pair_forces(MLIAPData *, double *) except +
LOADED_MODEL = None
# @property sans getter
def write_only_property(fset):
return property(fget=None, fset=fset)
# Cython implementation of MLIAPData
# Automatically converts between C arrays and numpy when needed
cdef class MLIAPDataPy:
cdef MLIAPData * data
def __cinit__(self):
self.data = NULL
def update_pair_energy(self, eij):
cdef double[:] eij_arr = eij
update_pair_energy(self.data, &eij_arr[0])
def update_pair_forces(self, fij):
cdef double[:, ::1] fij_arr = fij
update_pair_forces(self.data, &fij_arr[0][0])
@property
def f(self):
if self.data.f is NULL:
return None
return np.asarray(<double[:self.ntotal, :3]> &self.data.f[0][0])
@property
def size_gradforce(self):
return self.data.size_gradforce
@write_only_property
def gradforce(self, value):
if self.data.gradforce is NULL:
raise ValueError("attempt to set NULL gradforce")
cdef double[:, :] gradforce_view = <double[:self.ntotal, :self.size_gradforce]> &self.data.gradforce[0][0]
cdef double[:, :] value_view = value
gradforce_view[:] = value_view
@write_only_property
def betas(self, value):
if self.data.betas is NULL:
raise ValueError("attempt to set NULL betas")
cdef double[:, :] betas_view = <double[:self.nlistatoms, :self.ndescriptors]> &self.data.betas[0][0]
cdef double[:, :] value_view = value
betas_view[:] = value_view
@write_only_property
def descriptors(self, value):
if self.data.descriptors is NULL:
raise ValueError("attempt to set NULL descriptors")
cdef double[:, :] descriptors_view = <double[:self.nlistatoms, :self.ndescriptors]> &self.data.descriptors[0][0]
cdef double[:, :] value_view = value
descriptors_view[:] = value_view
@write_only_property
def eatoms(self, value):
if self.data.eatoms is NULL:
raise ValueError("attempt to set NULL eatoms")
cdef double[:] eatoms_view = <double[:self.nlistatoms]> &self.data.eatoms[0]
cdef double[:] value_view = value
eatoms_view[:] = value_view
@write_only_property
def energy(self, value):
self.data.energy = <double>value
@property
def ndescriptors(self):
return self.data.ndescriptors
@property
def nparams(self):
return self.data.nparams
@property
def nelements(self):
return self.data.nelements
# data structures for grad-grad list (gamma)
@property
def gamma_nnz(self):
return self.data.gamma_nnz
@property
def gamma(self):
if self.data.gamma is NULL:
return None
return np.asarray(<double[:self.nlistatoms, :self.gama_nnz]> &self.data.gamma[0][0])
@property
def gamma_row_index(self):
if self.data.gamma_row_index is NULL:
return None
return np.asarray(<int[:self.nlistatoms, :self.gamma_nnz]> &self.data.gamma_row_index[0][0])
@property
def gamma_col_index(self):
if self.data.gamma_col_index is NULL:
return None
return np.asarray(<int[:self.nlistatoms, :self.gamma_nnz]> &self.data.gamma_col_index[0][0])
@property
def egradient(self):
if self.data.egradient is NULL:
return None
return np.asarray(<double[:self.nelements*self.nparams]> &self.data.egradient[0])
# data structures for mliap neighbor list
# only neighbors strictly inside descriptor cutoff
@property
def ntotal(self):
return self.data.ntotal
@property
def elems(self):
if self.data.elems is NULL:
return None
return np.asarray(<int[:self.ntotal]> &self.data.elems[0])
@property
def nlistatoms(self):
return self.data.nlistatoms
@property
def natomneigh(self):
return self.data.natomneigh
@property
def numneighs(self):
if self.data.numneighs is NULL:
return None
return np.asarray(<int[:self.natomneigh]> &self.data.numneighs[0])
@property
def iatoms(self):
if self.data.iatoms is NULL:
return None
return np.asarray(<int[:self.natomneigh]> &self.data.iatoms[0])
@property
def ielems(self):
if self.data.ielems is NULL:
return None
return np.asarray(<int[:self.natomneigh]> &self.data.ielems[0])
@property
def npairs(self):
return self.data.npairs
@property
def pair_i(self):
if self.data.pair_i is NULL:
return None
return np.asarray(<int[:self.npairs]> &self.data.pair_i[0])
@property
def pair_j(self):
return self.jatoms
@property
def jatoms(self):
if self.data.jatoms is NULL:
return None
return np.asarray(<int[:self.npairs]> &self.data.jatoms[0])
@property
def jelems(self):
if self.data.jelems is NULL:
return None
return np.asarray(<int[:self.npairs]> &self.data.jelems[0])
@property
def rij(self):
if self.data.rij is NULL:
return None
return np.asarray(<double[:self.npairs, :3]> &self.data.rij[0][0])
@write_only_property
def graddesc(self, value):
if self.data.graddesc is NULL:
raise ValueError("attempt to set NULL graddesc")
cdef double[:, :, :] graddesc_view = <double[:self.npairs, :self.ndescriptors, :3]> &self.data.graddesc[0][0][0]
cdef double[:, :, :] value_view = value
graddesc_view[:] = value_view
@property
def eflag(self):
return self.data.eflag
@property
def vflag(self):
return self.data.vflag
# Interface between C and Python compute functions
cdef class MLIAPUnifiedInterface:
cdef MLIAPDummyModel * model
cdef MLIAPDummyDescriptor * descriptor
cdef unified_impl
def __init__(self, unified_impl):
self.model = NULL
self.descriptor = NULL
self.unified_impl = unified_impl
def compute_gradients(self, data):
self.unified_impl.compute_gradients(data)
def compute_descriptors(self, data):
self.unified_impl.compute_descriptors(data)
def compute_forces(self, data):
self.unified_impl.compute_forces(data)
cdef public void compute_gradients_python(unified_int, MLIAPData *data) except * with gil:
pydata = MLIAPDataPy()
pydata.data = data
unified_int.compute_gradients(pydata)
cdef public void compute_descriptors_python(unified_int, MLIAPData *data) except * with gil:
pydata = MLIAPDataPy()
pydata.data = data
unified_int.compute_descriptors(pydata)
cdef public void compute_forces_python(unified_int, MLIAPData *data) except * with gil:
pydata = MLIAPDataPy()
pydata.data = data
unified_int.compute_forces(pydata)
# Create a MLIAPUnifiedInterface and connect it to the dummy model, descriptor
cdef public object mliap_unified_connect(char *fname, MLIAPDummyModel * model,
MLIAPDummyDescriptor * descriptor) with gil:
str_fname = fname.decode('utf-8')
if str_fname == 'EXISTS':
if LOADED_MODEL is None:
raise ValueError("No unified model loaded")
unified = LOADED_MODEL
elif str_fname.endswith(".pt") or str_fname.endswith('.pth'):
import torch
unified = torch.load(str_fname)
else:
with open(str_fname, 'rb') as pfile:
unified = pickle.load(pfile)
unified_int = MLIAPUnifiedInterface(unified)
unified_int.model = model
unified_int.descriptor = descriptor
unified.interface = unified_int
if unified.ndescriptors is None:
raise ValueError("no descriptors set")
unified_int.descriptor.ndescriptors = <int>unified.ndescriptors
unified_int.descriptor.rcutfac = <double>unified.rcutfac
unified_int.model.ndescriptors = <int>unified.ndescriptors
unified_int.model.nparams = <int>unified.nparams
if unified.element_types is None:
raise ValueError("no element type set")
cdef int nelements = <int>len(unified.element_types)
cdef char **elements = <char**>malloc(nelements * sizeof(char*))
if not elements:
raise MemoryError("failed to allocate memory for element names")
cdef char *elem_name
for i, elem in enumerate(unified.element_types):
elem_name_bytes = elem.encode('UTF-8')
elem_name = elem_name_bytes
elements[i] = &elem_name[0]
unified_int.descriptor.set_elements(elements, nelements)
unified_int.model.nelements = nelements
free(elements)
return unified_int
# For pre-loading a Python model
def load_from_python(unified):
global LOADED_MODEL
LOADED_MODEL = unified

81
src/ML-IAP/pair_mliap.cpp
View File

@ -25,14 +25,17 @@
#include "mliap_model_nn.h"
#include "mliap_model_quadratic.h"
#ifdef MLIAP_PYTHON
#include "mliap_unified.h"
#include "mliap_model_python.h"
#endif
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "memory.h"
#include "neighbor.h"
#include "neigh_request.h"
#include <cmath>
#include <cstring>
@ -64,6 +67,7 @@ PairMLIAP::~PairMLIAP()
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(cutghost);
memory->destroy(map);
}
}
@ -78,10 +82,12 @@ void PairMLIAP::compute(int eflag, int vflag)
// consistency checks
if (data->ndescriptors != model->ndescriptors)
error->all(FLERR, "Incompatible model and descriptor descriptor count");
error->all(FLERR, "Inconsistent model and descriptor descriptor count: {} vs {}",
model->ndescriptors, data->ndescriptors);
if (data->nelements != model->nelements)
error->all(FLERR, "Incompatible model and descriptor element count");
error->all(FLERR, "Inconsistent model and descriptor element count: {} vs {}",
model->nelements, data->nelements);
ev_init(eflag, vflag);
data->generate_neighdata(list, eflag, vflag);
@ -93,11 +99,11 @@ void PairMLIAP::compute(int eflag, int vflag)
// compute E_i and beta_i = dE_i/dB_i for all i in list
model->compute_gradients(data);
e_tally(data);
// calculate force contributions beta_i*dB_i/dR_j
descriptor->compute_forces(data);
e_tally(data);
// calculate stress
@ -115,6 +121,7 @@ void PairMLIAP::allocate()
memory->create(setflag,n+1,n+1,"pair:setflag");
memory->create(cutsq,n+1,n+1,"pair:cutsq");
memory->create(cutghost,n+1,n+1,"pair:cutghost");
memory->create(map,n+1,"pair:map");
}
@ -124,8 +131,7 @@ void PairMLIAP::allocate()
void PairMLIAP::settings(int narg, char ** arg)
{
if (narg < 4)
error->all(FLERR,"Illegal pair_style command");
if (narg < 2) utils::missing_cmd_args(FLERR, "pair_style mliap", error);
// set flags for required keywords
@ -140,47 +146,70 @@ void PairMLIAP::settings(int narg, char ** arg)
while (iarg < narg) {
if (strcmp(arg[iarg],"model") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal pair_style mliap command");
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "pair_style mliap model", error);
if (modelflag) error->all(FLERR,"Illegal multiple pair_style mliap model definition");
if (strcmp(arg[iarg+1],"linear") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal pair_style mliap command");
if (iarg+3 > narg) utils::missing_cmd_args(FLERR, "pair_style mliap model linear", error);
model = new MLIAPModelLinear(lmp,arg[iarg+2]);
iarg += 3;
} else if (strcmp(arg[iarg+1],"quadratic") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal pair_style mliap command");
if (iarg+3 > narg) utils::missing_cmd_args(FLERR, "pair_style mliap model quadratic", error);
model = new MLIAPModelQuadratic(lmp,arg[iarg+2]);
iarg += 3;
} else if (strcmp(arg[iarg+1],"nn") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal pair_style mliap command");
if (iarg+3 > narg) utils::missing_cmd_args(FLERR, "pair_style mliap model nn", error);
model = new MLIAPModelNN(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");
#ifdef MLIAP_PYTHON
if (iarg+3 > narg) utils::missing_cmd_args(FLERR, "pair_style mliap mliappy", error);
model = new MLIAPModelPython(lmp,arg[iarg+2]);
iarg += 3;
#else
error->all(FLERR,"Using pair_style mliap model mliappy requires ML-IAP with python support");
#endif
} else error->all(FLERR,"Illegal pair_style mliap command");
} else error->all(FLERR,"Unknown pair_style mliap model keyword: {}", arg[iarg]);
modelflag = 1;
} else if (strcmp(arg[iarg],"descriptor") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal pair_style mliap command");
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "pair_style mliap descriptor", error);
if (descriptorflag) error->all(FLERR,"Illegal multiple pair_style mliap descriptor definition");
if (strcmp(arg[iarg+1],"sna") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal pair_style mliap command");
if (iarg+3 > narg) utils::missing_cmd_args(FLERR, "pair_style mliap descriptor sna", error);
descriptor = new MLIAPDescriptorSNAP(lmp,arg[iarg+2]);
iarg += 3;
} else if (strcmp(arg[iarg+1],"so3") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal pair_style mliap command");
if (iarg+3 > narg) utils::missing_cmd_args(FLERR, "pair_style mliap descriptor so3", error);
descriptor = new MLIAPDescriptorSO3(lmp,arg[iarg+2]);
iarg += 3;
} else error->all(FLERR,"Illegal pair_style mliap command");
descriptorflag = 1;
} else if (strcmp(arg[iarg], "unified") == 0) {
#ifdef MLIAP_PYTHON
if (modelflag) error->all(FLERR,"Illegal multiple pair_style mliap model definitions");
if (descriptorflag) error->all(FLERR,"Illegal multiple pair_style mliap descriptor definitions");
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "pair_style mliap unified", error);
MLIAPBuildUnified_t build = build_unified(arg[iarg+1], data, lmp);
if (iarg+3 > narg) {
ghostneigh = 0;
} else {
ghostneigh = utils::logical(FLERR, arg[iarg+2], false, lmp);
}
iarg += 3;
model = build.model;
descriptor = build.descriptor;
modelflag = 1;
descriptorflag = 1;
#else
error->all(FLERR,"Using pair_style mliap unified requires ML-IAP with python support");
#endif
} else
error->all(FLERR,"Illegal pair_style mliap command");
error->all(FLERR,"Unknown pair_style mliap keyword: {}", arg[iarg]);
}
if (modelflag == 0 || descriptorflag == 0)
error->all(FLERR,"Illegal pair_style command");
error->all(FLERR,"Incomplete pair_style mliap setup: need model and descriptor, or unified");
}
/* ----------------------------------------------------------------------
@ -196,11 +225,6 @@ void PairMLIAP::coeff(int narg, char **arg)
char* type2 = arg[1];
char** elemtypes = &arg[2];
// insure I,J args are * *
if (strcmp(type1,"*") != 0 || strcmp(type2,"*") != 0)
error->all(FLERR,"Incorrect args for pair coefficients");
// read args that map atom types to elements
// map[i] = which element the Ith atom type is, -1 if not mapped
// map[0] is not used
@ -313,7 +337,11 @@ void PairMLIAP::init_style()
// need a full neighbor list
neighbor->add_request(this, NeighConst::REQ_FULL);
if (ghostneigh == 1) {
neighbor->add_request(this, NeighConst::REQ_FULL | NeighConst::REQ_GHOST);
} else {
neighbor->add_request(this, NeighConst::REQ_FULL);
}
}
@ -324,7 +352,9 @@ void PairMLIAP::init_style()
double PairMLIAP::init_one(int i, int j)
{
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
return sqrt(descriptor->cutsq[map[i]][map[j]]);
double cutmax = sqrt(descriptor->cutsq[map[i]][map[j]]);
cutghost[i][j] = cutghost[j][i] = 2.0 * cutmax + neighbor->skin;
return cutmax;
}
/* ----------------------------------------------------------------------
@ -338,6 +368,7 @@ double PairMLIAP::memory_usage()
int n = atom->ntypes+1;
bytes += (double)n*n*sizeof(int); // setflag
bytes += (double)n*n*sizeof(int); // cutsq
bytes += (double)n*n*sizeof(int); // cutghost
bytes += (double)n*sizeof(int); // map
bytes += descriptor->memory_usage(); // Descriptor object
bytes += model->memory_usage(); // Model object

5
src/PYTHON/python_impl.cpp
View File

@ -29,10 +29,12 @@
#ifdef MLIAP_PYTHON
#include "mliap_model_python.h"
#include "mliap_unified.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"
#include "mliap_unified_couple.h"
#endif
using namespace LAMMPS_NS;
@ -66,6 +68,9 @@ PythonImpl::PythonImpl(LAMMPS *lmp) : Pointers(lmp)
// 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.");
err = PyImport_AppendInittab("mliap_unified_couple", PyInit_mliap_unified_couple);
if (err) error->all(FLERR, "Could not register MLIAPPY unified embedded python module.");
#endif
Py_Initialize();

27
src/library.cpp
View File

@ -57,6 +57,10 @@
#include "exceptions.h"
#endif
#if defined(LMP_PYTHON)
#include <Python.h>
#endif
using namespace LAMMPS_NS;
// for printing the non-null pointer argument warning only once
@ -5654,6 +5658,29 @@ int lammps_get_last_error_message(void *handle, char *buffer, int buf_size) {
return 0;
}
/* ---------------------------------------------------------------------- */
/** Return API version of embedded Python interpreter
\verbatim embed:rst
This function is used by the ML-IAP python code (mliappy) to verify
the API version of the embedded python interpreter of the PYTHON
package. It returns -1 if the PYTHON package is not enabled.
.. versionadded:: TBD
\endverbatim
*
* \return PYTHON_API_VERSION constant of the python interpreter or -1 */
int lammps_python_api_version() {
#if defined(LMP_PYTHON)
return PYTHON_API_VERSION;
#else
return -1;
#endif
}
// Local Variables:
// fill-column: 72
// End:

2
src/library.h
View File

@ -270,6 +270,8 @@ void lammps_force_timeout(void *handle);
int lammps_has_error(void *handle);
int lammps_get_last_error_message(void *handle, char *buffer, int buf_size);
int lammps_python_api_version();
#ifdef __cplusplus
}
#endif

15
unittest/force-styles/CMakeLists.txt
View File

@ -79,6 +79,7 @@ else()
set(FORCE_TEST_ENVIRONMENT PYTHONPATH=${TEST_INPUT_FOLDER}:$ENV{PYTHONPATH})
endif()
list(APPEND FORCE_TEST_ENVIRONMENT "PYTHONUNBUFFERED=1")
list(APPEND FORCE_TEST_ENVIRONMENT "PYTHONDONTWRITEBYTECODE=1")
list(APPEND FORCE_TEST_ENVIRONMENT "OMP_PROC_BIND=false")
list(APPEND FORCE_TEST_ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR}")
get_property(BUILD_IS_MULTI_CONFIG GLOBAL PROPERTY GENERATOR_IS_MULTI_CONFIG)
@ -213,7 +214,7 @@ foreach(TEST ${FIX_TIMESTEP_TESTS})
if(WIN32)
set_tests_properties(${TNAME} PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR};PYTHONPATH=${TEST_INPUT_FOLDER}\\\;${LAMMPS_PYTHON_DIR}")
else()
set_tests_properties(${TNAME} PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR};PYTHONPATH=${TEST_INPUT_FOLDER}:${LAMMPS_PYTHON_DIR}:$ENV{PYTHONPATH}")
set_tests_properties(${TNAME} PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR};PYTHONPATH=${TEST_INPUT_FOLDER}:${LAMMPS_PYTHON_DIR}:$ENV{PYTHONPATH};PYTHONDONTWRITEBYTECODE=1")
endif()
set_tests_properties(${TNAME} PROPERTIES LABELS "${TEST_TAGS}")
endforeach()
@ -231,7 +232,7 @@ foreach(TEST ${DIHEDRAL_TESTS})
continue()
endif()
add_test(NAME ${TNAME} COMMAND test_dihedral_style ${TEST})
set_tests_properties(${TNAME} PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR};PYTHONPATH=${TEST_INPUT_FOLDER}:$ENV{PYTHONPATH}")
set_tests_properties(${TNAME} PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR};PYTHONPATH=${TEST_INPUT_FOLDER}:$ENV{PYTHONPATH};PYTHONDONTWRITEBYTECODE=1")
set_tests_properties(${TNAME} PROPERTIES LABELS "${TEST_TAGS}")
endforeach()
@ -248,6 +249,14 @@ foreach(TEST ${IMPROPER_TESTS})
continue()
endif()
add_test(NAME ${TNAME} COMMAND test_improper_style ${TEST})
set_tests_properties(${TNAME} PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR};PYTHONPATH=${TEST_INPUT_FOLDER}:$ENV{PYTHONPATH}")
set_tests_properties(${TNAME} PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR};PYTHONPATH=${TEST_INPUT_FOLDER}:$ENV{PYTHONPATH};PYTHONDONTWRITEBYTECODE=1")
set_tests_properties(${TNAME} PROPERTIES LABELS "${TEST_TAGS}")
endforeach()
if(MLIAP_ENABLE_PYTHON AND (NOT WIN32))
add_executable(test_mliappy_unified test_mliappy_unified.cpp)
target_link_libraries(test_mliappy_unified PRIVATE lammps GTest::GMockMain)
add_test(NAME TestMliapPyUnified COMMAND test_mliappy_unified)
set_tests_properties(TestMliapPyUnified PROPERTIES ENVIRONMENT "PYTHONPATH=${LAMMPS_PYTHON_DIR};PYTHONDONTWRITEBYTECODE=1")
endif()

116
unittest/force-styles/test_mliappy_unified.cpp Normal file
View File

@ -0,0 +1,116 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, 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.
------------------------------------------------------------------------- */
#include "library.h"
#include "gmock/gmock.h"
#include "gtest/gtest.h"
const char pickle[] = "python create_pickle here \"\"\"\n"
"import lammps\n"
"import lammps.mliap\n"
"from lammps.mliap.mliap_unified_lj import MLIAPUnifiedLJ\n"
"def create_pickle():\n"
" unified = MLIAPUnifiedLJ(['Ar'])\n"
" unified.pickle('mliap_unified_lj_Ar.pkl')\n"
"\"\"\"\n";
const char first[] = "units lj\n"
"atom_style atomic\n"
"lattice fcc 0.8442\n"
"region box block 0 2 0 2 0 2\n"
"create_box 1 box\n"
"create_atoms 1 box\n"
"mass 1 1.0\n"
"velocity all create 3.0 87287 loop geom\n";
const char second[] = "neighbor 0.3 bin\n"
"neigh_modify every 20 delay 0 check no\n"
"fix 1 all nve\n"
"run 2 post no\n";
namespace LAMMPS_NS {
TEST(MliapUnified, VersusLJMelt)
{
const char *lmpargv[] = {"melt", "-log", "none", "-nocite"};
int lmpargc = sizeof(lmpargv) / sizeof(const char *);
void *ljmelt = lammps_open_no_mpi(lmpargc, (char **)lmpargv, nullptr);
void *mliap = lammps_open_no_mpi(lmpargc, (char **)lmpargv, nullptr);
lammps_commands_string(ljmelt, first);
lammps_command(ljmelt, "pair_style lj/cut 2.5");
lammps_command(ljmelt, "pair_coeff * * 1.0 1.0");
lammps_commands_string(ljmelt, second);
lammps_command(mliap, pickle);
lammps_command(mliap, "python create_pickle invoke");
lammps_commands_string(mliap, first);
lammps_command(mliap, "pair_style mliap unified mliap_unified_lj_Ar.pkl 0");
lammps_command(mliap, "pair_coeff * * Ar");
lammps_commands_string(mliap, second);
double lj_pe = lammps_get_thermo(ljmelt, "pe");
double ml_pe = lammps_get_thermo(mliap, "pe");
EXPECT_NEAR(lj_pe, ml_pe, 1.0e-14);
double lj_ke = lammps_get_thermo(ljmelt, "ke");
double ml_ke = lammps_get_thermo(mliap, "ke");
EXPECT_NEAR(lj_ke, ml_ke, 1.0e-14);
double lj_press = lammps_get_thermo(ljmelt, "press");
double ml_press = lammps_get_thermo(mliap, "press");
EXPECT_NEAR(lj_press, ml_press, 1.0e-14);
lammps_command(mliap, "shell rm mliap_unified_lj_Ar.pkl");
lammps_close(ljmelt);
lammps_close(mliap);
}
TEST(MliapUnified, VersusLJMeltGhost)
{
const char *lmpargv[] = {"melt", "-log", "none", "-nocite"};
int lmpargc = sizeof(lmpargv) / sizeof(const char *);
void *ljmelt = lammps_open_no_mpi(lmpargc, (char **)lmpargv, nullptr);
void *mliap = lammps_open_no_mpi(lmpargc, (char **)lmpargv, nullptr);
lammps_commands_string(ljmelt, first);
lammps_command(ljmelt, "pair_style lj/cut 2.5");
lammps_command(ljmelt, "pair_coeff * * 1.0 1.0");
lammps_commands_string(ljmelt, second);
lammps_command(mliap, pickle);
lammps_command(mliap, "python create_pickle invoke");
lammps_commands_string(mliap, first);
lammps_command(mliap, "pair_style mliap unified mliap_unified_lj_Ar.pkl 1");
lammps_command(mliap, "pair_coeff * * Ar");
lammps_commands_string(mliap, second);
double lj_pe = lammps_get_thermo(ljmelt, "pe");
double ml_pe = lammps_get_thermo(mliap, "pe");
EXPECT_NEAR(lj_pe, ml_pe, 1.0e-14);
double lj_ke = lammps_get_thermo(ljmelt, "ke");
double ml_ke = lammps_get_thermo(mliap, "ke");
EXPECT_NEAR(lj_ke, ml_ke, 1.0e-14);
double lj_press = lammps_get_thermo(ljmelt, "press");
double ml_press = lammps_get_thermo(mliap, "press");
EXPECT_NEAR(lj_press, ml_press, 1.0e-14);
lammps_command(mliap, "shell rm mliap_unified_lj_Ar.pkl");
lammps_close(ljmelt);
lammps_close(mliap);
}
} // namespace LAMMPS_NS

1
unittest/python/CMakeLists.txt
View File

@ -44,6 +44,7 @@ else()
endif()
list(APPEND PYTHON_TEST_ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR}")
list(APPEND PYTHON_TEST_ENVIRONMENT "PYTHONUNBUFFERED=1")
list(APPEND PYTHON_TEST_ENVIRONMENT "PYTHONDONTWRITEBYTECODE=1")
if(APPLE)
list(APPEND PYTHON_TEST_ENVIRONMENT "DYLD_LIBRARY_PATH=${LAMMPS_LIB_PATH}:$ENV{DYLD_LIBRARY_PATH}")
elseif(WIN32)