diff --git a/src/SNAP/compute_snap.cpp b/src/SNAP/compute_snap.cpp new file mode 100644 index 0000000000..31b50e1902 --- /dev/null +++ b/src/SNAP/compute_snap.cpp @@ -0,0 +1,377 @@ +/* ---------------------------------------------------------------------- + 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 "compute_snap.h" +#include +#include +#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; + +ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg), cutsq(NULL), list(NULL), snap(NULL), + radelem(NULL), wjelem(NULL) +{ + 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; + yoffset = nperdim; + zoffset = 2*nperdim; + size_peratom_cols = 3*nperdim*atom->ntypes; + comm_reverse = size_peratom_cols; + + nmax = 0; + snap = NULL; +} + +/* ---------------------------------------------------------------------- */ + +ComputeSnap::~ComputeSnap() +{ + memory->destroy(snap); + 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(); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSnap::init_list(int /*id*/, NeighList *ptr) +{ + list = ptr; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSnap::compute() +{ + int ntotal = atom->nlocal + atom->nghost; + + invoked_peratom = update->ntimestep; + + // grow snap array if necessary + + if (atom->nmax > nmax) { + memory->destroy(snap); + nmax = atom->nmax; + memory->create(snap,nmax,size_peratom_cols, + "snap:snap"); + array = snap; + } + + // clear local array + + for (int i = 0; i < ntotal; i++) + for (int icoeff = 0; icoeff < size_peratom_cols; icoeff++) { + snap[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 = threencoeff*(atom->type[i]-1); + // const int quadraticoffset = threencoeff*atom->ntypes + + // threencoeffq*(atom->type[i]-1); + const int typeoffset = 3*nperdim*(atom->type[i]-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(); + if (quadraticflag) { + 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 *snapi = snap[i]+typeoffset; + double *snapj = snap[j]+typeoffset; + + for (int icoeff = 0; icoeff < ncoeff; icoeff++) { + snapi[icoeff] += snaptr->dblist[icoeff][0]; + snapi[icoeff+yoffset] += snaptr->dblist[icoeff][1]; + snapi[icoeff+zoffset] += snaptr->dblist[icoeff][2]; + snapj[icoeff] -= snaptr->dblist[icoeff][0]; + snapj[icoeff+yoffset] -= snaptr->dblist[icoeff][1]; + snapj[icoeff+zoffset] -= snaptr->dblist[icoeff][2]; + } + + if (quadraticflag) { + const int quadraticoffset = ncoeff; + snapi += quadraticoffset; + snapj += 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; + + snapi[ncount] += dbxtmp; + snapi[ncount+yoffset] += dbytmp; + snapi[ncount+zoffset] += dbztmp; + snapj[ncount] -= dbxtmp; + snapj[ncount+yoffset] -= dbytmp; + snapj[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]; + + snapi[ncount] += dbxtmp; + snapi[ncount+yoffset] += dbytmp; + snapi[ncount+zoffset] += dbztmp; + snapj[ncount] -= dbxtmp; + snapj[ncount+yoffset] -= dbytmp; + snapj[ncount+zoffset] -= dbztmp; + ncount++; + } + } + } + } + } + } + + // communicate snap contributions between neighbor procs + + comm->reverse_comm_compute(this); + +} + +/* ---------------------------------------------------------------------- */ + +int ComputeSnap::pack_reverse_comm(int n, int first, double *buf) +{ + int i,m,last,icoeff; + + m = 0; + last = first + n; + for (i = first; i < last; i++) + for (icoeff = 0; icoeff < size_peratom_cols; icoeff++) + buf[m++] = snap[i][icoeff]; + return m; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSnap::unpack_reverse_comm(int n, int *list, double *buf) +{ + int i,j,m,icoeff; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + for (icoeff = 0; icoeff < size_peratom_cols; icoeff++) + snap[j][icoeff] += buf[m++]; + } +} + +/* ---------------------------------------------------------------------- + memory usage +------------------------------------------------------------------------- */ + +double ComputeSnap::memory_usage() +{ + + double bytes = nmax*size_peratom_cols * sizeof(double); // snap + bytes += snaptr->memory_usage(); // SNA object + + return bytes; +} diff --git a/src/SNAP/compute_snap.h b/src/SNAP/compute_snap.h new file mode 100644 index 0000000000..9fa998981a --- /dev/null +++ b/src/SNAP/compute_snap.h @@ -0,0 +1,77 @@ +/* -*- 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(); + int pack_reverse_comm(int, int, double *); + void unpack_reverse_comm(int, int *, double *); + double memory_usage(); + + private: + int nmax; + int ncoeff, nperdim, yoffset, zoffset; + double **cutsq; + class NeighList *list; + double **snap; + double rcutfac; + double *radelem; + double *wjelem; + class SNA* snaptr; + double cutmax; + int quadraticflag; +}; + +} + +#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. + +*/