diff --git a/doc/src/Commands_compute.rst b/doc/src/Commands_compute.rst index a0a84f29f4..6841485f8f 100644 --- a/doc/src/Commands_compute.rst +++ b/doc/src/Commands_compute.rst @@ -80,6 +80,7 @@ KOKKOS, o = USER-OMP, t = OPT. * :doc:`ke/eff ` * :doc:`ke/rigid ` * :doc:`mesont ` + * :doc:`mliap ` * :doc:`momentum ` * :doc:`msd ` * :doc:`msd/chunk ` diff --git a/doc/src/compute.rst b/doc/src/compute.rst index da2c8d28eb..e31f4fde96 100644 --- a/doc/src/compute.rst +++ b/doc/src/compute.rst @@ -225,6 +225,7 @@ The individual style names on the :doc:`Commands compute ` doc * :doc:`ke/atom/eff ` - per-atom translational and radial kinetic energy in the electron force field model * :doc:`ke/eff ` - kinetic energy of a group of nuclei and electrons in the electron force field model * :doc:`ke/rigid ` - translational kinetic energy of rigid bodies +* :doc:`mliap ` - gradients of energy and forces w.r.t. model parameters and related quantities for training machine learning interatomic potentials * :doc:`momentum ` - translational momentum * :doc:`msd ` - mean-squared displacement of group of atoms * :doc:`msd/chunk ` - mean-squared displacement for each chunk @@ -274,7 +275,7 @@ The individual style names on the :doc:`Commands compute ` doc * :doc:`smd/ulsph/strain/rate ` - * :doc:`smd/ulsph/stress ` - per-particle Cauchy stress tensor and von Mises equivalent stress in Smooth Mach Dynamics * :doc:`smd/vol ` - per-particle volumes and their sum in Smooth Mach Dynamics -* :doc:`snap ` - bispectrum components and related quantities for a group of atoms +* :doc:`snap ` - gradients of SNAP energy and forces w.r.t. linear coefficients and related quantities for fitting SNAP potentials * :doc:`sna/atom ` - bispectrum components for each atom * :doc:`snad/atom ` - derivative of bispectrum components for each atom * :doc:`snav/atom ` - virial contribution from bispectrum components for each atom diff --git a/doc/src/compute_mliap.rst b/doc/src/compute_mliap.rst new file mode 100644 index 0000000000..91253b29c2 --- /dev/null +++ b/doc/src/compute_mliap.rst @@ -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 ` 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 ` 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 `, 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 ` 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 `. +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 `, 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 ` command can remove pairwise + interactions between atoms in the same bond, angle, or dihedral. This + is the default setting for the :doc:`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 ` + 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 ` 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 ` 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 `_. + +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 ` doc page for more info. + +Related commands +"""""""""""""""" + +:doc:`pair_style mliap ` + +**Default:** none diff --git a/doc/src/compute_sna_atom.rst b/doc/src/compute_sna_atom.rst index 4e42c4523b..cf9339a30d 100644 --- a/doc/src/compute_sna_atom.rst +++ b/doc/src/compute_sna_atom.rst @@ -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 ` 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 `_. Restrictions """""""""""" diff --git a/doc/src/pair_mliap.rst b/doc/src/pair_mliap.rst index ed7caea86f..d96e35a7f4 100644 --- a/doc/src/pair_mliap.rst +++ b/doc/src/pair_mliap.rst @@ -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 ` 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 ` coefficient file. Specifically, the line containing the element weight and radius is omitted, @@ -131,6 +146,6 @@ See the :doc:`Build package ` doc page for more info. Related commands """""""""""""""" -:doc:`pair_style snap `, +:doc:`pair_style snap `, :doc:`compute mliap ` **Default:** none diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index 726018a8bf..682969714c 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -2083,6 +2083,7 @@ Novint np Npair Npairs +nparams nparticle npernode nph diff --git a/examples/mliap/compute.mliap.descriptor b/examples/mliap/compute.mliap.descriptor new file mode 100644 index 0000000000..c3fc432886 --- /dev/null +++ b/examples/mliap/compute.mliap.descriptor @@ -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 diff --git a/examples/mliap/in.mliap.quadratic.compute b/examples/mliap/in.mliap.quadratic.compute new file mode 100644 index 0000000000..6b346392fa --- /dev/null +++ b/examples/mliap/in.mliap.quadratic.compute @@ -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} diff --git a/examples/mliap/in.mliap.snap.compute b/examples/mliap/in.mliap.snap.compute new file mode 100644 index 0000000000..4b88a3e52d --- /dev/null +++ b/examples/mliap/in.mliap.snap.compute @@ -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} diff --git a/examples/mliap/log.03Jul20.mliap.quadratic.compute.g++.1 b/examples/mliap/log.03Jul20.mliap.quadratic.compute.g++.1 new file mode 100644 index 0000000000..d5226cec89 --- /dev/null +++ b/examples/mliap/log.03Jul20.mliap.quadratic.compute.g++.1 @@ -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 diff --git a/examples/mliap/log.03Jul20.mliap.quadratic.compute.g++.4 b/examples/mliap/log.03Jul20.mliap.quadratic.compute.g++.4 new file mode 100644 index 0000000000..dfa6483b53 --- /dev/null +++ b/examples/mliap/log.03Jul20.mliap.quadratic.compute.g++.4 @@ -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 diff --git a/examples/mliap/log.03Jul20.mliap.snap.compute.g++.1 b/examples/mliap/log.03Jul20.mliap.snap.compute.g++.1 new file mode 100644 index 0000000000..083d451afa --- /dev/null +++ b/examples/mliap/log.03Jul20.mliap.snap.compute.g++.1 @@ -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 diff --git a/examples/mliap/log.03Jul20.mliap.snap.compute.g++.4 b/examples/mliap/log.03Jul20.mliap.snap.compute.g++.4 new file mode 100644 index 0000000000..ebedaee049 --- /dev/null +++ b/examples/mliap/log.03Jul20.mliap.snap.compute.g++.4 @@ -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 diff --git a/examples/snap/in.snap.compute b/examples/snap/in.snap.compute index e698a134a6..b0c7314882 100644 --- a/examples/snap/in.snap.compute +++ b/examples/snap/in.snap.compute @@ -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[*] diff --git a/examples/snap/in.snap.compute.quadratic b/examples/snap/in.snap.compute.quadratic index 00e46bd3a8..e03d4af3bf 100644 --- a/examples/snap/in.snap.compute.quadratic +++ b/examples/snap/in.snap.compute.quadratic @@ -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[*] diff --git a/examples/snap/log.27Nov18.snap.compute.g++.1 b/examples/snap/log.03Jul20.snap.compute.g++.1 similarity index 87% rename from examples/snap/log.27Nov18.snap.compute.g++.1 rename to examples/snap/log.03Jul20.snap.compute.g++.1 index b189c552da..8b9fe80a07 100644 --- a/examples/snap/log.27Nov18.snap.compute.g++.1 +++ b/examples/snap/log.03Jul20.snap.compute.g++.1 @@ -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 diff --git a/examples/snap/log.27Nov18.snap.compute.g++.4 b/examples/snap/log.03Jul20.snap.compute.g++.4 similarity index 86% rename from examples/snap/log.27Nov18.snap.compute.g++.4 rename to examples/snap/log.03Jul20.snap.compute.g++.4 index f979123b2a..ff68bf75d5 100644 --- a/examples/snap/log.27Nov18.snap.compute.g++.4 +++ b/examples/snap/log.03Jul20.snap.compute.g++.4 @@ -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 diff --git a/examples/snap/log.27Nov18.snap.compute.quadratic.g++.1 b/examples/snap/log.03Jul20.snap.compute.quadratic.g++.1 similarity index 87% rename from examples/snap/log.27Nov18.snap.compute.quadratic.g++.1 rename to examples/snap/log.03Jul20.snap.compute.quadratic.g++.1 index c2054023ed..a52810063c 100644 --- a/examples/snap/log.27Nov18.snap.compute.quadratic.g++.1 +++ b/examples/snap/log.03Jul20.snap.compute.quadratic.g++.1 @@ -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 diff --git a/examples/snap/log.27Nov18.snap.compute.quadratic.g++.4 b/examples/snap/log.03Jul20.snap.compute.quadratic.g++.4 similarity index 86% rename from examples/snap/log.27Nov18.snap.compute.quadratic.g++.4 rename to examples/snap/log.03Jul20.snap.compute.quadratic.g++.4 index 819397f745..ba674cd4fa 100644 --- a/examples/snap/log.27Nov18.snap.compute.quadratic.g++.4 +++ b/examples/snap/log.03Jul20.snap.compute.quadratic.g++.4 @@ -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 diff --git a/src/MLIAP/README b/src/MLIAP/README new file mode 100644 index 0000000000..77e38d3c62 --- /dev/null +++ b/src/MLIAP/README @@ -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>. diff --git a/src/MLIAP/compute_mliap.cpp b/src/MLIAP/compute_mliap.cpp new file mode 100644 index 0000000000..655e70c6f9 --- /dev/null +++ b/src/MLIAP/compute_mliap.cpp @@ -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 +#include +#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; +} diff --git a/src/MLIAP/compute_mliap.h b/src/MLIAP/compute_mliap.h new file mode 100644 index 0000000000..74b5fc1f93 --- /dev/null +++ b/src/MLIAP/compute_mliap.h @@ -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. + +*/ diff --git a/src/MLIAP/mliap_descriptor.h b/src/MLIAP/mliap_descriptor.h index cd42cb3be6..072a5b4729 100644 --- a/src/MLIAP/mliap_descriptor.h +++ b/src/MLIAP/mliap_descriptor.h @@ -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: }; diff --git a/src/MLIAP/mliap_descriptor_snap.cpp b/src/MLIAP/mliap_descriptor_snap.cpp index 5568ad5b4f..c2879e942a 100644 --- a/src/MLIAP/mliap_descriptor_snap.cpp +++ b/src/MLIAP/mliap_descriptor_snap.cpp @@ -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; qidxu_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; } /* ---------------------------------------------------------------------- diff --git a/src/MLIAP/mliap_descriptor_snap.h b/src/MLIAP/mliap_descriptor_snap.h index 15691fabfe..a929e01149 100644 --- a/src/MLIAP/mliap_descriptor_snap.h +++ b/src/MLIAP/mliap_descriptor_snap.h @@ -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 diff --git a/src/MLIAP/mliap_model.cpp b/src/MLIAP/mliap_model.cpp index 2facd65ee6..1f88465304 100644 --- a/src/MLIAP/mliap_model.cpp +++ b/src/MLIAP/mliap_model.cpp @@ -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); diff --git a/src/MLIAP/mliap_model.h b/src/MLIAP/mliap_model.h index aeb16cb299..57f6f76214 100644 --- a/src/MLIAP/mliap_model.h +++ b/src/MLIAP/mliap_model.h @@ -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 diff --git a/src/MLIAP/mliap_model_linear.cpp b/src/MLIAP/mliap_model_linear.cpp index 36c0a054f1..b8d06155ba 100644 --- a/src/MLIAP/mliap_model_linear.cpp +++ b/src/MLIAP/mliap_model_linear.cpp @@ -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; +} diff --git a/src/MLIAP/mliap_model_linear.h b/src/MLIAP/mliap_model_linear.h index 53c7c36cef..db4cb26514 100644 --- a/src/MLIAP/mliap_model_linear.h +++ b/src/MLIAP/mliap_model_linear.h @@ -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: }; diff --git a/src/MLIAP/mliap_model_quadratic.cpp b/src/MLIAP/mliap_model_quadratic.cpp index 57b89caacf..f5aae1d89d 100644 --- a/src/MLIAP/mliap_model_quadratic.cpp +++ b/src/MLIAP/mliap_model_quadratic.cpp @@ -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; +} + diff --git a/src/MLIAP/mliap_model_quadratic.h b/src/MLIAP/mliap_model_quadratic.h index 6ca0697919..b597fc7664 100644 --- a/src/MLIAP/mliap_model_quadratic.h +++ b/src/MLIAP/mliap_model_quadratic.h @@ -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: }; diff --git a/src/MLIAP/pair_mliap.cpp b/src/MLIAP/pair_mliap.cpp index fac25ae578..2dd70fc6c2 100644 --- a/src/MLIAP/pair_mliap.cpp +++ b/src/MLIAP/pair_mliap.cpp @@ -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 diff --git a/src/SNAP/compute_snap.cpp b/src/SNAP/compute_snap.cpp index c9e1006571..7eb216a63d 100644 --- a/src/SNAP/compute_snap.cpp +++ b/src/SNAP/compute_snap.cpp @@ -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]; diff --git a/src/thermo.cpp b/src/thermo.cpp index dcec00fafe..dddfb8b21c 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -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"); }