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

Merge pull request #2641 from pedroantoniosantosf/mliap_nn

Browse Source 
Add support for neural network (nn) models to MLIAP package to compute energies and forces without going through Python
This commit is contained in:
Axel Kohlmeyer
2021-03-08 20:27:01 -05:00
committed by GitHub
parent 37afa0e6d5 d1c79d54fa
commit 7ae2da0ace
36 changed files with 1452 additions and 308 deletions
Download Patch File Download Diff File Expand all files Collapse all files

25
doc/src/pair_mliap.rst
View File

@ -16,7 +16,7 @@ Syntax
.. parsed-literal:: .. parsed-literal::
*model* values = style filename *model* values = style filename
style = *linear* or *quadratic* or *mliappy* style = *linear* or *quadratic* or *nn* or *mliappy*
filename = name of file containing model definitions filename = name of file containing model definitions
*descriptor* values = style filename *descriptor* values = style filename
style = *sna* style = *sna*
@ -45,7 +45,7 @@ pair style currently supports just one descriptor style, but it is
is straightforward to add new descriptor styles. is straightforward to add new descriptor styles.
The SNAP descriptor style *sna* is the same as that used by :doc:`pair_style snap <pair_snap>`, The SNAP descriptor style *sna* is the same as that used by :doc:`pair_style snap <pair_snap>`,
including the linear, quadratic, and chem variants. including the linear, quadratic, and chem variants.
The available models are *linear*, *quadratic*, and *mliappy*. The available models are *linear*, *quadratic*, *nn*, and *mliappy*.
The *mliappy* style can be used to couple python models, The *mliappy* style can be used to couple python models,
e.g. PyTorch neural network energy models, and requires building e.g. PyTorch neural network energy models, and requires building
LAMMPS with the PYTHON package (see below). LAMMPS with the PYTHON package (see below).
@ -77,13 +77,32 @@ line must contain two integers:
* nelems = Number of elements * nelems = Number of elements
* nparams = Number of parameters * nparams = Number of parameters
This is followed by one block for each of the *nelem* elements. When the *model* keyword is *linear* or *quadratic*,
this is followed by one block for each of the *nelem* elements.
Each block consists of *nparams* parameters, one per line. Each block consists of *nparams* parameters, one per line.
Note that this format is similar, but not identical to that used Note that this format is similar, but not identical to that used
for the :doc:`pair_style snap <pair_snap>` coefficient file. for the :doc:`pair_style snap <pair_snap>` coefficient file.
Specifically, the line containing the element weight and radius is omitted, Specifically, the line containing the element weight and radius is omitted,
since these are handled by the *descriptor*. since these are handled by the *descriptor*.
When the *model* keyword is *nn* (neural networks), the model file can contain
blank and comment lines (start with #) anywhere. The second non-blank non-comment
line must contain the string NET, followed by two integers:
* ndescriptors = Number of descriptors
* nlayers = Number of layers (including the hidden layers and the output layer)
and followed by a sequence of a string and an integer for each layer:
* Activation function (linear, sigmoid, tanh or relu)
* nnodes = Number of nodes
This is followed by one block for each of the *nelem* elements. Each block consists
of *scale0* minimum value, *scale1* (maximum - minimum) value,
in order to normalize the descriptors, followed by *nparams* parameters,
including *bias* and *weights* of the model, starting with the first node of the first layer
and so on, with a maximum of 30 values per line.
Notes on mliappy models: Notes on mliappy models:
When the *model* keyword is *mliappy*, the filename should end in '.pt', 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 '.pth' for pytorch models, or be a pickle file. To load a model from

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

@ -2098,6 +2098,7 @@ ncol
ncorr ncorr
ncount ncount
nd nd
ndescriptors
ndihedrals ndihedrals
ndihedraltypes ndihedraltypes
Ndihedraltype Ndihedraltype
@ -2158,6 +2159,7 @@ niu
Nk Nk
nktv nktv
nl nl
nlayers
nlen nlen
Nlines Nlines
nlo nlo
@ -2173,6 +2175,7 @@ Nmin
nmin nmin
Nmols Nmols
nn nn
nnodes
Nocedal Nocedal
nocite nocite
nocoeff nocoeff
@ -2659,6 +2662,7 @@ relaxbox
relink relink
relres relres
relTol relTol
relu
remappings remappings
remd remd
Ren Ren
@ -2888,6 +2892,7 @@ si
SiC SiC
Siepmann Siepmann
Sievers Sievers
sigmoid
Sij Sij
Sikandar Sikandar
Silbert Silbert

1
examples/mliap/Cu.nn.mliap.model Symbolic link
View File

@ -0,0 +1 @@
../../potentials/Cu.nn.mliap.model

1
examples/mliap/Cu.snap.mliap.descriptor Symbolic link
View File

@ -0,0 +1 @@
../../potentials/Cu.snap.mliap.descriptor

9
examples/mliap/README
View File

@ -101,3 +101,12 @@ pickle data that will execute arbitrary code during unpickling. Never
load data that could have come from an untrusted source, or that load data that could have come from an untrusted source, or that
could have been tampered with. Only load data you trust. could have been tampered with. Only load data you trust.
in.mliap.nn.Ta06A
-------------------
Run linear SNAP using the "nn" model style, equivalent to examples/snap/in.snap.Ta06A
in.mliap.nn.cu
-------------------------
Run a neural network potential for Cu, a combination of SNAP descriptors and the "nn" model style

1
examples/mliap/Ta06A.nn.mliap Symbolic link
View File

@ -0,0 +1 @@
../../potentials/Ta06A.nn.mliap

1
examples/mliap/Ta06A.nn.mliap.model Symbolic link
View File

@ -0,0 +1 @@
../../potentials/Ta06A.nn.mliap.model

1
examples/mliap/Ta06A.pytorch.mliap Symbolic link
View File

@ -0,0 +1 @@
../../potentials/Ta06A.pytorch.mliap

53
examples/mliap/in.mliap.nn.Cu Normal file
View File

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

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

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

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

@ -24,7 +24,7 @@ mass 1 180.88
# choose potential # choose potential
include Ta06A.mliap.pytorch include Ta06A.pytorch.mliap
# Setup output # Setup output

134
examples/mliap/log.03Mar21.mliap.nn.Cu.g++.1 Normal file
View File

@ -0,0 +1,134 @@
LAMMPS (10 Feb 2021)
using 1 OpenMP thread(s) per MPI task
# Demonstrate MLIAP interface to SNAP nn Cu potential
# Initialize simulation
variable nsteps index 100
variable nrep equal 4
variable a equal 3.621262
units metal
# generate the box and atom positions using a FCC lattice
variable nx equal ${nrep}
variable nx equal 4
variable ny equal ${nrep}
variable ny equal 4
variable nz equal ${nrep}
variable nz equal 4
boundary p p p
lattice fcc $a
lattice fcc 3.621262
Lattice spacing in x,y,z = 3.6212620 3.6212620 3.6212620
region box block 0 ${nx} 0 ${ny} 0 ${nz}
region box block 0 4 0 ${ny} 0 ${nz}
region box block 0 4 0 4 0 ${nz}
region box block 0 4 0 4 0 4
create_box 1 box
Created orthogonal box = (0.0000000 0.0000000 0.0000000) to (14.485048 14.485048 14.485048)
1 by 1 by 1 MPI processor grid
create_atoms 1 box
Created 256 atoms
create_atoms CPU = 0.000 seconds
mass 1 63.546
# choose potential
pair_style mliap model nn Cu.nn.mliap.model descriptor sna Cu.snap.mliap.descriptor
Reading potential file Cu.nn.mliap.model with DATE: 2021-03-03
Reading potential file Cu.snap.mliap.descriptor with DATE: 2021-03-03
SNAP keyword rcutfac 5.0
SNAP keyword twojmax 6
SNAP keyword nelems 1
SNAP keyword elems Cu
SNAP keyword radelems 0.5
SNAP keyword welems 1.0
SNAP keyword rfac0 0.99363
SNAP keyword rmin0 0
SNAP keyword bzeroflag 0
SNAP keyword switchflag 1
pair_coeff * * Cu
# Setup output
compute eatom all pe/atom
compute energy all reduce sum c_eatom
compute satom all stress/atom NULL
compute str all reduce sum c_satom[1] c_satom[2] c_satom[3]
variable press equal (c_str[1]+c_str[2]+c_str[3])/(3*vol)
thermo_style custom step temp epair c_energy etotal press v_press
thermo 10
thermo_modify norm yes
# Set up NVE run
timestep 0.5e-3
neighbor 1.0 bin
neigh_modify once no every 1 delay 0 check yes
# Run MD
velocity all create 300.0 4928459 loop geom
fix 1 all nve
run ${nsteps}
run 100
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 6
ghost atom cutoff = 6
binsize = 3, bins = 5 5 5
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) = 31.17 | 31.17 | 31.17 Mbytes
Step Temp E_pair c_energy TotEng Press v_press
0 300 -4.0935581 -4.0935581 -4.0549315 -2916.7505 2916.7505
10 294.19698 -4.0927997 -4.0927997 -4.0549203 -2653.7504 2653.7504
20 277.15991 -4.0905816 -4.0905816 -4.0548958 -1875.7718 1875.7718
30 250.30462 -4.0870777 -4.0870777 -4.0548496 -618.45378 618.45378
40 215.94455 -4.0825861 -4.0825861 -4.0547821 1014.6151 -1014.6151
50 177.29307 -4.0775299 -4.0775299 -4.0547025 2846.7675 -2846.7675
60 138.26899 -4.072429 -4.072429 -4.0546261 4698.9654 -4698.9654
70 103.04321 -4.0678341 -4.0678341 -4.0545667 6401.8077 -6401.8077
80 75.426905 -4.0642423 -4.0642423 -4.0545307 7766.2217 -7766.2217
90 58.150738 -4.0620015 -4.0620015 -4.0545143 8668.0598 -8668.0598
100 52.301012 -4.0612408 -4.0612408 -4.0545067 9049.3141 -9049.3141
Loop time of 12.7584 on 1 procs for 100 steps with 256 atoms
Performance: 0.339 ns/day, 70.880 hours/ns, 7.838 timesteps/s
100.0% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 12.755 | 12.755 | 12.755 | 0.0 | 99.97
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.0014069 | 0.0014069 | 0.0014069 | 0.0 | 0.01
Output | 0.00085897 | 0.00085897 | 0.00085897 | 0.0 | 0.01
Modify | 0.00075486 | 0.00075486 | 0.00075486 | 0.0 | 0.01
Other | | 0.0003242 | | | 0.00
Nlocal: 256.000 ave 256 max 256 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 1431.00 ave 1431 max 1431 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 0.00000 ave 0 max 0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 19968.0 ave 19968 max 19968 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 19968
Ave neighs/atom = 78.000000
Neighbor list builds = 0
Dangerous builds = 0
Total wall time: 0:00:12

134
examples/mliap/log.03Mar21.mliap.nn.Cu.g++.4 Normal file
View File

