diff --git a/doc/src/compute_fabric.rst b/doc/src/compute_fabric.rst new file mode 100644 index 0000000000..1da39c248a --- /dev/null +++ b/doc/src/compute_fabric.rst @@ -0,0 +1,149 @@ +.. index:: compute fabric + +compute fabric command +============================== + +Syntax +"""""" + +.. parsed-literal:: + + compute ID group-ID fabric cutoff attribute1 attribute2 ... + +* ID, group-ID are documented in :doc:`compute ` command +* fabric = style name of this compute command +* cutoff = *type* or *radius* + + .. parsed-literal:: + + *type* = cutoffs determined for the types of the two atoms + *radius* = cutoffs determined based on finite size of atoms + +* one or more attributes may be appended + + .. parsed-literal:: + + *contact* = contact tensor + *branch* = branch tensor + *force/normal* = normal force tensor + *force/tangential* = tangential force tensor + +Examples +"""""""" + +.. code-block:: LAMMPS + + compute 1 all fabric type contact force/normal + compute 1 all fabric radius force/normal force/tangential + +Description +""""""""""" + +Define a computation that calculates various fabric tensors for pairwise +interactions. The *type* and *radius* settings are used to select whether +the compute will be used used with regular interactions with cutoffs +determined by atom types or with granular interactions with interaction +lengths determined by particle radii, respectively. + +Four fabric tensors can be computed: a contact, branch, normal force, or +tangential force tensor. The contact tensor is calculated as + +.. math:: + + C_{ab} = \frac{15}{2} (\phi_{ab} - \mathrm{tr}(\phi) \delta_{ab}) + +where :math:`a` and :math:`b` are the :math:`x`, :math:`y`, :math:`z` +directions, :math:`\delta_{ab}` is the Kronecker delta function, and +the tensor :math:`\phi` is defined as + +.. math:: + + \phi_{ab} = \sum_{n = 1}^{N_p} \frac{r_{a} r_{b}}{r^2} + +where :math:`N_p` is the number of pair interactions in the simulation, +:math:`r_{a}` is the :math:`a` component of the radial vector between the +two particles, and :math:`r` is the magnitude of the radial vector. + +The branch tensor is calculated as + +.. math:: + + B_{ab} = \frac{15}{6 tr(D)} (D_{ab} - tr(D) \delta_{ab}) + +where the tensor :math:`D` is defined as + +.. math:: + + D_{ab} = \sum_{n = 1}^{N_p} + \frac{1}{N_c (r^2 + C_{cd} r_c r_d)} + \frac{r_{a} r_{b}}{r} + +where :math:`N_c` is the total number of contacts in the system and the subscripts +:math:`c` and :math:`d` are summed according to Einstein notation. + +The normal force fabric tensor is calculated as + +.. math:: + + F^n_{ab} = \frac{15}{6 tr(N)} (N_{ab} - tr(N) \delta_{ab}) + +where the tensor :math:`N` is defined as + +.. math:: + + N_{ab} = \sum_{n = 1}^{N_p} + \frac{1}{N_c (r^2 + C_{cd} r_c r_d)} + \frac{r_{a} r_{b}}{r^2} f_n + +where :math:`f_n` is the magnitude of the normal, central-body force between the two atoms. + +Finally, the tangential force fabric tensor is only defined for pair styles that +apply tangential forces to particles, namely granular pair styles. It is calculated +as + +.. math:: + + F^t_{ab} = \frac{15}{9 tr(N)} (T_{ab} - tr(T) \delta_{ab}) + +where the tensor :math:`T` is defined as + +.. math:: + + T_{ab} = \sum_{n = 1}^{N_p} + \frac{1}{N_c (r^2 + C_{cd} r_c r_d)} + \frac{r_{a} r_{b}}{r^2} f_t + +where :math:`f_t` is the magnitude of the tagential force between the two atoms. + + +Output info +""""""""""" + +This compute calculates a local vector of doubles. The vector stores the +unique components of the first requested tensor in the order +xx, yy, zz, xy, xz, yz followed by the same components for all +subsequent tensors. The length of the vector is therefore six times +the number of requested tensors. + +Restrictions +"""""""""""" + +Currently, compute *fabric* does not support pair styles +with many-body interactions. It also does not +support models with long-range Coulombic or dispersion forces, +i.e. the kspace_style command in LAMMPS. It also does not support the +following fixes which add rigid-body constraints: :doc:`fix shake +`, :doc:`fix rattle `, :doc:`fix rigid +`, :doc:`fix rigid/small `. It does not support +granular pair styles that extend beyond the contact (e.g. JKR and DMT). + + +Related commands +"""""""""""""""" + +none + +Default +""""""" + +none diff --git a/src/.gitignore b/src/.gitignore index 72bb8a3004..97229711cb 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -379,6 +379,8 @@ /compute_erotate_rigid.h /compute_event_displace.cpp /compute_event_displace.h +/compute_fabric.cpp +/compute_fabric.h /compute_fep.cpp /compute_fep.h /compute_force_tally.cpp diff --git a/src/GRANULAR/compute_fabric.cpp b/src/GRANULAR/compute_fabric.cpp new file mode 100644 index 0000000000..d08ce0087b --- /dev/null +++ b/src/GRANULAR/compute_fabric.cpp @@ -0,0 +1,407 @@ +/* ---------------------------------------------------------------------- + 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_fabric.h" +#include +#include +#include +#include "atom.h" +#include "update.h" +#include "force.h" +#include "pair.h" +#include "modify.h" +#include "neighbor.h" +#include "neigh_request.h" +#include "neigh_list.h" +#include "error.h" + +using namespace LAMMPS_NS; + +enum{OTHER,GRANULAR}; +enum{TYPE,RADIUS}; +enum{CN,BR,FN,FT}; + +/* ---------------------------------------------------------------------- */ + +ComputeFabric::ComputeFabric(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg), + tensor_style(NULL) +{ + if (narg < 3) error->all(FLERR,"Illegal compute fabric command"); + + if (strcmp(arg[3],"type") == 0) cutstyle = TYPE; + else if (strcmp(arg[3],"radius") == 0) cutstyle = RADIUS; + else error->all(FLERR,"Illegal compute fabric command"); + + if (cutstyle == RADIUS && !atom->radius_flag) + error->all(FLERR,"Compute fabric radius style requires atom attribute radius"); + + ntensors = narg - 4; + tensor_style = new int[ntensors]; + + cn_flag = 0; + br_flag = 0; + fn_flag = 0; + ft_flag = 0; + + ntensors = 0; + int iarg = 4; + while (iarg < narg) { + if (strcmp(arg[iarg],"contact") == 0) { + cn_flag = 1; + tensor_style[ntensors++] = CN; + } else if (strcmp(arg[iarg],"branch") == 0) { + br_flag = 1; + tensor_style[ntensors++] = BR; + } else if (strcmp(arg[iarg],"force/normal") == 0) { + fn_flag = 1; + tensor_style[ntensors++] = FN; + } else if (strcmp(arg[iarg],"force/tangential") == 0) { + ft_flag = 1; + tensor_style[ntensors++] = FT; + } else error->all(FLERR,"Illegal compute fabric command"); + iarg++; + } + + vector_flag = 1; + size_vector = 1+ntensors*6; + extvector = 0; + + vector = new double[size_vector]; +} + +/* ---------------------------------------------------------------------- */ + +ComputeFabric::~ComputeFabric() +{ + delete [] vector; + delete [] tensor_style; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeFabric::init() +{ + if (force->pair == NULL) + error->all(FLERR,"No pair style is defined for compute fabric"); + if (force->pair->single_enable == 0 && (fn_flag || ft_flag)) + error->all(FLERR,"Pair style does not support compute fabric normal or tangential force"); + + // Find if granular or gran + pstyle = OTHER; + if (force->pair_match("granular",0) || + force->pair_match("gran/hooke",0) || + force->pair_match("gran/hertz",0) || + force->pair_match("gran/hooke/history",0) || + force->pair_match("gran/hertz/history",0)) pstyle = GRANULAR; + + if (pstyle != GRANULAR && ft_flag) + error->all(FLERR,"Pair style does not have tangential forces for compute fabric"); + + if(force->pair->beyond_contact) + error->all(FLERR, "Compute fabric does not support pair styles that extend beyond contact"); + + // need an occasional half neighbor list + // set size to same value as request made by force->pair + // this should enable it to always be a copy list (e.g. for granular pstyle) + + int irequest = neighbor->request(this,instance_me); + neighbor->requests[irequest]->pair = 0; + neighbor->requests[irequest]->compute = 1; + neighbor->requests[irequest]->occasional = 1; + NeighRequest *pairrequest = neighbor->find_request((void *) force->pair); + if (pairrequest) neighbor->requests[irequest]->size = pairrequest->size; + +} + +/* ---------------------------------------------------------------------- */ + +void ComputeFabric::init_list(int /*id*/, NeighList *ptr) +{ + list = ptr; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeFabric::compute_vector() +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + tagint itag,jtag; + double xtmp,ytmp,ztmp,delx,dely,delz; + double r,rinv,rsq,radsum,eng,fpair; + + double nx, ny, nz, fx, fy ,fz; + double ncinv, denom, fn, ft, prefactor; + double cn_tensor[6], br_tensor[6], ft_tensor[6], fn_tensor[6]; + double trace_phi, trace_D, trace_Xfn, trace_Xft; + double phi_ij[6] = {0.0}; + double Ac_ij[6] = {0.0}; + double D_ij[6] = {0.0}; + double Xfn_ij[6] = {0.0}; + double Xft_ij[6] = {0.0}; + int nc; + int temp_int; + double temp_dbl[6]; + + int *ilist,*jlist,*numneigh,**firstneigh; + + double **x = atom->x; + double *radius = atom->radius; + tagint *tag = atom->tag; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + int newton_pair = force->newton_pair; + + // invoke half neighbor list (will copy or build if necessary) + neighbor->build_one(list); + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + Pair *pair = force->pair; + double **cutsq = force->pair->cutsq; + + // First loop through and calculate contact fabric tensor + nc = 0; + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + if (!(mask[i] & groupbit)) continue; + + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itag = tag[i]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + if (!(mask[j] & groupbit)) continue; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + jtype = type[j]; + + if (cutstyle == TYPE) { + if (rsq >= cutsq[itype][jtype]) continue; + } else { + radsum = radius[i] + radius[j]; + if (rsq >= radsum*radsum) continue; + } + + nc += 1; + + r = sqrt(rsq); + rinv = 1.0/r; + nx = delx*rinv; + ny = dely*rinv; + nz = delz*rinv; + + + phi_ij[0] += nx*nx; + phi_ij[1] += ny*ny; + phi_ij[2] += nz*nz; + phi_ij[3] += nx*ny; + phi_ij[4] += nx*nz; + phi_ij[5] += ny*nz; + } + } + + //Count total contacts across processors + MPI_Allreduce(&nc,&temp_int,1,MPI_INT,MPI_SUM,world); + nc = temp_int; + + // If no contacts, everything will be zero + if (nc == 0) { + for (i = 0; i < size_vector; i ++) vector[i] = 0.0; + return; + } + + vector[ntensors*6] = nc; + ncinv = 1.0/nc; + + //Sum phi across processors + MPI_Allreduce(phi_ij,temp_dbl,6,MPI_DOUBLE,MPI_SUM,world); + for (i = 0; i < 6; i++) phi_ij[i] = temp_dbl[i]*ncinv; + + trace_phi = (1.0/3.0)*(phi_ij[0]+phi_ij[1]+phi_ij[2]); + + Ac_ij[0] = (15.0/2.0)*(phi_ij[0] - trace_phi); + Ac_ij[1] = (15.0/2.0)*(phi_ij[1] - trace_phi); + Ac_ij[2] = (15.0/2.0)*(phi_ij[2] - trace_phi); + Ac_ij[3] = (15.0/2.0)*(phi_ij[3]); + Ac_ij[4] = (15.0/2.0)*(phi_ij[4]); + Ac_ij[5] = (15.0/2.0)*(phi_ij[5]); + + // If needed, loop through and calculate other fabric tensors + if(br_flag || fn_flag || ft_flag) { + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + if (!(mask[i] & groupbit)) continue; + + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itag = tag[i]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + if (!(mask[j] & groupbit)) continue; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + jtype = type[j]; + + if (cutstyle == TYPE) { + if (rsq >= cutsq[itype][jtype]) continue; + } else { + radsum = radius[i] + radius[j]; + if (rsq >= radsum*radsum) continue; + } + + if (fn_flag || ft_flag) + eng = pair->single(i,j,itype,jtype,rsq,1.0,1.0,fpair); + + r = sqrt(rsq); + rinv = 1.0/r; + nx = delx*rinv; + ny = dely*rinv; + nz = delz*rinv; + + denom = 1 + Ac_ij[0]*nx*nx + Ac_ij[1]*ny*ny + Ac_ij[2]*nz*nz; + denom += 2*Ac_ij[3]*nx*ny + 2*Ac_ij[4]*nx*nz + 2*Ac_ij[5]*ny*nz; + prefactor = ncinv/denom; + + if (br_flag) { + D_ij[0] += prefactor*nx*nx*r; + D_ij[1] += prefactor*ny*ny*r; + D_ij[2] += prefactor*nz*nz*r; + D_ij[3] += prefactor*nx*ny*r; + D_ij[4] += prefactor*nx*nz*r; + D_ij[5] += prefactor*ny*nz*r; + } + + if (fn_flag || ft_flag) { + fn = r*fpair; + + Xfn_ij[0] += prefactor*nx*nx*fn; + Xfn_ij[1] += prefactor*ny*ny*fn; + Xfn_ij[2] += prefactor*nz*nz*fn; + Xfn_ij[3] += prefactor*nx*ny*fn; + Xfn_ij[4] += prefactor*nx*nz*fn; + Xfn_ij[5] += prefactor*ny*nz*fn; + + if (ft_flag) { + ft = force->pair->svector[3]; + + Xft_ij[0] += prefactor*nx*nx*ft; + Xft_ij[1] += prefactor*ny*ny*ft; + Xft_ij[2] += prefactor*nz*nz*ft; + Xft_ij[3] += prefactor*nx*ny*ft; + Xft_ij[4] += prefactor*nx*nz*ft; + Xft_ij[5] += prefactor*ny*nz*ft; + } + } + } + } + } + + // Output results + + if (cn_flag) { + for (i = 0; i < ntensors; i ++) { + if (tensor_style[i] == CN) { + for (j = 0; j < 6; j++) vector[6*i + j] = Ac_ij[j]; + } + } + } + + if (br_flag) { + MPI_Allreduce(D_ij,temp_dbl,6,MPI_DOUBLE,MPI_SUM,world); + for (i = 0; i < 6; i++) D_ij[i] = temp_dbl[i]; + + trace_D = (1.0/3.0)*(D_ij[0]+D_ij[1]+D_ij[2]); + + br_tensor[0] = (15.0/(6.0*trace_D))*(D_ij[0] - trace_D); + br_tensor[1] = (15.0/(6.0*trace_D))*(D_ij[1] - trace_D); + br_tensor[2] = (15.0/(6.0*trace_D))*(D_ij[2] - trace_D); + br_tensor[3] = (15.0/(6.0*trace_D))*(D_ij[3]); + br_tensor[4] = (15.0/(6.0*trace_D))*(D_ij[4]); + br_tensor[5] = (15.0/(6.0*trace_D))*(D_ij[5]); + + for (i = 0; i < ntensors; i ++) { + if (tensor_style[i] == BR) { + for (j = 0; j < 6; j++) vector[6*i + j] = br_tensor[j]; + } + } + } + + if (fn_flag || ft_flag) { + MPI_Allreduce(Xfn_ij,temp_dbl,6,MPI_DOUBLE,MPI_SUM,world); + for (i = 0; i < 6; i++) Xfn_ij[i] = temp_dbl[i]; + + trace_Xfn = (1.0/3.0)*(Xfn_ij[0]+Xfn_ij[1]+Xfn_ij[2]); + } + + if (fn_flag) { + + fn_tensor[0] = (15.0/(6.0*trace_Xfn))*(Xfn_ij[0] - trace_Xfn); + fn_tensor[1] = (15.0/(6.0*trace_Xfn))*(Xfn_ij[1] - trace_Xfn); + fn_tensor[2] = (15.0/(6.0*trace_Xfn))*(Xfn_ij[2] - trace_Xfn); + fn_tensor[3] = (15.0/(6.0*trace_Xfn))*(Xfn_ij[3]); + fn_tensor[4] = (15.0/(6.0*trace_Xfn))*(Xfn_ij[4]); + fn_tensor[5] = (15.0/(6.0*trace_Xfn))*(Xfn_ij[5]); + + for (i = 0; i < ntensors; i ++) { + if (tensor_style[i] == FN) { + for (j = 0; j < 6; j++) vector[6*i + j] = fn_tensor[j]; + } + } + } + + if (ft_flag) { + MPI_Allreduce(Xft_ij,temp_dbl,6,MPI_DOUBLE,MPI_SUM,world); + for (i = 0; i < 6; i++) Xft_ij[i] = temp_dbl[i]; + + trace_Xft = (1.0/3.0)*(Xft_ij[0]+Xft_ij[1]+Xft_ij[2]); + + ft_tensor[0] = (15.0/(9.0*trace_Xfn))*(Xft_ij[0] - trace_Xft); + ft_tensor[1] = (15.0/(9.0*trace_Xfn))*(Xft_ij[1] - trace_Xft); + ft_tensor[2] = (15.0/(9.0*trace_Xfn))*(Xft_ij[2] - trace_Xft); + ft_tensor[3] = (15.0/(9.0*trace_Xfn))*(Xft_ij[3]); + ft_tensor[4] = (15.0/(9.0*trace_Xfn))*(Xft_ij[4]); + ft_tensor[5] = (15.0/(9.0*trace_Xfn))*(Xft_ij[5]); + + for (i = 0; i < ntensors; i ++) { + if (tensor_style[i] == FT) { + for (j = 0; j < 6; j++) vector[6*i + j] = ft_tensor[j]; + } + } + } +} diff --git a/src/GRANULAR/compute_fabric.h b/src/GRANULAR/compute_fabric.h new file mode 100644 index 0000000000..a65fc50c61 --- /dev/null +++ b/src/GRANULAR/compute_fabric.h @@ -0,0 +1,76 @@ +/* -*- 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(fabric,ComputeFabric) + +#else + +#ifndef LMP_COMPUTE_FABRIC_H +#define LMP_COMPUTE_FABRIC_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeFabric : public Compute { + public: + ComputeFabric(class LAMMPS *, int, char **); + ~ComputeFabric(); + void init(); + void init_list(int, class NeighList *); + void compute_vector(); + + private: + int ntensors, pstyle, cutstyle; + int *tensor_style; + class NeighList *list; + + int cn_flag, br_flag, fn_flag, ft_flag; +}; + +} + +#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 fabric radius style requires atom attribute radius + +Self-explanatory. + +E: No pair style is defined for compute fabric + +Self-explanatory. + +E: Pair style does not support compute fabric normal or tangential force + +Pair style must be single enabled to calculate the normal or tangential force tensors + +E: Pair style does not have tangential forces for compute fabric + +The tangential force tensor can only be calculated for granular pair styles + +E: Compute fabric does not support pair styles that extend beyond contact + +Granular pair styles that extend beyond contact such as JKR or DMT are not supported + +*/