diff --git a/cmake/Modules/Packages/ML-IAP.cmake b/cmake/Modules/Packages/ML-IAP.cmake index 63f91ba8d3..0e711cde8f 100644 --- a/cmake/Modules/Packages/ML-IAP.cmake +++ b/cmake/Modules/Packages/ML-IAP.cmake @@ -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() diff --git a/doc/src/Library_utility.rst b/doc/src/Library_utility.rst index da64e3b8f0..9d16daae0b 100644 --- a/doc/src/Library_utility.rst +++ b/doc/src/Library_utility.rst @@ -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 + diff --git a/doc/src/pair_mliap.rst b/doc/src/pair_mliap.rst index 900a2c2e6a..dd021eb84c 100644 --- a/doc/src/pair_mliap.rst +++ b/doc/src/pair_mliap.rst @@ -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 `, including the linear, quadratic, and @@ -61,11 +68,11 @@ basis :ref:`(Bartok) ` and :ref:`(Zagaceta) `. 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) ` +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 """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index aba474d257..c719620fc6 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -1206,6 +1206,7 @@ gewald Gezelter Gflop gfortran +ghostneigh ghostwhite Giacomo gif diff --git a/examples/mliap/README b/examples/mliap/README index c16fa218e0..070ce86bfd 100644 --- a/examples/mliap/README +++ b/examples/mliap/README @@ -1,6 +1,6 @@ -This directory contains multiple examples of -machine-learning potentials defined using the -MLIAP package in LAMMPS. The input files +This directory contains multiple examples of +machine-learning potentials defined using the +ML-IAP package in LAMMPS. The input files are described below. in.mliap.snap.Ta06A @@ -21,14 +21,14 @@ 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 +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 +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 @@ -41,8 +41,8 @@ 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 +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 @@ -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. @@ -62,7 +62,7 @@ working by running ../examples/in.python. Before testing this, ensure that the previous method (running a LAMMPS executable) works. -You can run the example in serial: +You can run the example in serial: `python mliap_pytorch_Ta06A.py` @@ -80,25 +80,25 @@ 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. +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, +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. +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`. +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 +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. in.mliap.nn.Ta06A @@ -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 diff --git a/examples/mliap/in.mliap.unified.lj.Ar b/examples/mliap/in.mliap.unified.lj.Ar new file mode 100644 index 0000000000..190ea1288f --- /dev/null +++ b/examples/mliap/in.mliap.unified.lj.Ar @@ -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 diff --git a/examples/mliap/log.23Sep22.mliap.unified.lj.Ar.g++.1 b/examples/mliap/log.23Sep22.mliap.unified.lj.Ar.g++.1 new file mode 100644 index 0000000000..1059b65538 --- /dev/null +++ b/examples/mliap/log.23Sep22.mliap.unified.lj.Ar.g++.1 @@ -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 diff --git a/examples/mliap/log.23Sep22.mliap.unified.lj.Ar.g++.4 b/examples/mliap/log.23Sep22.mliap.unified.lj.Ar.g++.4 new file mode 100644 index 0000000000..236e415ae6 --- /dev/null +++ b/examples/mliap/log.23Sep22.mliap.unified.lj.Ar.g++.4 @@ -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 diff --git a/examples/mliap/mliap_numpy_Ta06A.py b/examples/mliap/mliap_numpy_Ta06A.py new file mode 100644 index 0000000000..25c928cd89 --- /dev/null +++ b/examples/mliap/mliap_numpy_Ta06A.py @@ -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) + diff --git a/examples/mliap/mliap_unified_lj_Ar.py b/examples/mliap/mliap_unified_lj_Ar.py new file mode 100644 index 0000000000..2b404bd98e --- /dev/null +++ b/examples/mliap/mliap_unified_lj_Ar.py @@ -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) diff --git a/python/lammps/mliap/__init__.py b/python/lammps/mliap/__init__.py index cfe3fb6b38..b888eb55ce 100644 --- a/python/lammps/mliap/__init__.py +++ b/python/lammps/mliap/__init__.py @@ -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 diff --git a/python/lammps/mliap/loader.py b/python/lammps/mliap/loader.py index dff791bfc1..e1e60c4f50 100644 --- a/python/lammps/mliap/loader.py +++ b/python/lammps/mliap/loader.py @@ -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) diff --git a/python/lammps/mliap/mliap_unified_abc.py b/python/lammps/mliap/mliap_unified_abc.py new file mode 100644 index 0000000000..3243df85c4 --- /dev/null +++ b/python/lammps/mliap/mliap_unified_abc.py @@ -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) diff --git a/python/lammps/mliap/mliap_unified_lj.py b/python/lammps/mliap/mliap_unified_lj.py new file mode 100644 index 0000000000..37e7170d4a --- /dev/null +++ b/python/lammps/mliap/mliap_unified_lj.py @@ -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 diff --git a/python/lammps/mliap/pytorch.py b/python/lammps/mliap/pytorch.py index d699c239b0..51b953fd61 100644 --- a/python/lammps/mliap/pytorch.py +++ b/python/lammps/mliap/pytorch.py @@ -80,10 +80,10 @@ class TorchWrapper(torch.nn.Module): n_params : torch.nn.Module (None) Number of NN model parameters - + device : torch.nn.Module (None) Accelerator device - + dtype : torch.dtype (torch.float64) Dtype to use on device """ @@ -325,6 +325,6 @@ class ElemwiseModels(torch.nn.Module): per_atom_attributes = torch.zeros(elems.size(dim=0), dtype=self.dtype) given_elems, elem_indices = torch.unique(elems, return_inverse=True) for i, elem in enumerate(given_elems): - self.subnets[elem].to(self.dtype) + self.subnets[elem].to(self.dtype) per_atom_attributes[elem_indices == i] = self.subnets[elem](descriptors[elem_indices == i]).flatten() return per_atom_attributes diff --git a/src/ML-IAP/Install.sh b/src/ML-IAP/Install.sh index 3a1a7493aa..27c9b19b01 100755 --- a/src/ML-IAP/Install.sh +++ b/src/ML-IAP/Install.sh @@ -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 diff --git a/src/ML-IAP/mliap_data.cpp b/src/ML-IAP/mliap_data.cpp index e96e87ba11..a91e66b211 100644 --- a/src/ML-IAP/mliap_data.cpp +++ b/src/ML-IAP/mliap_data.cpp @@ -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" diff --git a/src/ML-IAP/mliap_data.h b/src/ML-IAP/mliap_data.h index 258261fa1a..562dcd63cb 100644 --- a/src/ML-IAP/mliap_data.h +++ b/src/ML-IAP/mliap_data.h @@ -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 diff --git a/src/ML-IAP/mliap_descriptor.cpp b/src/ML-IAP/mliap_descriptor.cpp index ff8cfb43a9..ce8abf4e56 100644 --- a/src/ML-IAP/mliap_descriptor.cpp +++ b/src/ML-IAP/mliap_descriptor.cpp @@ -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); } diff --git a/src/ML-IAP/mliap_descriptor.h b/src/ML-IAP/mliap_descriptor.h index 156f8d1e92..6db872ccaf 100644 --- a/src/ML-IAP/mliap_descriptor.h +++ b/src/ML-IAP/mliap_descriptor.h @@ -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: }; diff --git a/src/ML-IAP/mliap_model_python.h b/src/ML-IAP/mliap_model_python.h index 4243c67332..915e79e53c 100644 --- a/src/ML-IAP/mliap_model_python.h +++ b/src/ML-IAP/mliap_model_python.h @@ -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: }; diff --git a/src/ML-IAP/mliap_unified.cpp b/src/ML-IAP/mliap_unified.cpp new file mode 100644 index 0000000000..6bcdc88289 --- /dev/null +++ b/src/ML-IAP/mliap_unified.cpp @@ -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 + +#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 diff --git a/src/ML-IAP/mliap_unified.h b/src/ML-IAP/mliap_unified.h new file mode 100644 index 0000000000..adc900d8c1 --- /dev/null +++ b/src/ML-IAP/mliap_unified.h @@ -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 + +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 diff --git a/src/ML-IAP/mliap_unified_couple.pyx b/src/ML-IAP/mliap_unified_couple.pyx new file mode 100644 index 0000000000..3fde99a25e --- /dev/null +++ b/src/ML-IAP/mliap_unified_couple.pyx @@ -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( &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 = &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 = &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 = &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 = &self.data.eatoms[0] + cdef double[:] value_view = value + eatoms_view[:] = value_view + + @write_only_property + def energy(self, value): + self.data.energy = 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( &self.data.gamma[0][0]) + + @property + def gamma_row_index(self): + if self.data.gamma_row_index is NULL: + return None + return np.asarray( &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( &self.data.gamma_col_index[0][0]) + + @property + def egradient(self): + if self.data.egradient is NULL: + return None + return np.asarray( &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( &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( &self.data.numneighs[0]) + + @property + def iatoms(self): + if self.data.iatoms is NULL: + return None + return np.asarray( &self.data.iatoms[0]) + + @property + def ielems(self): + if self.data.ielems is NULL: + return None + return np.asarray( &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( &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( &self.data.jatoms[0]) + + @property + def jelems(self): + if self.data.jelems is NULL: + return None + return np.asarray( &self.data.jelems[0]) + + @property + def rij(self): + if self.data.rij is NULL: + return None + return np.asarray( &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 = &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 = unified.ndescriptors + unified_int.descriptor.rcutfac = unified.rcutfac + unified_int.model.ndescriptors = unified.ndescriptors + unified_int.model.nparams = unified.nparams + + if unified.element_types is None: + raise ValueError("no element type set") + + cdef int nelements = len(unified.element_types) + cdef char **elements = 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 diff --git a/src/ML-IAP/pair_mliap.cpp b/src/ML-IAP/pair_mliap.cpp index 179649caff..af8768568e 100644 --- a/src/ML-IAP/pair_mliap.cpp +++ b/src/ML-IAP/pair_mliap.cpp @@ -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 #include @@ -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 diff --git a/src/PYTHON/python_impl.cpp b/src/PYTHON/python_impl.cpp index 5a7098814b..aba8e88e10 100644 --- a/src/PYTHON/python_impl.cpp +++ b/src/PYTHON/python_impl.cpp @@ -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(); diff --git a/src/library.cpp b/src/library.cpp index cc547106f2..16381a089d 100644 --- a/src/library.cpp +++ b/src/library.cpp @@ -57,6 +57,10 @@ #include "exceptions.h" #endif +#if defined(LMP_PYTHON) +#include +#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: diff --git a/src/library.h b/src/library.h index 4e50cbee84..1eec57898e 100644 --- a/src/library.h +++ b/src/library.h @@ -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 diff --git a/unittest/force-styles/CMakeLists.txt b/unittest/force-styles/CMakeLists.txt index d2345ceaa7..294ece3a1f 100644 --- a/unittest/force-styles/CMakeLists.txt +++ b/unittest/force-styles/CMakeLists.txt @@ -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() + diff --git a/unittest/force-styles/test_mliappy_unified.cpp b/unittest/force-styles/test_mliappy_unified.cpp new file mode 100644 index 0000000000..1eb800a954 --- /dev/null +++ b/unittest/force-styles/test_mliappy_unified.cpp @@ -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 diff --git a/unittest/python/CMakeLists.txt b/unittest/python/CMakeLists.txt index 0aebc6fc7c..78b0d18446 100644 --- a/unittest/python/CMakeLists.txt +++ b/unittest/python/CMakeLists.txt @@ -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)