@ -0,0 +1,134 @@
LAMMPS (10 Feb 2021)
using 1 OpenMP thread(s) per MPI task
# Demonstrate MLIAP interface to SNAP nn Cu potential
# Initialize simulation
variable nsteps index 100
variable nrep equal 4
variable a equal 3.621262
units metal
# generate the box and atom positions using a FCC lattice
variable nx equal ${nrep}
variable nx equal 4
variable ny equal ${nrep}
variable ny equal 4
variable nz equal ${nrep}
variable nz equal 4
boundary p p p
lattice fcc $a
lattice fcc 3.621262
Lattice spacing in x,y,z = 3.6212620 3.6212620 3.6212620
region box block 0 ${nx} 0 ${ny} 0 ${nz}
region box block 0 4 0 ${ny} 0 ${nz}
region box block 0 4 0 4 0 ${nz}
region box block 0 4 0 4 0 4
create_box 1 box
Created orthogonal box = (0.0000000 0.0000000 0.0000000) to (14.485048 14.485048 14.485048)
1 by 2 by 2 MPI processor grid
create_atoms 1 box
Created 256 atoms
create_atoms CPU = 0.000 seconds
mass 1 63.546
# choose potential
pair_style mliap model nn Cu.nn.mliap.model descriptor sna Cu.snap.mliap.descriptor
Reading potential file Cu.nn.mliap.model with DATE: 2021-03-03
Reading potential file Cu.snap.mliap.descriptor with DATE: 2021-03-03
SNAP keyword rcutfac 5.0
SNAP keyword twojmax 6
SNAP keyword nelems 1
SNAP keyword elems Cu
SNAP keyword radelems 0.5
SNAP keyword welems 1.0
SNAP keyword rfac0 0.99363
SNAP keyword rmin0 0
SNAP keyword bzeroflag 0
SNAP keyword switchflag 1
pair_coeff * * Cu
# Setup output
compute eatom all pe/atom
compute energy all reduce sum c_eatom
compute satom all stress/atom NULL
compute str all reduce sum c_satom[1] c_satom[2] c_satom[3]
variable press equal (c_str[1]+c_str[2]+c_str[3])/(3*vol)
thermo_style custom step temp epair c_energy etotal press v_press
thermo 10
thermo_modify norm yes
# Set up NVE run
timestep 0.5e-3
neighbor 1.0 bin
neigh_modify once no every 1 delay 0 check yes
# Run MD
velocity all create 300.0 4928459 loop geom
fix 1 all nve
run ${nsteps}
run 100
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 6
ghost atom cutoff = 6
binsize = 3, bins = 5 5 5
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) = 30.85 | 30.85 | 30.85 Mbytes
Step Temp E_pair c_energy TotEng Press v_press
0 300 -4.0935581 -4.0935581 -4.0549315 -2916.7505 2916.7505
10 294.19698 -4.0927997 -4.0927997 -4.0549203 -2653.7504 2653.7504
20 277.15991 -4.0905816 -4.0905816 -4.0548958 -1875.7718 1875.7718
30 250.30462 -4.0870777 -4.0870777 -4.0548496 -618.45378 618.45378
40 215.94455 -4.0825861 -4.0825861 -4.0547821 1014.6151 -1014.6151
50 177.29307 -4.0775299 -4.0775299 -4.0547025 2846.7675 -2846.7675
60 138.26899 -4.072429 -4.072429 -4.0546261 4698.9654 -4698.9654
70 103.04321 -4.0678341 -4.0678341 -4.0545667 6401.8077 -6401.8077
80 75.426905 -4.0642423 -4.0642423 -4.0545307 7766.2217 -7766.2217
90 58.150738 -4.0620015 -4.0620015 -4.0545143 8668.0598 -8668.0598
100 52.301012 -4.0612408 -4.0612408 -4.0545067 9049.3141 -9049.3141
Loop time of 3.22769 on 4 procs for 100 steps with 256 atoms
Performance: 1.338 ns/day, 17.932 hours/ns, 30.982 timesteps/s
100.0% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 3.19 | 3.2037 | 3.2196 | 0.7 | 99.26
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.0061415 | 0.022081 | 0.035855 | 8.9 | 0.68
Output | 0.00093221 | 0.00094824 | 0.00098861 | 0.0 | 0.03
Modify | 0.0002173 | 0.00022773 | 0.00024477 | 0.0 | 0.01
Other | | 0.0006805 | | | 0.02
Nlocal: 64.0000 ave 64 max 64 min
Histogram: 4 0 0 0 0 0 0 0 0 0
Nghost: 843.000 ave 843 max 843 min
Histogram: 4 0 0 0 0 0 0 0 0 0
Neighs: 0.00000 ave 0 max 0 min
Histogram: 4 0 0 0 0 0 0 0 0 0
FullNghs: 4992.00 ave 4992 max 4992 min
Histogram: 4 0 0 0 0 0 0 0 0 0
Total # of neighbors = 19968
Ave neighs/atom = 78.000000
Neighbor list builds = 0
Dangerous builds = 0
Total wall time: 0:00:03

157
examples/mliap/log.03Mar21.mliap.nn.Ta06A.g++.1 Normal file
View File

@ -0,0 +1,157 @@
LAMMPS (10 Feb 2021)
using 1 OpenMP thread(s) per MPI task
# Demonstrate MLIAP interface to linear SNAP potential
# Initialize simulation
variable nsteps index 100
variable nrep equal 4
variable a equal 3.316
units metal
# generate the box and atom positions using a BCC lattice
variable nx equal ${nrep}
variable nx equal 4
variable ny equal ${nrep}
variable ny equal 4
variable nz equal ${nrep}
variable nz equal 4
boundary p p p
lattice bcc $a
lattice bcc 3.316
Lattice spacing in x,y,z = 3.3160000 3.3160000 3.3160000
region box block 0 ${nx} 0 ${ny} 0 ${nz}
region box block 0 4 0 ${ny} 0 ${nz}
region box block 0 4 0 4 0 ${nz}
region box block 0 4 0 4 0 4
create_box 1 box
Created orthogonal box = (0.0000000 0.0000000 0.0000000) to (13.264000 13.264000 13.264000)
1 by 1 by 1 MPI processor grid
create_atoms 1 box
Created 128 atoms
create_atoms CPU = 0.000 seconds
mass 1 180.88
# choose potential
include Ta06A.nn.mliap
# 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 nn Ta06A.nn.mliap.model descriptor sna Ta06A.mliap.descriptor
pair_style hybrid/overlay zbl 4 ${zblcutouter} mliap model nn Ta06A.nn.mliap.model descriptor sna Ta06A.mliap.descriptor
pair_style hybrid/overlay zbl 4 4.8 mliap model nn Ta06A.nn.mliap.model descriptor sna Ta06A.mliap.descriptor
Reading potential file Ta06A.nn.mliap.model with DATE: 2021-03-03
Reading potential file Ta06A.mliap.descriptor with DATE: 2014-09-05
SNAP keyword rcutfac 4.67637
SNAP keyword twojmax 6
SNAP keyword nelems 1
SNAP keyword elems Ta
SNAP keyword radelems 0.5
SNAP keyword welems 1
SNAP keyword rfac0 0.99363
SNAP keyword rmin0 0
SNAP keyword bzeroflag 0
pair_coeff 1 1 zbl ${zblz} ${zblz}
pair_coeff 1 1 zbl 73 ${zblz}
pair_coeff 1 1 zbl 73 73
pair_coeff * * mliap Ta
# Setup output
compute eatom all pe/atom
compute energy all reduce sum c_eatom
compute satom all stress/atom NULL
compute str all reduce sum c_satom[1] c_satom[2] c_satom[3]
variable press equal (c_str[1]+c_str[2]+c_str[3])/(3*vol)
thermo_style custom step temp epair c_energy etotal press v_press
thermo 10
thermo_modify norm yes
# Set up NVE run
timestep 0.5e-3
neighbor 1.0 bin
neigh_modify once no every 1 delay 0 check yes
# Run MD
velocity all create 300.0 4928459 loop geom
fix 1 all nve
run ${nsteps}
run 100
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 5.8
ghost atom cutoff = 5.8
binsize = 2.9, bins = 5 5 5
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair zbl, perpetual, half/full from (2)
attributes: half, newton on
pair build: halffull/newton
stencil: none
bin: none
(2) pair mliap, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 18.15 | 18.15 | 18.15 Mbytes
Step Temp E_pair c_energy TotEng Press v_press
0 300 -11.85157 -11.85157 -11.813095 2717.1661 -2717.1661
10 296.01467 -11.851059 -11.851059 -11.813095 2697.4796 -2697.4796
20 284.53666 -11.849587 -11.849587 -11.813095 2289.1527 -2289.1527
30 266.51577 -11.847275 -11.847275 -11.813095 1851.7131 -1851.7131
40 243.05007 -11.844266 -11.844266 -11.813095 1570.684 -1570.684
50 215.51032 -11.840734 -11.840734 -11.813094 1468.1899 -1468.1899
60 185.48331 -11.836883 -11.836883 -11.813094 1524.8757 -1524.8757
70 154.6736 -11.832931 -11.832931 -11.813094 1698.3351 -1698.3351
80 124.79303 -11.829099 -11.829099 -11.813094 1947.0715 -1947.0715
90 97.448054 -11.825592 -11.825592 -11.813094 2231.9563 -2231.9563
100 74.035418 -11.822589 -11.822589 -11.813094 2515.8526 -2515.8526
Loop time of 3.05207 on 1 procs for 100 steps with 128 atoms
Performance: 1.415 ns/day, 16.956 hours/ns, 32.765 timesteps/s
100.0% 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.0503 | 3.0503 | 3.0503 | 0.0 | 99.94
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.00072395 | 0.00072395 | 0.00072395 | 0.0 | 0.02
Output | 0.00050893 | 0.00050893 | 0.00050893 | 0.0 | 0.02
Modify | 0.00037375 | 0.00037375 | 0.00037375 | 0.0 | 0.01
Other | | 0.0001878 | | | 0.01
Nlocal: 128.000 ave 128 max 128 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 727.000 ave 727 max 727 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 3712.00 ave 3712 max 3712 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 7424.00 ave 7424 max 7424 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 7424
Ave neighs/atom = 58.000000
Neighbor list builds = 0
Dangerous builds = 0
Total wall time: 0:00:03

157
examples/mliap/log.03Mar21.mliap.nn.Ta06A.g++.4 Normal file
View File

