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

Update example and docs.

Browse Source
This commit is contained in:
rohskopf
2022-06-22 09:32:41 -06:00
parent 0e6bbf8dff
commit 2396c16026
7 changed files with 139 additions and 726 deletions
Download Patch File Download Diff File Expand all files Collapse all files

43
doc/src/compute_sna_atom.rst
View File

@ -33,7 +33,7 @@ Syntax
* R_1, R_2,... = list of cutoff radii, one for each type (distance units)
* w_1, w_2,... = list of neighbor weights, one for each type
* zero or more keyword/value pairs may be appended
* keyword = *rmin0* or *switchflag* or *bzeroflag* or *quadraticflag* or *chem* or *bnormflag* or *wselfallflag* or *bikflag* or *switchinnerflag* or *sinner* or *dinner*
* keyword = *rmin0* or *switchflag* or *bzeroflag* or *quadraticflag* or *chem* or *bnormflag* or *wselfallflag* or *bikflag* or *switchinnerflag* or *sinner* or *dinner* or *dgradflag*
.. parsed-literal::
@ -66,6 +66,9 @@ Syntax
*sinnerlist* = *ntypes* values of *Sinner* (distance units)
*dinner* values = *dinnerlist*
*dinnerlist* = *ntypes* values of *Dinner* (distance units)
*dgradflag* value = *0* or *1*
*0* = bispectrum descriptor gradients are summed over neighbors
*1* = bispectrum descriptor gradients are not summed over neighbors
Examples
""""""""
@ -340,6 +343,14 @@ When the central atom and the neighbor atom have different types,
the values of :math:`S_{inner}` and :math:`D_{inner}` are
the arithmetic means of the values for both types.
The keyword *dgradflag* determines whether or not to sum the bispectrum descriptor gradients over neighboring atoms *i'*
as explained with *snad/atom* above. If *dgradflag* is set to 1 then the descriptor gradient rows of the global snap array
are not summed over atoms *i'*. Instead, each row corresponds to a single term :math:`\frac{\partial {B_{i,k} }}{\partial {r}^a_j}`
where :math:`a` is the Cartesian direction for the gradient. This also changes
the number of columns to be equal to the number of bispectrum components, with 3
additional columns representing the indices :math:`i`, :math:`j`, and :math:`a`,
as explained more in the Output info section below. The option *dgradflag=1* must be used with *bikflag=1*.
.. note::
If you have a bonded system, then the settings of :doc:`special_bonds
@ -435,6 +446,36 @@ components. For the purposes of handling contributions to force, virial,
and quadratic combinations, these :math:`N_{elem}^3` sub-blocks are
treated as a single block of :math:`K N_{elem}^3` columns.
If the *bik* keyword is set to 1, then the first :math:`N` rows of the snap array
correspond to :math:`B_{i,k}` instead of the sum over atoms :math:`i`. In this case, the entries in the final column for these rows
are set to zero.
If the *dgradflag* keyword is set to 1, this changes the structure of the snap array completely.
Here the *snad/atom* quantities are replaced with rows corresponding to descriptor
gradient components
.. math::
\frac{\partial {B_{i,k} }}{\partial {r}^a_j}
where :math:`a` is the Cartesian direction for the gradient. The rows are organized in chunks, where each chunk corresponds to
an atom :math:`j` in the system of :math:`N` atoms. The rows in an atom :math:`j` chunk correspond to neighbors :math:`i` of :math:`j`.
The number of rows in the atom :math:`j` chunk is therefore equal to the number of neighbors :math:`N_{neighs}[j]` within the SNAP
potential cutoff radius of atom :math:`j`, times 3 for each Cartesian direction.
The total number of rows for these descriptor gradients is therefore
.. math::
3 \sum_j^{N} N_{neighs}[j].
For *dgradflag=1*, the number of columns is equal to the number of bispectrum components,
plus 3 additional columns representing the indices :math:`i`, :math:`j`, and :math:`a` which
identify the atoms :math:`i` and :math:`j`, and Cartesian direction :math:`a` for which
a particular gradient :math:`\frac{\partial {B_{i,k} }}{\partial {r}^a_j}` belongs to. The reference energy and forces are also located in different parts of the array.
The last 3 columns of the first :math:`N` rows belong to the reference potential force components.
The first column of the last row, after the first :math:`N + 3 \sum_j^{N} N_{neighs}[j]` rows,
contains the reference potential energy. The virial components are not used with this option.
These values can be accessed by any command that uses per-atom values
from a compute as input. See the :doc:`Howto output <Howto_output>` doc
page for an overview of LAMMPS output options. To see how this command

10
examples/snap/README.md
View File

@ -1,9 +1 @@
See `compute_snap_dgrad.py` for a test that compares the dBi/dRj from compute snap (`dbirjflag=1`) to the sum of dBi/dRj from usual compute snap (`dbirjflag=0`).
The format of the global array from `dbirjflag=1` is as follows.
The first N rows belong to bispectrum components for each atom, if we use `bikflag=1`. The first `K` columns correspond to each bispectrum coefficient. The final 3 columns contain reference force components for each atom.
The rows after the first N rows contain dBi/dRj values for all pairs. These values are arranged in row chunks for each atom `j`, where all the rows in a chunk are associated with the neighbors `i` of `j`, as well as the self-terms where `i=j`. So for atom `j`, the number of rows is equal to the number of atoms within the SNAP cutoff, plus 1 for the `i=j` terms, times 3 for each Cartesian component. The total number of dBi/dRj rows is therefore equal to `N*(Nneigh+1)*3`, and `Nneigh` may be different for each atom. To facilitate with determining which row belong to which atom pair `ij`, the last 3 columns contain indices; the 3rd to last column contains global indices of atoms `i` (the neighbors), the 2nd to last column contains global indices of atoms `j`, and the last column contains an index 0,1,2 for the Cartesian component. Like the `bik` rows, the first `K` columns correspond to each bispectrum coefficient.
Finally, the first column of the last row contains the reference energy.
See `compute_snap_dgrad.py` for a test that compares the dBi/dRj from compute snap (`dgradflag=1`) to the sum of dBi/dRj from usual compute snap (`dgradflag=0`).

24
examples/snap/compute_snap_dgrad.py
View File

@ -1,8 +1,8 @@
"""
compute_snap_dgrad.py
Purpose: Demonstrate extraction of descriptor gradient (dB/dR) array from compute snap.
Show that dBi/dRj components summed over neighbors i yields same output as regular compute snap with dbirjflag=0.
This shows that the dBi/dRj components extracted with dbirjflag=1 are correct.
Show that dBi/dRj components summed over neighbors i yields same output as regular compute snap with dgradflag=0.
This shows that the dBi/dRj components extracted with dgradflag=1 are correct.
Serial syntax:
python compute_snap_dgrad.py
Parallel syntax:
@ -24,7 +24,7 @@ from lammps import lammps, LMP_TYPE_ARRAY, LMP_STYLE_GLOBAL
cmds = ["-screen", "none", "-log", "none"]
lmp = lammps(cmdargs=cmds)
def run_lammps(dbirjflag):
def run_lammps(dgradflag):
lmp.command("clear")
lmp.command("units metal")
lmp.command("boundary p p p")
@ -36,7 +36,7 @@ def run_lammps(dbirjflag):
lmp.command("mass * 180.88")
lmp.command("displace_atoms all random 0.01 0.01 0.01 123456")
# Pair style
snap_options=f'{rcutfac} {rfac0} {twojmax} {radelem1} {radelem2} {wj1} {wj2} rmin0 {rmin0} quadraticflag {quadratic} bzeroflag {bzero} switchflag {switch} bikflag {bikflag} dbirjflag {dbirjflag}'
snap_options=f'{rcutfac} {rfac0} {twojmax} {radelem1} {radelem2} {wj1} {wj2} rmin0 {rmin0} quadraticflag {quadratic} bzeroflag {bzero} switchflag {switch} bikflag {bikflag} dgradflag {dgradflag}'
lmp.command(f"pair_style zero {rcutfac}")
lmp.command(f"pair_coeff * *")
lmp.command(f"pair_style zbl {zblcutinner} {zblcutouter}")
@ -70,7 +70,7 @@ quadratic=0
bzero=0
switch=0
bikflag=1
dbirjflag=1
dgradflag=1
# Declare reference potential variables
zblcutinner=4.0
@ -84,7 +84,8 @@ else:
nd = int(m*(m+1)*(m+2)/3)
print(f"Number of descriptors based on twojmax : {nd}")
# Run lammps with dbirjflag on
# Run lammps with dgradflag on
print("Running with dgradflag on")
run_lammps(1)
# Get global snap array
@ -93,12 +94,12 @@ lmp_snap = lmp.numpy.extract_compute("snap",0, 2)
# Extract dBj/dRi (includes dBi/dRi)
natoms = lmp.get_natoms()
fref1 = lmp_snap[0:natoms,-3:].flatten()
eref1 = lmp_snap[-6,0]
dbdr_length = np.shape(lmp_snap)[0]-(natoms)-6 # Length of neighborlist pruned dbdr array
eref1 = lmp_snap[-1,0] #lmp_snap[-6,0]
dbdr_length = np.shape(lmp_snap)[0]-(natoms) - 1 #-6 # Length of neighborlist pruned dbdr array
dBdR = lmp_snap[natoms:(natoms+dbdr_length),:]
force_indices = lmp_snap[natoms:(natoms+dbdr_length),-3:].astype(np.int32)
# Sum over neighbors j for each atom i, like dbirjflag=0 does.
# Sum over neighbors j for each atom i, like dgradflag=0 does.
array1 = np.zeros((3*natoms,nd))
a = 0
for k in range(0,nd):
@ -111,7 +112,8 @@ for k in range(0,nd):
if (a>2):
a=0
# Run lammps with dbirjflag off
# Run lammps with dgradflag off
print("Running with dgradflag off")
run_lammps(0)
# Get global snap array
@ -121,7 +123,7 @@ fref2 = lmp_snap[natoms:(natoms+3*natoms),-1]
eref2 = lmp_snap[0,-1]
array2 = lmp_snap[natoms:natoms+(3*natoms), nd:-1]
# Sum the arrays obtained from dbirjflag on and off.
# Sum the arrays obtained from dgradflag on and off.
summ = array1 + array2
#np.savetxt("sum.dat", summ)

138
src/ML-SNAP/compute_snap.cpp
View File

@ -56,7 +56,7 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) :
bzeroflag = 1;
quadraticflag = 0;
bikflag = 0;
dbirjflag = 0;
dgradflag = 0;
chemflag = 0;
bnormflag = 0;
wselfallflag = 0;
@ -148,10 +148,10 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) :
error->all(FLERR,"Illegal compute snap command");
bikflag = atoi(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"dbirjflag") == 0) {
} else if (strcmp(arg[iarg],"dgradflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snap command");
dbirjflag = atoi(arg[iarg+1]);
dgradflag = atoi(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"switchinnerflag") == 0) {
if (iarg+2 > narg)
@ -185,6 +185,9 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) :
if (!switchinnerflag && (sinnerflag || dinnerflag))
error->all(FLERR,"Illegal compute snap command: switchinnerflag = 0, unexpected sinner/dinner keyword");
if (dgradflag && !bikflag)
error->all(FLERR,"Illegal compute snap command: dgradflag=1 requires bikflag=1");
snaptr = new SNA(lmp, rfac0, twojmax,
rmin0, switchflag, bzeroflag,
chemflag, bnormflag, wselfallflag,
@ -200,9 +203,9 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) :
natoms = atom->natoms;
bik_rows = 1;
if (bikflag) bik_rows = natoms;
dbirj_rows = ndims_force*natoms;
size_array_rows = bik_rows+dbirj_rows+ndims_virial;
if (dbirjflag) size_array_cols = nperdim+3; // plus 3 for tag[i], tag[j], and cartesian index
dgrad_rows = ndims_force*natoms;
size_array_rows = bik_rows+dgrad_rows+ndims_virial;
if (dgradflag) size_array_cols = nperdim+3; // plus 3 for tag[i], tag[j], and cartesian index
else size_array_cols = nperdim*atom->ntypes+1;
lastcol = size_array_cols-1;
@ -230,8 +233,8 @@ ComputeSnap::~ComputeSnap()
memory->destroy(sinnerelem);
memory->destroy(dinnerelem);
}
if (dbirjflag){
memory->destroy(dbirj);
if (dgradflag){
memory->destroy(dgrad);
memory->destroy(nneighs);
memory->destroy(neighsum);
memory->destroy(icounter);
@ -260,10 +263,19 @@ void ComputeSnap::init()
// allocate memory for global array
if (dgradflag){
// Initially allocate natoms^2 rows, will prune with neighborlist later
memory->create(snap,natoms*natoms,ncoeff+3,
"snap:snap");
memory->create(snapall,natoms*natoms,ncoeff+3,
"snap:snapall");
}
else{
memory->create(snap,size_array_rows,size_array_cols,
"snap:snap");
memory->create(snapall,size_array_rows,size_array_cols,
"snap:snapall");
}
array = snapall;
// find compute for reference energy
@ -299,8 +311,8 @@ void ComputeSnap::init_list(int /*id*/, NeighList *ptr)
void ComputeSnap::compute_array()
{
if (dbirjflag){
get_dbirj_length();
if (dgradflag){
get_dgrad_length();
}
int ntotal = atom->nlocal + atom->nghost;
@ -345,7 +357,7 @@ void ComputeSnap::compute_array()
const int* const mask = atom->mask;
int ninside;
int numneigh_sum = 0;
int dbirj_row_indx;
int dgrad_row_indx;
for (int ii = 0; ii < inum; ii++) {
int irow = 0;
if (bikflag) irow = atom->tag[ilist[ii] & NEIGHMASK]-1;
@ -416,8 +428,8 @@ void ComputeSnap::compute_array()
for (int jj = 0; jj < ninside; jj++) {
const int j = snaptr->inside[jj];
if (dbirjflag){
dbirj_row_indx = 3*neighsum[atom->tag[j]-1] + 3*icounter[atom->tag[j]-1] ;
if (dgradflag){
dgrad_row_indx = 3*neighsum[atom->tag[j]-1] + 3*icounter[atom->tag[j]-1] ;
icounter[atom->tag[j]-1] += 1;
}
@ -440,20 +452,20 @@ void ComputeSnap::compute_array()
snadj[icoeff+zoffset] -= snaptr->dblist[icoeff][2];
if (dbirjflag){
dbirj[dbirj_row_indx+0][icoeff] = snaptr->dblist[icoeff][0];
dbirj[dbirj_row_indx+1][icoeff] = snaptr->dblist[icoeff][1];
dbirj[dbirj_row_indx+2][icoeff] = snaptr->dblist[icoeff][2];
if (dgradflag){
dgrad[dgrad_row_indx+0][icoeff] = snaptr->dblist[icoeff][0];
dgrad[dgrad_row_indx+1][icoeff] = snaptr->dblist[icoeff][1];
dgrad[dgrad_row_indx+2][icoeff] = snaptr->dblist[icoeff][2];
if (icoeff==(ncoeff-1)){
dbirj[dbirj_row_indx+0][ncoeff] = atom->tag[i]-1;
dbirj[dbirj_row_indx+0][ncoeff+1] = atom->tag[j]-1;
dbirj[dbirj_row_indx+0][ncoeff+2] = 0;
dbirj[dbirj_row_indx+1][ncoeff] = atom->tag[i]-1;
dbirj[dbirj_row_indx+1][ncoeff+1] = atom->tag[j]-1;
dbirj[dbirj_row_indx+1][ncoeff+2] = 1;
dbirj[dbirj_row_indx+2][ncoeff] = atom->tag[i]-1;
dbirj[dbirj_row_indx+2][ncoeff+1] = atom->tag[j]-1;
dbirj[dbirj_row_indx+2][ncoeff+2] = 2;
dgrad[dgrad_row_indx+0][ncoeff] = atom->tag[i]-1;
dgrad[dgrad_row_indx+0][ncoeff+1] = atom->tag[j]-1;
dgrad[dgrad_row_indx+0][ncoeff+2] = 0;
dgrad[dgrad_row_indx+1][ncoeff] = atom->tag[i]-1;
dgrad[dgrad_row_indx+1][ncoeff+1] = atom->tag[j]-1;
dgrad[dgrad_row_indx+1][ncoeff+2] = 1;
dgrad[dgrad_row_indx+2][ncoeff] = atom->tag[i]-1;
dgrad[dgrad_row_indx+2][ncoeff+1] = atom->tag[j]-1;
dgrad[dgrad_row_indx+2][ncoeff+2] = 2;
}
// Accumulate dBi/dRi = sum (-dBi/dRj) for neighbors j of if i
dbiri[3*(atom->tag[i]-1)+0][icoeff] -= snaptr->dblist[icoeff][0];
@ -530,7 +542,7 @@ void ComputeSnap::compute_array()
// linear contributions
int k;
if (dbirjflag) k = 0;
if (dgradflag) k = 0;
else k = typeoffset_global;
for (int icoeff = 0; icoeff < ncoeff; icoeff++){
snap[irow][k++] += snaptr->blist[icoeff];
@ -555,7 +567,7 @@ void ComputeSnap::compute_array()
// Accumulate contributions to global array
if (dbirjflag){
if (dgradflag){
int dbiri_indx;
int irow;
@ -568,34 +580,34 @@ void ComputeSnap::compute_array()
for (int jj=0; jj<nneighs[i]; jj++){
int dbirj_row_indx = 3*neighsum[atom->tag[i]-1] + 3*jj;
int dgrad_row_indx = 3*neighsum[atom->tag[i]-1] + 3*jj;
int snap_row_indx = 3*neighsum[atom->tag[i]-1] + 3*(atom->tag[i]-1) + 3*jj;
irow = snap_row_indx + bik_rows;
// x-coordinate
snap[irow][icoeff+typeoffset_global] += dbirj[dbirj_row_indx+0][icoeff];
snap[irow][icoeff+typeoffset_global] += dgrad[dgrad_row_indx+0][icoeff];
if (icoeff==(ncoeff-1)){
snap[irow][ncoeff] += dbirj[dbirj_row_indx+0][ncoeff];
snap[irow][ncoeff+1] += dbirj[dbirj_row_indx+0][ncoeff+1];
snap[irow][ncoeff+2] += dbirj[dbirj_row_indx+0][ncoeff+2];
snap[irow][ncoeff] += dgrad[dgrad_row_indx+0][ncoeff];
snap[irow][ncoeff+1] += dgrad[dgrad_row_indx+0][ncoeff+1];
snap[irow][ncoeff+2] += dgrad[dgrad_row_indx+0][ncoeff+2];
}
irow++;
// y-coordinate
snap[irow][icoeff+typeoffset_global] += dbirj[dbirj_row_indx+1][icoeff];
snap[irow][icoeff+typeoffset_global] += dgrad[dgrad_row_indx+1][icoeff];
if (icoeff==(ncoeff-1)){
snap[irow][ncoeff] += dbirj[dbirj_row_indx+1][ncoeff];
snap[irow][ncoeff+1] += dbirj[dbirj_row_indx+1][ncoeff+1];
snap[irow][ncoeff+2] += dbirj[dbirj_row_indx+1][ncoeff+2];
snap[irow][ncoeff] += dgrad[dgrad_row_indx+1][ncoeff];
snap[irow][ncoeff+1] += dgrad[dgrad_row_indx+1][ncoeff+1];
snap[irow][ncoeff+2] += dgrad[dgrad_row_indx+1][ncoeff+2];
}
irow++;
// z-coordinate
snap[irow][icoeff+typeoffset_global] += dbirj[dbirj_row_indx+2][icoeff];
snap[irow][icoeff+typeoffset_global] += dgrad[dgrad_row_indx+2][icoeff];
if (icoeff==(ncoeff-1)){
snap[irow][ncoeff] += dbirj[dbirj_row_indx+2][ncoeff];
snap[irow][ncoeff+1] += dbirj[dbirj_row_indx+2][ncoeff+1];
snap[irow][ncoeff+2] += dbirj[dbirj_row_indx+2][ncoeff+2];
snap[irow][ncoeff] += dgrad[dgrad_row_indx+2][ncoeff];
snap[irow][ncoeff+1] += dgrad[dgrad_row_indx+2][ncoeff+1];
snap[irow][ncoeff+2] += dgrad[dgrad_row_indx+2][ncoeff+2];
}
dbiri_indx = dbiri_indx+3;
}
@ -659,8 +671,8 @@ void ComputeSnap::compute_array()
// accumulate forces to global array
if (dbirjflag){
// for dbirjflag=1, put forces at last 3 columns of bik rows
if (dgradflag){
// for dgradflag=1, put forces at last 3 columns of bik rows
for (int i=0; i<atom->nlocal; i++){
int iglobal = atom->tag[i];
snap[iglobal-1][ncoeff+0] = atom->f[i][0];
@ -681,8 +693,8 @@ void ComputeSnap::compute_array()
// accumulate bispectrum virial contributions to global array
if (dbirjflag){
// no virial terms for dbirj yet
if (dgradflag){
// no virial terms for dgrad yet
}
else{
dbdotr_compute();
@ -691,10 +703,10 @@ void ComputeSnap::compute_array()
// 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
if (dbirjflag){
// Assign reference energy right after the dbirj rows, first column.
if (dgradflag){
// Assign reference energy right after the dgrad rows, first column.
// Add 3N for the dBi/dRi rows.
int irow = bik_rows + dbirj_rows + 3*natoms;
int irow = bik_rows + dgrad_rows + 3*natoms;
double reference_energy = c_pe->compute_scalar();
snapall[irow][0] = reference_energy;
}
@ -709,10 +721,11 @@ void ComputeSnap::compute_array()
// switch to Voigt notation
c_virial->compute_vector();
if (dbirjflag){
// no virial terms for dbirj yet
if (dgradflag){
// no virial terms for dgrad yet
}
else{
c_virial->compute_vector();
int irow = 3*natoms+bik_rows;
snapall[irow++][lastcol] = c_virial->vector[0];
snapall[irow++][lastcol] = c_virial->vector[1];
@ -760,10 +773,10 @@ void ComputeSnap::dbdotr_compute()
}
/* ----------------------------------------------------------------------
compute dbirj length
compute dgrad length
------------------------------------------------------------------------- */
void ComputeSnap::get_dbirj_length()
void ComputeSnap::get_dgrad_length()
{
int rank = universe->me; // for MPI debugging
@ -772,7 +785,7 @@ void ComputeSnap::get_dbirj_length()
memory->destroy(snapall);
// invoke full neighbor list
neighbor->build_one(list);
dbirj_rows = 0;
dgrad_rows = 0;
const int inum = list->inum;
const int* const ilist = list->ilist;
const int* const numneigh = list->numneigh;
@ -791,7 +804,7 @@ void ComputeSnap::get_dbirj_length()
memory->create(icounter, natoms, "snap:icounter");
memory->create(dbiri, 3*natoms,ncoeff+3, "snap:dbiri");
if (atom->nlocal != natoms){
error->all(FLERR,"Compute snap dbirjflag=1 does not support parallelism.");
error->all(FLERR,"Compute snap dgradflag=1 does not support parallelism.");
}
for (int ii=0; ii<3*natoms; ii++){
for (int icoeff=0; icoeff<ncoeff; icoeff++){
@ -822,7 +835,7 @@ void ComputeSnap::get_dbirj_length()
int jtype = type[j];
if (rsq < cutsq[itype][jtype]&&rsq>1e-20) {
dbirj_rows += 1; //jnum + 1;
dgrad_rows += 1; //jnum + 1;
jnum_cutoff += 1;
nneighs[i]+=1;
}
@ -830,7 +843,7 @@ void ComputeSnap::get_dbirj_length()
}
}
dbirj_rows *= ndims_force;
dgrad_rows *= ndims_force;
// Loop over all atoms again to calculate neighsum.
for (int ii = 0; ii < inum; ii++) {
@ -850,14 +863,15 @@ void ComputeSnap::get_dbirj_length()
}
}
memory->create(dbirj, dbirj_rows, ncoeff+3, "snap:dbirj");
for (int i=0; i<dbirj_rows; i++){
memory->create(dgrad, dgrad_rows, ncoeff+3, "snap:dgrad");
for (int i=0; i<dgrad_rows; i++){
for (int j=0; j<ncoeff+3; j++){
dbirj[i][j]=0.0;
dgrad[i][j]=0.0;
}
}
// Set size array rows which now depends on dbirj_rows.
size_array_rows = bik_rows+dbirj_rows+ndims_virial+3*atom->nlocal; // Add 3*N for dBi/dRi
// Set size array rows which now depends on dgrad_rows.
//size_array_rows = bik_rows+dgrad_rows+ndims_virial+3*atom->nlocal; // Add 3*N for dBi/dRi
size_array_rows = bik_rows + dgrad_rows + 3*atom->nlocal + 1; // Add 3*N for dBi/dRi. and add 1 for reference energy
memory->create(snap,size_array_rows,size_array_cols, "snap:snap");
memory->create(snapall,size_array_rows,size_array_cols, "snap:snapall");
@ -871,7 +885,7 @@ void ComputeSnap::get_dbirj_length()
double ComputeSnap::compute_scalar()
{
if (dbirjflag) get_dbirj_length();
if (dgradflag) get_dgrad_length();
return size_array_rows;
}

6
src/ML-SNAP/compute_snap.h
View File

@ -56,8 +56,8 @@ class ComputeSnap : public Compute {
int quadraticflag;
//int bikflag;
//int bik_rows;
int bikflag, bik_rows, dbirjflag, dbirj_rows;
double **dbirj;
int bikflag, bik_rows, dgradflag, dgrad_rows;
double **dgrad;
double **dbiri; // dBi/dRi = sum(-dBi/dRj) over neighbors j
int *nneighs; // number of neighs inside the snap cutoff.
int *neighsum;
@ -67,7 +67,7 @@ class ComputeSnap : public Compute {
Compute *c_virial;
void dbdotr_compute();
void get_dbirj_length();
void get_dgrad_length();
};
} // namespace LAMMPS_NS

563
src/ML-SNAP/compute_snapneigh.cpp
View File

@ -1,563 +0,0 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "compute_snapneigh.h"
#include "sna.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "force.h"
#include "pair.h"
#include "comm.h"
#include "memory.h"
#include "error.h"
#include <cstring>
using namespace LAMMPS_NS;
enum{SCALAR,VECTOR,ARRAY};
ComputeSnapneigh::ComputeSnapneigh(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), cutsq(nullptr), list(nullptr), radelem(nullptr), wjelem(nullptr),
sinnerelem(nullptr), dinnerelem(nullptr)
{
array_flag = 1;
//vector_flag = 1;
extarray = 0;
double rfac0, rmin0;
int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag;
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;
bikflag = 0;
dbirjflag = 0;
chemflag = 0;
bnormflag = 0;
wselfallflag = 0;
switchinnerflag = 0;
nelements = 1;
// process required arguments
memory->create(radelem,ntypes+1,"snapneigh:radelem"); // offset by 1 to match up with types
memory->create(wjelem,ntypes+1,"snapneigh: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,"snapneigh:cutsq");
for (int i = 1; i <= ntypes; i++) {
cut = 2.0*radelem[i]*rcutfac;
//printf("cut: %f\n", cut);
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;
}
}
// set local input checks
int sinnerflag = 0;
int dinnerflag = 0;
// 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 if (strcmp(arg[iarg],"chem") == 0) {
if (iarg+2+ntypes > narg)
error->all(FLERR,"Illegal compute snap command");
chemflag = 1;
memory->create(map,ntypes+1,"compute_snapneigh:map");
nelements = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
for (int i = 0; i < ntypes; i++) {
int jelem = utils::inumeric(FLERR,arg[iarg+2+i],false,lmp);
if (jelem < 0 || jelem >= nelements)
error->all(FLERR,"Illegal compute snap command");
map[i+1] = jelem;
}
iarg += 2+ntypes;
} else if (strcmp(arg[iarg],"bnormflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snap command");
bnormflag = atoi(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"wselfallflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snap command");
wselfallflag = atoi(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"bikflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snap command");
bikflag = atoi(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"dbirjflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snap command");
dbirjflag = atoi(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"switchinnerflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snap command");
switchinnerflag = atoi(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"sinner") == 0) {
iarg++;
if (iarg+ntypes > narg)
error->all(FLERR,"Illegal compute snap command");
memory->create(sinnerelem,ntypes+1,"snapneigh:sinnerelem");
for (int i = 0; i < ntypes; i++)
sinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp);
sinnerflag = 1;
iarg += ntypes;
} else if (strcmp(arg[iarg],"dinner") == 0) {
iarg++;
if (iarg+ntypes > narg)
error->all(FLERR,"Illegal compute snap command");
memory->create(dinnerelem,ntypes+1,"snapneigh:dinnerelem");
for (int i = 0; i < ntypes; i++)
dinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp);
dinnerflag = 1;
iarg += ntypes;
} else error->all(FLERR,"Illegal compute snap command");
}
if (switchinnerflag && !(sinnerflag && dinnerflag))
error->all(FLERR,"Illegal compute snap command: switchinnerflag = 1, missing sinner/dinner keyword");
if (!switchinnerflag && (sinnerflag || dinnerflag))
error->all(FLERR,"Illegal compute snap command: switchinnerflag = 0, unexpected sinner/dinner keyword");
/*
snaptr = new SNA(lmp, rfac0, twojmax,
rmin0, switchflag, bzeroflag,
chemflag, bnormflag, wselfallflag,
nelements, switchinnerflag);
*/
//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;
bik_rows = 1;
if (bikflag) bik_rows = natoms;
//size_array_rows = bik_rows+ndims_force*natoms+ndims_virial;
dbirj_rows = ndims_force*natoms;
size_array_rows = bik_rows+dbirj_rows+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;
}
/* ---------------------------------------------------------------------- */
ComputeSnapneigh::~ComputeSnapneigh()
{
memory->destroy(neighs);
memory->destroy(radelem);
memory->destroy(wjelem);
memory->destroy(cutsq);
//delete snaptr;
if (chemflag) memory->destroy(map);
if (switchinnerflag) {
memory->destroy(sinnerelem);
memory->destroy(dinnerelem);
}
if (dbirjflag){
//printf("dbirj_rows: %d\n", dbirj_rows);
//memory->destroy(dbirj);
memory->destroy(nneighs);
memory->destroy(neighsum);
memory->destroy(icounter);
//memory->destroy(dbiri);
}
}
/* ---------------------------------------------------------------------- */
void ComputeSnapneigh::init()
{
if (force->pair == nullptr)
error->all(FLERR,"Compute snap requires a pair style be defined");
if (cutmax > force->pair->cutforce){
//printf("----- cutmax cutforce: %f %f\n", cutmax, force->pair->cutforce);
error->all(FLERR,"Compute snap cutoff is longer than pairwise cutoff");
}
// need an occasional full neighbor list
neighbor->add_request(this, NeighConst::REQ_FULL | NeighConst::REQ_OCCASIONAL);
}
/* ---------------------------------------------------------------------- */
void ComputeSnapneigh::init_list(int /*id*/, NeighList *ptr)
{
list = ptr;
}
/* ---------------------------------------------------------------------- */
void ComputeSnapneigh::compute_array()
{
if (dbirjflag){
//printf("----- dbirjflag true.\n");
get_dbirj_length();
//printf("----- got dbirj_length\n");
}
//printf("----- cutmax cutforce: %f %f\n", cutmax, force->pair->cutforce);
//else{
int ntotal = atom->nlocal + atom->nghost;
invoked_array = update->ntimestep;
// 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;
//printf("----- inum: %d\n", inum);
//printf("----- NEIGHMASK: %d\n", NEIGHMASK);
int ninside;
int numneigh_sum = 0;
int dbirj_row_indx;
int dbiri_indx=0;
for (int ii = 0; ii < inum; ii++) {
int irow = 0;
if (bikflag) irow = atom->tag[ilist[ii] & NEIGHMASK]-1;
//printf("----- i, itag: %d %d\n", ilist[ii] & NEIGHMASK, atom->tag[ilist[ii]]);
const int i = ilist[ii];
//printf("----- ii, i: %d %d\n", ii, i);
//printf("----- mask[i] groupbit: %d %d\n", mask[i], groupbit);
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];
int ielem = 0;
if (chemflag)
ielem = map[itype];
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
/*
This loop assigns quantities in snaptr.
snaptr is a SNA class instance, see sna.h
*/
//int ninside = 0;
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];
int jelem = 0;
if (chemflag)
jelem = map[jtype];
if (rsq < cutsq[itype][jtype]&&rsq>1e-20) {
if (dbirjflag){
//dbirj_row_indx = 3*neighsum[atom->tag[j]-1] + 3*icounter[atom->tag[j]-1] ; // THIS IS WRONG, SEE NEXT VAR.
//dbirj_row_indx = 3*neighsum[atom->tag[i]-1] + 3*(atom->tag[i]-1) + 3*icounter[atom->tag[j]-1]; // 3*tagi is to leave space for dBi/dRi
dbirj_row_indx = 3*neighsum[atom->tag[j]-1] + 3*icounter[atom->tag[j]-1] + 3*(atom->tag[j]-1); // THIS IS WRONG, SEE NEXT VAR.
//printf("--- %d %d %d %d\n", dbirj_row_indx, 3*neighsum[atom->tag[i]-1], 3*(atom->tag[i]-1), jj);
//printf("jtag, icounter, dbirj_row_indx: %d, %d, %d %d %d\n", atom->tag[j], icounter[atom->tag[j]-1], dbirj_row_indx+0, dbirj_row_indx+1, dbirj_row_indx+2);
icounter[atom->tag[j]-1] += 1;
neighs[dbirj_row_indx+0][0] = atom->tag[i];
neighs[dbirj_row_indx+1][0] = atom->tag[i];
neighs[dbirj_row_indx+2][0] = atom->tag[i];
neighs[dbirj_row_indx+0][1] = atom->tag[j];
neighs[dbirj_row_indx+1][1] = atom->tag[j];
neighs[dbirj_row_indx+2][1] = atom->tag[j];
neighs[dbirj_row_indx+0][2] = 0;
neighs[dbirj_row_indx+1][2] = 1;
neighs[dbirj_row_indx+2][2] = 2;
dbiri_indx = dbiri_indx+3;
}
}
}
//printf("--- dbiri_indx: %d\n", dbiri_indx);
// Put dBi/dRi in
neighs[dbiri_indx+0][0] = atom->tag[i];
neighs[dbiri_indx+1][0] = atom->tag[i];
neighs[dbiri_indx+2][0] = atom->tag[i];
neighs[dbiri_indx+0][1] = atom->tag[i];
neighs[dbiri_indx+1][1] = atom->tag[i];
neighs[dbiri_indx+2][1] = atom->tag[i];
neighs[dbiri_indx+0][2] = 0;
neighs[dbiri_indx+1][2] = 1;
neighs[dbiri_indx+2][2] = 2;
dbiri_indx = dbiri_indx+3;
}
numneigh_sum += ninside;
} // for (int ii = 0; ii < inum; ii++)
// Check icounter.
/*
for (int ii = 0; ii < inum; ii++) {
const int i = ilist[ii];
if (mask[i] & groupbit) {
printf("icounter[i]: %d\n", icounter[i]);
}
}
*/
// sum up over all processes
// I'll need to do something like this...
/*
MPI_Allreduce(&snap[0][0],&snapall[0][0],size_array_rows*size_array_cols,MPI_DOUBLE,MPI_SUM,world);
// assign energy to last column
for (int i = 0; i < bik_rows; i++) snapall[i][lastcol] = 0;
int irow = 0;
double reference_energy = c_pe->compute_scalar();
snapall[irow][lastcol] = reference_energy;
*/
/*
for (int i=0; i<dbirj_rows; i++){
printf("----- %d: %f %f\n", i, array[i][0], array[i][1]);
}
//printf("vector[0]: %d\n", vector[0]);
printf("----- End of compute snapneigh.\n");
*/
}
/* ----------------------------------------------------------------------
compute dbirj length
------------------------------------------------------------------------- */
void ComputeSnapneigh::get_dbirj_length()
{
// invoke full neighbor list (will copy or build if necessary)
neighbor->build_one(list);
dbirj_rows = 0;
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;
const int* const mask = atom->mask;
double** const x = atom->x;
//printf("----- inum: %d\n", inum);
memory->create(neighsum, atom->nlocal, "snapneigh:neighsum");
memory->create(nneighs, atom->nlocal, "snapneigh:nneighs");
memory->create(icounter, atom->nlocal, "snapneigh:icounter");
//memory->create(dbiri, 3*inum,ncoeff, "snapneigh:dbiri");
/*
for (int ii=0; ii<inum; ii++){
for (int icoeff=0; icoeff<ncoeff; icoeff++){
dbiri[ii][icoeff]=0.0;
}
}
*/
for (int ii = 0; ii < inum; ii++) {
const int i = ilist[ii];
if (mask[i] & groupbit) {
icounter[i]=0;
neighsum[i] = 0;
nneighs[i] = 0;
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* const jlist = firstneigh[i];
const int jnum = numneigh[i];
//printf("----- jnum: %d\n", jnum);
int jnum_cutoff = 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) {
dbirj_rows += 1; //jnum + 1;
jnum_cutoff += 1;
nneighs[i]+=1;
}
}
//printf("----- jnum_cutoff: %d\n", jnum_cutoff);
}
}
dbirj_rows *= ndims_force;
// Loop over all atoms again to calculate neighsum.
for (int ii = 0; ii < inum; ii++) {
const int i = ilist[ii];
if (mask[i] & groupbit) {
//printf("nneighs[i]: %d\n", nneighs[i]);
//neighsum[i] = 0;
//printf("i nneighs[i]: %d %d\n", i, nneighs[i]);
if (i==0){
neighsum[i]=0;
}
else{
for (int jj=0; jj < ii; jj++){
const int j = ilist[jj];
if (mask[j] & groupbit) {
//printf(" j nneighs[j-1]: %d %d\n", j, nneighs[j]);
neighsum[i] += nneighs[j];
}
}
//neighsum[i] += 1; // Add 1 for the self term dBi/dRi
}
}
//printf("%d\n", neighsum[i]);
}
size_array_rows = dbirj_rows+(3*atom->nlocal);
size_array_cols = 3;
//memory->create(dbirj, dbirj_rows, ncoeff, "snapneigh:dbirj");
memory->create(neighs, size_array_rows, size_array_cols, "snapneigh:neighs");
//vector = neighs;
array = neighs;
// Set size array rows which now depends on dbirj_rows.
//size_array_rows = bik_rows+dbirj_rows+ndims_virial;
//printf("----- dbirj_rows: %d\n", dbirj_rows);
//printf("----- end of dbirj length.\n");
}
/* ----------------------------------------------------------------------
compute array length
------------------------------------------------------------------------- */
double ComputeSnapneigh::compute_scalar()
{
if (dbirjflag) get_dbirj_length();
return size_array_rows;
}
/* ----------------------------------------------------------------------
memory usage
------------------------------------------------------------------------- */
double ComputeSnapneigh::memory_usage()
{
double bytes = (double)size_array_rows*size_array_cols *
sizeof(double); // array
return bytes;
}

73
src/ML-SNAP/compute_snapneigh.h
View File

@ -1,73 +0,0 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef COMPUTE_CLASS
// clang-format off
ComputeStyle(snapneigh,ComputeSnapneigh);
// clang-format on
#else
#ifndef LMP_COMPUTE_SNAPNEIGH_H
#define LMP_COMPUTE_SNAPNEIGH_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeSnapneigh : public Compute {
public:
ComputeSnapneigh(class LAMMPS *, int, char **);
~ComputeSnapneigh() override;
void init() override;
void init_list(int, class NeighList *) override;
void compute_array() override;
double compute_scalar() override;
double memory_usage() override;
private:
FILE * fh_d;
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;
int *map; // map types to [0,nelements)
int nelements, chemflag;
int switchinnerflag;
double *sinnerelem;
double *dinnerelem;
//class SNA *snaptr;
double cutmax;
int quadraticflag;
//int bikflag;
//int bik_rows;
int bikflag, bik_rows, dbirjflag, dbirj_rows;
double **dbirj;
double **dbiri; // dBi/dRi = sum(-dBi/dRj) over neighbors j
int *nneighs; // number of neighs inside the snap cutoff.
int *neighsum;
int *icounter; // counting atoms i for each j.
double **neighs; // neighborlist for neural networks
void get_dbirj_length();
};
} // namespace LAMMPS_NS
#endif
#endif