diff --git a/doc/src/compute_sna_atom.rst b/doc/src/compute_sna_atom.rst index f156133521..7573a2031a 100644 --- a/doc/src/compute_sna_atom.rst +++ b/doc/src/compute_sna_atom.rst @@ -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* +* keyword = *rmin0* or *switchflag* or *bzeroflag* or *quadraticflag* or *chem* or *bnormflag* or *wselfallflag* or *bikflag* .. parsed-literal:: @@ -56,6 +56,9 @@ Syntax *wselfallflag* value = *0* or *1* *0* = self-contribution only for element of central atom *1* = self-contribution for all elements + *bikflag* value = *0* or *1* (only implemented for compute snap) + *0* = per-atom bispectrum descriptors are summed over atoms + *1* = per-atom bispectrum descriptors are not summed over atoms Examples """""""" @@ -296,6 +299,15 @@ This option is typically used in conjunction with the *chem* keyword, and LAMMPS will generate a warning if both *chem* and *bnormflag* are not both set or not both unset. +The keyword *bikflag* determines whether or not to expand the bispectrum +rows of the global array returned by compute snap. If *bikflag* is set +to *1* then the bispectrum row, which is typically the per-atom bispectrum +descriptors :math:`B_{i,k}` summed over all atoms *i* to produce +:math:`B_k`, becomes bispectrum rows equal to the number of atoms. Thus, +the resulting bispectrum rows are :math:`B_{i,k}` instead of just +:math:`B_k`. In this case, the entries in the final column for these rows +are set to zero. + .. note:: If you have a bonded system, then the settings of :doc:`special_bonds diff --git a/src/ML-SNAP/compute_snap.cpp b/src/ML-SNAP/compute_snap.cpp index 63deff9f8f..4d24616aa5 100644 --- a/src/ML-SNAP/compute_snap.cpp +++ b/src/ML-SNAP/compute_snap.cpp @@ -57,6 +57,7 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) : switchflag = 1; bzeroflag = 1; quadraticflag = 0; + bikflag = 0; chemflag = 0; bnormflag = 0; wselfallflag = 0; @@ -137,6 +138,11 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) : 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 error->all(FLERR,"Illegal compute snap command"); } @@ -152,7 +158,9 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) : yoffset = nperdim; zoffset = 2*nperdim; natoms = atom->natoms; - size_array_rows = 1+ndims_force*natoms+ndims_virial; + bik_rows = 1; + if (bikflag) bik_rows = natoms; + size_array_rows = bik_rows+ndims_force*natoms+ndims_virial; size_array_cols = nperdim*atom->ntypes+1; lastcol = size_array_cols-1; @@ -287,6 +295,8 @@ void ComputeSnap::compute_array() const int* const mask = atom->mask; for (int ii = 0; ii < inum; ii++) { + int irow = 0; + if (bikflag) irow = atom->tag[ilist[ii] & NEIGHMASK]-1; const int i = ilist[ii]; if (mask[i] & groupbit) { @@ -417,17 +427,17 @@ void ComputeSnap::compute_array() int k = typeoffset_global; for (int icoeff = 0; icoeff < ncoeff; icoeff++) - snap[0][k++] += snaptr->blist[icoeff]; + snap[irow][k++] += snaptr->blist[icoeff]; // quadratic contributions if (quadraticflag) { for (int icoeff = 0; icoeff < ncoeff; icoeff++) { double bveci = snaptr->blist[icoeff]; - snap[0][k++] += 0.5*bveci*bveci; + snap[irow][k++] += 0.5*bveci*bveci; for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) { double bvecj = snaptr->blist[jcoeff]; - snap[0][k++] += bveci*bvecj; + snap[irow][k++] += bveci*bvecj; } } } @@ -443,10 +453,10 @@ void ComputeSnap::compute_array() 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]; + int irow = 3*(iglobal-1)+bik_rows; + snap[irow++][icoeff+typeoffset_global] += snadi[icoeff]; + snap[irow++][icoeff+typeoffset_global] += snadi[icoeff+yoffset]; + snap[irow][icoeff+typeoffset_global] += snadi[icoeff+zoffset]; } } } @@ -455,10 +465,10 @@ void ComputeSnap::compute_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]; + int irow = 3*(iglobal-1)+bik_rows; + snap[irow++][lastcol] = atom->f[i][0]; + snap[irow++][lastcol] = atom->f[i][1]; + snap[irow][lastcol] = atom->f[i][2]; } // accumulate bispectrum virial contributions to global array @@ -470,22 +480,22 @@ void ComputeSnap::compute_array() 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; + snapall[irow][lastcol] = reference_energy; // assign virial stress to last column // switch to Voigt notation c_virial->compute_vector(); - irow += 3*natoms; + irow += 3*natoms+bik_rows; 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]; + snapall[irow][lastcol] = c_virial->vector[3]; } @@ -497,7 +507,7 @@ void ComputeSnap::compute_array() void ComputeSnap::dbdotr_compute() { double **x = atom->x; - int irow0 = 1+ndims_force*natoms; + int irow0 = bik_rows+ndims_force*natoms; // sum over bispectrum contributions to forces // on all particles including ghosts @@ -518,7 +528,7 @@ void ComputeSnap::dbdotr_compute() 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]; + snap[irow][icoeff+typeoffset_global] += dbdy*x[i][0]; } } } diff --git a/src/ML-SNAP/compute_snap.h b/src/ML-SNAP/compute_snap.h index f092ede729..ccdc0e107e 100644 --- a/src/ML-SNAP/compute_snap.h +++ b/src/ML-SNAP/compute_snap.h @@ -49,6 +49,8 @@ class ComputeSnap : public Compute { class SNA *snaptr; double cutmax; int quadraticflag; + int bikflag; + int bik_rows; Compute *c_pe; Compute *c_virial;