@ -0,0 +1,157 @@
LAMMPS (10 Feb 2021)
using 1 OpenMP thread(s) per MPI task
# Demonstrate MLIAP interface to linear SNAP potential
# Initialize simulation
variable nsteps index 100
variable nrep equal 4
variable a equal 3.316
units metal
# generate the box and atom positions using a BCC lattice
variable nx equal ${nrep}
variable nx equal 4
variable ny equal ${nrep}
variable ny equal 4
variable nz equal ${nrep}
variable nz equal 4
boundary p p p
lattice bcc $a
lattice bcc 3.316
Lattice spacing in x,y,z = 3.3160000 3.3160000 3.3160000
region box block 0 ${nx} 0 ${ny} 0 ${nz}
region box block 0 4 0 ${ny} 0 ${nz}
region box block 0 4 0 4 0 ${nz}
region box block 0 4 0 4 0 4
create_box 1 box
Created orthogonal box = (0.0000000 0.0000000 0.0000000) to (13.264000 13.264000 13.264000)
1 by 2 by 2 MPI processor grid
create_atoms 1 box
Created 128 atoms
create_atoms CPU = 0.000 seconds
mass 1 180.88
# choose potential
include Ta06A.nn.mliap
# 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 nn Ta06A.nn.mliap.model descriptor sna Ta06A.mliap.descriptor
pair_style hybrid/overlay zbl 4 ${zblcutouter} mliap model nn Ta06A.nn.mliap.model descriptor sna Ta06A.mliap.descriptor
pair_style hybrid/overlay zbl 4 4.8 mliap model nn Ta06A.nn.mliap.model descriptor sna Ta06A.mliap.descriptor
Reading potential file Ta06A.nn.mliap.model with DATE: 2021-03-03
Reading potential file Ta06A.mliap.descriptor with DATE: 2014-09-05
SNAP keyword rcutfac 4.67637
SNAP keyword twojmax 6
SNAP keyword nelems 1
SNAP keyword elems Ta
SNAP keyword radelems 0.5
SNAP keyword welems 1
SNAP keyword rfac0 0.99363
SNAP keyword rmin0 0
SNAP keyword bzeroflag 0
pair_coeff 1 1 zbl ${zblz} ${zblz}
pair_coeff 1 1 zbl 73 ${zblz}
pair_coeff 1 1 zbl 73 73
pair_coeff * * mliap Ta
# Setup output
compute eatom all pe/atom
compute energy all reduce sum c_eatom
compute satom all stress/atom NULL
compute str all reduce sum c_satom[1] c_satom[2] c_satom[3]
variable press equal (c_str[1]+c_str[2]+c_str[3])/(3*vol)
thermo_style custom step temp epair c_energy etotal press v_press
thermo 10
thermo_modify norm yes
# Set up NVE run
timestep 0.5e-3
neighbor 1.0 bin
neigh_modify once no every 1 delay 0 check yes
# Run MD
velocity all create 300.0 4928459 loop geom
fix 1 all nve
run ${nsteps}
run 100
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 5.8
ghost atom cutoff = 5.8
binsize = 2.9, bins = 5 5 5
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair zbl, perpetual, half/full from (2)
attributes: half, newton on
pair build: halffull/newton
stencil: none
bin: none
(2) pair mliap, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 18.06 | 18.06 | 18.06 Mbytes
Step Temp E_pair c_energy TotEng Press v_press
0 300 -11.85157 -11.85157 -11.813095 2717.1661 -2717.1661
10 296.01467 -11.851059 -11.851059 -11.813095 2697.4796 -2697.4796
20 284.53666 -11.849587 -11.849587 -11.813095 2289.1527 -2289.1527
30 266.51577 -11.847275 -11.847275 -11.813095 1851.7131 -1851.7131
40 243.05007 -11.844266 -11.844266 -11.813095 1570.684 -1570.684
50 215.51032 -11.840734 -11.840734 -11.813094 1468.1899 -1468.1899
60 185.48331 -11.836883 -11.836883 -11.813094 1524.8757 -1524.8757
70 154.6736 -11.832931 -11.832931 -11.813094 1698.3351 -1698.3351
80 124.79303 -11.829099 -11.829099 -11.813094 1947.0715 -1947.0715
90 97.448054 -11.825592 -11.825592 -11.813094 2231.9563 -2231.9563
100 74.035418 -11.822589 -11.822589 -11.813094 2515.8526 -2515.8526
Loop time of 0.858301 on 4 procs for 100 steps with 128 atoms
Performance: 5.033 ns/day, 4.768 hours/ns, 116.509 timesteps/s
100.2% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.75351 | 0.79999 | 0.85429 | 4.2 | 93.21
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.0026857 | 0.057019 | 0.10354 | 15.8 | 6.64
Output | 0.00078796 | 0.00080864 | 0.00085671 | 0.0 | 0.09
Modify | 0.00011437 | 0.00012275 | 0.00013339 | 0.0 | 0.01
Other | | 0.0003564 | | | 0.04
Nlocal: 32.0000 ave 32 max 32 min
Histogram: 4 0 0 0 0 0 0 0 0 0
Nghost: 431.000 ave 431 max 431 min
Histogram: 4 0 0 0 0 0 0 0 0 0
Neighs: 928.000 ave 928 max 928 min
Histogram: 4 0 0 0 0 0 0 0 0 0
FullNghs: 1856.00 ave 1856 max 1856 min
Histogram: 4 0 0 0 0 0 0 0 0 0
Total # of neighbors = 7424
Ave neighs/atom = 58.000000
Neighbor list builds = 0
Dangerous builds = 0
Total wall time: 0:00:00

32
examples/mliap/log.04Dec20.mliap.pytorch.Ta06A.g++.1
View File

