forked from lijiext/lammps
Merge pull request #2211 from athomps/compute-mliap
Add compute style mliap to MLIAP package
This commit is contained in:
Axel Kohlmeyer
GitHub
commit
2977a8aa15
|
|
@ -80,6 +80,7 @@ KOKKOS, o = USER-OMP, t = OPT.
|
|||
* :doc:`ke/eff <compute_ke_eff>`
|
||||
* :doc:`ke/rigid <compute_ke_rigid>`
|
||||
* :doc:`mesont <compute_mesont>`
|
||||
* :doc:`mliap <compute_mliap>`
|
||||
* :doc:`momentum <compute_momentum>`
|
||||
* :doc:`msd <compute_msd>`
|
||||
* :doc:`msd/chunk <compute_msd_chunk>`
|
||||
|
|
|
|||
|
|
@ -225,6 +225,7 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` doc
|
|||
* :doc:`ke/atom/eff <compute_ke_atom_eff>` - per-atom translational and radial kinetic energy in the electron force field model
|
||||
* :doc:`ke/eff <compute_ke_eff>` - kinetic energy of a group of nuclei and electrons in the electron force field model
|
||||
* :doc:`ke/rigid <compute_ke_rigid>` - translational kinetic energy of rigid bodies
|
||||
* :doc:`mliap <compute_mliap>` - gradients of energy and forces w.r.t. model parameters and related quantities for training machine learning interatomic potentials
|
||||
* :doc:`momentum <compute_momentum>` - translational momentum
|
||||
* :doc:`msd <compute_msd>` - mean-squared displacement of group of atoms
|
||||
* :doc:`msd/chunk <compute_msd_chunk>` - mean-squared displacement for each chunk
|
||||
|
|
@ -274,7 +275,7 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` doc
|
|||
* :doc:`smd/ulsph/strain/rate <compute_smd_ulsph_strain_rate>` -
|
||||
* :doc:`smd/ulsph/stress <compute_smd_ulsph_stress>` - per-particle Cauchy stress tensor and von Mises equivalent stress in Smooth Mach Dynamics
|
||||
* :doc:`smd/vol <compute_smd_vol>` - per-particle volumes and their sum in Smooth Mach Dynamics
|
||||
* :doc:`snap <compute_sna_atom>` - bispectrum components and related quantities for a group of atoms
|
||||
* :doc:`snap <compute_sna_atom>` - gradients of SNAP energy and forces w.r.t. linear coefficients and related quantities for fitting SNAP potentials
|
||||
* :doc:`sna/atom <compute_sna_atom>` - bispectrum components for each atom
|
||||
* :doc:`snad/atom <compute_sna_atom>` - derivative of bispectrum components for each atom
|
||||
* :doc:`snav/atom <compute_sna_atom>` - virial contribution from bispectrum components for each atom
|
||||
|
|
|
|||
|
|
@ -0,0 +1,170 @@
|
|||
.. index:: compute mliap
|
||||
|
||||
compute mliap command
|
||||
=====================
|
||||
|
||||
Syntax
|
||||
""""""
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
compute ID group-ID mliap ... keyword values ...
|
||||
|
||||
* ID, group-ID are documented in :doc:`compute <compute>` command
|
||||
* mliap = style name of this compute command
|
||||
* two keyword/value pairs must be appended
|
||||
* keyword = *model* or *descriptor*
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
*model* values = style Nelems Nparams
|
||||
style = *linear* or *quadratic*
|
||||
Nelems = number of elements
|
||||
Nparams = number of parameters per element
|
||||
*descriptor* values = style filename
|
||||
style = *sna*
|
||||
filename = name of file containing descriptor definitions
|
||||
|
||||
Examples
|
||||
""""""""
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
compute mliap model linear 2 31 descriptor sna Ta06A.mliap.descriptor
|
||||
|
||||
Description
|
||||
"""""""""""
|
||||
|
||||
Compute style *mliap* provides a general interface to the gradient
|
||||
of machine-learning interatomic potentials w.r.t. model parameters.
|
||||
It is used primarily for calculating the gradient of energy, force, and
|
||||
stress components w.r.t. model parameters, which is useful when training
|
||||
:doc:`mliap pair_style <pair_mliap>` models to match target data.
|
||||
It provides 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. Currently, the
|
||||
compute supports just two models, *linear* and *quadratic*,
|
||||
and one descriptor, *sna*, the SNAP descriptor used by
|
||||
:doc:`pair_style snap <pair_snap>`, including the linear, quadratic,
|
||||
and chem variants. Work is currently underway to extend
|
||||
the interface to handle neural network energy models,
|
||||
and it is also straightforward to add new descriptor styles.
|
||||
|
||||
The compute *mliap* command must be followed by two keywords
|
||||
*model* and *descriptor* in either order.
|
||||
|
||||
The *model* keyword is followed by a model style, currently limited to
|
||||
either *linear* or *quadratic*. In both cases,
|
||||
this is followed by two arguments. *nelems* is the number of elements.
|
||||
It must be equal to the number of LAMMPS atom types. *nparams*
|
||||
is the number of parameters per element for this model i.e.
|
||||
the number of parameter gradients for each element. Note these definitions
|
||||
are identical to those of *nelems* and *nparams* in the
|
||||
:doc:`pair_style mliap <pair_mliap>` model file.
|
||||
|
||||
The *descriptor* keyword is followed by a descriptor style, and additional arguments.
|
||||
Currently the only descriptor style is *sna*, indicating the bispectrum component
|
||||
descriptors used by the Spectral Neighbor Analysis Potential (SNAP) potentials of
|
||||
:doc:`pair_style snap <pair_snap>`.
|
||||
The \'p\' in SNAP is dropped, because keywords that match pair_styles are silently stripped
|
||||
out by the LAMMPS command parser. A single additional argument specifies the descriptor filename
|
||||
containing the parameters and setting used by the SNAP descriptor.
|
||||
The descriptor filename usually ends in the *.mliap.descriptor* extension.
|
||||
The format of this file is identical to the descriptor file in the
|
||||
:doc:`pair_style mliap <pair_mliap>`, and is described in detail
|
||||
there.
|
||||
|
||||
.. note::
|
||||
|
||||
The number of LAMMPS atom types (and the value of *nelems* in the model)
|
||||
must match the value of *nelems* in the descriptor file.
|
||||
|
||||
Compute *mliap* calculates a global array containing gradient information.
|
||||
The number of columns in the array is :math:`nelems \times nparams + 1`.
|
||||
The first row of the array contain the derivative of potential energy w.r.t. to
|
||||
each parameter and each element. The last six rows
|
||||
of the array contain the corresponding derivatives of the
|
||||
virial stress tensor, listed in Voigt notation: *pxx*, *pyy*, *pzz*,
|
||||
*pyz*, *pxz*, *pxy*. In between the energy and stress rows are
|
||||
the 3\*\ *N* rows containing the derivatives of the force components.
|
||||
See section below on output for a detailed description of how
|
||||
rows and columns are ordered.
|
||||
|
||||
The element in the last column of each row contains
|
||||
the potential energy, force, or stress, according to the row.
|
||||
These quantities correspond to the user-specified reference potential
|
||||
that must be subtracted from the target data when training a model.
|
||||
The potential energy calculation uses the built in compute *thermo_pe*.
|
||||
The stress calculation uses a compute called *mliap_press* that is
|
||||
automatically created behind the scenes, according to the following
|
||||
command:
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
compute mliap_press all pressure NULL virial
|
||||
|
||||
See section below on output for a detailed explanation of the data
|
||||
layout in the global array.
|
||||
|
||||
Atoms not in the group do not contribute to this compute.
|
||||
Neighbor atoms not in the group do not contribute to this compute.
|
||||
The neighbor list needed to compute this quantity is constructed each
|
||||
time the calculation is performed (i.e. each time a snapshot of atoms
|
||||
is dumped). Thus it can be inefficient to compute/dump this quantity
|
||||
too frequently.
|
||||
|
||||
.. note::
|
||||
|
||||
If the user-specified reference potentials includes bonded and
|
||||
non-bonded pairwise interactions, then the settings of
|
||||
:doc:`special_bonds <special_bonds>` command can remove pairwise
|
||||
interactions between atoms in the same bond, angle, or dihedral. This
|
||||
is the default setting for the :doc:`special_bonds <special_bonds>`
|
||||
command, and means those pairwise interactions do not appear in the
|
||||
neighbor list. Because this fix uses the neighbor list, it also means
|
||||
those pairs will not be included in the calculation. The :doc:`rerun <rerun>`
|
||||
command is not an option here, since the reference potential is required
|
||||
for the last column of the global array. A work-around is to prevent
|
||||
pairwise interactions from being removed by explicitly adding a
|
||||
*tiny* positive value for every pairwise interaction that would otherwise be
|
||||
set to zero in the :doc:`special_bonds <special_bonds>` command.
|
||||
|
||||
----------
|
||||
|
||||
**Output info:**
|
||||
|
||||
Compute *mliap* evaluates a global array.
|
||||
The columns are arranged into
|
||||
*nelems* blocks, listed in order of element *I*\ . Each block
|
||||
contains one column for each of the *nparams* model parameters.
|
||||
A final column contains the corresponding energy, force component
|
||||
on an atom, or virial stress component. The rows of the array appear
|
||||
in the following order:
|
||||
|
||||
* 1 row: Derivatives of potential energy w.r.t. each parameter of each element.
|
||||
* 3\*\ *N* rows: Derivatives of force components. x, y, and z components of force on atom *i* appearing in consecutive rows. The atoms are sorted based on atom ID.
|
||||
* 6 rows: Derivatives of virial stress tensor w.r.t. each parameter of each element. The ordering of the rows follows Voigt notation: *pxx*, *pyy*, *pzz*, *pyz*, *pxz*, *pxy*.
|
||||
|
||||
These values can be accessed by any command that uses a global array
|
||||
from a compute as input. See the :doc:`Howto output <Howto_output>` doc
|
||||
page for an overview of LAMMPS output options. 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>`_.
|
||||
|
||||
Restrictions
|
||||
""""""""""""
|
||||
|
||||
This compute is part of the MLIAP package. It is only enabled if
|
||||
LAMMPS was built with that package. In addition, building LAMMPS with the MLIAP package
|
||||
requires building LAMMPS with the SNAP package.
|
||||
See the :doc:`Build package <Build_package>` doc page for more info.
|
||||
|
||||
Related commands
|
||||
""""""""""""""""
|
||||
|
||||
:doc:`pair_style mliap <pair_mliap>`
|
||||
|
||||
**Default:** none
|
||||
|
|
@ -391,7 +391,9 @@ of :math:`K N_{elem}^3` columns.
|
|||
|
||||
These values can be accessed by any command that uses per-atom values
|
||||
from a compute as input. See the :doc:`Howto output <Howto_output>` doc
|
||||
page for an overview of LAMMPS output options.
|
||||
page for an overview of LAMMPS output options. To see how this command
|
||||
can be used within a Python workflow to train SNAP potentials,
|
||||
see the examples in `FitSNAP <https://github.com/FitSNAP/FitSNAP>`_.
|
||||
|
||||
Restrictions
|
||||
""""""""""""
|
||||
|
|
|
|||
|
|
@ -8,7 +8,19 @@ Syntax
|
|||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
pair_style mliap
|
||||
pair_style mliap ... keyword values ...
|
||||
|
||||
* two keyword/value pairs must be appended
|
||||
* keyword = *model* or *descriptor*
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
*model* values = style filename
|
||||
style = *linear* or *quadratic*
|
||||
filename = name of file containing model definitions
|
||||
*descriptor* values = style filename
|
||||
style = *sna*
|
||||
filename = name of file containing descriptor definitions
|
||||
|
||||
Examples
|
||||
""""""""
|
||||
|
|
@ -22,8 +34,8 @@ Examples
|
|||
Description
|
||||
"""""""""""
|
||||
|
||||
Pair style *mliap* provides a general interface to families of
|
||||
machine-learning interatomic potentials. It provides separate
|
||||
Pair style *mliap* provides a general interface to families of
|
||||
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,
|
||||
|
|
@ -34,6 +46,9 @@ and one descriptor, *sna*, the SNAP descriptor used by :doc:`pair_style snap <pa
|
|||
and chem variants. Work is currently underway to extend
|
||||
the interface to handle neural network energy models,
|
||||
and it is also straightforward to add new descriptor styles.
|
||||
In order to train a model, it is useful to know the gradient or derivative
|
||||
of energy, force, and stress w.r.t. model parameters. This information
|
||||
can be accessed using the related :doc:`compute mliap <compute_mliap>` command.
|
||||
|
||||
The pair_style *mliap* command must be followed by two keywords
|
||||
*model* and *descriptor* in either order. A single
|
||||
|
|
@ -46,10 +61,10 @@ where N is the number of LAMMPS atom types.
|
|||
|
||||
The *model* keyword is followed by a model style, currently limited to
|
||||
either *linear* or *quadratic*. In both cases,
|
||||
this is followed by a single argument specifying the model filename containing the
|
||||
linear or quadratic coefficients for a set of elements.
|
||||
this is followed by a single argument specifying the model filename containing the
|
||||
parameters for a set of elements.
|
||||
The model filename usually ends in the *.mliap.model* extension.
|
||||
It may contain coefficients for many elements. The only requirement is that it
|
||||
It may contain parameters for many elements. The only requirement is that it
|
||||
contain at least those element names appearing in the
|
||||
*pair_coeff* command.
|
||||
|
||||
|
|
@ -58,10 +73,10 @@ but follows a strict format after that. The first non-blank non-comment
|
|||
line must contain two integers:
|
||||
|
||||
* nelems = Number of elements
|
||||
* ncoeff = Number of coefficients
|
||||
* nparams = Number of parameters
|
||||
|
||||
This is followed by one block for each of the *nelem* elements.
|
||||
Each block consists of *ncoeff* coefficients, one per line.
|
||||
Each block consists of *nparams* parameters, one per line.
|
||||
Note that this format is similar, but not identical to that used
|
||||
for the :doc:`pair_style snap <pair_snap>` coefficient file.
|
||||
Specifically, the line containing the element weight and radius is omitted,
|
||||
|
|
@ -131,6 +146,6 @@ See the :doc:`Build package <Build_package>` doc page for more info.
|
|||
Related commands
|
||||
""""""""""""""""
|
||||
|
||||
:doc:`pair_style snap <pair_snap>`,
|
||||
:doc:`pair_style snap <pair_snap>`, :doc:`compute mliap <compute_mliap>`
|
||||
|
||||
**Default:** none
|
||||
|
|
|
|||
|
|
@ -2083,6 +2083,7 @@ Novint
|
|||
np
|
||||
Npair
|
||||
Npairs
|
||||
nparams
|
||||
nparticle
|
||||
npernode
|
||||
nph
|
||||
|
|
|
|||
|
|
@ -0,0 +1,19 @@
|
|||
# LAMMPS SNAP parameters for MLIAP compute example
|
||||
|
||||
# required
|
||||
rcutfac 1.0
|
||||
twojmax 2
|
||||
|
||||
# elements
|
||||
|
||||
nelems 2
|
||||
elems A B
|
||||
radelems 2.3 2.0
|
||||
welems 1.0 0.96
|
||||
|
||||
# optional
|
||||
|
||||
rfac0 0.99363
|
||||
rmin0 0
|
||||
bzeroflag 0
|
||||
switchflag 0
|
||||
|
|
@ -0,0 +1,95 @@
|
|||
# Demonstrate MLIAP quadratic compute
|
||||
|
||||
# initialize simulation
|
||||
|
||||
variable nsteps index 0
|
||||
variable nrep equal 1
|
||||
variable a equal 2.0
|
||||
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
|
||||
|
||||
atom_modify map hash
|
||||
lattice bcc $a
|
||||
region box block 0 ${nx} 0 ${ny} 0 ${nz}
|
||||
create_box 2 box
|
||||
create_atoms 2 box
|
||||
|
||||
mass * 180.88
|
||||
|
||||
displace_atoms all random 0.1 0.1 0.1 123456
|
||||
|
||||
# set up reference potential
|
||||
|
||||
variable zblcutinner equal 4
|
||||
variable zblcutouter equal 4.8
|
||||
variable zblz equal 73
|
||||
pair_style zbl ${zblcutinner} ${zblcutouter}
|
||||
pair_coeff * * ${zblz} ${zblz}
|
||||
|
||||
# choose SNA parameters
|
||||
|
||||
variable twojmax equal 2
|
||||
variable rcutfac equal 1.0
|
||||
variable rfac0 equal 0.99363
|
||||
variable rmin0 equal 0
|
||||
variable radelem1 equal 2.3
|
||||
variable radelem2 equal 2.0
|
||||
variable wj1 equal 1.0
|
||||
variable wj2 equal 0.96
|
||||
variable quadratic equal 1
|
||||
variable bzero equal 0
|
||||
variable switch equal 0
|
||||
variable snap_options string &
|
||||
"${rcutfac} ${rfac0} ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}"
|
||||
|
||||
# set up per-atom computes
|
||||
|
||||
compute b all sna/atom ${snap_options}
|
||||
compute vb all snav/atom ${snap_options}
|
||||
compute db all snad/atom ${snap_options}
|
||||
|
||||
# perform sums over atoms
|
||||
|
||||
group snapgroup1 type 1
|
||||
group snapgroup2 type 2
|
||||
compute bsum1 snapgroup1 reduce sum c_b[*]
|
||||
compute bsum2 snapgroup2 reduce sum c_b[*]
|
||||
# fix bsum1 all ave/time 1 1 1 c_bsum1 file bsum1.dat mode vector
|
||||
# fix bsum2 all ave/time 1 1 1 c_bsum2 file bsum2.dat mode vector
|
||||
compute vbsum all reduce sum c_vb[*]
|
||||
# fix vbsum all ave/time 1 1 1 c_vbsum file vbsum.dat mode vector
|
||||
variable db_2_100 equal c_db[2][100]
|
||||
|
||||
# set up compute snap generating global array
|
||||
|
||||
compute snap all mliap descriptor sna compute.mliap.descriptor model quadratic 2 21
|
||||
fix snap all ave/time 1 1 1 c_snap[*] file compute.quadratic.dat mode vector
|
||||
|
||||
thermo 100
|
||||
|
||||
# test output: 1: total potential energy
|
||||
# 2: xy component of stress tensor
|
||||
# 3: Sum(0.5*(B_{222}^i)^2, all i of type 2)
|
||||
# 4: xz component of Sum(Sum(r_j*(0.5*(dB_{222}^i)^2/dR[j]), all i of type 2), all j)
|
||||
# 5: y component of -Sum(d(0.5*(B_{222}^i)^2/dR[2]), all i of type 2)
|
||||
#
|
||||
# followed by 5 counterparts from compute snap
|
||||
|
||||
thermo_style custom &
|
||||
pe pxy c_bsum2[20] c_vbsum[220] v_db_2_100 &
|
||||
c_snap[1][43] c_snap[13][43] c_snap[1][42] c_snap[12][42] c_snap[6][42]
|
||||
thermo_modify norm no
|
||||
|
||||
# dump mydump_db all custom 1000 dump_db id c_db[*]
|
||||
# dump_modify mydump_db sort id
|
||||
|
||||
# Run MD
|
||||
|
||||
run ${nsteps}
|
||||
|
|
@ -0,0 +1,95 @@
|
|||
# Demonstrate bispectrum computes
|
||||
|
||||
# initialize simulation
|
||||
|
||||
variable nsteps index 0
|
||||
variable nrep equal 1
|
||||
variable a equal 2.0
|
||||
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
|
||||
|
||||
atom_modify map hash
|
||||
lattice bcc $a
|
||||
region box block 0 ${nx} 0 ${ny} 0 ${nz}
|
||||
create_box 2 box
|
||||
create_atoms 2 box
|
||||
|
||||
mass * 180.88
|
||||
|
||||
displace_atoms all random 0.1 0.1 0.1 123456
|
||||
|
||||
# set up reference potential
|
||||
|
||||
variable zblcutinner equal 4
|
||||
variable zblcutouter equal 4.8
|
||||
variable zblz equal 73
|
||||
pair_style zbl ${zblcutinner} ${zblcutouter}
|
||||
pair_coeff * * ${zblz} ${zblz}
|
||||
|
||||
# choose SNA parameters
|
||||
|
||||
variable twojmax equal 2
|
||||
variable rcutfac equal 1.0
|
||||
variable rfac0 equal 0.99363
|
||||
variable rmin0 equal 0
|
||||
variable radelem1 equal 2.3
|
||||
variable radelem2 equal 2.0
|
||||
variable wj1 equal 1.0
|
||||
variable wj2 equal 0.96
|
||||
variable quadratic equal 0
|
||||
variable bzero equal 0
|
||||
variable switch equal 0
|
||||
variable snap_options string &
|
||||
"${rcutfac} ${rfac0} ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}"
|
||||
|
||||
# set up per-atom computes
|
||||
|
||||
compute b all sna/atom ${snap_options}
|
||||
compute vb all snav/atom ${snap_options}
|
||||
compute db all snad/atom ${snap_options}
|
||||
|
||||
# perform sums over atoms
|
||||
|
||||
group snapgroup1 type 1
|
||||
group snapgroup2 type 2
|
||||
compute bsum1 snapgroup1 reduce sum c_b[*]
|
||||
compute bsum2 snapgroup2 reduce sum c_b[*]
|
||||
# fix bsum1 all ave/time 1 1 1 c_bsum1 file bsum1.dat mode vector
|
||||
# fix bsum2 all ave/time 1 1 1 c_bsum2 file bsum2.dat mode vector
|
||||
compute vbsum all reduce sum c_vb[*]
|
||||
# fix vbsum all ave/time 1 1 1 c_vbsum file vbsum.dat mode vector
|
||||
variable db_2_25 equal c_db[2][25]
|
||||
|
||||
# set up compute snap generating global array
|
||||
|
||||
compute snap all mliap descriptor sna compute.mliap.descriptor model linear 2 6
|
||||
fix snap all ave/time 1 1 1 c_snap[*] file compute.snap.dat mode vector
|
||||
|
||||
thermo 100
|
||||
|
||||
# test output: 1: total potential energy
|
||||
# 2: xy component of stress tensor
|
||||
# 3: Sum(B_{000}^i, all i of type 2)
|
||||
# 4: xz component of Sum(Sum(r_j*dB_{222}^i/dR[j]), all i of type 2), all j)
|
||||
# 5: y component of -Sum(d(B_{222}^i)/dR[2]), all i of type 2)
|
||||
#
|
||||
# followed by 5 counterparts from compute snap
|
||||
|
||||
thermo_style custom &
|
||||
pe pxy c_bsum2[1] c_vbsum[55] v_db_2_25 &
|
||||
c_snap[1][13] c_snap[13][13] c_snap[1][8] c_snap[12][12] c_snap[6][12]
|
||||
thermo_modify norm no
|
||||
|
||||
# dump mydump_db all custom 1000 dump_db id c_db[*]
|
||||
# dump_modify mydump_db sort id
|
||||
|
||||
# Run MD
|
||||
|
||||
run ${nsteps}
|
||||
|
|
@ -0,0 +1,199 @@
|
|||
LAMMPS (15 Jun 2020)
|
||||
# Demonstrate MLIAP quadratic compute
|
||||
|
||||
# initialize simulation
|
||||
|
||||
variable nsteps index 0
|
||||
variable nrep equal 1
|
||||
variable a equal 2.0
|
||||
units metal
|
||||
|
||||
# generate the box and atom positions using a BCC lattice
|
||||
|
||||
variable nx equal ${nrep}
|
||||
variable nx equal 1
|
||||
variable ny equal ${nrep}
|
||||
variable ny equal 1
|
||||
variable nz equal ${nrep}
|
||||
variable nz equal 1
|
||||
|
||||
boundary p p p
|
||||
|
||||
atom_modify map hash
|
||||
lattice bcc $a
|
||||
lattice bcc 2
|
||||
Lattice spacing in x,y,z = 2 2 2
|
||||
region box block 0 ${nx} 0 ${ny} 0 ${nz}
|
||||
region box block 0 1 0 ${ny} 0 ${nz}
|
||||
region box block 0 1 0 1 0 ${nz}
|
||||
region box block 0 1 0 1 0 1
|
||||
create_box 2 box
|
||||
Created orthogonal box = (0.0 0.0 0.0) to (2.0 2.0 2.0)
|
||||
1 by 1 by 1 MPI processor grid
|
||||
create_atoms 2 box
|
||||
Created 2 atoms
|
||||
create_atoms CPU = 0.000 seconds
|
||||
|
||||
mass * 180.88
|
||||
|
||||
displace_atoms all random 0.1 0.1 0.1 123456
|
||||
|
||||
# set up reference potential
|
||||
|
||||
variable zblcutinner equal 4
|
||||
variable zblcutouter equal 4.8
|
||||
variable zblz equal 73
|
||||
pair_style zbl ${zblcutinner} ${zblcutouter}
|
||||
pair_style zbl 4 ${zblcutouter}
|
||||
pair_style zbl 4 4.8
|
||||
pair_coeff * * ${zblz} ${zblz}
|
||||
pair_coeff * * 73 ${zblz}
|
||||
pair_coeff * * 73 73
|
||||
|
||||
# choose SNA parameters
|
||||
|
||||
variable twojmax equal 2
|
||||
variable rcutfac equal 1.0
|
||||
variable rfac0 equal 0.99363
|
||||
variable rmin0 equal 0
|
||||
variable radelem1 equal 2.3
|
||||
variable radelem2 equal 2.0
|
||||
variable wj1 equal 1.0
|
||||
variable wj2 equal 0.96
|
||||
variable quadratic equal 1
|
||||
variable bzero equal 0
|
||||
variable switch equal 0
|
||||
variable snap_options string "${rcutfac} ${rfac0} ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}"
|
||||
1 ${rfac0} ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 2 ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 2 2.3 ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 2 2.3 2 ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 2 2.3 2 1 ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 2 2.3 2 1 0.96 rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 1 bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 1 bzeroflag 0 switchflag ${switch}
|
||||
1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 1 bzeroflag 0 switchflag 0
|
||||
|
||||
# set up per-atom computes
|
||||
|
||||
compute b all sna/atom ${snap_options}
|
||||
compute b all sna/atom 1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 1 bzeroflag 0 switchflag 0
|
||||
compute vb all snav/atom ${snap_options}
|
||||
compute vb all snav/atom 1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 1 bzeroflag 0 switchflag 0
|
||||
compute db all snad/atom ${snap_options}
|
||||
compute db all snad/atom 1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 1 bzeroflag 0 switchflag 0
|
||||
|
||||
# perform sums over atoms
|
||||
|
||||
group snapgroup1 type 1
|
||||
0 atoms in group snapgroup1
|
||||
group snapgroup2 type 2
|
||||
2 atoms in group snapgroup2
|
||||
compute bsum1 snapgroup1 reduce sum c_b[*]
|
||||
compute bsum2 snapgroup2 reduce sum c_b[*]
|
||||
# fix bsum1 all ave/time 1 1 1 c_bsum1 file bsum1.dat mode vector
|
||||
# fix bsum2 all ave/time 1 1 1 c_bsum2 file bsum2.dat mode vector
|
||||
compute vbsum all reduce sum c_vb[*]
|
||||
# fix vbsum all ave/time 1 1 1 c_vbsum file vbsum.dat mode vector
|
||||
variable db_2_100 equal c_db[2][100]
|
||||
|
||||
# set up compute snap generating global array
|
||||
|
||||
compute snap all mliap descriptor sna compute.mliap.descriptor model quadratic 2 21
|
||||
SNAP keyword rcutfac 1.0
|
||||
SNAP keyword twojmax 2
|
||||
SNAP keyword nelems 2
|
||||
SNAP keyword elems A
|
||||
SNAP keyword radelems 2.3
|
||||
SNAP keyword welems 1.0
|
||||
SNAP keyword rfac0 0.99363
|
||||
SNAP keyword rmin0 0
|
||||
SNAP keyword bzeroflag 0
|
||||
SNAP keyword switchflag 0
|
||||
fix snap all ave/time 1 1 1 c_snap[*] file compute.quadratic.dat mode vector
|
||||
|
||||
thermo 100
|
||||
|
||||
# test output: 1: total potential energy
|
||||
# 2: xy component of stress tensor
|
||||
# 3: Sum(0.5*(B_{222}^i)^2, all i of type 2)
|
||||
# 4: xz component of Sum(Sum(r_j*(0.5*(dB_{222}^i)^2/dR[j]), all i of type 2), all j)
|
||||
# 5: y component of -Sum(d(0.5*(B_{222}^i)^2/dR[2]), all i of type 2)
|
||||
#
|
||||
# followed by 5 counterparts from compute snap
|
||||
|
||||
thermo_style custom pe pxy c_bsum2[20] c_vbsum[220] v_db_2_100 c_snap[1][43] c_snap[13][43] c_snap[1][42] c_snap[12][42] c_snap[6][42]
|
||||
thermo_modify norm no
|
||||
|
||||
# dump mydump_db all custom 1000 dump_db id c_db[*]
|
||||
# dump_modify mydump_db sort id
|
||||
|
||||
# Run MD
|
||||
|
||||
run ${nsteps}
|
||||
run 0
|
||||
Neighbor list info ...
|
||||
update every 1 steps, delay 10 steps, check yes
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
master list distance cutoff = 6.8
|
||||
ghost atom cutoff = 6.8
|
||||
binsize = 3.4, bins = 1 1 1
|
||||
5 neighbor lists, perpetual/occasional/extra = 1 4 0
|
||||
(1) pair zbl, perpetual
|
||||
attributes: half, newton on
|
||||
pair build: half/bin/atomonly/newton
|
||||
stencil: half/bin/3d/newton
|
||||
bin: standard
|
||||
(2) compute sna/atom, occasional
|
||||
attributes: full, newton on
|
||||
pair build: full/bin/atomonly
|
||||
stencil: full/bin/3d
|
||||
bin: standard
|
||||
(3) compute snav/atom, occasional
|
||||
attributes: full, newton on
|
||||
pair build: full/bin/atomonly
|
||||
stencil: full/bin/3d
|
||||
bin: standard
|
||||
(4) compute snad/atom, occasional
|
||||
attributes: full, newton on
|
||||
pair build: full/bin/atomonly
|
||||
stencil: full/bin/3d
|
||||
bin: standard
|
||||
(5) compute mliap, occasional
|
||||
attributes: full, newton on
|
||||
pair build: full/bin/atomonly
|
||||
stencil: full/bin/3d
|
||||
bin: standard
|
||||
Per MPI rank memory allocation (min/avg/max) = 22.45 | 22.45 | 22.45 Mbytes
|
||||
PotEng Pxy c_bsum2[20] c_vbsum[220] v_db_2_100 c_snap[1][43] c_snap[13][43] c_snap[1][42] c_snap[12][42] c_snap[6][42]
|
||||
322.86952 1505558.1 4.2492771e+08 -4952473 28484035 322.86952 1505558.1 4.2492771e+08 -4952473 28484035
|
||||
Loop time of 1e-06 on 1 procs for 0 steps with 2 atoms
|
||||
|
||||
100.0% 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 | 0 | 0 | 0.0 | 0.00
|
||||
Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Comm | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Output | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Modify | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Other | | 1e-06 | | |100.00
|
||||
|
||||
Nlocal: 2 ave 2 max 2 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 853 ave 853 max 853 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 330 ave 330 max 330 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
FullNghs: 660 ave 660 max 660 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 660
|
||||
Ave neighs/atom = 330
|
||||
Neighbor list builds = 0
|
||||
Dangerous builds = 0
|
||||
Total wall time: 0:00:00
|
||||
|
|
@ -0,0 +1,200 @@
|
|||
LAMMPS (15 Jun 2020)
|
||||
# Demonstrate MLIAP quadratic compute
|
||||
|
||||
# initialize simulation
|
||||
|
||||
variable nsteps index 0
|
||||
variable nrep equal 1
|
||||
variable a equal 2.0
|
||||
units metal
|
||||
|
||||
# generate the box and atom positions using a BCC lattice
|
||||
|
||||
variable nx equal ${nrep}
|
||||
variable nx equal 1
|
||||
variable ny equal ${nrep}
|
||||
variable ny equal 1
|
||||
variable nz equal ${nrep}
|
||||
variable nz equal 1
|
||||
|
||||
boundary p p p
|
||||
|
||||
atom_modify map hash
|
||||
lattice bcc $a
|
||||
lattice bcc 2
|
||||
Lattice spacing in x,y,z = 2 2 2
|
||||
region box block 0 ${nx} 0 ${ny} 0 ${nz}
|
||||
region box block 0 1 0 ${ny} 0 ${nz}
|
||||
region box block 0 1 0 1 0 ${nz}
|
||||
region box block 0 1 0 1 0 1
|
||||
create_box 2 box
|
||||
Created orthogonal box = (0.0 0.0 0.0) to (2.0 2.0 2.0)
|
||||
1 by 2 by 2 MPI processor grid
|
||||
create_atoms 2 box
|
||||
Created 2 atoms
|
||||
create_atoms CPU = 0.001 seconds
|
||||
|
||||
mass * 180.88
|
||||
|
||||
displace_atoms all random 0.1 0.1 0.1 123456
|
||||
|
||||
# set up reference potential
|
||||
|
||||
variable zblcutinner equal 4
|
||||
variable zblcutouter equal 4.8
|
||||
variable zblz equal 73
|
||||
pair_style zbl ${zblcutinner} ${zblcutouter}
|
||||
pair_style zbl 4 ${zblcutouter}
|
||||
pair_style zbl 4 4.8
|
||||
pair_coeff * * ${zblz} ${zblz}
|
||||
pair_coeff * * 73 ${zblz}
|
||||
pair_coeff * * 73 73
|
||||
|
||||
# choose SNA parameters
|
||||
|
||||
variable twojmax equal 2
|
||||
variable rcutfac equal 1.0
|
||||
variable rfac0 equal 0.99363
|
||||
variable rmin0 equal 0
|
||||
variable radelem1 equal 2.3
|
||||
variable radelem2 equal 2.0
|
||||
variable wj1 equal 1.0
|
||||
variable wj2 equal 0.96
|
||||
variable quadratic equal 1
|
||||
variable bzero equal 0
|
||||
variable switch equal 0
|
||||
variable snap_options string "${rcutfac} ${rfac0} ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}"
|
||||
1 ${rfac0} ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 2 ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 2 2.3 ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 2 2.3 2 ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 2 2.3 2 1 ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 2 2.3 2 1 0.96 rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 1 bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 1 bzeroflag 0 switchflag ${switch}
|
||||
1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 1 bzeroflag 0 switchflag 0
|
||||
|
||||
# set up per-atom computes
|
||||
|
||||
compute b all sna/atom ${snap_options}
|
||||
compute b all sna/atom 1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 1 bzeroflag 0 switchflag 0
|
||||
compute vb all snav/atom ${snap_options}
|
||||
compute vb all snav/atom 1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 1 bzeroflag 0 switchflag 0
|
||||
compute db all snad/atom ${snap_options}
|
||||
compute db all snad/atom 1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 1 bzeroflag 0 switchflag 0
|
||||
|
||||
# perform sums over atoms
|
||||
|
||||
group snapgroup1 type 1
|
||||
0 atoms in group snapgroup1
|
||||
group snapgroup2 type 2
|
||||
2 atoms in group snapgroup2
|
||||
compute bsum1 snapgroup1 reduce sum c_b[*]
|
||||
compute bsum2 snapgroup2 reduce sum c_b[*]
|
||||
# fix bsum1 all ave/time 1 1 1 c_bsum1 file bsum1.dat mode vector
|
||||
# fix bsum2 all ave/time 1 1 1 c_bsum2 file bsum2.dat mode vector
|
||||
compute vbsum all reduce sum c_vb[*]
|
||||
# fix vbsum all ave/time 1 1 1 c_vbsum file vbsum.dat mode vector
|
||||
variable db_2_100 equal c_db[2][100]
|
||||
|
||||
# set up compute snap generating global array
|
||||
|
||||
compute snap all mliap descriptor sna compute.mliap.descriptor model quadratic 2 21
|
||||
SNAP keyword rcutfac 1.0
|
||||
SNAP keyword twojmax 2
|
||||
SNAP keyword nelems 2
|
||||
SNAP keyword elems A
|
||||
SNAP keyword radelems 2.3
|
||||
SNAP keyword welems 1.0
|
||||
SNAP keyword rfac0 0.99363
|
||||
SNAP keyword rmin0 0
|
||||
SNAP keyword bzeroflag 0
|
||||
SNAP keyword switchflag 0
|
||||
fix snap all ave/time 1 1 1 c_snap[*] file compute.quadratic.dat mode vector
|
||||
|
||||
thermo 100
|
||||
|
||||
# test output: 1: total potential energy
|
||||
# 2: xy component of stress tensor
|
||||
# 3: Sum(0.5*(B_{222}^i)^2, all i of type 2)
|
||||
# 4: xz component of Sum(Sum(r_j*(0.5*(dB_{222}^i)^2/dR[j]), all i of type 2), all j)
|
||||
# 5: y component of -Sum(d(0.5*(B_{222}^i)^2/dR[2]), all i of type 2)
|
||||
#
|
||||
# followed by 5 counterparts from compute snap
|
||||
|
||||
thermo_style custom pe pxy c_bsum2[20] c_vbsum[220] v_db_2_100 c_snap[1][43] c_snap[13][43] c_snap[1][42] c_snap[12][42] c_snap[6][42]
|
||||
thermo_modify norm no
|
||||
|
||||
# dump mydump_db all custom 1000 dump_db id c_db[*]
|
||||
# dump_modify mydump_db sort id
|
||||
|
||||
# Run MD
|
||||
|
||||
run ${nsteps}
|
||||
run 0
|
||||
Neighbor list info ...
|
||||
update every 1 steps, delay 10 steps, check yes
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
master list distance cutoff = 6.8
|
||||
ghost atom cutoff = 6.8
|
||||
binsize = 3.4, bins = 1 1 1
|
||||
5 neighbor lists, perpetual/occasional/extra = 1 4 0
|
||||
(1) pair zbl, perpetual
|
||||
attributes: half, newton on
|
||||
pair build: half/bin/atomonly/newton
|
||||
stencil: half/bin/3d/newton
|
||||
bin: standard
|
||||
(2) compute sna/atom, occasional
|
||||
attributes: full, newton on
|
||||
pair build: full/bin/atomonly
|
||||
stencil: full/bin/3d
|
||||
bin: standard
|
||||
(3) compute snav/atom, occasional
|
||||
attributes: full, newton on
|
||||
pair build: full/bin/atomonly
|
||||
stencil: full/bin/3d
|
||||
bin: standard
|
||||
(4) compute snad/atom, occasional
|
||||
attributes: full, newton on
|
||||
pair build: full/bin/atomonly
|
||||
stencil: full/bin/3d
|
||||
bin: standard
|
||||
(5) compute mliap, occasional
|
||||
attributes: full, newton on
|
||||
pair build: full/bin/atomonly
|
||||
stencil: full/bin/3d
|
||||
bin: standard
|
||||
WARNING: Proc sub-domain size < neighbor skin, could lead to lost atoms (../domain.cpp:964)
|
||||
Per MPI rank memory allocation (min/avg/max) = 22.18 | 22.18 | 22.18 Mbytes
|
||||
PotEng Pxy c_bsum2[20] c_vbsum[220] v_db_2_100 c_snap[1][43] c_snap[13][43] c_snap[1][42] c_snap[12][42] c_snap[6][42]
|
||||
322.86952 1505558.1 4.2492771e+08 -4952473 28484035 322.86952 1505558.1 4.2492771e+08 -4952473 28484035
|
||||
Loop time of 2e-06 on 4 procs for 0 steps with 2 atoms
|
||||
|
||||
100.0% CPU use with 4 MPI tasks x no OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Comm | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Output | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Modify | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Other | | 2e-06 | | |100.00
|
||||
|
||||
Nlocal: 0.5 ave 1 max 0 min
|
||||
Histogram: 2 0 0 0 0 0 0 0 0 2
|
||||
Nghost: 734.5 ave 735 max 734 min
|
||||
Histogram: 2 0 0 0 0 0 0 0 0 2
|
||||
Neighs: 82.5 ave 177 max 0 min
|
||||
Histogram: 2 0 0 0 0 0 0 0 1 1
|
||||
FullNghs: 165 ave 330 max 0 min
|
||||
Histogram: 2 0 0 0 0 0 0 0 0 2
|
||||
|
||||
Total # of neighbors = 660
|
||||
Ave neighs/atom = 330
|
||||
Neighbor list builds = 0
|
||||
Dangerous builds = 0
|
||||
Total wall time: 0:00:00
|
||||
|
|
@ -0,0 +1,199 @@
|
|||
LAMMPS (15 Jun 2020)
|
||||
# Demonstrate bispectrum computes
|
||||
|
||||
# initialize simulation
|
||||
|
||||
variable nsteps index 0
|
||||
variable nrep equal 1
|
||||
variable a equal 2.0
|
||||
units metal
|
||||
|
||||
# generate the box and atom positions using a BCC lattice
|
||||
|
||||
variable nx equal ${nrep}
|
||||
variable nx equal 1
|
||||
variable ny equal ${nrep}
|
||||
variable ny equal 1
|
||||
variable nz equal ${nrep}
|
||||
variable nz equal 1
|
||||
|
||||
boundary p p p
|
||||
|
||||
atom_modify map hash
|
||||
lattice bcc $a
|
||||
lattice bcc 2
|
||||
Lattice spacing in x,y,z = 2 2 2
|
||||
region box block 0 ${nx} 0 ${ny} 0 ${nz}
|
||||
region box block 0 1 0 ${ny} 0 ${nz}
|
||||
region box block 0 1 0 1 0 ${nz}
|
||||
region box block 0 1 0 1 0 1
|
||||
create_box 2 box
|
||||
Created orthogonal box = (0.0 0.0 0.0) to (2.0 2.0 2.0)
|
||||
1 by 1 by 1 MPI processor grid
|
||||
create_atoms 2 box
|
||||
Created 2 atoms
|
||||
create_atoms CPU = 0.001 seconds
|
||||
|
||||
mass * 180.88
|
||||
|
||||
displace_atoms all random 0.1 0.1 0.1 123456
|
||||
|
||||
# set up reference potential
|
||||
|
||||
variable zblcutinner equal 4
|
||||
variable zblcutouter equal 4.8
|
||||
variable zblz equal 73
|
||||
pair_style zbl ${zblcutinner} ${zblcutouter}
|
||||
pair_style zbl 4 ${zblcutouter}
|
||||
pair_style zbl 4 4.8
|
||||
pair_coeff * * ${zblz} ${zblz}
|
||||
pair_coeff * * 73 ${zblz}
|
||||
pair_coeff * * 73 73
|
||||
|
||||
# choose SNA parameters
|
||||
|
||||
variable twojmax equal 2
|
||||
variable rcutfac equal 1.0
|
||||
variable rfac0 equal 0.99363
|
||||
variable rmin0 equal 0
|
||||
variable radelem1 equal 2.3
|
||||
variable radelem2 equal 2.0
|
||||
variable wj1 equal 1.0
|
||||
variable wj2 equal 0.96
|
||||
variable quadratic equal 0
|
||||
variable bzero equal 0
|
||||
variable switch equal 0
|
||||
variable snap_options string "${rcutfac} ${rfac0} ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}"
|
||||
1 ${rfac0} ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 2 ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 2 2.3 ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 2 2.3 2 ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 2 2.3 2 1 ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 2 2.3 2 1 0.96 rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 0 bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag ${switch}
|
||||
1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 0
|
||||
|
||||
# set up per-atom computes
|
||||
|
||||
compute b all sna/atom ${snap_options}
|
||||
compute b all sna/atom 1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 0
|
||||
compute vb all snav/atom ${snap_options}
|
||||
compute vb all snav/atom 1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 0
|
||||
compute db all snad/atom ${snap_options}
|
||||
compute db all snad/atom 1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 0
|
||||
|
||||
# perform sums over atoms
|
||||
|
||||
group snapgroup1 type 1
|
||||
0 atoms in group snapgroup1
|
||||
group snapgroup2 type 2
|
||||
2 atoms in group snapgroup2
|
||||
compute bsum1 snapgroup1 reduce sum c_b[*]
|
||||
compute bsum2 snapgroup2 reduce sum c_b[*]
|
||||
# fix bsum1 all ave/time 1 1 1 c_bsum1 file bsum1.dat mode vector
|
||||
# fix bsum2 all ave/time 1 1 1 c_bsum2 file bsum2.dat mode vector
|
||||
compute vbsum all reduce sum c_vb[*]
|
||||
# fix vbsum all ave/time 1 1 1 c_vbsum file vbsum.dat mode vector
|
||||
variable db_2_25 equal c_db[2][25]
|
||||
|
||||
# set up compute snap generating global array
|
||||
|
||||
compute snap all mliap descriptor sna compute.mliap.descriptor model linear 2 6
|
||||
SNAP keyword rcutfac 1.0
|
||||
SNAP keyword twojmax 2
|
||||
SNAP keyword nelems 2
|
||||
SNAP keyword elems A
|
||||
SNAP keyword radelems 2.3
|
||||
SNAP keyword welems 1.0
|
||||
SNAP keyword rfac0 0.99363
|
||||
SNAP keyword rmin0 0
|
||||
SNAP keyword bzeroflag 0
|
||||
SNAP keyword switchflag 0
|
||||
fix snap all ave/time 1 1 1 c_snap[*] file compute.snap.dat mode vector
|
||||
|
||||
thermo 100
|
||||
|
||||
# test output: 1: total potential energy
|
||||
# 2: xy component of stress tensor
|
||||
# 3: Sum(B_{000}^i, all i of type 2)
|
||||
# 4: xz component of Sum(Sum(r_j*dB_{222}^i/dR[j]), all i of type 2), all j)
|
||||
# 5: y component of -Sum(d(B_{222}^i)/dR[2]), all i of type 2)
|
||||
#
|
||||
# followed by 5 counterparts from compute snap
|
||||
|
||||
thermo_style custom pe pxy c_bsum2[1] c_vbsum[55] v_db_2_25 c_snap[1][13] c_snap[13][13] c_snap[1][8] c_snap[12][12] c_snap[6][12]
|
||||
thermo_modify norm no
|
||||
|
||||
# dump mydump_db all custom 1000 dump_db id c_db[*]
|
||||
# dump_modify mydump_db sort id
|
||||
|
||||
# Run MD
|
||||
|
||||
run ${nsteps}
|
||||
run 0
|
||||
Neighbor list info ...
|
||||
update every 1 steps, delay 10 steps, check yes
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
master list distance cutoff = 6.8
|
||||
ghost atom cutoff = 6.8
|
||||
binsize = 3.4, bins = 1 1 1
|
||||
5 neighbor lists, perpetual/occasional/extra = 1 4 0
|
||||
(1) pair zbl, perpetual
|
||||
attributes: half, newton on
|
||||
pair build: half/bin/atomonly/newton
|
||||
stencil: half/bin/3d/newton
|
||||
bin: standard
|
||||
(2) compute sna/atom, occasional
|
||||
attributes: full, newton on
|
||||
pair build: full/bin/atomonly
|
||||
stencil: full/bin/3d
|
||||
bin: standard
|
||||
(3) compute snav/atom, occasional
|
||||
attributes: full, newton on
|
||||
pair build: full/bin/atomonly
|
||||
stencil: full/bin/3d
|
||||
bin: standard
|
||||
(4) compute snad/atom, occasional
|
||||
attributes: full, newton on
|
||||
pair build: full/bin/atomonly
|
||||
stencil: full/bin/3d
|
||||
bin: standard
|
||||
(5) compute mliap, occasional
|
||||
attributes: full, newton on
|
||||
pair build: full/bin/atomonly
|
||||
stencil: full/bin/3d
|
||||
bin: standard
|
||||
Per MPI rank memory allocation (min/avg/max) = 10.73 | 10.73 | 10.73 Mbytes
|
||||
PotEng Pxy c_bsum2[1] c_vbsum[55] v_db_2_25 c_snap[1][13] c_snap[13][13] c_snap[1][8] c_snap[12][12] c_snap[6][12]
|
||||
322.86952 1505558.1 364182.88 -240.25066 1381.7961 322.86952 1505558.1 364182.88 -240.25066 1381.7961
|
||||
Loop time of 0 on 1 procs for 0 steps with 2 atoms
|
||||
|
||||
0.0% 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 | 0 | 0 | 0.0 | 0.00
|
||||
Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Comm | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Output | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Modify | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Other | | 0 | | | 0.00
|
||||
|
||||
Nlocal: 2 ave 2 max 2 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 853 ave 853 max 853 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 330 ave 330 max 330 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
FullNghs: 660 ave 660 max 660 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 660
|
||||
Ave neighs/atom = 330
|
||||
Neighbor list builds = 0
|
||||
Dangerous builds = 0
|
||||
Total wall time: 0:00:00
|
||||
|
|
@ -0,0 +1,200 @@
|
|||
LAMMPS (15 Jun 2020)
|
||||
# Demonstrate bispectrum computes
|
||||
|
||||
# initialize simulation
|
||||
|
||||
variable nsteps index 0
|
||||
variable nrep equal 1
|
||||
variable a equal 2.0
|
||||
units metal
|
||||
|
||||
# generate the box and atom positions using a BCC lattice
|
||||
|
||||
variable nx equal ${nrep}
|
||||
variable nx equal 1
|
||||
variable ny equal ${nrep}
|
||||
variable ny equal 1
|
||||
variable nz equal ${nrep}
|
||||
variable nz equal 1
|
||||
|
||||
boundary p p p
|
||||
|
||||
atom_modify map hash
|
||||
lattice bcc $a
|
||||
lattice bcc 2
|
||||
Lattice spacing in x,y,z = 2 2 2
|
||||
region box block 0 ${nx} 0 ${ny} 0 ${nz}
|
||||
region box block 0 1 0 ${ny} 0 ${nz}
|
||||
region box block 0 1 0 1 0 ${nz}
|
||||
region box block 0 1 0 1 0 1
|
||||
create_box 2 box
|
||||
Created orthogonal box = (0.0 0.0 0.0) to (2.0 2.0 2.0)
|
||||
1 by 2 by 2 MPI processor grid
|
||||
create_atoms 2 box
|
||||
Created 2 atoms
|
||||
create_atoms CPU = 0.001 seconds
|
||||
|
||||
mass * 180.88
|
||||
|
||||
displace_atoms all random 0.1 0.1 0.1 123456
|
||||
|
||||
# set up reference potential
|
||||
|
||||
variable zblcutinner equal 4
|
||||
variable zblcutouter equal 4.8
|
||||
variable zblz equal 73
|
||||
pair_style zbl ${zblcutinner} ${zblcutouter}
|
||||
pair_style zbl 4 ${zblcutouter}
|
||||
pair_style zbl 4 4.8
|
||||
pair_coeff * * ${zblz} ${zblz}
|
||||
pair_coeff * * 73 ${zblz}
|
||||
pair_coeff * * 73 73
|
||||
|
||||
# choose SNA parameters
|
||||
|
||||
variable twojmax equal 2
|
||||
variable rcutfac equal 1.0
|
||||
variable rfac0 equal 0.99363
|
||||
variable rmin0 equal 0
|
||||
variable radelem1 equal 2.3
|
||||
variable radelem2 equal 2.0
|
||||
variable wj1 equal 1.0
|
||||
variable wj2 equal 0.96
|
||||
variable quadratic equal 0
|
||||
variable bzero equal 0
|
||||
variable switch equal 0
|
||||
variable snap_options string "${rcutfac} ${rfac0} ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}"
|
||||
1 ${rfac0} ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 2 ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 2 2.3 ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 2 2.3 2 ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 2 2.3 2 1 ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 2 2.3 2 1 0.96 rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 0 bzeroflag ${bzero} switchflag ${switch}
|
||||
1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag ${switch}
|
||||
1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 0
|
||||
|
||||
# set up per-atom computes
|
||||
|
||||
compute b all sna/atom ${snap_options}
|
||||
compute b all sna/atom 1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 0
|
||||
compute vb all snav/atom ${snap_options}
|
||||
compute vb all snav/atom 1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 0
|
||||
compute db all snad/atom ${snap_options}
|
||||
compute db all snad/atom 1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 0
|
||||
|
||||
# perform sums over atoms
|
||||
|
||||
group snapgroup1 type 1
|
||||
0 atoms in group snapgroup1
|
||||
group snapgroup2 type 2
|
||||
2 atoms in group snapgroup2
|
||||
compute bsum1 snapgroup1 reduce sum c_b[*]
|
||||
compute bsum2 snapgroup2 reduce sum c_b[*]
|
||||
# fix bsum1 all ave/time 1 1 1 c_bsum1 file bsum1.dat mode vector
|
||||
# fix bsum2 all ave/time 1 1 1 c_bsum2 file bsum2.dat mode vector
|
||||
compute vbsum all reduce sum c_vb[*]
|
||||
# fix vbsum all ave/time 1 1 1 c_vbsum file vbsum.dat mode vector
|
||||
variable db_2_25 equal c_db[2][25]
|
||||
|
||||
# set up compute snap generating global array
|
||||
|
||||
compute snap all mliap descriptor sna compute.mliap.descriptor model linear 2 6
|
||||
SNAP keyword rcutfac 1.0
|
||||
SNAP keyword twojmax 2
|
||||
SNAP keyword nelems 2
|
||||
SNAP keyword elems A
|
||||
SNAP keyword radelems 2.3
|
||||
SNAP keyword welems 1.0
|
||||
SNAP keyword rfac0 0.99363
|
||||
SNAP keyword rmin0 0
|
||||
SNAP keyword bzeroflag 0
|
||||
SNAP keyword switchflag 0
|
||||
fix snap all ave/time 1 1 1 c_snap[*] file compute.snap.dat mode vector
|
||||
|
||||
thermo 100
|
||||
|
||||
# test output: 1: total potential energy
|
||||
# 2: xy component of stress tensor
|
||||
# 3: Sum(B_{000}^i, all i of type 2)
|
||||
# 4: xz component of Sum(Sum(r_j*dB_{222}^i/dR[j]), all i of type 2), all j)
|
||||
# 5: y component of -Sum(d(B_{222}^i)/dR[2]), all i of type 2)
|
||||
#
|
||||
# followed by 5 counterparts from compute snap
|
||||
|
||||
thermo_style custom pe pxy c_bsum2[1] c_vbsum[55] v_db_2_25 c_snap[1][13] c_snap[13][13] c_snap[1][8] c_snap[12][12] c_snap[6][12]
|
||||
thermo_modify norm no
|
||||
|
||||
# dump mydump_db all custom 1000 dump_db id c_db[*]
|
||||
# dump_modify mydump_db sort id
|
||||
|
||||
# Run MD
|
||||
|
||||
run ${nsteps}
|
||||
run 0
|
||||
Neighbor list info ...
|
||||
update every 1 steps, delay 10 steps, check yes
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
master list distance cutoff = 6.8
|
||||
ghost atom cutoff = 6.8
|
||||
binsize = 3.4, bins = 1 1 1
|
||||
5 neighbor lists, perpetual/occasional/extra = 1 4 0
|
||||
(1) pair zbl, perpetual
|
||||
attributes: half, newton on
|
||||
pair build: half/bin/atomonly/newton
|
||||
stencil: half/bin/3d/newton
|
||||
bin: standard
|
||||
(2) compute sna/atom, occasional
|
||||
attributes: full, newton on
|
||||
pair build: full/bin/atomonly
|
||||
stencil: full/bin/3d
|
||||
bin: standard
|
||||
(3) compute snav/atom, occasional
|
||||
attributes: full, newton on
|
||||
pair build: full/bin/atomonly
|
||||
stencil: full/bin/3d
|
||||
bin: standard
|
||||
(4) compute snad/atom, occasional
|
||||
attributes: full, newton on
|
||||
pair build: full/bin/atomonly
|
||||
stencil: full/bin/3d
|
||||
bin: standard
|
||||
(5) compute mliap, occasional
|
||||
attributes: full, newton on
|
||||
pair build: full/bin/atomonly
|
||||
stencil: full/bin/3d
|
||||
bin: standard
|
||||
WARNING: Proc sub-domain size < neighbor skin, could lead to lost atoms (../domain.cpp:964)
|
||||
Per MPI rank memory allocation (min/avg/max) = 10.69 | 10.69 | 10.7 Mbytes
|
||||
PotEng Pxy c_bsum2[1] c_vbsum[55] v_db_2_25 c_snap[1][13] c_snap[13][13] c_snap[1][8] c_snap[12][12] c_snap[6][12]
|
||||
322.86952 1505558.1 364182.88 -240.25066 1381.7961 322.86952 1505558.1 364182.88 -240.25066 1381.7961
|
||||
Loop time of 3.25e-06 on 4 procs for 0 steps with 2 atoms
|
||||
|
||||
115.4% CPU use with 4 MPI tasks x no OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Comm | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Output | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Modify | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Other | | 3.25e-06 | | |100.00
|
||||
|
||||
Nlocal: 0.5 ave 1 max 0 min
|
||||
Histogram: 2 0 0 0 0 0 0 0 0 2
|
||||
Nghost: 734.5 ave 735 max 734 min
|
||||
Histogram: 2 0 0 0 0 0 0 0 0 2
|
||||
Neighs: 82.5 ave 177 max 0 min
|
||||
Histogram: 2 0 0 0 0 0 0 0 1 1
|
||||
FullNghs: 165 ave 330 max 0 min
|
||||
Histogram: 2 0 0 0 0 0 0 0 0 2
|
||||
|
||||
Total # of neighbors = 660
|
||||
Ave neighs/atom = 330
|
||||
Neighbor list builds = 0
|
||||
Dangerous builds = 0
|
||||
Total wall time: 0:00:00
|
||||
|
|
@ -70,7 +70,7 @@ compute bsum2 snapgroup2 reduce sum c_b[*]
|
|||
# fix bsum2 all ave/time 1 1 1 c_bsum2 file bsum2.dat mode vector
|
||||
compute vbsum all reduce sum c_vb[*]
|
||||
# fix vbsum all ave/time 1 1 1 c_vbsum file vbsum.dat mode vector
|
||||
variable db_2_30 equal c_db[2][30]
|
||||
variable db_2_25 equal c_db[2][25]
|
||||
|
||||
# set up compute snap generating global array
|
||||
|
||||
|
|
@ -82,14 +82,14 @@ thermo 100
|
|||
# test output: 1: total potential energy
|
||||
# 2: xy component of stress tensor
|
||||
# 3: Sum(B_{000}^i, all i of type 2)
|
||||
# 4: xy component of Sum(Sum(r_j*dB_{222}^i/dR[j]), all i of type 2), all j)
|
||||
# 5: z component of -Sum(d(B_{222}^i)/dR[2]), all i of type 2)
|
||||
# 4: xz component of Sum(Sum(r_j*dB_{222}^i/dR[j]), all i of type 2), all j)
|
||||
# 5: y component of -Sum(d(B_{222}^i)/dR[2]), all i of type 2)
|
||||
#
|
||||
# followed by 5 counterparts from compute snap
|
||||
|
||||
thermo_style custom &
|
||||
pe pxy c_bsum2[1] c_vbsum[60] v_db_2_30 &
|
||||
c_snap[1][11] c_snap[13][11] c_snap[1][6] c_snap[13][10] c_snap[7][10]
|
||||
pe pxy c_bsum2[1] c_vbsum[55] v_db_2_25 &
|
||||
c_snap[1][11] c_snap[13][11] c_snap[1][6] c_snap[12][10] c_snap[6][10]
|
||||
thermo_modify norm no
|
||||
|
||||
# dump mydump_db all custom 1000 dump_db id c_db[*]
|
||||
|
|
|
|||
|
|
@ -70,7 +70,7 @@ compute bsum2 snapgroup2 reduce sum c_b[*]
|
|||
# fix bsum2 all ave/time 1 1 1 c_bsum2 file bsum2.dat mode vector
|
||||
compute vbsum all reduce sum c_vb[*]
|
||||
# fix vbsum all ave/time 1 1 1 c_vbsum file vbsum.dat mode vector
|
||||
variable db_2_120 equal c_db[2][120]
|
||||
variable db_2_100 equal c_db[2][100]
|
||||
|
||||
# set up compute snap generating global array
|
||||
|
||||
|
|
@ -82,14 +82,14 @@ thermo 100
|
|||
# test output: 1: total potential energy
|
||||
# 2: xy component of stress tensor
|
||||
# 3: Sum(0.5*(B_{222}^i)^2, all i of type 2)
|
||||
# 4: xy component of Sum(Sum(r_j*(0.5*(dB_{222}^i)^2/dR[j]), all i of type 2), all j)
|
||||
# 5: z component of -Sum(d(0.5*(B_{222}^i)^2/dR[2]), all i of type 2)
|
||||
# 4: xz component of Sum(Sum(r_j*(0.5*(dB_{222}^i)^2/dR[j]), all i of type 2), all j)
|
||||
# 5: y component of -Sum(d(0.5*(B_{222}^i)^2/dR[2]), all i of type 2)
|
||||
#
|
||||
# followed by 5 counterparts from compute snap
|
||||
|
||||
thermo_style custom &
|
||||
pe pxy c_bsum2[20] c_vbsum[240] v_db_2_120 &
|
||||
c_snap[1][41] c_snap[13][41] c_snap[1][40] c_snap[13][40] c_snap[7][40]
|
||||
pe pxy c_bsum2[20] c_vbsum[220] v_db_2_100 &
|
||||
c_snap[1][41] c_snap[13][41] c_snap[1][40] c_snap[12][40] c_snap[6][40]
|
||||
thermo_modify norm no
|
||||
|
||||
# dump mydump_db all custom 1000 dump_db id c_db[*]
|
||||
|
|
|
|||
|
|
@ -1,6 +1,4 @@
|
|||
LAMMPS (20 Nov 2019)
|
||||
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:93)
|
||||
using 1 OpenMP thread(s) per MPI task
|
||||
LAMMPS (15 Jun 2020)
|
||||
# Demonstrate bispectrum computes
|
||||
|
||||
# initialize simulation
|
||||
|
|
@ -30,11 +28,11 @@ region box block 0 1 0 ${ny} 0 ${nz}
|
|||
region box block 0 1 0 1 0 ${nz}
|
||||
region box block 0 1 0 1 0 1
|
||||
create_box 2 box
|
||||
Created orthogonal box = (0 0 0) to (2 2 2)
|
||||
Created orthogonal box = (0.0 0.0 0.0) to (2.0 2.0 2.0)
|
||||
1 by 1 by 1 MPI processor grid
|
||||
create_atoms 2 box
|
||||
Created 2 atoms
|
||||
create_atoms CPU = 0.000478029 secs
|
||||
create_atoms CPU = 0.000 seconds
|
||||
|
||||
mass * 180.88
|
||||
|
||||
|
|
@ -105,7 +103,7 @@ compute bsum2 snapgroup2 reduce sum c_b[*]
|
|||
# fix bsum2 all ave/time 1 1 1 c_bsum2 file bsum2.dat mode vector
|
||||
compute vbsum all reduce sum c_vb[*]
|
||||
# fix vbsum all ave/time 1 1 1 c_vbsum file vbsum.dat mode vector
|
||||
variable db_2_30 equal c_db[2][30]
|
||||
variable db_2_25 equal c_db[2][25]
|
||||
|
||||
# set up compute snap generating global array
|
||||
|
||||
|
|
@ -118,12 +116,12 @@ thermo 100
|
|||
# test output: 1: total potential energy
|
||||
# 2: xy component of stress tensor
|
||||
# 3: Sum(B_{000}^i, all i of type 2)
|
||||
# 4: xy component of Sum(Sum(r_j*dB_{222}^i/dR[j]), all i of type 2), all j)
|
||||
# 5: z component of -Sum(d(B_{222}^i)/dR[2]), all i of type 2)
|
||||
# 4: xz component of Sum(Sum(r_j*dB_{222}^i/dR[j]), all i of type 2), all j)
|
||||
# 5: y component of -Sum(d(B_{222}^i)/dR[2]), all i of type 2)
|
||||
#
|
||||
# followed by 5 counterparts from compute snap
|
||||
|
||||
thermo_style custom pe pxy c_bsum2[1] c_vbsum[60] v_db_2_30 c_snap[1][11] c_snap[13][11] c_snap[1][6] c_snap[13][10] c_snap[7][10]
|
||||
thermo_style custom pe pxy c_bsum2[1] c_vbsum[55] v_db_2_25 c_snap[1][11] c_snap[13][11] c_snap[1][6] c_snap[12][10] c_snap[6][10]
|
||||
thermo_modify norm no
|
||||
|
||||
# dump mydump_db all custom 1000 dump_db id c_db[*]
|
||||
|
|
@ -166,11 +164,11 @@ Neighbor list info ...
|
|||
stencil: full/bin/3d
|
||||
bin: standard
|
||||
Per MPI rank memory allocation (min/avg/max) = 10.06 | 10.06 | 10.06 Mbytes
|
||||
PotEng Pxy c_bsum2[1] c_vbsum[60] v_db_2_30 c_snap[1][11] c_snap[13][11] c_snap[1][6] c_snap[13][10] c_snap[7][10]
|
||||
322.86952 1505558.1 364182.88 381.32218 -855.04473 322.86952 1505558.1 364182.88 381.32218 -855.04473
|
||||
Loop time of 9.53674e-07 on 1 procs for 0 steps with 2 atoms
|
||||
PotEng Pxy c_bsum2[1] c_vbsum[55] v_db_2_25 c_snap[1][11] c_snap[13][11] c_snap[1][6] c_snap[12][10] c_snap[6][10]
|
||||
322.86952 1505558.1 364182.88 -240.25066 1381.7961 322.86952 1505558.1 364182.88 -240.25066 1381.7961
|
||||
Loop time of 1e-06 on 1 procs for 0 steps with 2 atoms
|
||||
|
||||
104.9% CPU use with 1 MPI tasks x 1 OpenMP threads
|
||||
200.0% CPU use with 1 MPI tasks x no OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
|
|
@ -180,7 +178,7 @@ Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
|||
Comm | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Output | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Modify | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Other | | 9.537e-07 | | |100.00
|
||||
Other | | 1e-06 | | |100.00
|
||||
|
||||
Nlocal: 2 ave 2 max 2 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
|
@ -1,6 +1,4 @@
|
|||
LAMMPS (20 Nov 2019)
|
||||
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:93)
|
||||
using 1 OpenMP thread(s) per MPI task
|
||||
LAMMPS (15 Jun 2020)
|
||||
# Demonstrate bispectrum computes
|
||||
|
||||
# initialize simulation
|
||||
|
|
@ -30,11 +28,11 @@ region box block 0 1 0 ${ny} 0 ${nz}
|
|||
region box block 0 1 0 1 0 ${nz}
|
||||
region box block 0 1 0 1 0 1
|
||||
create_box 2 box
|
||||
Created orthogonal box = (0 0 0) to (2 2 2)
|
||||
Created orthogonal box = (0.0 0.0 0.0) to (2.0 2.0 2.0)
|
||||
1 by 2 by 2 MPI processor grid
|
||||
create_atoms 2 box
|
||||
Created 2 atoms
|
||||
create_atoms CPU = 0.000610113 secs
|
||||
create_atoms CPU = 0.001 seconds
|
||||
|
||||
mass * 180.88
|
||||
|
||||
|
|
@ -105,7 +103,7 @@ compute bsum2 snapgroup2 reduce sum c_b[*]
|
|||
# fix bsum2 all ave/time 1 1 1 c_bsum2 file bsum2.dat mode vector
|
||||
compute vbsum all reduce sum c_vb[*]
|
||||
# fix vbsum all ave/time 1 1 1 c_vbsum file vbsum.dat mode vector
|
||||
variable db_2_30 equal c_db[2][30]
|
||||
variable db_2_25 equal c_db[2][25]
|
||||
|
||||
# set up compute snap generating global array
|
||||
|
||||
|
|
@ -118,12 +116,12 @@ thermo 100
|
|||
# test output: 1: total potential energy
|
||||
# 2: xy component of stress tensor
|
||||
# 3: Sum(B_{000}^i, all i of type 2)
|
||||
# 4: xy component of Sum(Sum(r_j*dB_{222}^i/dR[j]), all i of type 2), all j)
|
||||
# 5: z component of -Sum(d(B_{222}^i)/dR[2]), all i of type 2)
|
||||
# 4: xz component of Sum(Sum(r_j*dB_{222}^i/dR[j]), all i of type 2), all j)
|
||||
# 5: y component of -Sum(d(B_{222}^i)/dR[2]), all i of type 2)
|
||||
#
|
||||
# followed by 5 counterparts from compute snap
|
||||
|
||||
thermo_style custom pe pxy c_bsum2[1] c_vbsum[60] v_db_2_30 c_snap[1][11] c_snap[13][11] c_snap[1][6] c_snap[13][10] c_snap[7][10]
|
||||
thermo_style custom pe pxy c_bsum2[1] c_vbsum[55] v_db_2_25 c_snap[1][11] c_snap[13][11] c_snap[1][6] c_snap[12][10] c_snap[6][10]
|
||||
thermo_modify norm no
|
||||
|
||||
# dump mydump_db all custom 1000 dump_db id c_db[*]
|
||||
|
|
@ -165,13 +163,13 @@ Neighbor list info ...
|
|||
pair build: full/bin/atomonly
|
||||
stencil: full/bin/3d
|
||||
bin: standard
|
||||
WARNING: Proc sub-domain size < neighbor skin, could lead to lost atoms (../domain.cpp:936)
|
||||
Per MPI rank memory allocation (min/avg/max) = 8.211 | 8.254 | 8.295 Mbytes
|
||||
PotEng Pxy c_bsum2[1] c_vbsum[60] v_db_2_30 c_snap[1][11] c_snap[13][11] c_snap[1][6] c_snap[13][10] c_snap[7][10]
|
||||
322.86952 1505558.1 364182.88 381.32218 -855.04473 322.86952 1505558.1 364182.88 381.32218 -855.04473
|
||||
Loop time of 2.38419e-06 on 4 procs for 0 steps with 2 atoms
|
||||
WARNING: Proc sub-domain size < neighbor skin, could lead to lost atoms (../domain.cpp:964)
|
||||
Per MPI rank memory allocation (min/avg/max) = 9.945 | 9.988 | 10.03 Mbytes
|
||||
PotEng Pxy c_bsum2[1] c_vbsum[55] v_db_2_25 c_snap[1][11] c_snap[13][11] c_snap[1][6] c_snap[12][10] c_snap[6][10]
|
||||
322.86952 1505558.1 364182.88 -240.25066 1381.7961 322.86952 1505558.1 364182.88 -240.25066 1381.7961
|
||||
Loop time of 2.75e-06 on 4 procs for 0 steps with 2 atoms
|
||||
|
||||
104.9% CPU use with 4 MPI tasks x 1 OpenMP threads
|
||||
90.9% CPU use with 4 MPI tasks x no OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
|
|
@ -181,7 +179,7 @@ Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
|||
Comm | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Output | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Modify | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Other | | 2.384e-06 | | |100.00
|
||||
Other | | 2.75e-06 | | |100.00
|
||||
|
||||
Nlocal: 0.5 ave 1 max 0 min
|
||||
Histogram: 2 0 0 0 0 0 0 0 0 2
|
||||
|
|
@ -1,6 +1,4 @@
|
|||
LAMMPS (20 Nov 2019)
|
||||
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:93)
|
||||
using 1 OpenMP thread(s) per MPI task
|
||||
LAMMPS (15 Jun 2020)
|
||||
# Demonstrate bispectrum computes
|
||||
|
||||
# initialize simulation
|
||||
|
|
@ -30,11 +28,11 @@ region box block 0 1 0 ${ny} 0 ${nz}
|
|||
region box block 0 1 0 1 0 ${nz}
|
||||
region box block 0 1 0 1 0 1
|
||||
create_box 2 box
|
||||
Created orthogonal box = (0 0 0) to (2 2 2)
|
||||
Created orthogonal box = (0.0 0.0 0.0) to (2.0 2.0 2.0)
|
||||
1 by 1 by 1 MPI processor grid
|
||||
create_atoms 2 box
|
||||
Created 2 atoms
|
||||
create_atoms CPU = 0.000473976 secs
|
||||
create_atoms CPU = 0.001 seconds
|
||||
|
||||
mass * 180.88
|
||||
|
||||
|
|
@ -105,7 +103,7 @@ compute bsum2 snapgroup2 reduce sum c_b[*]
|
|||
# fix bsum2 all ave/time 1 1 1 c_bsum2 file bsum2.dat mode vector
|
||||
compute vbsum all reduce sum c_vb[*]
|
||||
# fix vbsum all ave/time 1 1 1 c_vbsum file vbsum.dat mode vector
|
||||
variable db_2_120 equal c_db[2][120]
|
||||
variable db_2_100 equal c_db[2][100]
|
||||
|
||||
# set up compute snap generating global array
|
||||
|
||||
|
|
@ -118,12 +116,12 @@ thermo 100
|
|||
# test output: 1: total potential energy
|
||||
# 2: xy component of stress tensor
|
||||
# 3: Sum(0.5*(B_{222}^i)^2, all i of type 2)
|
||||
# 4: xy component of Sum(Sum(r_j*(0.5*(dB_{222}^i)^2/dR[j]), all i of type 2), all j)
|
||||
# 5: z component of -Sum(d(0.5*(B_{222}^i)^2/dR[2]), all i of type 2)
|
||||
# 4: xz component of Sum(Sum(r_j*(0.5*(dB_{222}^i)^2/dR[j]), all i of type 2), all j)
|
||||
# 5: y component of -Sum(d(0.5*(B_{222}^i)^2/dR[2]), all i of type 2)
|
||||
#
|
||||
# followed by 5 counterparts from compute snap
|
||||
|
||||
thermo_style custom pe pxy c_bsum2[20] c_vbsum[240] v_db_2_120 c_snap[1][41] c_snap[13][41] c_snap[1][40] c_snap[13][40] c_snap[7][40]
|
||||
thermo_style custom pe pxy c_bsum2[20] c_vbsum[220] v_db_2_100 c_snap[1][41] c_snap[13][41] c_snap[1][40] c_snap[12][40] c_snap[6][40]
|
||||
thermo_modify norm no
|
||||
|
||||
# dump mydump_db all custom 1000 dump_db id c_db[*]
|
||||
|
|
@ -166,11 +164,11 @@ Neighbor list info ...
|
|||
stencil: full/bin/3d
|
||||
bin: standard
|
||||
Per MPI rank memory allocation (min/avg/max) = 21.78 | 21.78 | 21.78 Mbytes
|
||||
PotEng Pxy c_bsum2[20] c_vbsum[240] v_db_2_120 c_snap[1][41] c_snap[13][41] c_snap[1][40] c_snap[13][40] c_snap[7][40]
|
||||
322.86952 1505558.1 4.2492771e+08 7860489.6 -17625699 322.86952 1505558.1 4.2492771e+08 7860489.6 -17625699
|
||||
Loop time of 2.14577e-06 on 1 procs for 0 steps with 2 atoms
|
||||
PotEng Pxy c_bsum2[20] c_vbsum[220] v_db_2_100 c_snap[1][41] c_snap[13][41] c_snap[1][40] c_snap[12][40] c_snap[6][40]
|
||||
322.86952 1505558.1 4.2492771e+08 -4952473 28484035 322.86952 1505558.1 4.2492771e+08 -4952473 28484035
|
||||
Loop time of 1e-06 on 1 procs for 0 steps with 2 atoms
|
||||
|
||||
93.2% CPU use with 1 MPI tasks x 1 OpenMP threads
|
||||
200.0% CPU use with 1 MPI tasks x no OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
|
|
@ -180,7 +178,7 @@ Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
|||
Comm | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Output | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Modify | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Other | | 2.146e-06 | | |100.00
|
||||
Other | | 1e-06 | | |100.00
|
||||
|
||||
Nlocal: 2 ave 2 max 2 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
|
@ -1,6 +1,4 @@
|
|||
LAMMPS (20 Nov 2019)
|
||||
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:93)
|
||||
using 1 OpenMP thread(s) per MPI task
|
||||
LAMMPS (15 Jun 2020)
|
||||
# Demonstrate bispectrum computes
|
||||
|
||||
# initialize simulation
|
||||
|
|
@ -30,11 +28,11 @@ region box block 0 1 0 ${ny} 0 ${nz}
|
|||
region box block 0 1 0 1 0 ${nz}
|
||||
region box block 0 1 0 1 0 1
|
||||
create_box 2 box
|
||||
Created orthogonal box = (0 0 0) to (2 2 2)
|
||||
Created orthogonal box = (0.0 0.0 0.0) to (2.0 2.0 2.0)
|
||||
1 by 2 by 2 MPI processor grid
|
||||
create_atoms 2 box
|
||||
Created 2 atoms
|
||||
create_atoms CPU = 0.000118971 secs
|
||||
create_atoms CPU = 0.001 seconds
|
||||
|
||||
mass * 180.88
|
||||
|
||||
|
|
@ -105,7 +103,7 @@ compute bsum2 snapgroup2 reduce sum c_b[*]
|
|||
# fix bsum2 all ave/time 1 1 1 c_bsum2 file bsum2.dat mode vector
|
||||
compute vbsum all reduce sum c_vb[*]
|
||||
# fix vbsum all ave/time 1 1 1 c_vbsum file vbsum.dat mode vector
|
||||
variable db_2_120 equal c_db[2][120]
|
||||
variable db_2_100 equal c_db[2][100]
|
||||
|
||||
# set up compute snap generating global array
|
||||
|
||||
|
|
@ -118,12 +116,12 @@ thermo 100
|
|||
# test output: 1: total potential energy
|
||||
# 2: xy component of stress tensor
|
||||
# 3: Sum(0.5*(B_{222}^i)^2, all i of type 2)
|
||||
# 4: xy component of Sum(Sum(r_j*(0.5*(dB_{222}^i)^2/dR[j]), all i of type 2), all j)
|
||||
# 5: z component of -Sum(d(0.5*(B_{222}^i)^2/dR[2]), all i of type 2)
|
||||
# 4: xz component of Sum(Sum(r_j*(0.5*(dB_{222}^i)^2/dR[j]), all i of type 2), all j)
|
||||
# 5: y component of -Sum(d(0.5*(B_{222}^i)^2/dR[2]), all i of type 2)
|
||||
#
|
||||
# followed by 5 counterparts from compute snap
|
||||
|
||||
thermo_style custom pe pxy c_bsum2[20] c_vbsum[240] v_db_2_120 c_snap[1][41] c_snap[13][41] c_snap[1][40] c_snap[13][40] c_snap[7][40]
|
||||
thermo_style custom pe pxy c_bsum2[20] c_vbsum[220] v_db_2_100 c_snap[1][41] c_snap[13][41] c_snap[1][40] c_snap[12][40] c_snap[6][40]
|
||||
thermo_modify norm no
|
||||
|
||||
# dump mydump_db all custom 1000 dump_db id c_db[*]
|
||||
|
|
@ -165,13 +163,13 @@ Neighbor list info ...
|
|||
pair build: full/bin/atomonly
|
||||
stencil: full/bin/3d
|
||||
bin: standard
|
||||
WARNING: Proc sub-domain size < neighbor skin, could lead to lost atoms (../domain.cpp:936)
|
||||
Per MPI rank memory allocation (min/avg/max) = 19.7 | 19.74 | 19.78 Mbytes
|
||||
PotEng Pxy c_bsum2[20] c_vbsum[240] v_db_2_120 c_snap[1][41] c_snap[13][41] c_snap[1][40] c_snap[13][40] c_snap[7][40]
|
||||
322.86952 1505558.1 4.2492771e+08 7860489.6 -17625699 322.86952 1505558.1 4.2492771e+08 7860489.6 -17625699
|
||||
Loop time of 2.80142e-06 on 4 procs for 0 steps with 2 atoms
|
||||
WARNING: Proc sub-domain size < neighbor skin, could lead to lost atoms (../domain.cpp:964)
|
||||
Per MPI rank memory allocation (min/avg/max) = 21.43 | 21.47 | 21.52 Mbytes
|
||||
PotEng Pxy c_bsum2[20] c_vbsum[220] v_db_2_100 c_snap[1][41] c_snap[13][41] c_snap[1][40] c_snap[12][40] c_snap[6][40]
|
||||
322.86952 1505558.1 4.2492771e+08 -4952473 28484035 322.86952 1505558.1 4.2492771e+08 -4952473 28484035
|
||||
Loop time of 2.25e-06 on 4 procs for 0 steps with 2 atoms
|
||||
|
||||
107.1% CPU use with 4 MPI tasks x 1 OpenMP threads
|
||||
111.1% CPU use with 4 MPI tasks x no OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
|
|
@ -181,7 +179,7 @@ Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
|||
Comm | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Output | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Modify | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Other | | 2.801e-06 | | |100.00
|
||||
Other | | 2.25e-06 | | |100.00
|
||||
|
||||
Nlocal: 0.5 ave 1 max 0 min
|
||||
Histogram: 2 0 0 0 0 0 0 0 0 2
|
||||
|
|
@ -0,0 +1,30 @@
|
|||
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. Currently, the pair_style
|
||||
supports just two models, *linear* and *quadratic*,
|
||||
and one descriptor, *sna*, the SNAP descriptor, including the
|
||||
linear, quadratic, and chem variants. Work is currently underway to extend
|
||||
the interface to handle neural network energy models,
|
||||
and it is also straightforward to add new descriptor styles.
|
||||
|
||||
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>.
|
||||
|
|
@ -0,0 +1,390 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include <cstring>
|
||||
#include <cstdlib>
|
||||
#include "mliap_model_linear.h"
|
||||
#include "mliap_model_quadratic.h"
|
||||
#include "mliap_descriptor_snap.h"
|
||||
#include "compute_mliap.h"
|
||||
#include "atom.h"
|
||||
#include "update.h"
|
||||
#include "modify.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "force.h"
|
||||
#include "pair.h"
|
||||
#include "comm.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
enum{SCALAR,VECTOR,ARRAY};
|
||||
|
||||
ComputeMLIAP::ComputeMLIAP(LAMMPS *lmp, int narg, char **arg) :
|
||||
Compute(lmp, narg, arg), list(NULL), mliap(NULL),
|
||||
gradforce(NULL), mliapall(NULL), map(NULL),
|
||||
descriptors(NULL), gamma_row_index(NULL), gamma_col_index(NULL),
|
||||
gamma(NULL), egradient(NULL), model(NULL), descriptor(NULL)
|
||||
{
|
||||
array_flag = 1;
|
||||
extarray = 0;
|
||||
|
||||
if (narg < 4)
|
||||
error->all(FLERR,"Illegal compute mliap command");
|
||||
|
||||
// set flags for required keywords
|
||||
|
||||
int modelflag = 0;
|
||||
int descriptorflag = 0;
|
||||
|
||||
// process keywords
|
||||
|
||||
int iarg = 3;
|
||||
|
||||
while (iarg < narg) {
|
||||
if (strcmp(arg[iarg],"model") == 0) {
|
||||
if (iarg+2 > narg) error->all(FLERR,"Illegal compute mliap command");
|
||||
if (strcmp(arg[iarg+1],"linear") == 0) {
|
||||
if (iarg+4 > narg) error->all(FLERR,"Illegal compute mliap command");
|
||||
int ntmp1 = atoi(arg[iarg+2]);
|
||||
int ntmp2 = atoi(arg[iarg+3]);
|
||||
model = new MLIAPModelLinear(lmp,ntmp1,ntmp2);
|
||||
iarg += 4;
|
||||
} else if (strcmp(arg[iarg+1],"quadratic") == 0) {
|
||||
if (iarg+4 > narg) error->all(FLERR,"Illegal compute mliap command");
|
||||
int ntmp1 = atoi(arg[iarg+2]);
|
||||
int ntmp2 = atoi(arg[iarg+3]);
|
||||
model = new MLIAPModelQuadratic(lmp,ntmp1,ntmp2);
|
||||
iarg += 4;
|
||||
} else error->all(FLERR,"Illegal compute mliap command");
|
||||
modelflag = 1;
|
||||
} else if (strcmp(arg[iarg],"descriptor") == 0) {
|
||||
if (iarg+2 > narg) error->all(FLERR,"Illegal compute mliap command");
|
||||
if (strcmp(arg[iarg+1],"sna") == 0) {
|
||||
if (iarg+3 > narg) error->all(FLERR,"Illegal compute mliap command");
|
||||
descriptor = new MLIAPDescriptorSNAP(lmp,arg[iarg+2]);
|
||||
iarg += 3;
|
||||
} else error->all(FLERR,"Illegal compute mliap command");
|
||||
descriptorflag = 1;
|
||||
} else
|
||||
error->all(FLERR,"Illegal compute mliap command");
|
||||
}
|
||||
|
||||
if (modelflag == 0 || descriptorflag == 0)
|
||||
error->all(FLERR,"Illegal compute_style command");
|
||||
|
||||
ndescriptors = descriptor->ndescriptors;
|
||||
nparams = model->nparams;
|
||||
nelements = model->nelements;
|
||||
gamma_nnz = model->get_gamma_nnz();
|
||||
ndims_force = 3;
|
||||
ndims_virial = 6;
|
||||
yoffset = nparams*nelements;
|
||||
zoffset = 2*yoffset;
|
||||
natoms = atom->natoms;
|
||||
size_array_rows = 1+ndims_force*natoms+ndims_virial;
|
||||
size_array_cols = nparams*nelements+1;
|
||||
lastcol = size_array_cols-1;
|
||||
|
||||
size_gradforce = ndims_force*nparams*nelements;
|
||||
|
||||
nmax = 0;
|
||||
gamma_max = 0;
|
||||
|
||||
// create a minimal map, placeholder for more general map
|
||||
|
||||
memory->create(map,atom->ntypes+1,"compute_mliap:map");
|
||||
|
||||
for (int i = 1; i <= atom->ntypes; i++)
|
||||
map[i] = i-1;
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ComputeMLIAP::~ComputeMLIAP()
|
||||
{
|
||||
memory->destroy(mliap);
|
||||
memory->destroy(mliapall);
|
||||
memory->destroy(gradforce);
|
||||
|
||||
memory->destroy(map);
|
||||
|
||||
memory->destroy(descriptors);
|
||||
memory->destroy(gamma_row_index);
|
||||
memory->destroy(gamma_col_index);
|
||||
memory->destroy(gamma);
|
||||
memory->destroy(egradient);
|
||||
|
||||
delete model;
|
||||
delete descriptor;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ComputeMLIAP::init()
|
||||
{
|
||||
if (force->pair == NULL)
|
||||
error->all(FLERR,"Compute mliap requires a pair style be defined");
|
||||
|
||||
if (descriptor->cutmax > force->pair->cutforce)
|
||||
error->all(FLERR,"Compute mliap cutoff is longer than pairwise cutoff");
|
||||
|
||||
// need an occasional full neighbor list
|
||||
|
||||
int irequest = neighbor->request(this,instance_me);
|
||||
neighbor->requests[irequest]->pair = 0;
|
||||
neighbor->requests[irequest]->compute = 1;
|
||||
neighbor->requests[irequest]->half = 0;
|
||||
neighbor->requests[irequest]->full = 1;
|
||||
neighbor->requests[irequest]->occasional = 1;
|
||||
|
||||
int count = 0;
|
||||
for (int i = 0; i < modify->ncompute; i++)
|
||||
if (strcmp(modify->compute[i]->style,"mliap") == 0) count++;
|
||||
if (count > 1 && comm->me == 0)
|
||||
error->warning(FLERR,"More than one compute mliap");
|
||||
|
||||
// initialize model and descriptor
|
||||
|
||||
model->init();
|
||||
descriptor->init();
|
||||
|
||||
// consistency checks
|
||||
|
||||
if (descriptor->ndescriptors != model->ndescriptors)
|
||||
error->all(FLERR,"Incompatible model and descriptor definitions");
|
||||
if (descriptor->nelements != model->nelements)
|
||||
error->all(FLERR,"Incompatible model and descriptor definitions");
|
||||
if (nelements != atom->ntypes)
|
||||
error->all(FLERR,"nelements must equal ntypes");
|
||||
|
||||
// allocate memory for global array
|
||||
|
||||
memory->create(mliap,size_array_rows,size_array_cols,
|
||||
"mliap:mliap");
|
||||
memory->create(mliapall,size_array_rows,size_array_cols,
|
||||
"mliap:mliapall");
|
||||
array = mliapall;
|
||||
|
||||
memory->create(egradient,nelements*nparams,"ComputeMLIAP:egradient");
|
||||
|
||||
// find compute for reference energy
|
||||
|
||||
std::string id_pe = std::string("thermo_pe");
|
||||
int ipe = modify->find_compute(id_pe);
|
||||
if (ipe == -1)
|
||||
error->all(FLERR,"compute thermo_pe does not exist.");
|
||||
c_pe = modify->compute[ipe];
|
||||
|
||||
// add compute for reference virial tensor
|
||||
|
||||
std::string id_virial = std::string("mliap_press");
|
||||
std::string pcmd = id_virial + " all pressure NULL virial";
|
||||
modify->add_compute(pcmd);
|
||||
|
||||
int ivirial = modify->find_compute(id_virial);
|
||||
if (ivirial == -1)
|
||||
error->all(FLERR,"compute mliap_press does not exist.");
|
||||
c_virial = modify->compute[ivirial];
|
||||
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ComputeMLIAP::init_list(int /*id*/, NeighList *ptr)
|
||||
{
|
||||
list = ptr;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ComputeMLIAP::compute_array()
|
||||
{
|
||||
int ntotal = atom->nlocal + atom->nghost;
|
||||
|
||||
invoked_array = update->ntimestep;
|
||||
|
||||
// grow gradforce array if necessary
|
||||
|
||||
if (atom->nmax > nmax) {
|
||||
memory->destroy(gradforce);
|
||||
nmax = atom->nmax;
|
||||
memory->create(gradforce,nmax,size_gradforce,
|
||||
"mliap:gradforce");
|
||||
}
|
||||
|
||||
// clear gradforce array
|
||||
|
||||
for (int i = 0; i < ntotal; i++)
|
||||
for (int j = 0; j < size_gradforce; j++) {
|
||||
gradforce[i][j] = 0.0;
|
||||
}
|
||||
|
||||
// clear global array
|
||||
|
||||
for (int irow = 0; irow < size_array_rows; irow++)
|
||||
for (int jcol = 0; jcol < size_array_cols; jcol++)
|
||||
mliap[irow][jcol] = 0.0;
|
||||
|
||||
// invoke full neighbor list (will copy or build if necessary)
|
||||
|
||||
neighbor->build_one(list);
|
||||
|
||||
if (gamma_max < list->inum) {
|
||||
memory->grow(descriptors,list->inum,ndescriptors,"ComputeMLIAP:descriptors");
|
||||
memory->grow(gamma_row_index,list->inum,gamma_nnz,"ComputeMLIAP:gamma_row_index");
|
||||
memory->grow(gamma_col_index,list->inum,gamma_nnz,"ComputeMLIAP:gamma_col_index");
|
||||
memory->grow(gamma,list->inum,gamma_nnz,"ComputeMLIAP:gamma");
|
||||
gamma_max = list->inum;
|
||||
}
|
||||
|
||||
// compute descriptors
|
||||
|
||||
descriptor->compute_descriptors(map, list, descriptors);
|
||||
|
||||
// calculate descriptor contributions to parameter gradients
|
||||
// and gamma = double gradient w.r.t. parameters and descriptors
|
||||
|
||||
// i.e. gamma = d2E/dsigma_l.dB_k
|
||||
// sigma_l is a parameter and B_k is a descriptor of atom i
|
||||
// for SNAP, this is a sparse natoms*nparams*ndescriptors matrix,
|
||||
// but in general it could be fully dense.
|
||||
|
||||
model->param_gradient(map, list, descriptors, gamma_row_index,
|
||||
gamma_col_index, gamma, egradient);
|
||||
|
||||
|
||||
// calculate descriptor gradient contributions to parameter gradients
|
||||
|
||||
descriptor->compute_gradients(map, list, gamma_nnz, gamma_row_index,
|
||||
gamma_col_index, gamma, gradforce,
|
||||
yoffset, zoffset);
|
||||
|
||||
// accumulate descriptor gradient contributions to global array
|
||||
|
||||
for (int ielem = 0; ielem < nelements; ielem++) {
|
||||
const int elemoffset = nparams*ielem;
|
||||
for (int jparam = 0; jparam < nparams; jparam++) {
|
||||
int irow = 1;
|
||||
for (int i = 0; i < ntotal; i++) {
|
||||
double *gradforcei = gradforce[i]+elemoffset;
|
||||
int iglobal = atom->tag[i];
|
||||
int irow = 3*(iglobal-1)+1;
|
||||
mliap[irow][jparam+elemoffset] += gradforcei[jparam];
|
||||
mliap[irow+1][jparam+elemoffset] += gradforcei[jparam+yoffset];
|
||||
mliap[irow+2][jparam+elemoffset] += gradforcei[jparam+zoffset];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// copy forces to global array
|
||||
|
||||
for (int i = 0; i < atom->nlocal; i++) {
|
||||
int iglobal = atom->tag[i];
|
||||
int irow = 3*(iglobal-1)+1;
|
||||
mliap[irow][lastcol] = atom->f[i][0];
|
||||
mliap[irow+1][lastcol] = atom->f[i][1];
|
||||
mliap[irow+2][lastcol] = atom->f[i][2];
|
||||
}
|
||||
|
||||
// accumulate bispectrum virial contributions to global array
|
||||
|
||||
dbdotr_compute();
|
||||
|
||||
// copy energy gradient contributions to global array
|
||||
|
||||
for (int ielem = 0; ielem < nelements; ielem++) {
|
||||
const int elemoffset = nparams*ielem;
|
||||
for (int jparam = 0; jparam < nparams; jparam++)
|
||||
mliap[0][jparam+elemoffset] = egradient[jparam+elemoffset];
|
||||
}
|
||||
|
||||
// sum up over all processes
|
||||
|
||||
MPI_Allreduce(&mliap[0][0],&mliapall[0][0],size_array_rows*size_array_cols,MPI_DOUBLE,MPI_SUM,world);
|
||||
|
||||
// copy energy to last column
|
||||
|
||||
int irow = 0;
|
||||
double reference_energy = c_pe->compute_scalar();
|
||||
mliapall[irow++][lastcol] = reference_energy;
|
||||
|
||||
// copy virial stress to last column
|
||||
// switch to Voigt notation
|
||||
|
||||
c_virial->compute_vector();
|
||||
irow += 3*natoms;
|
||||
mliapall[irow++][lastcol] = c_virial->vector[0];
|
||||
mliapall[irow++][lastcol] = c_virial->vector[1];
|
||||
mliapall[irow++][lastcol] = c_virial->vector[2];
|
||||
mliapall[irow++][lastcol] = c_virial->vector[5];
|
||||
mliapall[irow++][lastcol] = c_virial->vector[4];
|
||||
mliapall[irow++][lastcol] = c_virial->vector[3];
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
compute global virial contributions via summing r_i.dB^j/dr_i over
|
||||
own & ghost atoms
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void ComputeMLIAP::dbdotr_compute()
|
||||
{
|
||||
double **x = atom->x;
|
||||
int irow0 = 1+ndims_force*natoms;
|
||||
|
||||
// sum over bispectrum contributions to forces
|
||||
// on all particles including ghosts
|
||||
|
||||
int nall = atom->nlocal + atom->nghost;
|
||||
for (int i = 0; i < nall; i++)
|
||||
for (int ielem = 0; ielem < nelements; ielem++) {
|
||||
const int elemoffset = nparams*ielem;
|
||||
double *gradforcei = gradforce[i]+elemoffset;
|
||||
for (int jparam = 0; jparam < nparams; jparam++) {
|
||||
double dbdx = gradforcei[jparam];
|
||||
double dbdy = gradforcei[jparam+yoffset];
|
||||
double dbdz = gradforcei[jparam+zoffset];
|
||||
int irow = irow0;
|
||||
mliap[irow++][jparam+elemoffset] += dbdx*x[i][0];
|
||||
mliap[irow++][jparam+elemoffset] += dbdy*x[i][1];
|
||||
mliap[irow++][jparam+elemoffset] += dbdz*x[i][2];
|
||||
mliap[irow++][jparam+elemoffset] += dbdz*x[i][1];
|
||||
mliap[irow++][jparam+elemoffset] += dbdz*x[i][0];
|
||||
mliap[irow++][jparam+elemoffset] += dbdy*x[i][0];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
memory usage
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double ComputeMLIAP::memory_usage()
|
||||
{
|
||||
|
||||
double bytes = size_array_rows*size_array_cols *
|
||||
sizeof(double); // mliap
|
||||
bytes += size_array_rows*size_array_cols *
|
||||
sizeof(double); // mliapall
|
||||
bytes += nmax*size_gradforce * sizeof(double); // gradforce
|
||||
int n = atom->ntypes+1;
|
||||
bytes += n*sizeof(int); // map
|
||||
|
||||
return bytes;
|
||||
}
|
||||
|
|
@ -0,0 +1,90 @@
|
|||
/* -*- 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.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef COMPUTE_CLASS
|
||||
|
||||
ComputeStyle(mliap,ComputeMLIAP)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_COMPUTE_MLIAP_H
|
||||
#define LMP_COMPUTE_MLIAP_H
|
||||
|
||||
#include "compute.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class ComputeMLIAP : public Compute {
|
||||
public:
|
||||
ComputeMLIAP(class LAMMPS *, int, char **);
|
||||
~ComputeMLIAP();
|
||||
void init();
|
||||
void init_list(int, class NeighList *);
|
||||
void compute_array();
|
||||
double memory_usage();
|
||||
|
||||
private:
|
||||
int natoms, nmax, size_gradforce, lastcol;
|
||||
int yoffset, zoffset;
|
||||
int ndims_force, ndims_virial;
|
||||
class NeighList *list;
|
||||
double **mliap, **mliapall;
|
||||
double **gradforce;
|
||||
int *map; // map types to [0,nelements)
|
||||
int nelements;
|
||||
|
||||
double** descriptors; // descriptors for all atoms in list
|
||||
int ndescriptors; // number of descriptors
|
||||
int gamma_max; // number of atoms allocated for beta, descriptors
|
||||
int nparams; // number of model paramters per element
|
||||
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
|
||||
|
||||
class MLIAPModel* model;
|
||||
class MLIAPDescriptor* descriptor;
|
||||
|
||||
Compute *c_pe;
|
||||
Compute *c_virial;
|
||||
|
||||
void dbdotr_compute();
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
/* ERROR/WARNING messages:
|
||||
|
||||
E: Illegal ... command
|
||||
|
||||
Self-explanatory. Check the input script syntax and compare to the
|
||||
documentation for the command. You can use -echo screen as a
|
||||
command-line option when running LAMMPS to see the offending line.
|
||||
|
||||
E: Compute snap requires a pair style be defined
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Compute snap cutoff is longer than pairwise cutoff
|
||||
|
||||
UNDOCUMENTED
|
||||
|
||||
W: More than one compute snad/atom
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
*/
|
||||
|
|
@ -22,17 +22,20 @@ class MLIAPDescriptor : protected Pointers {
|
|||
public:
|
||||
MLIAPDescriptor(LAMMPS*);
|
||||
~MLIAPDescriptor();
|
||||
virtual void forward(int*, class NeighList*, double**)=0;
|
||||
virtual void backward(class PairMLIAP*, class NeighList*, double**, int)=0;
|
||||
virtual void compute_descriptors(int*, class NeighList*, double**)=0;
|
||||
virtual void compute_forces(class PairMLIAP*, class NeighList*, double**, int)=0;
|
||||
virtual void compute_gradients(int*, class NeighList*, int, int**, int**, double**,
|
||||
double**, int, int)=0;
|
||||
virtual void compute_descriptor_gradients(int*, class NeighList*, int, int**, int**, double**,
|
||||
double**, int, int)=0;
|
||||
virtual void init()=0;
|
||||
virtual double get_cutoff(int, int)=0;
|
||||
virtual double get_cutmax()=0;
|
||||
virtual double memory_usage()=0;
|
||||
|
||||
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
|
||||
protected:
|
||||
|
||||
};
|
||||
|
|
|
|||
|
|
@ -45,6 +45,13 @@ MLIAPDescriptorSNAP::MLIAPDescriptorSNAP(LAMMPS *lmp, char *paramfilename):
|
|||
wjelem = NULL;
|
||||
snaptr = NULL;
|
||||
read_paramfile(paramfilename);
|
||||
|
||||
snaptr = new SNA(lmp, rfac0, twojmax,
|
||||
rmin0, switchflag, bzeroflag,
|
||||
chemflag, bnormflag, wselfallflag, nelements);
|
||||
|
||||
ndescriptors = snaptr->ncoeff;
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
|
@ -58,6 +65,7 @@ MLIAPDescriptorSNAP::~MLIAPDescriptorSNAP()
|
|||
delete[] elements;
|
||||
memory->destroy(radelem);
|
||||
memory->destroy(wjelem);
|
||||
memory->destroy(cutsq);
|
||||
}
|
||||
|
||||
delete snaptr;
|
||||
|
|
@ -68,7 +76,7 @@ MLIAPDescriptorSNAP::~MLIAPDescriptorSNAP()
|
|||
compute descriptors for each atom
|
||||
---------------------------------------------------------------------- */
|
||||
|
||||
void MLIAPDescriptorSNAP::forward(int* map, NeighList* list, double **descriptors)
|
||||
void MLIAPDescriptorSNAP::compute_descriptors(int* map, NeighList* list, double **descriptors)
|
||||
{
|
||||
int i,j,jnum,ninside;
|
||||
double delx,dely,delz,rsq;
|
||||
|
|
@ -85,7 +93,6 @@ void MLIAPDescriptorSNAP::forward(int* map, NeighList* list, double **descriptor
|
|||
const double ztmp = x[i][2];
|
||||
const int itype = type[i];
|
||||
const int ielem = map[itype];
|
||||
const double radi = radelem[ielem];
|
||||
|
||||
jlist = list->firstneigh[i];
|
||||
jnum = list->numneigh[i];
|
||||
|
|
@ -109,19 +116,16 @@ void MLIAPDescriptorSNAP::forward(int* map, NeighList* list, double **descriptor
|
|||
delz = x[j][2] - ztmp;
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
int jtype = type[j];
|
||||
int jelem = map[jtype];
|
||||
const int jelem = map[jtype];
|
||||
|
||||
// printf("i = %d j = %d itype = %d jtype = %d cutsq[i][j] = %g rsq = %g\n",i,j,itype,jtype,cutsq[itype][jtype],rsq);
|
||||
|
||||
double rcutsqtmp = get_cutoff(ielem, jelem);
|
||||
if (rsq < rcutsqtmp*rcutsqtmp) {
|
||||
if (rsq < cutsq[ielem][jelem]) {
|
||||
snaptr->rij[ninside][0] = delx;
|
||||
snaptr->rij[ninside][1] = dely;
|
||||
snaptr->rij[ninside][2] = delz;
|
||||
snaptr->inside[ninside] = j;
|
||||
snaptr->wj[ninside] = wjelem[jelem];
|
||||
snaptr->rcutij[ninside] = (radi + radelem[jelem])*rcutfac;
|
||||
snaptr->element[ninside] = jelem; // element index for chem snap
|
||||
snaptr->wj[ninside] = wjelem[jelem];
|
||||
snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]);
|
||||
snaptr->element[ninside] = jelem; // element index for chem snap
|
||||
ninside++;
|
||||
}
|
||||
}
|
||||
|
|
@ -131,6 +135,7 @@ void MLIAPDescriptorSNAP::forward(int* map, NeighList* list, double **descriptor
|
|||
else
|
||||
snaptr->compute_ui(ninside, 0);
|
||||
snaptr->compute_zi();
|
||||
|
||||
if (chemflag)
|
||||
snaptr->compute_bi(ielem);
|
||||
else
|
||||
|
|
@ -146,7 +151,7 @@ void MLIAPDescriptorSNAP::forward(int* map, NeighList* list, double **descriptor
|
|||
compute forces for each atom
|
||||
---------------------------------------------------------------------- */
|
||||
|
||||
void MLIAPDescriptorSNAP::backward(PairMLIAP* pairmliap, NeighList* list, double **beta, int vflag)
|
||||
void MLIAPDescriptorSNAP::compute_forces(PairMLIAP* pairmliap, NeighList* list, double **beta, int vflag)
|
||||
{
|
||||
int i,j,jnum,ninside;
|
||||
double delx,dely,delz,rsq;
|
||||
|
|
@ -168,7 +173,6 @@ void MLIAPDescriptorSNAP::backward(PairMLIAP* pairmliap, NeighList* list, double
|
|||
const double ztmp = x[i][2];
|
||||
const int itype = type[i];
|
||||
const int ielem = pairmliap->map[itype];
|
||||
const double radi = radelem[ielem];
|
||||
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
|
|
@ -194,13 +198,13 @@ void MLIAPDescriptorSNAP::backward(PairMLIAP* pairmliap, NeighList* list, double
|
|||
int jtype = type[j];
|
||||
int jelem = pairmliap->map[jtype];
|
||||
|
||||
if (rsq < pairmliap->cutsq[itype][jtype]&&rsq>1e-20) {
|
||||
if (rsq < cutsq[ielem][jelem]) {
|
||||
snaptr->rij[ninside][0] = delx;
|
||||
snaptr->rij[ninside][1] = dely;
|
||||
snaptr->rij[ninside][2] = delz;
|
||||
snaptr->inside[ninside] = j;
|
||||
snaptr->wj[ninside] = wjelem[jelem];
|
||||
snaptr->rcutij[ninside] = (radi + radelem[jelem])*rcutfac;
|
||||
snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]);
|
||||
snaptr->element[ninside] = jelem; // element index for chem snap
|
||||
ninside++;
|
||||
}
|
||||
|
|
@ -218,9 +222,6 @@ void MLIAPDescriptorSNAP::backward(PairMLIAP* pairmliap, NeighList* list, double
|
|||
// add to Fi, subtract from Fj
|
||||
|
||||
snaptr->compute_yi(beta[ii]);
|
||||
//for (int q=0; q<snaptr->idxu_max*2; q++){
|
||||
// fprintf(screen, "%i %f\n",q, snaptr->ylist_r[q]);
|
||||
//}
|
||||
|
||||
for (int jj = 0; jj < ninside; jj++) {
|
||||
int j = snaptr->inside[jj];
|
||||
|
|
@ -240,7 +241,7 @@ void MLIAPDescriptorSNAP::backward(PairMLIAP* pairmliap, NeighList* list, double
|
|||
f[j][1] -= fij[1];
|
||||
f[j][2] -= fij[2];
|
||||
|
||||
// add in gloabl and per-atom virial contributions
|
||||
// add in global and per-atom virial contributions
|
||||
// this is optional and has no effect on force calculation
|
||||
|
||||
if (vflag)
|
||||
|
|
@ -254,21 +255,228 @@ void MLIAPDescriptorSNAP::backward(PairMLIAP* pairmliap, NeighList* list, double
|
|||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
compute force gradient for each atom
|
||||
---------------------------------------------------------------------- */
|
||||
|
||||
void MLIAPDescriptorSNAP::compute_gradients(int *map, NeighList* list,
|
||||
int gamma_nnz, int **gamma_row_index,
|
||||
int **gamma_col_index, double **gamma, double **gradforce,
|
||||
int yoffset, int zoffset)
|
||||
{
|
||||
int i,j,jnum,ninside;
|
||||
double delx,dely,delz,evdwl,rsq;
|
||||
double fij[3];
|
||||
int *jlist,*numneigh,**firstneigh;
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
int *type = atom->type;
|
||||
int nlocal = atom->nlocal;
|
||||
int newton_pair = force->newton_pair;
|
||||
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
for (int ii = 0; ii < list->inum; ii++) {
|
||||
i = list->ilist[ii];
|
||||
|
||||
const double xtmp = x[i][0];
|
||||
const double ytmp = x[i][1];
|
||||
const double ztmp = x[i][2];
|
||||
const int itype = type[i];
|
||||
const int ielem = map[itype];
|
||||
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
|
||||
// insure rij, inside, wj, and rcutij are of size jnum
|
||||
|
||||
snaptr->grow_rij(jnum);
|
||||
|
||||
// rij[][3] = displacements between atom I and those neighbors
|
||||
// inside = indices of neighbors of I within cutoff
|
||||
// wj = weights for neighbors of I within cutoff
|
||||
// rcutij = cutoffs for neighbors of I within cutoff
|
||||
// note Rij sign convention => dU/dRij = dU/dRj = -dU/dRi
|
||||
|
||||
ninside = 0;
|
||||
for (int jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
j &= NEIGHMASK;
|
||||
delx = x[j][0] - xtmp;
|
||||
dely = x[j][1] - ytmp;
|
||||
delz = x[j][2] - ztmp;
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
int jtype = type[j];
|
||||
const int jelem = map[jtype];
|
||||
|
||||
if (rsq < cutsq[ielem][jelem]) {
|
||||
snaptr->rij[ninside][0] = delx;
|
||||
snaptr->rij[ninside][1] = dely;
|
||||
snaptr->rij[ninside][2] = delz;
|
||||
snaptr->inside[ninside] = j;
|
||||
snaptr->wj[ninside] = wjelem[jelem];
|
||||
snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]);
|
||||
snaptr->element[ninside] = jelem; // element index for chem snap
|
||||
ninside++;
|
||||
}
|
||||
}
|
||||
|
||||
if (chemflag)
|
||||
snaptr->compute_ui(ninside, ielem);
|
||||
else
|
||||
snaptr->compute_ui(ninside, 0);
|
||||
|
||||
snaptr->compute_zi();
|
||||
if (chemflag)
|
||||
snaptr->compute_bi(ielem);
|
||||
else
|
||||
snaptr->compute_bi(0);
|
||||
|
||||
for (int jj = 0; jj < ninside; jj++) {
|
||||
const int j = snaptr->inside[jj];
|
||||
|
||||
if(chemflag)
|
||||
snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj],
|
||||
snaptr->rcutij[jj],jj, snaptr->element[jj]);
|
||||
else
|
||||
snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj],
|
||||
snaptr->rcutij[jj],jj, 0);
|
||||
|
||||
snaptr->compute_dbidrj();
|
||||
|
||||
// Accumulate gamma_lk*dB_k/dRi, -gamma_lk**dB_k/dRj
|
||||
|
||||
for (int inz = 0; inz < gamma_nnz; inz++) {
|
||||
const int l = gamma_row_index[ii][inz];
|
||||
const int k = gamma_col_index[ii][inz];
|
||||
gradforce[i][l] += gamma[ii][inz]*snaptr->dblist[k][0];
|
||||
gradforce[i][l+yoffset] += gamma[ii][inz]*snaptr->dblist[k][1];
|
||||
gradforce[i][l+zoffset] += gamma[ii][inz]*snaptr->dblist[k][2];
|
||||
gradforce[j][l] -= gamma[ii][inz]*snaptr->dblist[k][0];
|
||||
gradforce[j][l+yoffset] -= gamma[ii][inz]*snaptr->dblist[k][1];
|
||||
gradforce[j][l+zoffset] -= gamma[ii][inz]*snaptr->dblist[k][2];
|
||||
}
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
compute descriptor gradients for each neighbor atom
|
||||
---------------------------------------------------------------------- */
|
||||
|
||||
void MLIAPDescriptorSNAP::compute_descriptor_gradients(int *map, NeighList* list,
|
||||
int gamma_nnz, int **gamma_row_index,
|
||||
int **gamma_col_index, double **gamma, double **graddesc,
|
||||
int yoffset, int zoffset)
|
||||
{
|
||||
int i,j,jnum,ninside;
|
||||
double delx,dely,delz,evdwl,rsq;
|
||||
double fij[3];
|
||||
int *jlist,*numneigh,**firstneigh;
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
int *type = atom->type;
|
||||
int nlocal = atom->nlocal;
|
||||
int newton_pair = force->newton_pair;
|
||||
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
for (int ii = 0; ii < list->inum; ii++) {
|
||||
i = list->ilist[ii];
|
||||
|
||||
const double xtmp = x[i][0];
|
||||
const double ytmp = x[i][1];
|
||||
const double ztmp = x[i][2];
|
||||
const int itype = type[i];
|
||||
const int ielem = map[itype];
|
||||
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
|
||||
// insure rij, inside, wj, and rcutij are of size jnum
|
||||
|
||||
snaptr->grow_rij(jnum);
|
||||
|
||||
// rij[][3] = displacements between atom I and those neighbors
|
||||
// inside = indices of neighbors of I within cutoff
|
||||
// wj = weights for neighbors of I within cutoff
|
||||
// rcutij = cutoffs for neighbors of I within cutoff
|
||||
// note Rij sign convention => dU/dRij = dU/dRj = -dU/dRi
|
||||
|
||||
ninside = 0;
|
||||
for (int jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
j &= NEIGHMASK;
|
||||
delx = x[j][0] - xtmp;
|
||||
dely = x[j][1] - ytmp;
|
||||
delz = x[j][2] - ztmp;
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
int jtype = type[j];
|
||||
const int jelem = map[jtype];
|
||||
|
||||
if (rsq < cutsq[ielem][jelem]) {
|
||||
snaptr->rij[ninside][0] = delx;
|
||||
snaptr->rij[ninside][1] = dely;
|
||||
snaptr->rij[ninside][2] = delz;
|
||||
snaptr->inside[ninside] = j;
|
||||
snaptr->wj[ninside] = wjelem[jelem];
|
||||
snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]);
|
||||
snaptr->element[ninside] = jelem; // element index for chem snap
|
||||
ninside++;
|
||||
}
|
||||
}
|
||||
|
||||
if (chemflag)
|
||||
snaptr->compute_ui(ninside, ielem);
|
||||
else
|
||||
snaptr->compute_ui(ninside, 0);
|
||||
|
||||
snaptr->compute_zi();
|
||||
if (chemflag)
|
||||
snaptr->compute_bi(ielem);
|
||||
else
|
||||
snaptr->compute_bi(0);
|
||||
|
||||
for (int jj = 0; jj < ninside; jj++) {
|
||||
const int j = snaptr->inside[jj];
|
||||
|
||||
if(chemflag)
|
||||
snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj],
|
||||
snaptr->rcutij[jj],jj, snaptr->element[jj]);
|
||||
else
|
||||
snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj],
|
||||
snaptr->rcutij[jj],jj, 0);
|
||||
|
||||
snaptr->compute_dbidrj();
|
||||
|
||||
// Accumulate dB_k^i/dRi, dB_k^i/dRj
|
||||
|
||||
for (int k = 0; k < ndescriptors; k++) {
|
||||
graddesc[i][k] = snaptr->dblist[k][0];
|
||||
graddesc[i][k] = snaptr->dblist[k][1];
|
||||
graddesc[i][k] = snaptr->dblist[k][2];
|
||||
graddesc[j][k] = -snaptr->dblist[k][0];
|
||||
graddesc[j][k] = -snaptr->dblist[k][1];
|
||||
graddesc[j][k] = -snaptr->dblist[k][2];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
set coeffs for one or more type pairs
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void MLIAPDescriptorSNAP::init()
|
||||
{
|
||||
|
||||
snaptr = new SNA(lmp, rfac0, twojmax,
|
||||
rmin0, switchflag, bzeroflag,
|
||||
chemflag, bnormflag, wselfallflag, nelements);
|
||||
|
||||
snaptr->init();
|
||||
|
||||
ndescriptors = snaptr->ncoeff;
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
|
@ -413,31 +621,20 @@ void MLIAPDescriptorSNAP::read_paramfile(char *paramfilename)
|
|||
!elementsflag || !radelemflag || !wjelemflag)
|
||||
error->all(FLERR,"Incorrect SNAP parameter file");
|
||||
|
||||
}
|
||||
// construct cutsq
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
provide cutoff distance for two elements
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double MLIAPDescriptorSNAP::get_cutoff(int ielem, int jelem)
|
||||
{
|
||||
return (radelem[ielem] + radelem[jelem])*rcutfac;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
calculate maximum cutoff distance
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double MLIAPDescriptorSNAP::get_cutmax()
|
||||
{
|
||||
double cut;
|
||||
double cutmax = 0.0;
|
||||
for(int ielem = 0; ielem <= nelements; ielem++) {
|
||||
cutmax = 0.0;
|
||||
memory->create(cutsq,nelements,nelements,"mliap/descriptor/snap:cutsq");
|
||||
for (int ielem = 0; ielem < nelements; ielem++) {
|
||||
cut = 2.0*radelem[ielem]*rcutfac;
|
||||
if (cut > cutmax) cutmax = cut;
|
||||
return cutmax;
|
||||
cutsq[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;
|
||||
}
|
||||
}
|
||||
return cutmax;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
|
|
|
|||
|
|
@ -22,11 +22,13 @@ class MLIAPDescriptorSNAP : public MLIAPDescriptor {
|
|||
public:
|
||||
MLIAPDescriptorSNAP(LAMMPS*, char*);
|
||||
~MLIAPDescriptorSNAP();
|
||||
virtual void forward(int*, class NeighList*, double**);
|
||||
virtual void backward(class PairMLIAP*, class NeighList*, double**, int);
|
||||
virtual void compute_descriptors(int*, class NeighList*, double**);
|
||||
virtual void compute_forces(class PairMLIAP*, class NeighList*, double**, int);
|
||||
virtual void compute_gradients(int*, class NeighList*, int, int**, int**, double**,
|
||||
double**, int, int);
|
||||
virtual void compute_descriptor_gradients(int*, class NeighList*, int, int**, int**, double**,
|
||||
double**, int, int);
|
||||
virtual void init();
|
||||
virtual double get_cutoff(int, int);
|
||||
virtual double get_cutmax();
|
||||
virtual double memory_usage();
|
||||
|
||||
double rcutfac; // declared public to workaround gcc 4.9
|
||||
|
|
|
|||
|
|
@ -41,6 +41,16 @@ MLIAPModel::MLIAPModel(LAMMPS* lmp, char* coefffilename) : Pointers(lmp)
|
|||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
MLIAPModel::MLIAPModel(LAMMPS* lmp, int nelements_in, int nparams_in) : Pointers(lmp)
|
||||
{
|
||||
nelements = nelements_in;
|
||||
nparams = nparams_in;
|
||||
coeffelem = NULL;
|
||||
nonlinearflag = 0;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
MLIAPModel::~MLIAPModel()
|
||||
{
|
||||
memory->destroy(coeffelem);
|
||||
|
|
|
|||
|
|
@ -21,8 +21,11 @@ namespace LAMMPS_NS {
|
|||
class MLIAPModel : protected Pointers {
|
||||
public:
|
||||
MLIAPModel(LAMMPS*, char*);
|
||||
MLIAPModel(LAMMPS*, int, int);
|
||||
~MLIAPModel();
|
||||
virtual void gradient(class PairMLIAP*, class NeighList*, double**, double**, int)=0;
|
||||
virtual void param_gradient(int*, class NeighList*, double**, int**, int**, double**, double*)=0;
|
||||
virtual int get_gamma_nnz()=0;
|
||||
virtual void init();
|
||||
virtual double memory_usage();
|
||||
int nelements; // # of unique elements
|
||||
|
|
|
|||
|
|
@ -31,7 +31,14 @@ using namespace LAMMPS_NS;
|
|||
MLIAPModelLinear::MLIAPModelLinear(LAMMPS* lmp, char* coefffilename) :
|
||||
MLIAPModel(lmp, coefffilename)
|
||||
{
|
||||
nonlinearflag = 0;
|
||||
ndescriptors = nparams - 1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
MLIAPModelLinear::MLIAPModelLinear(LAMMPS* lmp, int nelements_in, int nparams_in) :
|
||||
MLIAPModel(lmp, nelements_in, nparams_in)
|
||||
{
|
||||
ndescriptors = nparams - 1;
|
||||
}
|
||||
|
||||
|
|
@ -40,7 +47,8 @@ MLIAPModelLinear::MLIAPModelLinear(LAMMPS* lmp, char* coefffilename) :
|
|||
MLIAPModelLinear::~MLIAPModelLinear(){}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Calculate model gradients w.r.t descriptors for each atom dE(B_i)/dB_i
|
||||
Calculate model gradients w.r.t descriptors
|
||||
for each atom beta_i = dE(B_i)/dB_i
|
||||
---------------------------------------------------------------------- */
|
||||
|
||||
void MLIAPModelLinear::gradient(PairMLIAP* pairmliap, NeighList* list, double **descriptors, double **beta, int eflag)
|
||||
|
|
@ -77,3 +85,64 @@ void MLIAPModelLinear::gradient(PairMLIAP* pairmliap, NeighList* list, double **
|
|||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
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 MLIAPModelLinear::param_gradient(int *map, NeighList* list,
|
||||
double **descriptors,
|
||||
int **gamma_row_index, int **gamma_col_index,
|
||||
double **gamma, double *egradient)
|
||||
{
|
||||
int i;
|
||||
int *type = atom->type;
|
||||
|
||||
// zero out energy gradients
|
||||
|
||||
for (int l = 0; l < nelements*nparams; l++)
|
||||
egradient[l] = 0.0;
|
||||
|
||||
for (int ii = 0; ii < list->inum; ii++) {
|
||||
|
||||
i = list->ilist[ii];
|
||||
const int itype = type[i];
|
||||
const int ielem = map[itype];
|
||||
const int elemoffset = nparams*ielem;
|
||||
|
||||
int l = elemoffset+1;
|
||||
for (int icoeff = 0; icoeff < ndescriptors; icoeff++) {
|
||||
gamma[ii][icoeff] = 1.0;
|
||||
gamma_row_index[ii][icoeff] = l++;
|
||||
gamma_col_index[ii][icoeff] = icoeff;
|
||||
}
|
||||
|
||||
// gradient of energy of atom I w.r.t. parameters
|
||||
|
||||
l = elemoffset;
|
||||
egradient[l++] += 1.0;
|
||||
for (int icoeff = 0; icoeff < ndescriptors; icoeff++)
|
||||
egradient[l++] += descriptors[ii][icoeff];
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
/* ----------------------------------------------------------------------
|
||||
count the number of non-zero entries in gamma matrix
|
||||
---------------------------------------------------------------------- */
|
||||
|
||||
int MLIAPModelLinear::get_gamma_nnz()
|
||||
{
|
||||
int inz = ndescriptors;
|
||||
return inz;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -21,8 +21,11 @@ namespace LAMMPS_NS {
|
|||
class MLIAPModelLinear : public MLIAPModel {
|
||||
public:
|
||||
MLIAPModelLinear(LAMMPS*, char*);
|
||||
MLIAPModelLinear(LAMMPS*, int, int);
|
||||
~MLIAPModelLinear();
|
||||
virtual void gradient(class PairMLIAP*, class NeighList*, double**, double**, int);
|
||||
virtual void param_gradient(int*, class NeighList*, double**, int**, int**, double**, double*);
|
||||
virtual int get_gamma_nnz();
|
||||
|
||||
protected:
|
||||
};
|
||||
|
|
|
|||
|
|
@ -37,6 +37,15 @@ MLIAPModelQuadratic::MLIAPModelQuadratic(LAMMPS* lmp, char* coefffilename) :
|
|||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
MLIAPModelQuadratic::MLIAPModelQuadratic(LAMMPS* lmp, int nelements_in, int nparams_in) :
|
||||
MLIAPModel(lmp, nelements_in, nparams_in)
|
||||
{
|
||||
nonlinearflag = 1;
|
||||
ndescriptors = sqrt(2*nparams)-1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
MLIAPModelQuadratic::~MLIAPModelQuadratic(){}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
|
|
@ -101,3 +110,110 @@ void MLIAPModelQuadratic::gradient(PairMLIAP* pairmliap, NeighList* list, double
|
|||
}
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
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 MLIAPModelQuadratic::param_gradient(int *map, NeighList* list,
|
||||
double **descriptors,
|
||||
int **gamma_row_index, int **gamma_col_index,
|
||||
double **gamma, double *egradient)
|
||||
{
|
||||
int i;
|
||||
int *type = atom->type;
|
||||
|
||||
// zero out energy gradients
|
||||
|
||||
for (int l = 0; l < nelements*nparams; l++)
|
||||
egradient[l] = 0.0;
|
||||
|
||||
for (int ii = 0; ii < list->inum; ii++) {
|
||||
|
||||
i = list->ilist[ii];
|
||||
const int itype = type[i];
|
||||
const int ielem = map[itype];
|
||||
const int elemoffset = nparams*ielem;
|
||||
|
||||
// linear contributions
|
||||
|
||||
int l = elemoffset+1;
|
||||
for (int icoeff = 0; icoeff < ndescriptors; icoeff++) {
|
||||
gamma[ii][icoeff] = 1.0;
|
||||
gamma_row_index[ii][icoeff] = l++;
|
||||
gamma_col_index[ii][icoeff] = icoeff;
|
||||
}
|
||||
|
||||
// quadratic contributions
|
||||
|
||||
int inz = ndescriptors;
|
||||
for (int icoeff = 0; icoeff < ndescriptors; icoeff++) {
|
||||
double bveci = descriptors[ii][icoeff];
|
||||
gamma[ii][inz] = bveci;
|
||||
gamma_row_index[ii][inz] = l++;
|
||||
gamma_col_index[ii][inz] = icoeff;
|
||||
inz++;
|
||||
for (int jcoeff = icoeff+1; jcoeff < ndescriptors; jcoeff++) {
|
||||
double bvecj = descriptors[ii][jcoeff];
|
||||
gamma[ii][inz] = bvecj; // derivative w.r.t. B[icoeff]
|
||||
gamma_row_index[ii][inz] = l;
|
||||
gamma_col_index[ii][inz] = icoeff;
|
||||
inz++;
|
||||
gamma[ii][inz] = bveci; // derivative w.r.t. B[jcoeff]
|
||||
gamma_row_index[ii][inz] = l;
|
||||
gamma_col_index[ii][inz] = jcoeff;
|
||||
inz++;
|
||||
l++;
|
||||
}
|
||||
}
|
||||
|
||||
// gradient of energy of atom I w.r.t. parameters
|
||||
|
||||
l = elemoffset;
|
||||
egradient[l++] += 1.0;
|
||||
for (int icoeff = 0; icoeff < ndescriptors; icoeff++)
|
||||
egradient[l++] += descriptors[ii][icoeff];
|
||||
|
||||
// quadratic contributions
|
||||
|
||||
for (int icoeff = 0; icoeff < ndescriptors; icoeff++) {
|
||||
double bveci = descriptors[ii][icoeff];
|
||||
egradient[l++] += 0.5*bveci*bveci;
|
||||
for (int jcoeff = icoeff+1; jcoeff < ndescriptors; jcoeff++) {
|
||||
double bvecj = descriptors[ii][jcoeff];
|
||||
egradient[l++] += bveci*bvecj;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
count the number of non-zero entries in gamma matrix
|
||||
---------------------------------------------------------------------- */
|
||||
|
||||
int MLIAPModelQuadratic::get_gamma_nnz()
|
||||
{
|
||||
int inz = ndescriptors;
|
||||
for (int icoeff = 0; icoeff < ndescriptors; icoeff++) {
|
||||
inz++;
|
||||
for (int jcoeff = icoeff+1; jcoeff < ndescriptors; jcoeff++) {
|
||||
inz++;
|
||||
inz++;
|
||||
}
|
||||
}
|
||||
|
||||
return inz;
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -21,8 +21,11 @@ namespace LAMMPS_NS {
|
|||
class MLIAPModelQuadratic : public MLIAPModel {
|
||||
public:
|
||||
MLIAPModelQuadratic(LAMMPS*, char*);
|
||||
MLIAPModelQuadratic(LAMMPS*, int, int);
|
||||
~MLIAPModelQuadratic();
|
||||
virtual void gradient(class PairMLIAP*, class NeighList*, double**, double**, int);
|
||||
virtual void param_gradient(int*, class NeighList*, double**, int**, int**, double**, double*);
|
||||
virtual int get_gamma_nnz();
|
||||
|
||||
protected:
|
||||
};
|
||||
|
|
|
|||
|
|
@ -86,7 +86,7 @@ void PairMLIAP::compute(int eflag, int vflag)
|
|||
// compute descriptors, if needed
|
||||
|
||||
if (model->nonlinearflag || eflag)
|
||||
descriptor->forward(map, list, descriptors);
|
||||
descriptor->compute_descriptors(map, list, descriptors);
|
||||
|
||||
// compute E_i and beta_i = dE_i/dB_i for all i in list
|
||||
|
||||
|
|
@ -94,7 +94,7 @@ void PairMLIAP::compute(int eflag, int vflag)
|
|||
|
||||
// calculate force contributions beta_i*dB_i/dR_j
|
||||
|
||||
descriptor->backward(this, list, beta, vflag);
|
||||
descriptor->compute_forces(this, list, beta, vflag);
|
||||
|
||||
// calculate stress
|
||||
|
||||
|
|
@ -310,7 +310,7 @@ 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 descriptor->get_cutoff(map[i],map[j]);
|
||||
return sqrt(descriptor->cutsq[map[i]][map[j]]);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
|
|
@ -323,7 +323,6 @@ double PairMLIAP::memory_usage()
|
|||
|
||||
int n = atom->ntypes+1;
|
||||
bytes += n*n*sizeof(int); // setflag
|
||||
bytes += n*n*sizeof(double); // cutsq
|
||||
bytes += beta_max*ndescriptors*sizeof(double); // descriptors
|
||||
bytes += beta_max*ndescriptors*sizeof(double); // beta
|
||||
|
||||
|
|
|
|||
|
|
@ -210,32 +210,18 @@ void ComputeSnap::init()
|
|||
"snap:snapall");
|
||||
array = snapall;
|
||||
|
||||
// INCOMPLETE: modify->find_compute()
|
||||
// was called 223960 times by snappy Ta example
|
||||
// that is over 600 times per config?
|
||||
// how is this possible???
|
||||
|
||||
// find compute for reference energy
|
||||
|
||||
char *id_pe = (char *) "thermo_pe";
|
||||
int ipe = modify->find_compute(id_pe);
|
||||
int ipe = modify->find_compute("thermo_pe");
|
||||
if (ipe == -1)
|
||||
error->all(FLERR,"compute thermo_pe does not exist.");
|
||||
c_pe = modify->compute[ipe];
|
||||
|
||||
// add compute for reference virial tensor
|
||||
|
||||
char *id_virial = (char *) "snap_press";
|
||||
char **newarg = new char*[5];
|
||||
newarg[0] = id_virial;
|
||||
newarg[1] = (char *) "all";
|
||||
newarg[2] = (char *) "pressure";
|
||||
newarg[3] = (char *) "NULL";
|
||||
newarg[4] = (char *) "virial";
|
||||
modify->add_compute(5,newarg);
|
||||
delete [] newarg;
|
||||
modify->add_compute("snap_press all pressure NULL virial");
|
||||
|
||||
int ivirial = modify->find_compute(id_virial);
|
||||
int ivirial = modify->find_compute("snap_press");
|
||||
if (ivirial == -1)
|
||||
error->all(FLERR,"compute snap_press does not exist.");
|
||||
c_virial = modify->compute[ivirial];
|
||||
|
|
|
|||
|
|
@ -316,6 +316,7 @@ void Thermo::header()
|
|||
|
||||
std::string hdr;
|
||||
for (int i = 0; i < nfield; i++) hdr += keyword[i] + std::string(" ");
|
||||
|
||||
if (me == 0) utils::logmesg(lmp,hdr+"\n");
|
||||
}
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue