Merge pull request #1779 from athomps/compute-snap

Compute snap
This commit is contained in:
Richard Berger 2019-12-01 12:52:22 -07:00 committed by GitHub
commit 3ce020eab2
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
9 changed files with 816 additions and 20 deletions

View File

@ -49,19 +49,19 @@ KOKKOS, o = USER-OMP, t = OPT.
+------------------------------------------------------------+--------------------------------------------------------+------------------------------------------------------------------+--------------------------------------------------------------+----------------------------------------------------------+--------------------------------------------------------------+
| :doc:`smd/tlsph/num/neighs <compute_smd_tlsph_num_neighs>` | :doc:`smd/tlsph/shape <compute_smd_tlsph_shape>` | :doc:`smd/tlsph/strain <compute_smd_tlsph_strain>` | :doc:`smd/tlsph/strain/rate <compute_smd_tlsph_strain_rate>` | :doc:`smd/tlsph/stress <compute_smd_tlsph_stress>` | :doc:`smd/triangle/vertices <compute_smd_triangle_vertices>` |
+------------------------------------------------------------+--------------------------------------------------------+------------------------------------------------------------------+--------------------------------------------------------------+----------------------------------------------------------+--------------------------------------------------------------+
| :doc:`smd/ulsph/num/neighs <compute_smd_ulsph_num_neighs>` | :doc:`smd/ulsph/strain <compute_smd_ulsph_strain>` | :doc:`smd/ulsph/strain/rate <compute_smd_ulsph_strain_rate>` | :doc:`smd/ulsph/stress <compute_smd_ulsph_stress>` | :doc:`smd/vol <compute_smd_vol>` | :doc:`sna/atom <compute_sna_atom>` |
| :doc:`smd/ulsph/num/neighs <compute_smd_ulsph_num_neighs>` | :doc:`smd/ulsph/strain <compute_smd_ulsph_strain>` | :doc:`smd/ulsph/strain/rate <compute_smd_ulsph_strain_rate>` | :doc:`smd/ulsph/stress <compute_smd_ulsph_stress>` | :doc:`smd/vol <compute_smd_vol>` | :doc:`snap <compute_sna_atom>` |
+------------------------------------------------------------+--------------------------------------------------------+------------------------------------------------------------------+--------------------------------------------------------------+----------------------------------------------------------+--------------------------------------------------------------+
| :doc:`snad/atom <compute_sna_atom>` | :doc:`snav/atom <compute_sna_atom>` | :doc:`spin <compute_spin>` | :doc:`stress/atom <compute_stress_atom>` | :doc:`stress/mop <compute_stress_mop>` | :doc:`stress/mop/profile <compute_stress_mop>` |
| :doc:`sna/atom <compute_sna_atom>` | :doc:`snad/atom <compute_sna_atom>` | :doc:`snav/atom <compute_sna_atom>` | :doc:`spin <compute_spin>` | :doc:`stress/atom <compute_stress_atom>` | :doc:`stress/mop <compute_stress_mop>` |
+------------------------------------------------------------+--------------------------------------------------------+------------------------------------------------------------------+--------------------------------------------------------------+----------------------------------------------------------+--------------------------------------------------------------+
| :doc:`stress/tally <compute_tally>` | :doc:`tdpd/cc/atom <compute_tdpd_cc_atom>` | :doc:`temp (k) <compute_temp>` | :doc:`temp/asphere <compute_temp_asphere>` | :doc:`temp/body <compute_temp_body>` | :doc:`temp/chunk <compute_temp_chunk>` |
| :doc:`stress/mop/profile <compute_stress_mop>` | :doc:`stress/tally <compute_tally>` | :doc:`tdpd/cc/atom <compute_tdpd_cc_atom>` | :doc:`temp (k) <compute_temp>` | :doc:`temp/asphere <compute_temp_asphere>` | :doc:`temp/body <compute_temp_body>` |
+------------------------------------------------------------+--------------------------------------------------------+------------------------------------------------------------------+--------------------------------------------------------------+----------------------------------------------------------+--------------------------------------------------------------+
| :doc:`temp/com <compute_temp_com>` | :doc:`temp/cs <compute_temp_cs>` | :doc:`temp/deform <compute_temp_deform>` | :doc:`temp/deform/eff <compute_temp_deform_eff>` | :doc:`temp/drude <compute_temp_drude>` | :doc:`temp/eff <compute_temp_eff>` |
| :doc:`temp/chunk <compute_temp_chunk>` | :doc:`temp/com <compute_temp_com>` | :doc:`temp/cs <compute_temp_cs>` | :doc:`temp/deform <compute_temp_deform>` | :doc:`temp/deform/eff <compute_temp_deform_eff>` | :doc:`temp/drude <compute_temp_drude>` |
+------------------------------------------------------------+--------------------------------------------------------+------------------------------------------------------------------+--------------------------------------------------------------+----------------------------------------------------------+--------------------------------------------------------------+
| :doc:`temp/partial <compute_temp_partial>` | :doc:`temp/profile <compute_temp_profile>` | :doc:`temp/ramp <compute_temp_ramp>` | :doc:`temp/region <compute_temp_region>` | :doc:`temp/region/eff <compute_temp_region_eff>` | :doc:`temp/rotate <compute_temp_rotate>` |
| :doc:`temp/eff <compute_temp_eff>` | :doc:`temp/partial <compute_temp_partial>` | :doc:`temp/profile <compute_temp_profile>` | :doc:`temp/ramp <compute_temp_ramp>` | :doc:`temp/region <compute_temp_region>` | :doc:`temp/region/eff <compute_temp_region_eff>` |
+------------------------------------------------------------+--------------------------------------------------------+------------------------------------------------------------------+--------------------------------------------------------------+----------------------------------------------------------+--------------------------------------------------------------+
| :doc:`temp/sphere <compute_temp_sphere>` | :doc:`temp/uef <compute_temp_uef>` | :doc:`ti <compute_ti>` | :doc:`torque/chunk <compute_torque_chunk>` | :doc:`vacf <compute_vacf>` | :doc:`vcm/chunk <compute_vcm_chunk>` |
| :doc:`temp/rotate <compute_temp_rotate>` | :doc:`temp/sphere <compute_temp_sphere>` | :doc:`temp/uef <compute_temp_uef>` | :doc:`ti <compute_ti>` | :doc:`torque/chunk <compute_torque_chunk>` | :doc:`vacf <compute_vacf>` |
+------------------------------------------------------------+--------------------------------------------------------+------------------------------------------------------------------+--------------------------------------------------------------+----------------------------------------------------------+--------------------------------------------------------------+
| :doc:`voronoi/atom <compute_voronoi_atom>` | :doc:`xrd <compute_xrd>` | | | | |
| :doc:`vcm/chunk <compute_vcm_chunk>` | :doc:`voronoi/atom <compute_voronoi_atom>` | :doc:`xrd <compute_xrd>` | | | |
+------------------------------------------------------------+--------------------------------------------------------+------------------------------------------------------------------+--------------------------------------------------------------+----------------------------------------------------------+--------------------------------------------------------------+