@ -1,5 +1,5 @@
LAMMPS (30 Nov 2020) LAMMPS (10 Feb 2021)
using 48 OpenMP thread(s) per MPI task using 1 OpenMP thread(s) per MPI task
# Demonstrate MLIAP/PyTorch interface to linear SNAP potential # Demonstrate MLIAP/PyTorch interface to linear SNAP potential
# Initialize simulation # Initialize simulation
@ -32,13 +32,13 @@ Created orthogonal box = (0.0000000 0.0000000 0.0000000) to (13.264000 13.264000
1 by 1 by 1 MPI processor grid 1 by 1 by 1 MPI processor grid
create_atoms 1 box create_atoms 1 box
Created 128 atoms Created 128 atoms
create_atoms CPU = 0.002 seconds create_atoms CPU = 0.000 seconds
mass 1 180.88 mass 1 180.88
# choose potential # choose potential
include Ta06A.mliap.pytorch include Ta06A.pytorch.mliap
# DATE: 2014-09-05 UNITS: metal CONTRIBUTOR: Aidan Thompson athomps@sandia.gov CITATION: Thompson, Swiler, Trott, Foiles and Tucker, arxiv.org, 1409.3880 (2014) # 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 # Definition of SNAP potential Ta_Cand06A
@ -50,9 +50,9 @@ variable zblz equal 73
# Specify hybrid with SNAP, ZBL # Specify hybrid with SNAP, ZBL
pair_style hybrid/overlay zbl ${zblcutinner} ${zblcutouter} mliap model mliappy Ta06A.mliap.pytorch.model.pkl descriptor sna Ta06A.mliap.descriptor pair_style hybrid/overlay zbl ${zblcutinner} ${zblcutouter} mliap model mliappy Ta06A.mliap.pytorch.model.pt descriptor sna Ta06A.mliap.descriptor
pair_style hybrid/overlay zbl 4 ${zblcutouter} mliap model mliappy Ta06A.mliap.pytorch.model.pkl descriptor sna Ta06A.mliap.descriptor pair_style hybrid/overlay zbl 4 ${zblcutouter} mliap model mliappy Ta06A.mliap.pytorch.model.pt descriptor sna Ta06A.mliap.descriptor
pair_style hybrid/overlay zbl 4 4.8 mliap model mliappy Ta06A.mliap.pytorch.model.pkl descriptor sna Ta06A.mliap.descriptor pair_style hybrid/overlay zbl 4 4.8 mliap model mliappy Ta06A.mliap.pytorch.model.pt descriptor sna Ta06A.mliap.descriptor
Loading python model complete. Loading python model complete.
Reading potential file Ta06A.mliap.descriptor with DATE: 2014-09-05 Reading potential file Ta06A.mliap.descriptor with DATE: 2014-09-05
SNAP keyword rcutfac 4.67637 SNAP keyword rcutfac 4.67637
@ -112,7 +112,7 @@ Neighbor list info ...
pair build: full/bin/atomonly pair build: full/bin/atomonly
stencil: full/bin/3d stencil: full/bin/3d
bin: standard bin: standard
Per MPI rank memory allocation (min/avg/max) = 159.8 | 159.8 | 159.8 Mbytes Per MPI rank memory allocation (min/avg/max) = 18.15 | 18.15 | 18.15 Mbytes
Step Temp E_pair c_energy TotEng Press v_press Step Temp E_pair c_energy TotEng Press v_press
0 300 -11.85157 -11.85157 -11.813095 2717.1661 -2717.1661 0 300 -11.85157 -11.85157 -11.813095 2717.1661 -2717.1661
10 296.01467 -11.851059 -11.851059 -11.813095 2697.4796 -2697.4796 10 296.01467 -11.851059 -11.851059 -11.813095 2697.4796 -2697.4796
@ -125,20 +125,20 @@ Step Temp E_pair c_energy TotEng Press v_press
80 124.79303 -11.829099 -11.829099 -11.813094 1947.0715 -1947.0715 80 124.79303 -11.829099 -11.829099 -11.813094 1947.0715 -1947.0715
90 97.448054 -11.825592 -11.825592 -11.813094 2231.9563 -2231.9563 90 97.448054 -11.825592 -11.825592 -11.813094 2231.9563 -2231.9563
100 74.035418 -11.822589 -11.822589 -11.813094 2515.8526 -2515.8526 100 74.035418 -11.822589 -11.822589 -11.813094 2515.8526 -2515.8526
Loop time of 2.00236 on 48 procs for 100 steps with 128 atoms Loop time of 2.0472 on 1 procs for 100 steps with 128 atoms
Performance: 2.157 ns/day, 11.124 hours/ns, 49.941 timesteps/s Performance: 2.110 ns/day, 11.373 hours/ns, 48.847 timesteps/s
288.8% CPU use with 1 MPI tasks x 48 OpenMP threads 100.2% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown: MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total Section | min time | avg time | max time |%varavg| %total
--------------------------------------------------------------- ---------------------------------------------------------------
Pair | 1.9998 | 1.9998 | 1.9998 | 0.0 | 99.87 Pair | 2.0443 | 2.0443 | 2.0443 | 0.0 | 99.86
Neigh | 0 | 0 | 0 | 0.0 | 0.00 Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.0011814 | 0.0011814 | 0.0011814 | 0.0 | 0.06 Comm | 0.001157 | 0.001157 | 0.001157 | 0.0 | 0.06
Output | 0.00059724 | 0.00059724 | 0.00059724 | 0.0 | 0.03 Output | 0.00080857 | 0.00080857 | 0.00080857 | 0.0 | 0.04
Modify | 0.00047352 | 0.00047352 | 0.00047352 | 0.0 | 0.02 Modify | 0.00053188 | 0.00053188 | 0.00053188 | 0.0 | 0.03
Other | | 0.0003468 | | | 0.02 Other | | 0.0004304 | | | 0.02
Nlocal: 128.000 ave 128 max 128 min Nlocal: 128.000 ave 128 max 128 min
Histogram: 1 0 0 0 0 0 0 0 0 0 Histogram: 1 0 0 0 0 0 0 0 0 0

34
examples/mliap/log.04Dec20.mliap.pytorch.Ta06A.g++.4
View File

@ -1,5 +1,5 @@
LAMMPS (30 Nov 2020) LAMMPS (10 Feb 2021)
using 48 OpenMP thread(s) per MPI task using 1 OpenMP thread(s) per MPI task
# Demonstrate MLIAP/PyTorch interface to linear SNAP potential # Demonstrate MLIAP/PyTorch interface to linear SNAP potential
# Initialize simulation # Initialize simulation
@ -32,13 +32,13 @@ Created orthogonal box = (0.0000000 0.0000000 0.0000000) to (13.264000 13.264000
1 by 2 by 2 MPI processor grid 1 by 2 by 2 MPI processor grid
create_atoms 1 box create_atoms 1 box
Created 128 atoms Created 128 atoms
create_atoms CPU = 0.002 seconds create_atoms CPU = 0.000 seconds
mass 1 180.88 mass 1 180.88
# choose potential # choose potential
include Ta06A.mliap.pytorch include Ta06A.pytorch.mliap
# DATE: 2014-09-05 UNITS: metal CONTRIBUTOR: Aidan Thompson athomps@sandia.gov CITATION: Thompson, Swiler, Trott, Foiles and Tucker, arxiv.org, 1409.3880 (2014) # 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 # Definition of SNAP potential Ta_Cand06A
@ -50,9 +50,9 @@ variable zblz equal 73
# Specify hybrid with SNAP, ZBL # Specify hybrid with SNAP, ZBL
pair_style hybrid/overlay zbl ${zblcutinner} ${zblcutouter} mliap model mliappy Ta06A.mliap.pytorch.model.pkl descriptor sna Ta06A.mliap.descriptor pair_style hybrid/overlay zbl ${zblcutinner} ${zblcutouter} mliap model mliappy Ta06A.mliap.pytorch.model.pt descriptor sna Ta06A.mliap.descriptor
pair_style hybrid/overlay zbl 4 ${zblcutouter} mliap model mliappy Ta06A.mliap.pytorch.model.pkl descriptor sna Ta06A.mliap.descriptor pair_style hybrid/overlay zbl 4 ${zblcutouter} mliap model mliappy Ta06A.mliap.pytorch.model.pt descriptor sna Ta06A.mliap.descriptor
pair_style hybrid/overlay zbl 4 4.8 mliap model mliappy Ta06A.mliap.pytorch.model.pkl descriptor sna Ta06A.mliap.descriptor pair_style hybrid/overlay zbl 4 4.8 mliap model mliappy Ta06A.mliap.pytorch.model.pt descriptor sna Ta06A.mliap.descriptor
Loading python model complete. Loading python model complete.
Reading potential file Ta06A.mliap.descriptor with DATE: 2014-09-05 Reading potential file Ta06A.mliap.descriptor with DATE: 2014-09-05
SNAP keyword rcutfac 4.67637 SNAP keyword rcutfac 4.67637
@ -112,7 +112,7 @@ Neighbor list info ...
pair build: full/bin/atomonly pair build: full/bin/atomonly
stencil: full/bin/3d stencil: full/bin/3d
bin: standard bin: standard
Per MPI rank memory allocation (min/avg/max) = 159.7 | 159.7 | 159.7 Mbytes Per MPI rank memory allocation (min/avg/max) = 18.06 | 18.06 | 18.06 Mbytes
Step Temp E_pair c_energy TotEng Press v_press Step Temp E_pair c_energy TotEng Press v_press
0 300 -11.85157 -11.85157 -11.813095 2717.1661 -2717.1661 0 300 -11.85157 -11.85157 -11.813095 2717.1661 -2717.1661
10 296.01467 -11.851059 -11.851059 -11.813095 2697.4796 -2697.4796 10 296.01467 -11.851059 -11.851059 -11.813095 2697.4796 -2697.4796
@ -125,20 +125,20 @@ Step Temp E_pair c_energy TotEng Press v_press
80 124.79303 -11.829099 -11.829099 -11.813094 1947.0715 -1947.0715 80 124.79303 -11.829099 -11.829099 -11.813094 1947.0715 -1947.0715
90 97.448054 -11.825592 -11.825592 -11.813094 2231.9563 -2231.9563 90 97.448054 -11.825592 -11.825592 -11.813094 2231.9563 -2231.9563
100 74.035418 -11.822589 -11.822589 -11.813094 2515.8526 -2515.8526 100 74.035418 -11.822589 -11.822589 -11.813094 2515.8526 -2515.8526
Loop time of 0.562802 on 192 procs for 100 steps with 128 atoms Loop time of 0.572344 on 4 procs for 100 steps with 128 atoms
Performance: 7.676 ns/day, 3.127 hours/ns, 177.682 timesteps/s Performance: 7.548 ns/day, 3.180 hours/ns, 174.720 timesteps/s
99.7% CPU use with 4 MPI tasks x 48 OpenMP threads 100.6% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown: MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total Section | min time | avg time | max time |%varavg| %total
--------------------------------------------------------------- ---------------------------------------------------------------
Pair | 0.53583 | 0.54622 | 0.55401 | 0.9 | 97.05 Pair | 0.51796 | 0.5401 | 0.56608 | 2.4 | 94.37
Neigh | 0 | 0 | 0 | 0.0 | 0.00 Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.0071442 | 0.01491 | 0.025289 | 5.4 | 2.65 Comm | 0.0044844 | 0.030393 | 0.052464 | 10.2 | 5.31
Output | 0.00092525 | 0.00095771 | 0.0010166 | 0.0 | 0.17 Output | 0.001 | 0.0010351 | 0.0011083 | 0.1 | 0.18
Modify | 0.00014479 | 0.00015043 | 0.00015893 | 0.0 | 0.03 Modify | 0.0001468 | 0.00016417 | 0.00017977 | 0.0 | 0.03
Other | | 0.0005624 | | | 0.10 Other | | 0.0006567 | | | 0.11
Nlocal: 32.0000 ave 32 max 32 min Nlocal: 32.0000 ave 32 max 32 min
Histogram: 4 0 0 0 0 0 0 0 0 0 Histogram: 4 0 0 0 0 0 0 0 0 0
@ -154,4 +154,4 @@ Ave neighs/atom = 58.000000
Neighbor list builds = 0 Neighbor list builds = 0
Dangerous builds = 0 Dangerous builds = 0
Total wall time: 0:00:02 Total wall time: 0:00:05

154
examples/mliap/log.snap.Ta06A.ref
View File

@ -1,154 +0,0 @@
LAMMPS (19 Mar 2020)
# Demonstrate SNAP Ta potential
# Initialize simulation
variable nsteps index 100
variable nrep equal 4
variable a equal 3.316
units metal
# generate the box and atom positions using a BCC lattice
variable nx equal ${nrep}
variable nx equal 4
variable ny equal ${nrep}
variable ny equal 4
variable nz equal ${nrep}
variable nz equal 4
boundary p p p
lattice bcc $a
lattice bcc 3.316
Lattice spacing in x,y,z = 3.316 3.316 3.316
region box block 0 ${nx} 0 ${ny} 0 ${nz}
region box block 0 4 0 ${ny} 0 ${nz}
region box block 0 4 0 4 0 ${nz}
region box block 0 4 0 4 0 4
create_box 1 box
Created orthogonal box = (0 0 0) to (13.264 13.264 13.264)
1 by 1 by 1 MPI processor grid
create_atoms 1 box
Created 128 atoms
create_atoms CPU = 0.000254 secs
mass 1 180.88
# choose potential
include Ta06A.snap
# DATE: 2014-09-05 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} snap
pair_style hybrid/overlay zbl 4 ${zblcutouter} snap
pair_style hybrid/overlay zbl 4 4.8 snap
pair_coeff 1 1 zbl ${zblz} ${zblz}
pair_coeff 1 1 zbl 73 ${zblz}
pair_coeff 1 1 zbl 73 73
pair_coeff * * snap Ta06A.snapcoeff Ta06A.snapparam Ta
Reading potential file Ta06A.snapcoeff with DATE: 2014-09-05
SNAP Element = Ta, Radius 0.5, Weight 1
Reading potential file Ta06A.snapparam with DATE: 2014-09-05
SNAP keyword rcutfac 4.67637
SNAP keyword twojmax 6
SNAP keyword rfac0 0.99363
SNAP keyword rmin0 0
SNAP keyword bzeroflag 0
SNAP keyword quadraticflag 0
# 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
fix 1 all nve
run ${nsteps}
run 100
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 5.8
ghost atom cutoff = 5.8
binsize = 2.9, bins = 5 5 5
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair zbl, perpetual, half/full from (2)
attributes: half, newton on
pair build: halffull/newton
stencil: none
bin: none
(2) pair snap, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 6.591 | 6.591 | 6.591 Mbytes
Step Temp E_pair c_energy TotEng Press v_press
0 300 -11.85157 -11.85157 -11.813095 2717.1661 -2717.1661
10 295.96579 -11.851053 -11.851053 -11.813095 2696.1559 -2696.1559
20 284.32535 -11.84956 -11.84956 -11.813095 2301.3713 -2301.3713
30 266.04602 -11.847215 -11.847215 -11.813095 1832.1745 -1832.1745
40 242.2862 -11.844168 -11.844168 -11.813095 1492.6765 -1492.6765
50 214.48968 -11.840603 -11.840603 -11.813094 1312.8908 -1312.8908
60 184.32523 -11.836734 -11.836734 -11.813094 1284.582 -1284.582
70 153.58055 -11.832791 -11.832791 -11.813094 1374.4457 -1374.4457
80 124.04276 -11.829003 -11.829003 -11.813094 1537.703 -1537.703
90 97.37622 -11.825582 -11.825582 -11.813094 1734.9662 -1734.9662
100 75.007873 -11.822714 -11.822714 -11.813094 1930.8005 -1930.8005
Loop time of 0.995328 on 1 procs for 100 steps with 128 atoms
Performance: 4.340 ns/day, 5.530 hours/ns, 100.469 timesteps/s
99.5% CPU use with 1 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.99426 | 0.99426 | 0.99426 | 0.0 | 99.89
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.000305 | 0.000305 | 0.000305 | 0.0 | 0.03
Output | 0.000413 | 0.000413 | 0.000413 | 0.0 | 0.04
Modify | 0.000159 | 0.000159 | 0.000159 | 0.0 | 0.02
Other | | 0.000191 | | | 0.02
Nlocal: 128 ave 128 max 128 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 727 ave 727 max 727 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 3712 ave 3712 max 3712 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 7424 ave 7424 max 7424 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 7424
Ave neighs/atom = 58
Neighbor list builds = 0
Dangerous builds = 0
Total wall time: 0:00:01

41
potentials/Cu.nn.mliap.model Normal file
View File

@ -0,0 +1,41 @@
# DATE: 2021-03-03 UNITS: metal CONTRIBUTOR: Pedro Antonio Santos Florez pedroantoniosantosf@gmail.com CITATION: none
# Neural networks weights generated in PyXtal_FF
# total_species nparams
1 71
# NET ndescriptors nlayers activation_func for layer 1, number of nodes for layer 1, ..., ..., ....,
NET 30 3 tanh 2 tanh 2 linear 1
# Scale0
98.5564567301067 0.1817576370006 -5.0513436999162 0.0453581866073 -4.4382278670569
-2.6985391112026 -4.9445068940273 7.8840372041149 -12.8729456412560 -2.9748781222888
-1.7702149043285 -4.6911979511282 -1.8167140506334 67.7067032503211 -12.5675074252146
-10.3889411514447 -2.6194029974308 1.8750503488944 -5.0934919399232 4.2037760919263
20.6488797181668 -8.2381178902428 -3.1311436510232 -1.1570117080256 2.1652328079269
0.0854872341687 73.5970916724043 -5.3230176040510 -5.3409170596296 -6.4491186818574
# Scale1
783.0204659783215 33.6453931227955 7.9068346594951 32.2509031640463 5.6624925894670
5.6511565657016 8.2175114613285 108.9434600935519 16.3162547782998 6.4485086884047
9.5767453287024 24.1576638495153 25.1723098340896 202.6741688060115 26.3843102888906
19.7937707236334 13.9016635198379 14.1399097533048 72.7078614093291 39.8358080634714
370.1745750951991 21.6864533954642 14.6554778261065 18.6456437491491 43.2946309680355
68.9168418942695 481.9819429675007 14.5585849595560 61.9026061535013 112.5968384702330
# Coefficients
-0.2503221860648 0.4607794046378 0.4080552184108 0.0043328332206 0.2778722163950
-0.0050555739235 0.0406748042003 -0.1188215011883 0.4669700426246 0.4050949681583
0.3004118858813 -0.3937378155782 0.5165632847373 0.1381866719908 -0.3996533578004
0.1330571447244 0.1124773980921 0.1548275555377 -0.6221152970188 -0.7813275942231
-0.4320969824072 0.2045533448767 -0.1687644566507 -0.0424066840175 -0.1981208883936
-0.3406471991751 0.1272788573915 -0.1236220514264 0.6482374783971 -0.0590428211363
0.1326060235944 0.1890472518561 -0.5428894028577 0.0229420637813 0.0750842032804
-0.2443049940974 0.5265826743321 -0.1802773940492 0.0309737044359 -0.4987806786020
-0.6233577810656 -0.9910225261357 0.3430680860751 0.2918559145139 0.0666915892882
-0.2558026088950 0.0417769401952 -0.2582828540822 -0.3685977011673 0.4793473883461
1.0290706767456 0.0400544448686 -0.1474888333296 0.2403156153173 0.1050109023846
0.1413378531949 0.0907250056918 -0.1185878739508 0.1471942864979 -0.2829701766758
0.1137318187684 0.0476123438306 0.8040349098084 -0.5432518550512 0.2312003570779
0.6886872913686 -0.1247049894693 0.6848986783339 -2.6353147047309 -0.7173579584864
-1.0543148168289

21
potentials/Cu.snap.mliap.descriptor Normal file
View File

@ -0,0 +1,21 @@
# DATE: 2021-03-03 UNITS: metal CONTRIBUTOR: Pedro Antonio Santos Florez pedroantoniosantosf@gmail.com CITATION: none
# LAMMPS SNAP parameters for Cu
# required
rcutfac 5.0
twojmax 6
# elements
nelems 1
elems Cu
radelems 0.5
welems 1.0
# optional
rfac0 0.99363
rmin0 0
bzeroflag 0
switchflag 1

3
potentials/README
View File

@ -100,6 +100,9 @@ meam.spline modified EAM (MEAM) spline potential
meam.sw.spline modified EAM (MEAM) Stillinger-Weber spline potential meam.sw.spline modified EAM (MEAM) Stillinger-Weber spline potential
mesocnt mesoscopic carbon nanotube (CNT) potential mesocnt mesoscopic carbon nanotube (CNT) potential
mgpt model generalized pseudopotential theory (MGPT) potential mgpt model generalized pseudopotential theory (MGPT) potential
mliap MLIAP potential
mliap.descriptor MLIAP potential descriptor
mliap.model MLIAP potential model
nb3b.harmonic nonbonded 3-body harmonic potential nb3b.harmonic nonbonded 3-body harmonic potential
poly polymorphic 3-body potential poly polymorphic 3-body potential
reax ReaxFF potential (see README.reax for more info) reax ReaxFF potential (see README.reax for more info)

18
potentials/Ta06A.nn.mliap Normal file
View File

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

15
potentials/Ta06A.nn.mliap.model Normal file
View File

@ -0,0 +1,15 @@
# DATE: 2021-03-03 UNITS: metal CONTRIBUTOR: Aidan Thompson athomps@sandia.gov CITATION: Thompson, Swiler, Trott, Foiles and Tucker, arxiv.org, 1409.3880 (2014)
# LAMMPS SNAP_nn coefficients for Ta_Cand06A
# nelements nparams
1 31
# NET #ndescriptors #nlayers #activation_func for layer 1 #number of nodes for layer 1, ..., ..., ....,
NET 30 1 linear 1
#scale0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
#scale1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
#coeff
-2.92477
-0.01137 -0.00775 -0.04907 -0.15047 0.09157 0.05590 0.05785 -0.11615 -0.17122 -0.10583 0.03941 -0.11284 0.03939 -0.07331 -0.06582 -0.09341 -0.10587 -0.15497 0.04820 0.00205 0.00060 -0.04898 -0.05084 -0.03371 -0.01441 -0.01501 -0.00599 -0.06373 0.03965 0.01072

0
examples/mliap/Ta06A.mliap.pytorch → potentials/Ta06A.pytorch.mliap
View File

57
src/MLIAP/README
View File

@ -1,57 +0,0 @@
This package provides a general interface to families of
machine-learning interatomic potentials (MLIAPs). This interface consists
of a mliap pair style and a mliap compute.
The mliap pair style is used when running simulations with energies and
forces calculated by an MLIAP. The interface 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
supports the following models: *linear*, *quadratic*, and
*mliappy* (general Python interface to things like PyTorch, see below
for build instructions).
It currently supports only one class of descriptors,
*sna*, the SNAP descriptors, including the
linear, quadratic, and chem variants.
It is straightforward to add new descriptor and model
styles.
The mliap compute style provides gradients of the energy, force,
and stress tensor w.r.t. model parameters.
These are useful when training MLIAPs to match target data.
Any *model or *descriptor* that has been implemented for the
*mliap* pair style can also be accessed by the *mliap* compute.
In addition to the energy, force, and stress gradients, w.r.t.
each *model* parameter, the compute also calculates the energy,
force, and stress contributions from a user-specified
reference potential.
To see how this command
can be used within a Python workflow to train machine-learning interatomic
potentials, see the examples in FitSNAP https://github.com/FitSNAP/FitSNAP.
*Additional instructions for building MLIAP with Python support enabled*
The *mliappy* energy model requires that the MLIAP package
be compiled with Python support enabled.
This extension, written by Nick Lubbers (LANL),
provides a coupling to PyTorch and other Python modules.
This should be automatically
enabled by default if the prerequisite software is installed. It can be
enforced during CMake configuration by setting the variable
MLIAP_ENABLE_PYTHON=yes or for conventional build by adding -DMLIAP_PYTHON
to the LMP_INC variable in your makefile and running the cythonize script
on the .pyx file(s) copied to the src folder.
This requires to also install the PYTHON package and have the cython
(https://cython.org) software installed. During configuration/compilation
the cythonize script will be used to convert the provided .pyx file(s)
to C++ code. Please do not run the cythonize script in the src/MLIAP folder.
If you have done so by accident, you need to delete the generated .cpp and .h
file(s) in the src/MLIAP folder or there may be problems during compilation.
More information on building LAMMPS with this package is here:
https://lammps.sandia.gov/doc/Build_extras.html#mliap

41
src/MLIAP/README.md Normal file
View File

@ -0,0 +1,41 @@
This package provides a general interface to families of machine-learning interatomic potentials (MLIAPs).
This interface consists of a `mliap pair style` and a `mliap compute`.
The `mliap pair style` is used when running simulations with energies and
forces calculated by an MLIAP. The interface 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
supports the following models:
- `linear`,
- `quadratic`,
- `nn` (neural networks)
- `mliappy` (general Python interface to things like PyTorch).
It currently supports only one class of descriptors, `sna`, the SNAP descriptors.
It is straightforward to add new descriptor and model styles.
The `mliap compute` style provides gradients of the energy, force,
and stress tensor w.r.t. model parameters.
These are useful when training MLIAPs to match target data.
Any `model` or `descriptor` that has been implemented for the
`mliap` pair style can also be accessed by the `mliap` compute.
In addition to the energy, force, and stress gradients, w.r.t.
each `model` parameter, the compute also calculates the energy,
force, and stress contributions from a user-specified
reference potential.
## Generating the model files from the third-party packages
- To train the `linear` and `quadratic` models with the SNAP descriptors, see the examples in [FitSNAP](https://github.com/FitSNAP/FitSNAP).
- To train the `nn` model with the SNAP descriptors, check the examples in [PyXtal\_FF](https://github.com/qzhu2017/PyXtal_FF).
## Building MLIAP with Python support
The `mliappy` energy model requires that the MLIAP package be compiled with Python support enabled. This extension, written by Nick Lubbers (LANL), provides a coupling to PyTorch and other Python modules. This should be automatically enabled by default if the prerequisite software is installed. It can be enforced during CMake configuration by setting the variable `MLIAP_ENABLE_PYTHON=yes` or for conventional build by adding `-DMLIAP_PYTHON` to the `LMP_INC` variable in your makefile and running the cythonize script on the .pyx file(s) copied to the src folder.
This requires to also install the PYTHON package and have the [cython](https://cython.org) software installed. During configuration/compilation
the cythonize script will be used to convert the provided .pyx file(s) to C++ code. Please do not run the cythonize script in the src/MLIAP folder. If you have done so by accident, you need to delete the generated .cpp and .h file(s) in the src/MLIAP folder or there may be problems during compilation.
More information on building LAMMPS with this package is [here](https://lammps.sandia.gov/doc/Build_extras.html#mliap).

8
src/MLIAP/compute_mliap.cpp
View File

@ -215,7 +215,7 @@ void ComputeMLIAP::init_list(int /*id*/, NeighList *ptr)
void ComputeMLIAP::compute_array() void ComputeMLIAP::compute_array()
{ {
int ntotal = atom->nlocal + atom->nghost; int nall = atom->nlocal + atom->nghost;
invoked_array = update->ntimestep; invoked_array = update->ntimestep;
// clear global array // clear global array
@ -261,7 +261,7 @@ void ComputeMLIAP::compute_array()
for (int ielem = 0; ielem < data->nelements; ielem++) { for (int ielem = 0; ielem < data->nelements; ielem++) {
const int elemoffset = data->nparams*ielem; const int elemoffset = data->nparams*ielem;
for (int jparam = 0; jparam < data->nparams; jparam++) { for (int jparam = 0; jparam < data->nparams; jparam++) {
for (int i = 0; i < ntotal; i++) { for (int i = 0; i < nall; i++) {
double *gradforcei = data->gradforce[i]+elemoffset; double *gradforcei = data->gradforce[i]+elemoffset;
tagint irow = 3*(atom->tag[i]-1)+1; tagint irow = 3*(atom->tag[i]-1)+1;
mliaparray[irow][jparam+elemoffset] += gradforcei[jparam]; mliaparray[irow][jparam+elemoffset] += gradforcei[jparam];
@ -307,7 +307,7 @@ void ComputeMLIAP::compute_array()
// switch to Voigt notation // switch to Voigt notation
c_virial->compute_vector(); c_virial->compute_vector();
irow += 3*data->natoms_array; irow += 3*data->natoms;
mliaparrayall[irow++][lastcol] = c_virial->vector[0]; mliaparrayall[irow++][lastcol] = c_virial->vector[0];
mliaparrayall[irow++][lastcol] = c_virial->vector[1]; mliaparrayall[irow++][lastcol] = c_virial->vector[1];
mliaparrayall[irow++][lastcol] = c_virial->vector[2]; mliaparrayall[irow++][lastcol] = c_virial->vector[2];
@ -325,7 +325,7 @@ void ComputeMLIAP::compute_array()
void ComputeMLIAP::dbdotr_compute() void ComputeMLIAP::dbdotr_compute()
{ {
double **x = atom->x; double **x = atom->x;
int irow0 = 1+data->ndims_force*data->natoms_array; int irow0 = 1+data->ndims_force*data->natoms;
// sum over bispectrum contributions to forces // sum over bispectrum contributions to forces
// on all particles including ghosts // on all particles including ghosts

81
src/MLIAP/mliap_data.cpp
View File

@ -18,6 +18,7 @@
#include "mliap_data.h" #include "mliap_data.h"
#include "atom.h" #include "atom.h"
#include "error.h"
#include "memory.h" #include "memory.h"
#include "mliap_descriptor.h" #include "mliap_descriptor.h"
#include "mliap_model.h" #include "mliap_model.h"
@ -50,12 +51,18 @@ MLIAPData::MLIAPData(LAMMPS *lmp, int gradgradflag_in, int *map_in,
ndims_virial = 6; ndims_virial = 6;
yoffset = nparams*nelements; yoffset = nparams*nelements;
zoffset = 2*yoffset; zoffset = 2*yoffset;
natoms_array = atom->natoms; natoms = atom->natoms;
size_array_rows = 1+ndims_force*natoms_array+ndims_virial;
// must check before assigning bigint expression to regular int
if (1+ndims_force*natoms+ndims_virial > MAXSMALLINT)
error->all(FLERR,"Too many atoms for MLIAP package");
size_array_rows = 1+ndims_force*natoms+ndims_virial;
size_array_cols = nparams*nelements+1; size_array_cols = nparams*nelements+1;
size_gradforce = ndims_force*nparams*nelements; size_gradforce = ndims_force*nparams*nelements;
natoms_max = 0; nlistatoms_max = 0;
natomneigh_max = 0; natomneigh_max = 0;
nneigh_max = 0; nneigh_max = 0;
nmax = 0; nmax = 0;
@ -102,6 +109,7 @@ void MLIAPData::generate_neighdata(NeighList* list_in, int eflag_in, int vflag_i
double **x = atom->x; double **x = atom->x;
int *type = atom->type; int *type = atom->type;
int *ilist = list->ilist;
int *numneigh = list->numneigh; int *numneigh = list->numneigh;
int **firstneigh = list->firstneigh; int **firstneigh = list->firstneigh;
@ -115,39 +123,38 @@ void MLIAPData::generate_neighdata(NeighList* list_in, int eflag_in, int vflag_i
// clear gradforce array // clear gradforce array
int ntotal = atom->nlocal + atom->nghost; int nall = atom->nlocal + atom->nghost;
for (int i = 0; i < ntotal; i++) for (int i = 0; i < nall; i++)
for (int j = 0; j < size_gradforce; j++) { for (int j = 0; j < size_gradforce; j++) {
gradforce[i][j] = 0.0; gradforce[i][j] = 0.0;
} }
// grow arrays if necessary
nlistatoms = list->inum;
if (nlistatoms_max < nlistatoms) {
memory->grow(betas,nlistatoms,ndescriptors,"MLIAPData:betas");
memory->grow(descriptors,nlistatoms,ndescriptors,"MLIAPData:descriptors");
memory->grow(eatoms,nlistatoms,"MLIAPData:eatoms");
nlistatoms_max = nlistatoms;
}
// grow gamma arrays if necessary // grow gamma arrays if necessary
if (gradgradflag == 1) { if (gradgradflag == 1) {
const int natomgamma = list->inum; if (natomgamma_max < nlistatoms) {
if (natomgamma_max < natomgamma) { memory->grow(gamma_row_index,nlistatoms,gamma_nnz,"MLIAPData:gamma_row_index");
memory->grow(gamma_row_index,natomgamma,gamma_nnz,"MLIAPData:gamma_row_index"); memory->grow(gamma_col_index,nlistatoms,gamma_nnz,"MLIAPData:gamma_col_index");
memory->grow(gamma_col_index,natomgamma,gamma_nnz,"MLIAPData:gamma_col_index"); memory->grow(gamma,nlistatoms,gamma_nnz,"MLIAPData:gamma");
memory->grow(gamma,natomgamma,gamma_nnz,"MLIAPData:gamma"); natomgamma_max = nlistatoms;
natomgamma_max = natomgamma;
} }
} }
// grow arrays if necessary
natoms = list->inum;
if (natoms_max < natoms) {
memory->grow(betas,natoms,ndescriptors,"MLIAPData:betas");
memory->grow(descriptors,natoms,ndescriptors,"MLIAPData:descriptors");
memory->grow(eatoms,natoms,"MLIAPData:eatoms");
natoms_max = natoms;
}
grow_neigharrays(); grow_neigharrays();
int ij = 0; int ij = 0;
for (int ii = 0; ii < list->inum; ii++) { for (int ii = 0; ii < nlistatoms; ii++) {
const int i = list->ilist[ii]; const int i = ilist[ii];
const double xtmp = x[i][0]; const double xtmp = x[i][0];
const double ytmp = x[i][1]; const double ytmp = x[i][1];
@ -197,30 +204,24 @@ void MLIAPData::grow_neigharrays()
// grow neighbor atom arrays if necessary // grow neighbor atom arrays if necessary
const int natomneigh = list->inum; if (natomneigh_max < nlistatoms) {
if (natomneigh_max < natomneigh) { memory->grow(iatoms,nlistatoms,"MLIAPData:iatoms");
memory->grow(iatoms,natomneigh,"MLIAPData:iatoms"); memory->grow(ielems,nlistatoms,"MLIAPData:ielems");
memory->grow(ielems,natomneigh,"MLIAPData:ielems"); memory->grow(numneighs,nlistatoms,"MLIAPData:numneighs");
memory->grow(numneighs,natomneigh,"MLIAPData:numneighs"); natomneigh_max = nlistatoms;
natomneigh_max = natomneigh;
} }
// grow neighbor arrays if necessary // grow neighbor arrays if necessary
int *ilist = list->ilist;
int *numneigh = list->numneigh; int *numneigh = list->numneigh;
int **firstneigh = list->firstneigh; int **firstneigh = list->firstneigh;
int iilast = list->inum-1;
int ilast = list->ilist[iilast];
int upperbound = firstneigh[ilast] - firstneigh[0] + numneigh[ilast];
if (nneigh_max >= upperbound) return;
double **x = atom->x; double **x = atom->x;
int *type = atom->type; int *type = atom->type;
int nneigh = 0; int nneigh = 0;
for (int ii = 0; ii < list->inum; ii++) { for (int ii = 0; ii < nlistatoms; ii++) {
const int i = list->ilist[ii]; const int i = ilist[ii];
const double xtmp = x[i][0]; const double xtmp = x[i][0];
const double ytmp = x[i][1]; const double ytmp = x[i][1];
@ -272,9 +273,9 @@ double MLIAPData::memory_usage()
gamma_nnz*sizeof(double); // gamma gamma_nnz*sizeof(double); // gamma
} }
bytes += (double)natoms*ndescriptors*sizeof(int); // betas bytes += (double)nlistatoms*ndescriptors*sizeof(int); // betas
bytes += (double)natoms*ndescriptors*sizeof(int); // descriptors bytes += (double)nlistatoms*ndescriptors*sizeof(int); // descriptors
bytes += (double)natoms*sizeof(double); // eatoms bytes += (double)nlistatoms*sizeof(double); // eatoms
bytes += (double)natomneigh_max*sizeof(int); // iatoms bytes += (double)natomneigh_max*sizeof(int); // iatoms
bytes += (double)natomneigh_max*sizeof(int); // ielems bytes += (double)natomneigh_max*sizeof(int); // ielems

7
src/MLIAP/mliap_data.h
View File

@ -30,7 +30,8 @@ class MLIAPData : protected Pointers {
double memory_usage(); double memory_usage();
int size_array_rows, size_array_cols; int size_array_rows, size_array_cols;
int natoms_array, size_gradforce; int natoms;
int size_gradforce;
int yoffset, zoffset; int yoffset, zoffset;
int ndims_force, ndims_virial; int ndims_force, ndims_virial;
double **gradforce; double **gradforce;
@ -54,8 +55,8 @@ class MLIAPData : protected Pointers {
// data structures for mliap neighbor list // data structures for mliap neighbor list
// only neighbors strictly inside descriptor cutoff // only neighbors strictly inside descriptor cutoff
int natoms; // current number of atoms int nlistatoms; // current number of atoms in neighborlist
int natoms_max; // allocated size of descriptor array int nlistatoms_max; // allocated size of descriptor array
int natomneigh_max; // allocated size of atom neighbor arrays int natomneigh_max; // allocated size of atom neighbor arrays
int *numneighs; // neighbors count for each atom int *numneighs; // neighbors count for each atom
int *iatoms; // index of each atom int *iatoms; // index of each atom

8
src/MLIAP/mliap_descriptor_snap.cpp
View File

@ -78,7 +78,7 @@ MLIAPDescriptorSNAP::~MLIAPDescriptorSNAP()
void MLIAPDescriptorSNAP::compute_descriptors(class MLIAPData* data) void MLIAPDescriptorSNAP::compute_descriptors(class MLIAPData* data)
{ {
int ij = 0; int ij = 0;
for (int ii = 0; ii < data->natoms; ii++) { for (int ii = 0; ii < data->nlistatoms; ii++) {
const int ielem = data->ielems[ii]; const int ielem = data->ielems[ii];
// insure rij, inside, wj, and rcutij are of size jnum // insure rij, inside, wj, and rcutij are of size jnum
@ -130,7 +130,7 @@ void MLIAPDescriptorSNAP::compute_forces(class MLIAPData* data)
double **f = atom->f; double **f = atom->f;
int ij = 0; int ij = 0;
for (int ii = 0; ii < data->natoms; ii++) { for (int ii = 0; ii < data->nlistatoms; ii++) {
const int i = data->iatoms[ii]; const int i = data->iatoms[ii];
const int ielem = data->ielems[ii]; const int ielem = data->ielems[ii];
@ -204,7 +204,7 @@ void MLIAPDescriptorSNAP::compute_forces(class MLIAPData* data)
void MLIAPDescriptorSNAP::compute_force_gradients(class MLIAPData* data) void MLIAPDescriptorSNAP::compute_force_gradients(class MLIAPData* data)
{ {
int ij = 0; int ij = 0;
for (int ii = 0; ii < data->natoms; ii++) { for (int ii = 0; ii < data->nlistatoms; ii++) {
const int i = data->iatoms[ii]; const int i = data->iatoms[ii];
const int ielem = data->ielems[ii]; const int ielem = data->ielems[ii];
@ -279,7 +279,7 @@ void MLIAPDescriptorSNAP::compute_force_gradients(class MLIAPData* data)
void MLIAPDescriptorSNAP::compute_descriptor_gradients(class MLIAPData* data) void MLIAPDescriptorSNAP::compute_descriptor_gradients(class MLIAPData* data)
{ {
int ij = 0; int ij = 0;
for (int ii = 0; ii < data->natoms; ii++) { for (int ii = 0; ii < data->nlistatoms; ii++) {
const int ielem = data->ielems[ii]; const int ielem = data->ielems[ii];
// insure rij, inside, wj, and rcutij are of size jnum // insure rij, inside, wj, and rcutij are of size jnum

6
src/MLIAP/mliap_model_linear.cpp
View File

@ -57,7 +57,7 @@ void MLIAPModelLinear::compute_gradients(MLIAPData* data)
{ {
data->energy = 0.0; data->energy = 0.0;
for (int ii = 0; ii < data->natoms; ii++) { for (int ii = 0; ii < data->nlistatoms; ii++) {
const int ielem = data->ielems[ii]; const int ielem = data->ielems[ii];
double* coeffi = coeffelem[ielem]; double* coeffi = coeffelem[ielem];
@ -107,7 +107,7 @@ void MLIAPModelLinear::compute_gradgrads(class MLIAPData* data)
for (int l = 0; l < data->nelements*data->nparams; l++) for (int l = 0; l < data->nelements*data->nparams; l++)
data->egradient[l] = 0.0; data->egradient[l] = 0.0;
for (int ii = 0; ii < data->natoms; ii++) { for (int ii = 0; ii < data->nlistatoms; ii++) {
const int ielem = data->ielems[ii]; const int ielem = data->ielems[ii];
const int elemoffset = data->nparams*ielem; const int elemoffset = data->nparams*ielem;
@ -143,7 +143,7 @@ void MLIAPModelLinear::compute_force_gradients(class MLIAPData* data)
data->egradient[l] = 0.0; data->egradient[l] = 0.0;
int ij = 0; int ij = 0;
for (int ii = 0; ii < data->natoms; ii++) { for (int ii = 0; ii < data->nlistatoms; ii++) {
const int i = data->iatoms[ii]; const int i = data->iatoms[ii];
const int ielem = data->ielems[ii]; const int ielem = data->ielems[ii];
const int elemoffset = data->nparams*ielem; const int elemoffset = data->nparams*ielem;

415
src/MLIAP/mliap_model_nn.cpp Normal file
View File

@ -0,0 +1,415 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://lammps.sandia.gov/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Pedro Antonio Santos Flórez (UNLV)
------------------------------------------------------------------------- */
#include "mliap_model_nn.h"
#include "pair_mliap.h"
#include "mliap_data.h"
#include "error.h"
#include "comm.h"
#include "memory.h"
#include "tokenizer.h"
#include <cmath>
using namespace LAMMPS_NS;
#define MAXLINE 1024
/* ---------------------------------------------------------------------- */
MLIAPModelNN::MLIAPModelNN(LAMMPS* lmp, char* coefffilename) :
MLIAPModel(lmp, coefffilename)
{
coeffelem = nullptr;
nnodes = nullptr;
activation = nullptr;
scale = nullptr;
if (coefffilename) read_coeffs(coefffilename);
}
/* ---------------------------------------------------------------------- */
MLIAPModelNN::~MLIAPModelNN()
{
memory->destroy(coeffelem);
memory->destroy(nnodes);
memory->destroy(activation);
memory->destroy(scale);
}
/* ----------------------------------------------------------------------
get number of parameters
---------------------------------------------------------------------- */
int MLIAPModelNN::get_nparams()
{
if (nparams == 0)
if (ndescriptors == 0) error->all(FLERR,"ndescriptors not defined");
return nparams;
}
void MLIAPModelNN::read_coeffs(char *coefffilename)
{
// open coefficient file on proc 0
FILE *fpcoeff;
if (comm->me == 0) {
fpcoeff = utils::open_potential(coefffilename,lmp,nullptr);
if (fpcoeff == nullptr)
error->one(FLERR,fmt::format("Cannot open MLIAPModel coeff file {}: {}",
coefffilename,utils::getsyserror()));
}
char line[MAXLINE], *ptr, *tstr;
int eof = 0;
int n;
int nwords = 0;
while (nwords == 0) {
if (comm->me == 0) {
ptr = fgets(line,MAXLINE,fpcoeff);
if (ptr == nullptr) {
eof = 1;
fclose(fpcoeff);
} else n = strlen(line) + 1;
}
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof) break;
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
// strip comment, skip line if blank
if ((ptr = strchr(line,'#'))) *ptr = '\0';
nwords = utils::count_words(line);
}
if (nwords != 2)
error->all(FLERR,"Incorrect format in MLIAPModel coefficient file");
// words = ptrs to all words in line
// strip single and double quotes from words
try {
ValueTokenizer coeffs(line);
nelements = coeffs.next_int();
nparams = coeffs.next_int();
} catch (TokenizerException &e) {
error->all(FLERR,fmt::format("Incorrect format in MLIAPModel coefficient "
"file: {}",e.what()));
}
// set up coeff lists
memory->create(coeffelem,nelements,nparams,"mliap_snap_model:coeffelem");
int stats = 0;
int ielem = 0;
int l = 0;
while (1) {
if (comm->me == 0) {
ptr = fgets(line,MAXLINE,fpcoeff);
if (ptr == nullptr) {
eof = 1;
fclose(fpcoeff);
} else n = strlen(line) + 1;
}
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof) break;
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
// strip comment, skip line if blank
if ((ptr = strchr(line,'#'))) *ptr = '\0';
nwords = utils::trim_and_count_words(line);
if (nwords == 0) continue;
if (stats == 0) { // Header NET
tstr = strtok(line,"' \t\n\r\f");
if (strncmp(tstr, "NET", 3) != 0) error->all(FLERR,"Incorrect format in MLIAPModel coefficient file");
ndescriptors = atoi(strtok(nullptr,"' \t\n\r\f"));
nlayers = atoi(strtok(nullptr,"' \t\n\r\f"));
memory->create(activation,nlayers,"mliap_model:activation");
memory->create(nnodes,nlayers,"mliap_model:nnodes");
memory->create(scale,nelements,2,ndescriptors,"mliap_model:scale");
for (int ilayer = 0; ilayer < nlayers; ilayer++) {
tstr = strtok(NULL,"' \t\n\r\f");
nnodes[ilayer] = atoi(strtok(NULL,"' \t\n\r\f"));
if (strncmp(tstr, "linear", 6) == 0) activation[ilayer] = 0;
else if (strncmp(tstr, "sigmoid", 7) == 0) activation[ilayer] = 1;
else if (strncmp(tstr, "tanh", 4) == 0) activation[ilayer] = 2;
else if (strncmp(tstr, "relu", 4) == 0) activation[ilayer] = 3;
else activation[ilayer] = 4;
}
stats = 1;
} else if (stats == 1) {
scale[ielem][0][l] = atof(strtok(line,"' \t\n\r\f"));
for (int icoeff = 1; icoeff < nwords; icoeff++) {
scale[ielem][0][l+icoeff] = atof(strtok(nullptr,"' \t\n\r\f"));
}
l += nwords;
if (l == ndescriptors) {
stats = 2;
l = 0;
}
} else if (stats == 2) {
scale[ielem][1][l] = atof(strtok(line,"' \t\n\r\f"));
for (int icoeff = 1; icoeff < nwords; icoeff++) {
scale[ielem][1][l+icoeff] = atof(strtok(nullptr,"' \t\n\r\f"));
}
l += nwords;
if (l == ndescriptors) {
stats = 3;
l = 0;
}
// set up coeff lists
} else if (stats == 3) {
if (nwords > 30) error->all(FLERR,"Incorrect format in MLIAPModel coefficient file");
coeffelem[ielem][l] = atof(strtok(line,"' \t\n\r\f"));
for (int icoeff = 1; icoeff < nwords; icoeff++) {
coeffelem[ielem][l+icoeff] = atof(strtok(nullptr,"' \t\n\r\f"));
}
l += nwords;
if (l == nparams) {
stats = 1;
l = 0;
ielem++;
}
}
}
}
/* ----------------------------------------------------------------------
Calculate model gradients w.r.t descriptors
for each atom beta_i = dE(B_i)/dB_i
---------------------------------------------------------------------- */
void MLIAPModelNN::compute_gradients(MLIAPData* data)
{
data->energy = 0.0;
for (int ii = 0; ii < data->nlistatoms; ii++) {
const int ielem = data->ielems[ii];
const int nl = nlayers;
double* coeffi = coeffelem[ielem];
double** scalei = scale[ielem];
double **nodes, **dnodes, **bnodes;
nodes = new double*[nl];
dnodes = new double*[nl];
bnodes = new double*[nl];
for (int l=0; l<nl; ++l) {
nodes[l] = new double[nnodes[l]];
dnodes[l] = new double[nnodes[l]];
bnodes[l] = new double[nnodes[l]];
}
// forwardprop
// input - hidden1
for (int n=0; n < nnodes[0]; n++) {
nodes[0][n] = 0;
for (int icoeff = 0; icoeff < data->ndescriptors; icoeff++) {
nodes[0][n] += coeffi[n*((data->ndescriptors)+1)+icoeff+1] *
(data->descriptors[ii][icoeff] - scalei[0][icoeff]) /
scalei[1][icoeff];
}
if (activation[0] == 1) {
nodes[0][n] = sigm(nodes[0][n] +
coeffi[n*((data->ndescriptors)+1)],
dnodes[0][n]);
}
else if (activation[0] == 2) {
nodes[0][n] = tanh(nodes[0][n] +
coeffi[n*((data->ndescriptors)+1)],
dnodes[0][n]);
}
else if (activation[0] == 3) {
nodes[0][n] = relu(nodes[0][n] +
coeffi[n*((data->ndescriptors)+1)],
dnodes[0][n]);
}
else {
nodes[0][n] += coeffi[n*((data->ndescriptors)+1)];
dnodes[0][n] = 1;
}
}
// hidden~output
int k = 0;
if (nl > 1) {
k += ((data->ndescriptors)+1)*nnodes[0];
for (int l=1; l < nl; l++) {
for (int n=0; n < nnodes[l]; n++) {
nodes[l][n] = 0;
for (int j=0; j < nnodes[l-1]; j++) {
nodes[l][n] += coeffi[k+n*(nnodes[l-1]+1)+j+1] *
nodes[l-1][j];
}
if (activation[l] == 1) {
nodes[l][n] = sigm(nodes[l][n] +
coeffi[k+n*(nnodes[l-1]+1)],
dnodes[l][n]);
}
else if (activation[l] == 2) {
nodes[l][n] = tanh(nodes[l][n] +
coeffi[k+n*(nnodes[l-1]+1)],
dnodes[l][n]);
}
else if (activation[l] == 3) {
nodes[l][n] = relu(nodes[l][n] +
coeffi[k+n*(nnodes[l-1]+1)],
dnodes[l][n]);
}
else {
nodes[l][n] += coeffi[k+n*(nnodes[l-1]+1)];
dnodes[l][n] = 1;
}
}
k += (nnodes[l-1]+1)*nnodes[l];
}
}
// backwardprop
// output layer dnode initialized to 1.
for (int n=0; n<nnodes[nl-1]; n++) {
if (activation[nl-1] == 0) {
bnodes[nl-1][n] = 1;
}
else {
bnodes[nl-1][n] = dnodes[nl-1][n];
}
}
if (nl > 1) {
for (int l=nl-1; l>0; l--) {
k -= (nnodes[l-1]+1)*nnodes[l];
for (int n=0; n<nnodes[l-1]; n++) {
bnodes[l-1][n] = 0;
for (int j=0; j<nnodes[l]; j++) {
bnodes[l-1][n] += coeffi[k+j*(nnodes[l-1]+1)+n+1] *
bnodes[l][j];
}
if (activation[l-1] >= 1) {
bnodes[l-1][n] *= dnodes[l-1][n];
}
}
}
}
for (int icoeff = 0; icoeff < data->ndescriptors; icoeff++) {
data->betas[ii][icoeff] = 0;
for (int j=0; j<nnodes[0]; j++) {
data->betas[ii][icoeff] +=
coeffi[j*((data->ndescriptors)+1)+icoeff+1] *
bnodes[0][j];
}
data->betas[ii][icoeff] = data->betas[ii][icoeff]/scalei[1][icoeff];
}
if (data->eflag) {
// energy of atom I (E_i)
double etmp = nodes[nl-1][0];
data->energy += etmp;
data->eatoms[ii] = etmp;
}
// Deleting the variables
for (int n=0; n<nl; n++) {
delete nodes[n];
delete dnodes[n];
delete bnodes[n];
}
delete nodes;
delete dnodes;
delete bnodes;
}
}
/* ----------------------------------------------------------------------
Calculate model double gradients w.r.t descriptors and parameters
for each atom energy gamma_lk = d2E(B)/dB_k/dsigma_l,
where sigma_l is a parameter, B_k a descriptor,
and atom subscript i is omitted
gamma is in CSR format:
nnz = number of non-zero values
gamma_row_index[inz] = l indices, 0 <= l < nparams
gamma_col_indexiinz] = k indices, 0 <= k < ndescriptors
gamma[i][inz] = non-zero values, 0 <= inz < nnz
egradient is derivative of energy w.r.t. parameters
---------------------------------------------------------------------- */
void MLIAPModelNN::compute_gradgrads(class MLIAPData* data)
{
error->all(FLERR,"compute_gradgrads not implemented");
}
/* ----------------------------------------------------------------------
calculate gradients of forces w.r.t. parameters
egradient is derivative of energy w.r.t. parameters
---------------------------------------------------------------------- */
void MLIAPModelNN::compute_force_gradients(class MLIAPData* data)
{
error->all(FLERR,"compute_force_gradients not implemented");
}
/* ----------------------------------------------------------------------
count the number of non-zero entries in gamma matrix
---------------------------------------------------------------------- */
int MLIAPModelNN::get_gamma_nnz(class MLIAPData* data)
{
// todo: get_gamma_nnz
return 0;
}
double MLIAPModelNN::memory_usage()
{
double bytes = 0;
bytes += (double)nelements*nparams*sizeof(double); // coeffelem
bytes += (double)nelements*2*ndescriptors*sizeof(double); // scale
bytes += (int)nlayers*sizeof(int); // nnodes
bytes += (int)nlayers*sizeof(int); // activation
return bytes;
}

68
src/MLIAP/mliap_model_nn.h Normal file
View File

@ -0,0 +1,68 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifndef LMP_MLIAP_MODEL_NN_H
#define LMP_MLIAP_MODEL_NN_H
#include "mliap_model.h"
namespace LAMMPS_NS {
class MLIAPModelNN : public MLIAPModel {
public:
MLIAPModelNN(LAMMPS*, char* = nullptr);
~MLIAPModelNN();
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();
int nlayers; // number of layers per element
protected:
int *activation; // activation functions
int *nnodes; // number of nodes per layer
double ***scale; // element scale values
double **coeffelem; // element coefficients
virtual void read_coeffs(char *);
inline double sigm(double x, double &deriv) {
double expl = 1./(1.+exp(-x));
deriv = expl*(1-expl);
return expl;
}
inline double tanh(double x, double &deriv) {
double expl = 2./(1.+exp(-2.*x))-1;
deriv = 1.-expl*expl;
return expl;
}
inline double relu(double x, double &deriv) {
if (x > 0) {
deriv = 1.;
return x;
} else {
deriv = 0.;
return 0;
}
}
};
}
#endif

4
src/MLIAP/mliap_model_python_couple.pyx
View File

@ -17,7 +17,7 @@ from libcpp.string cimport string
cdef extern from "mliap_data.h" namespace "LAMMPS_NS": cdef extern from "mliap_data.h" namespace "LAMMPS_NS":
cdef cppclass MLIAPData: cdef cppclass MLIAPData:
# Array shapes # Array shapes
int natoms int nlistatoms
int ndescriptors int ndescriptors
# Input data # Input data
@ -103,7 +103,7 @@ cdef public void MLIAPPY_compute_gradients(MLIAPModelPython * c_model, MLIAPData
model = retrieve(c_model) model = retrieve(c_model)
n_d = data.ndescriptors n_d = data.ndescriptors
n_a = data.natoms n_a = data.nlistatoms
# Make numpy arrays from pointers # Make numpy arrays from pointers
beta_np = np.asarray(<double[:n_a,:n_d] > &data.betas[0][0]) beta_np = np.asarray(<double[:n_a,:n_d] > &data.betas[0][0])

6
src/MLIAP/mliap_model_quadratic.cpp
View File

@ -59,7 +59,7 @@ void MLIAPModelQuadratic::compute_gradients(MLIAPData* data)
{ {
data->energy = 0.0; data->energy = 0.0;
for (int ii = 0; ii < data->natoms; ii++) { for (int ii = 0; ii < data->nlistatoms; ii++) {
const int ielem = data->ielems[ii]; const int ielem = data->ielems[ii];
double* coeffi = coeffelem[ielem]; double* coeffi = coeffelem[ielem];
@ -133,7 +133,7 @@ void MLIAPModelQuadratic::compute_gradgrads(class MLIAPData* data)
for (int l = 0; l < data->nelements*data->nparams; l++) for (int l = 0; l < data->nelements*data->nparams; l++)
data->egradient[l] = 0.0; data->egradient[l] = 0.0;
for (int ii = 0; ii < data->natoms; ii++) { for (int ii = 0; ii < data->nlistatoms; ii++) {
const int ielem = data->ielems[ii]; const int ielem = data->ielems[ii];
const int elemoffset = data->nparams*ielem; const int elemoffset = data->nparams*ielem;
@ -215,7 +215,7 @@ void MLIAPModelQuadratic::compute_force_gradients(class MLIAPData* data) {
data->egradient[l] = 0.0; data->egradient[l] = 0.0;
int ij = 0; int ij = 0;
for (int ii = 0; ii < data->natoms; ii++) { for (int ii = 0; ii < data->nlistatoms; ii++) {
const int i = data->iatoms[ii]; const int i = data->iatoms[ii];
const int ielem = data->ielems[ii]; const int ielem = data->ielems[ii];
const int elemoffset = data->nparams*ielem; const int elemoffset = data->nparams*ielem;

7
src/MLIAP/pair_mliap.cpp
View File

@ -20,6 +20,7 @@
#include "mliap_data.h" #include "mliap_data.h"
#include "mliap_model_linear.h" #include "mliap_model_linear.h"
#include "mliap_model_quadratic.h" #include "mliap_model_quadratic.h"
#include "mliap_model_nn.h"
#include "mliap_descriptor_snap.h" #include "mliap_descriptor_snap.h"
#ifdef MLIAP_PYTHON #ifdef MLIAP_PYTHON
#include "mliap_model_python.h" #include "mliap_model_python.h"
@ -152,6 +153,10 @@ void PairMLIAP::settings(int narg, char ** arg)
if (iarg+3 > narg) error->all(FLERR,"Illegal pair_style mliap command"); if (iarg+3 > narg) error->all(FLERR,"Illegal pair_style mliap command");
model = new MLIAPModelQuadratic(lmp,arg[iarg+2]); model = new MLIAPModelQuadratic(lmp,arg[iarg+2]);
iarg += 3; iarg += 3;
} else if (strcmp(arg[iarg+1],"nn") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal pair_style mliap command");
model = new MLIAPModelNN(lmp,arg[iarg+2]);
iarg += 3;
#ifdef MLIAP_PYTHON #ifdef MLIAP_PYTHON
} else if (strcmp(arg[iarg+1],"mliappy") == 0) { } else if (strcmp(arg[iarg+1],"mliappy") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal pair_style mliap command"); if (iarg+3 > narg) error->all(FLERR,"Illegal pair_style mliap command");
@ -250,7 +255,7 @@ void PairMLIAP::e_tally(MLIAPData* data)
{ {
if (eflag_global) eng_vdwl += data->energy; if (eflag_global) eng_vdwl += data->energy;
if (eflag_atom) if (eflag_atom)
for (int ii = 0; ii < data->natoms; ii++) { for (int ii = 0; ii < data->nlistatoms; ii++) {
const int i = data->iatoms[ii]; const int i = data->iatoms[ii];
eatom[i] += data->eatoms[ii]; eatom[i] += data->eatoms[ii];
} }