View File

@ -287,9 +287,10 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` doc
* :doc:`smd/ulsph/strain/rate <compute_smd_ulsph_strain_rate>` -
* :doc:`smd/ulsph/stress <compute_smd_ulsph_stress>` - per-particle Cauchy stress tensor and von Mises equivalent stress in Smooth Mach Dynamics
* :doc:`smd/vol <compute_smd_vol>` - per-particle volumes and their sum in Smooth Mach Dynamics
* :doc:`sna/atom <compute_sna_atom>` - calculate bispectrum coefficients for each atom
* :doc:`snad/atom <compute_sna_atom>` - derivative of bispectrum coefficients for each atom
* :doc:`snav/atom <compute_sna_atom>` - virial contribution from bispectrum coefficients for each atom
* :doc:`snap <compute_sna_atom>` - bispectrum components and related quantities for a group of atoms
* :doc:`sna/atom <compute_sna_atom>` - bispectrum components for each atom
* :doc:`snad/atom <compute_sna_atom>` - derivative of bispectrum components for each atom
* :doc:`snav/atom <compute_sna_atom>` - virial contribution from bispectrum components for each atom
* :doc:`spin <compute_spin>` - magnetic quantities for a system of atoms having spins
* :doc:`stress/atom <compute_stress_atom>` - stress tensor for each atom
* :doc:`stress/mop <compute_stress_mop>` - normal components of the local stress tensor using the method of planes

View File

@ -9,6 +9,9 @@ compute snad/atom command
compute snav/atom command
=========================
compute snap command
====================
Syntax
""""""
@ -17,7 +20,8 @@ Syntax
compute ID group-ID sna/atom rcutfac rfac0 twojmax R_1 R_2 ... w_1 w_2 ... keyword values ...
compute ID group-ID snad/atom rcutfac rfac0 twojmax R_1 R_2 ... w_1 w_2 ... keyword values ...
compute ID group-ID snav/atom rcutfac rfac0 twojmax R_1 R_2 ... w_1 w_2 ... keyword values ...
compute ID group-ID snav/atom rcutfac rfac0 twojmax R_1 R_2 ... w_1 w_2 ... keyword values ...
compute ID group-ID snap rcutfac rfac0 twojmax R_1 R_2 ... w_1 w_2 ... keyword values ...
* ID, group-ID are documented in :doc:`compute <compute>` command
* sna/atom = style name of this compute command
@ -53,12 +57,17 @@ Examples
compute b all sna/atom 1.4 0.99363 6 2.0 2.4 0.75 1.0 rmin0 0.0
compute db all sna/atom 1.4 0.95 6 2.0 1.0
compute vb all sna/atom 1.4 0.95 6 2.0 1.0
compute snap all snap 1.4 0.95 6 2.0 1.0
Description
"""""""""""
Define a computation that calculates a set of bispectrum components
for each atom in a group.
Define a computation that calculates a set of quantities related to the
bispectrum components of the atoms in a group. These computes are
used primarily for calculating the dependence of energy, force, and
stress components on the linear coefficients in the
:doc:`snap pair\_style <pair_snap>`, which is useful when training a
SNAP potential to match target data.
Bispectrum components of an atom are order parameters characterizing
the radial and angular distribution of neighbor atoms. The detailed
@ -148,6 +157,30 @@ Again, the sum is over all atoms *i'* of atom type *I*\ . For each atom
virial components, each atom type, and each bispectrum component. See
section below on output for a detailed explanation.
Compute *snap* calculates a global array contains information related
to all three of the above per-atom computes *sna/atom*\ , *snad/atom*\ ,
and *snav/atom*\ . The first row of the array contains the summation of
*sna/atom* over all atoms, but broken out by type. The last six rows
of the array contain the summation of *snav/atom* over all atoms, broken
out by type. In between these are 3\*\ *N* rows containing the same values
computed by *snad/atom* (these are already summed over all atoms and
broken out by type). 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 fitting SNAP.
The potential energy calculation uses the built in compute *thermo\_pe*.
The stress calculation uses a compute called *snap\_press* that is
automatically created behind the scenes, according to the following
command:
.. parsed-literal::
compute snap_press all pressure NULL virial
See section below on output for a detailed explanation of the data
layout in the global array.
The value of all bispectrum components will be zero for atoms not in
the group. Neighbor atoms not in the group do not contribute to the
bispectrum of atoms in the group.
@ -239,10 +272,25 @@ block contains six sub-blocks corresponding to the *xx*\ , *yy*\ , *zz*\ ,
notation. Each of these sub-blocks contains one column for each
bispectrum component, the same as for compute *sna/atom*
Compute *snap* evaluates a global array.
The columns are arranged into
*ntypes* blocks, listed in order of atom type *I*\ . Each block
contains one column for each bispectrum component, the same as for compute
*sna/atom*\ . 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: *sna/atom* quantities summed for all atoms of type *I*
* 3\*\ *N* rows: *snad/atom* quantities, with derivatives w.r.t. x, y, and z coordinate of atom *i* appearing in consecutive rows. The atoms are sorted based on atom ID.
* 6 rows: *snav/atom* quantities summed for all atoms of type *I*
For example, if *K* =30 and ntypes=1, the number of columns in the per-atom
arrays generated by *sna/atom*\ , *snad/atom*\ , and *snav/atom*
are 30, 90, and 180, respectively. With *quadratic* value=1,
the numbers of columns are 930, 2790, and 5580, respectively.
The number of columns in the global array generated by *snap*
are 31, and 931, respectively, while the number of rows is
1+3\*\ *N*\ +6, where *N* is the total number of atoms.
If the *quadratic* keyword value is set to 1, then additional
columns are generated, corresponding to

View File

@ -131,6 +131,7 @@ KOKKOS, o = USER-OMP, t = OPT.
"smd/ulsph/strain/rate"_compute_smd_ulsph_strain_rate.html,
"smd/ulsph/stress"_compute_smd_ulsph_stress.html,
"smd/vol"_compute_smd_vol.html,
"snap"_compute_sna_atom.html,
"sna/atom"_compute_sna_atom.html,
"snad/atom"_compute_sna_atom.html,
"snav/atom"_compute_sna_atom.html,

View File

@ -278,9 +278,10 @@ compute"_Commands_compute.html doc page are followed by one or more of
"smd/ulsph/strain/rate"_compute_smd_ulsph_strain_rate.html -
"smd/ulsph/stress"_compute_smd_ulsph_stress.html - per-particle Cauchy stress tensor and von Mises equivalent stress in Smooth Mach Dynamics
"smd/vol"_compute_smd_vol.html - per-particle volumes and their sum in Smooth Mach Dynamics
"sna/atom"_compute_sna_atom.html - calculate bispectrum coefficients for each atom
"snad/atom"_compute_sna_atom.html - derivative of bispectrum coefficients for each atom
"snav/atom"_compute_sna_atom.html - virial contribution from bispectrum coefficients for each atom
"snap"_compute_sna_atom.html - bispectrum components and related quantities for a group of atoms
"sna/atom"_compute_sna_atom.html - bispectrum components for each atom
"snad/atom"_compute_sna_atom.html - derivative of bispectrum components for each atom
"snav/atom"_compute_sna_atom.html - virial contribution from bispectrum components for each atom
"spin"_compute_spin.html - magnetic quantities for a system of atoms having spins
"stress/atom"_compute_stress_atom.html - stress tensor for each atom
"stress/mop"_compute_stress_mop.html - normal components of the local stress tensor using the method of planes

View File

@ -9,12 +9,14 @@
compute sna/atom command :h3
compute snad/atom command :h3
compute snav/atom command :h3
compute snap command :h3
[Syntax:]
compute ID group-ID sna/atom rcutfac rfac0 twojmax R_1 R_2 ... w_1 w_2 ... keyword values ...
compute ID group-ID snad/atom rcutfac rfac0 twojmax R_1 R_2 ... w_1 w_2 ... keyword values ...
compute ID group-ID snav/atom rcutfac rfac0 twojmax R_1 R_2 ... w_1 w_2 ... keyword values ... :pre
compute ID group-ID snav/atom rcutfac rfac0 twojmax R_1 R_2 ... w_1 w_2 ... keyword values ...
compute ID group-ID snap rcutfac rfac0 twojmax R_1 R_2 ... w_1 w_2 ... keyword values ... :pre
ID, group-ID are documented in "compute"_compute.html command :ulb,l
sna/atom = style name of this compute command :l
@ -41,12 +43,17 @@ keyword = {rmin0} or {switchflag} or {bzeroflag} or {quadraticflag} :l
compute b all sna/atom 1.4 0.99363 6 2.0 2.4 0.75 1.0 rmin0 0.0
compute db all sna/atom 1.4 0.95 6 2.0 1.0
compute vb all sna/atom 1.4 0.95 6 2.0 1.0 :pre
compute vb all sna/atom 1.4 0.95 6 2.0 1.0
compute snap all snap 1.4 0.95 6 2.0 1.0 :pre
[Description:]
Define a computation that calculates a set of bispectrum components
for each atom in a group.
Define a computation that calculates a set of quantities related to the
bispectrum components of the atoms in a group. These computes are
used primarily for calculating the dependence of energy, force, and
stress components on the linear coefficients in the
"snap pair_style"_pair_snap.html, which is useful when training a
SNAP potential to match target data.
Bispectrum components of an atom are order parameters characterizing
the radial and angular distribution of neighbor atoms. The detailed
@ -130,6 +137,27 @@ Again, the sum is over all atoms {i'} of atom type {I}. For each atom
virial components, each atom type, and each bispectrum component. See
section below on output for a detailed explanation.
Compute {snap} calculates a global array contains information related
to all three of the above per-atom computes {sna/atom}, {snad/atom},
and {snav/atom}. The first row of the array contains the summation of
{sna/atom} over all atoms, but broken out by type. The last six rows
of the array contain the summation of {snav/atom} over all atoms, broken
out by type. In between these are 3*{N} rows containing the same values
computed by {snad/atom} (these are already summed over all atoms and
broken out by type). 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 fitting SNAP.
The potential energy calculation uses the built in compute {thermo_pe}.
The stress calculation uses a compute called {snap_press} that is
automatically created behind the scenes, according to the following
command:
compute snap_press all pressure NULL virial :pre
See section below on output for a detailed explanation of the data
layout in the global array.
The value of all bispectrum components will be zero for atoms not in
the group. Neighbor atoms not in the group do not contribute to the
bispectrum of atoms in the group.
@ -214,10 +242,25 @@ block contains six sub-blocks corresponding to the {xx}, {yy}, {zz},
notation. Each of these sub-blocks contains one column for each
bispectrum component, the same as for compute {sna/atom}
Compute {snap} evaluates a global array.
The columns are arranged into
{ntypes} blocks, listed in order of atom type {I}. Each block
contains one column for each bispectrum component, the same as for compute
{sna/atom}. 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: {sna/atom} quantities summed for all atoms of type {I}
3*{N} rows: {snad/atom} quantities, with derivatives w.r.t. x, y, and z coordinate of atom {i} appearing in consecutive rows. The atoms are sorted based on atom ID.
6 rows: {snav/atom} quantities summed for all atoms of type {I} :ul
For example, if {K} =30 and ntypes=1, the number of columns in the per-atom
arrays generated by {sna/atom}, {snad/atom}, and {snav/atom}
are 30, 90, and 180, respectively. With {quadratic} value=1,
the numbers of columns are 930, 2790, and 5580, respectively.
The number of columns in the global array generated by {snap}
are 31, and 931, respectively, while the number of rows is
1+3*{N}+6, where {N} is the total number of atoms.
If the {quadratic} keyword value is set to 1, then additional
columns are generated, corresponding to

View File

@ -0,0 +1,94 @@
# Demonstrate bispectrum computes
# initialize simulation
variable nsteps index 0
variable nrep equal 1
#variable a equal 3.316
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
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
# 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 snap_options string &
"${rcutfac} ${rfac0} ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag 0 bzeroflag 0 switchflag 0"
# set up dummy potential to satisfy cutoff
pair_style zero ${rcutfac}
pair_coeff * *
# 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}
# 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
# set up compute snap generating global array
compute snap all snap ${snap_options}
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: xy component of Sum(Sum(r_j*dB_{222}^i/dr_j), all i of type 2), all j)
# followed by counterparts from compute snap
thermo_style custom &
pe pxy c_bsum2[1] c_vbsum[60] &
c_snap[1][11] c_snap[13][11] c_snap[1][6] c_snap[13][10]
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}

526
src/SNAP/compute_snap.cpp Normal file
View File

@ -0,0 +1,526 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
/* IDEAS
-DONE: Need to define a local peratom array for snad and snad on local and ghost atoms
-DONE: Reverse communicate local peratom array
-DONE: Copy peratom array into output array
-DONE: size_array_cols = nperdim (ncoeff [+quadratic])
-DONE: size_array_rows = 1 + total number of atoms + 6
-DONE: size_peratom = (3+6)*nperdim*ntypes
INCOMPLETE: Mappy from local to global
INCOMPLETE: modify->find_compute()
DONE: eliminate local peratom array for viral, replace with fdotr
*/
#include "compute_snap.h"
#include <cstring>
#include <cstdlib>
#include "sna.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};
ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), cutsq(NULL), list(NULL), snap(NULL),
radelem(NULL), wjelem(NULL), snap_peratom(NULL), snapall(NULL)
{
array_flag = 1;
extarray = 0;
double rfac0, rmin0;
int twojmax, switchflag, bzeroflag;
radelem = NULL;
wjelem = NULL;
int ntypes = atom->ntypes;
int nargmin = 6+2*ntypes;
if (narg < nargmin) error->all(FLERR,"Illegal compute snap command");
// default values
rmin0 = 0.0;
switchflag = 1;
bzeroflag = 1;
quadraticflag = 0;
// process required arguments
memory->create(radelem,ntypes+1,"snap:radelem"); // offset by 1 to match up with types
memory->create(wjelem,ntypes+1,"snap:wjelem");
rcutfac = atof(arg[3]);
rfac0 = atof(arg[4]);
twojmax = atoi(arg[5]);
for(int i = 0; i < ntypes; i++)
radelem[i+1] = atof(arg[6+i]);
for(int i = 0; i < ntypes; i++)
wjelem[i+1] = atof(arg[6+ntypes+i]);
// construct cutsq
double cut;
cutmax = 0.0;
memory->create(cutsq,ntypes+1,ntypes+1,"snap:cutsq");
for(int i = 1; i <= ntypes; i++) {
cut = 2.0*radelem[i]*rcutfac;
if (cut > cutmax) cutmax = cut;
cutsq[i][i] = cut*cut;
for(int j = i+1; j <= ntypes; j++) {
cut = (radelem[i]+radelem[j])*rcutfac;
cutsq[i][j] = cutsq[j][i] = cut*cut;
}
}
// process optional args
int iarg = nargmin;
while (iarg < narg) {
if (strcmp(arg[iarg],"rmin0") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snap command");
rmin0 = atof(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"bzeroflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snap command");
bzeroflag = atoi(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"switchflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snap command");
switchflag = atoi(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"quadraticflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snap command");
quadraticflag = atoi(arg[iarg+1]);
iarg += 2;
} else error->all(FLERR,"Illegal compute snap command");
}
snaptr = new SNA(lmp,rfac0,twojmax,
rmin0,switchflag,bzeroflag);
ncoeff = snaptr->ncoeff;
nperdim = ncoeff;
if (quadraticflag) nperdim += (ncoeff*(ncoeff+1))/2;
ndims_force = 3;
ndims_virial = 6;
yoffset = nperdim;
zoffset = 2*nperdim;
natoms = atom->natoms;
size_array_rows = 1+ndims_force*natoms+ndims_virial;
size_array_cols = nperdim*atom->ntypes+1;
lastcol = size_array_cols-1;
ndims_peratom = ndims_force;
size_peratom = ndims_peratom*nperdim*atom->ntypes;
nmax = 0;
}
/* ---------------------------------------------------------------------- */
ComputeSnap::~ComputeSnap()
{
memory->destroy(snap);
memory->destroy(snapall);
memory->destroy(snap_peratom);
memory->destroy(radelem);
memory->destroy(wjelem);
memory->destroy(cutsq);
delete snaptr;
}
/* ---------------------------------------------------------------------- */
void ComputeSnap::init()
{
if (force->pair == NULL)
error->all(FLERR,"Compute snap requires a pair style be defined");
if (cutmax > force->pair->cutforce)
error->all(FLERR,"Compute snap 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,"snap") == 0) count++;
if (count > 1 && comm->me == 0)
error->warning(FLERR,"More than one compute snap");
snaptr->init();
// allocate memory for global array
memory->create(snap,size_array_rows,size_array_cols,
"snap:snap");
memory->create(snapall,size_array_rows,size_array_cols,
"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);
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;
int ivirial = modify->find_compute(id_virial);
if (ivirial == -1)
error->all(FLERR,"compute snap_press does not exist.");
c_virial = modify->compute[ivirial];
}
/* ---------------------------------------------------------------------- */
void ComputeSnap::init_list(int /*id*/, NeighList *ptr)
{
list = ptr;
}
/* ---------------------------------------------------------------------- */
void ComputeSnap::compute_array()
{
int ntotal = atom->nlocal + atom->nghost;
invoked_array = update->ntimestep;
// grow snap_peratom array if necessary
if (atom->nmax > nmax) {
memory->destroy(snap_peratom);
nmax = atom->nmax;
memory->create(snap_peratom,nmax,size_peratom,
"snap:snap_peratom");
}
// clear global array
for (int irow = 0; irow < size_array_rows; irow++)
for (int icoeff = 0; icoeff < size_array_cols; icoeff++)
snap[irow][icoeff] = 0.0;
// clear local peratom array
for (int i = 0; i < ntotal; i++)
for (int icoeff = 0; icoeff < size_peratom; icoeff++) {
snap_peratom[i][icoeff] = 0.0;
}
// invoke full neighbor list (will copy or build if necessary)
neighbor->build_one(list);
const int inum = list->inum;
const int* const ilist = list->ilist;
const int* const numneigh = list->numneigh;
int** const firstneigh = list->firstneigh;
int * const type = atom->type;
// compute sna derivatives for each atom in group
// use full neighbor list to count atoms less than cutoff
double** const x = atom->x;
const int* const mask = atom->mask;
for (int ii = 0; ii < inum; ii++) {
const int i = ilist[ii];
if (mask[i] & groupbit) {
const double xtmp = x[i][0];
const double ytmp = x[i][1];
const double ztmp = x[i][2];
const int itype = type[i];
const double radi = radelem[itype];
const int* const jlist = firstneigh[i];
const int jnum = numneigh[i];
const int typeoffset_local = ndims_peratom*nperdim*(itype-1);
const int typeoffset_global = nperdim*(itype-1);
// insure rij, inside, and typej 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
// typej = types of neighbors of I within cutoff
// note Rij sign convention => dU/dRij = dU/dRj = -dU/dRi
int ninside = 0;
for (int jj = 0; jj < jnum; jj++) {
int j = jlist[jj];
j &= NEIGHMASK;
const double delx = x[j][0] - xtmp;
const double dely = x[j][1] - ytmp;
const double delz = x[j][2] - ztmp;
const double rsq = delx*delx + dely*dely + delz*delz;
int jtype = type[j];
if (rsq < cutsq[itype][jtype]&&rsq>1e-20) {
snaptr->rij[ninside][0] = delx;
snaptr->rij[ninside][1] = dely;
snaptr->rij[ninside][2] = delz;
snaptr->inside[ninside] = j;
snaptr->wj[ninside] = wjelem[jtype];
snaptr->rcutij[ninside] = (radi+radelem[jtype])*rcutfac;
ninside++;
}
}
snaptr->compute_ui(ninside);
snaptr->compute_zi();
snaptr->compute_bi();
for (int jj = 0; jj < ninside; jj++) {
const int j = snaptr->inside[jj];
snaptr->compute_duidrj(snaptr->rij[jj],
snaptr->wj[jj],
snaptr->rcutij[jj],jj);
snaptr->compute_dbidrj();
// Accumulate dBi/dRi, -dBi/dRj
double *snadi = snap_peratom[i]+typeoffset_local;
double *snadj = snap_peratom[j]+typeoffset_local;
for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
snadi[icoeff] += snaptr->dblist[icoeff][0];
snadi[icoeff+yoffset] += snaptr->dblist[icoeff][1];
snadi[icoeff+zoffset] += snaptr->dblist[icoeff][2];
snadj[icoeff] -= snaptr->dblist[icoeff][0];
snadj[icoeff+yoffset] -= snaptr->dblist[icoeff][1];
snadj[icoeff+zoffset] -= snaptr->dblist[icoeff][2];
}
if (quadraticflag) {
const int quadraticoffset = ncoeff;
snadi += quadraticoffset;
snadj += quadraticoffset;
int ncount = 0;
for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
double bi = snaptr->blist[icoeff];
double bix = snaptr->dblist[icoeff][0];
double biy = snaptr->dblist[icoeff][1];
double biz = snaptr->dblist[icoeff][2];
// diagonal elements of quadratic matrix
double dbxtmp = bi*bix;
double dbytmp = bi*biy;
double dbztmp = bi*biz;
snadi[ncount] += dbxtmp;
snadi[ncount+yoffset] += dbytmp;
snadi[ncount+zoffset] += dbztmp;
snadj[ncount] -= dbxtmp;
snadj[ncount+yoffset] -= dbytmp;
snadj[ncount+zoffset] -= dbztmp;
ncount++;
// upper-triangular elements of quadratic matrix
for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) {
double dbxtmp = bi*snaptr->dblist[jcoeff][0]
+ bix*snaptr->blist[jcoeff];
double dbytmp = bi*snaptr->dblist[jcoeff][1]
+ biy*snaptr->blist[jcoeff];
double dbztmp = bi*snaptr->dblist[jcoeff][2]
+ biz*snaptr->blist[jcoeff];
snadi[ncount] += dbxtmp;
snadi[ncount+yoffset] += dbytmp;
snadi[ncount+zoffset] += dbztmp;
snadj[ncount] -= dbxtmp;
snadj[ncount+yoffset] -= dbytmp;
snadj[ncount+zoffset] -= dbztmp;
ncount++;
}
}
}
}
// Accumulate Bi
// linear contributions
for (int icoeff = 0; icoeff < ncoeff; icoeff++)
snap[0][icoeff+typeoffset_global] += snaptr->blist[icoeff];
// quadratic contributions
if (quadraticflag) {
for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
double bveci = snaptr->blist[icoeff];
snap[0][icoeff+typeoffset_global] += 0.5*bveci*bveci;
for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) {
double bvecj = snaptr->blist[jcoeff];
snap[0][icoeff+typeoffset_global] += bveci*bvecj;
}
}
}
}
}
// accumulate bispectrum force contributions to global array
for (int itype = 0; itype < atom->ntypes; itype++) {
const int typeoffset_local = ndims_peratom*nperdim*itype;
const int typeoffset_global = nperdim*itype;
for (int icoeff = 0; icoeff < nperdim; icoeff++) {
int irow = 1;
for (int i = 0; i < ntotal; i++) {
double *snadi = snap_peratom[i]+typeoffset_local;
int iglobal = atom->tag[i];
int irow = 3*(iglobal-1)+1;
snap[irow][icoeff+typeoffset_global] += snadi[icoeff];
snap[irow+1][icoeff+typeoffset_global] += snadi[icoeff+yoffset];
snap[irow+2][icoeff+typeoffset_global] += snadi[icoeff+zoffset];
}
}
}
// accumulate forces to global array
for (int i = 0; i < atom->nlocal; i++) {
int iglobal = atom->tag[i];
int irow = 3*(iglobal-1)+1;
snap[irow][lastcol] = atom->f[i][0];
snap[irow+1][lastcol] = atom->f[i][1];
snap[irow+2][lastcol] = atom->f[i][2];
}
// accumulate bispectrum virial contributions to global array
dbdotr_compute();
// sum up over all processes
MPI_Allreduce(&snap[0][0],&snapall[0][0],size_array_rows*size_array_cols,MPI_DOUBLE,MPI_SUM,world);
// assign energy to last column
int irow = 0;
double reference_energy = c_pe->compute_scalar();
snapall[irow++][lastcol] = reference_energy;
// assign virial stress to last column
// switch to Voigt notation
c_virial->compute_vector();
irow += 3*natoms;
snapall[irow++][lastcol] = c_virial->vector[0];
snapall[irow++][lastcol] = c_virial->vector[1];
snapall[irow++][lastcol] = c_virial->vector[2];
snapall[irow++][lastcol] = c_virial->vector[5];
snapall[irow++][lastcol] = c_virial->vector[4];
snapall[irow++][lastcol] = c_virial->vector[3];
}
/* ----------------------------------------------------------------------
compute global virial contributions via summing r_i.dB^j/dr_i over
own & ghost atoms
------------------------------------------------------------------------- */
void ComputeSnap::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 itype = 0; itype < atom->ntypes; itype++) {
const int typeoffset_local = ndims_peratom*nperdim*itype;
const int typeoffset_global = nperdim*itype;
double *snadi = snap_peratom[i]+typeoffset_local;
for (int icoeff = 0; icoeff < nperdim; icoeff++) {
double dbdx = snadi[icoeff];
double dbdy = snadi[icoeff+yoffset];
double dbdz = snadi[icoeff+zoffset];
int irow = irow0;
snap[irow++][icoeff+typeoffset_global] += dbdx*x[i][0];
snap[irow++][icoeff+typeoffset_global] += dbdy*x[i][1];
snap[irow++][icoeff+typeoffset_global] += dbdz*x[i][2];
snap[irow++][icoeff+typeoffset_global] += dbdz*x[i][1];
snap[irow++][icoeff+typeoffset_global] += dbdz*x[i][0];
snap[irow++][icoeff+typeoffset_global] += dbdy*x[i][0];
}
}
}
/* ----------------------------------------------------------------------
memory usage
------------------------------------------------------------------------- */
double ComputeSnap::memory_usage()
{
double bytes = size_array_rows*size_array_cols *
sizeof(double); // snap
bytes += size_array_rows*size_array_cols *
sizeof(double); // snapall
bytes += nmax*size_peratom * sizeof(double); // snap_peratom
bytes += snaptr->memory_usage(); // SNA object
return bytes;
}

82
src/SNAP/compute_snap.h Normal file
View File

@ -0,0 +1,82 @@
/* -*- 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(snap,ComputeSnap)
#else
#ifndef LMP_COMPUTE_SNAP_H
#define LMP_COMPUTE_SNAP_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeSnap : public Compute {
public:
ComputeSnap(class LAMMPS *, int, char **);
~ComputeSnap();
void init();
void init_list(int, class NeighList *);
void compute_array();
double memory_usage();
private:
int natoms, nmax, size_peratom, lastcol;
int ncoeff, nperdim, yoffset, zoffset;
int ndims_peratom, ndims_force, ndims_virial;
double **cutsq;
class NeighList *list;
double **snap, **snapall;
double **snap_peratom;
double rcutfac;
double *radelem;
double *wjelem;
class SNA* snaptr;
double cutmax;
int quadraticflag;
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.
*/