diff --git a/src/PERI/Install.sh b/src/PERI/Install.sh index 766cd8758c..4d783b0a64 100644 --- a/src/PERI/Install.sh +++ b/src/PERI/Install.sh @@ -4,11 +4,13 @@ if (test $1 = 1) then cp atom_vec_peri.cpp .. cp pair_peri_pmb.cpp .. + cp pair_peri_lps.cpp .. cp fix_peri_neigh.cpp .. cp compute_damage_atom.cpp .. cp atom_vec_peri.h .. cp pair_peri_pmb.h .. + cp pair_peri_lps.h .. cp fix_peri_neigh.h .. cp compute_damage_atom.h .. @@ -16,11 +18,13 @@ elif (test $1 = 0) then rm -f ../atom_vec_peri.cpp rm -f ../pair_peri_pmb.cpp + rm -f ../pair_peri_lps.cpp rm -f ../fix_peri_neigh.cpp rm -f ../compute_damage_atom.cpp rm -f ../atom_vec_peri.h rm -f ../pair_peri_pmb.h + rm -f ../pair_peri_lps.h rm -f ../fix_peri_neigh.h rm -f ../compute_damage_atom.h diff --git a/src/PERI/fix_peri_neigh.cpp b/src/PERI/fix_peri_neigh.cpp index a37f83eebc..1e41f0c511 100644 --- a/src/PERI/fix_peri_neigh.cpp +++ b/src/PERI/fix_peri_neigh.cpp @@ -17,18 +17,21 @@ #include "math.h" #include "fix_peri_neigh.h" +#include "pair_peri_pmb.h" +#include "pair_peri_lps.h" #include "atom.h" +#include "domain.h" #include "force.h" -#include "pair.h" +#include "comm.h" +#include "update.h" #include "neighbor.h" #include "neigh_list.h" #include "neigh_request.h" +#include "pair.h" +#include "lattice.h" #include "memory.h" #include "error.h" -#include "comm.h" -#include "update.h" - using namespace LAMMPS_NS; #define MIN(A,B) ((A) < (B)) ? (A) : (B) @@ -39,6 +42,7 @@ using namespace LAMMPS_NS; FixPeriNeigh::FixPeriNeigh(LAMMPS *lmp,int narg, char **arg) : Fix(lmp, narg, arg) { + restart_global = 1; restart_peratom = 1; first = 1; @@ -51,6 +55,7 @@ FixPeriNeigh::FixPeriNeigh(LAMMPS *lmp,int narg, char **arg) : partner = NULL; r0 = NULL; vinter = NULL; + wvolume = NULL; grow_arrays(atom->nmax); atom->add_callback(0); @@ -60,6 +65,10 @@ FixPeriNeigh::FixPeriNeigh(LAMMPS *lmp,int narg, char **arg) : int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) npartner[i] = 0; + + // set comm sizes needed by this fix + + comm_forward = 1; } /* ---------------------------------------------------------------------- */ @@ -78,6 +87,7 @@ FixPeriNeigh::~FixPeriNeigh() memory->destroy_2d_int_array(partner); memory->destroy_2d_double_array(r0); memory->sfree(vinter); + memory->sfree(wvolume); } /* ---------------------------------------------------------------------- */ @@ -92,6 +102,8 @@ int FixPeriNeigh::setmask() void FixPeriNeigh::init() { + if (!first) return; + // need a full neighbor list once int irequest = neighbor->request((void *) this); @@ -127,7 +139,7 @@ void FixPeriNeigh::setup(int vflag) int *tag = atom->tag; int nlocal = atom->nlocal; - // only build list of bonds on very first run + // only build list of bonds on very first run if (!first) return; first = 0; @@ -190,6 +202,7 @@ void FixPeriNeigh::setup(int vflag) for (i = 0; i < nlocal; i++) { npartner[i] = 0; vinter[i] = 0.0; + wvolume[i] = 0.0; } for (ii = 0; ii < inum; ii++) { @@ -203,7 +216,7 @@ void FixPeriNeigh::setup(int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - + delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; @@ -211,14 +224,74 @@ void FixPeriNeigh::setup(int vflag) jtype = type[j]; if (rsq <= cutsq[itype][jtype]) { - partner[i][npartner[i]] = tag[j]; - r0[i][npartner[i]] = sqrt(rsq); - npartner[i]++; + partner[i][npartner[i]] = tag[j]; + r0[i][npartner[i]] = sqrt(rsq); + npartner[i]++; vinter[i] += vfrac[j]; } } } + // compute wvolume for each atom + + double **x0 = atom->x0; + double half_lc = 0.5*(domain->lattice->xlattice); + double vfrac_scale; + PairPeriLPS *pairlps = dynamic_cast(anypair); + PairPeriPMB *pairpmb = dynamic_cast(anypair); + for (i = 0; i < nlocal; i++) { + double xtmp0 = x0[i][0]; + double ytmp0 = x0[i][1]; + double ztmp0 = x0[i][2]; + jnum = npartner[i]; + itype = type[i]; + + // loop over partners of particle i + + for (jj = 0; jj < jnum; jj++) { + + // if bond already broken, skip this partner + + if (partner[i][jj] == 0) continue; + + // lookup local index of partner particle + + j = atom->map(partner[i][jj]); + + // skip if particle is "lost" + + if (j < 0) continue; + + double delx0 = xtmp0 - x0[j][0]; + double dely0 = ytmp0 - x0[j][1]; + double delz0 = ztmp0 - x0[j][2]; + double rsq0 = delx0*delx0 + dely0*dely0 + delz0*delz0; + + jtype = type[j]; + double delta = sqrt(cutsq[itype][jtype]); + + // scale vfrac[j] if particle j near the horizon + + if ((fabs(r0[i][jj] - delta)) <= half_lc) + vfrac_scale = (-1.0/(2*half_lc))*(r0[i][jj]) + + (1.0 + ((delta - half_lc)/(2*half_lc) ) ); + else vfrac_scale = 1.0; + + if (pairpmb != NULL) // define influence function to be 1.0 + wvolume[i] += 1.0 * rsq0 * vfrac[j] * vfrac_scale; + else if (pairlps != NULL) // call the PairPeriLPS influence function + wvolume[i] += pairlps->influence_function(delx0,dely0,delz0) * + rsq0 * vfrac[j] * vfrac_scale; + else + error->all("Unknown peridynamic pair style in FixPeriNeigh."); + + } + } + + // communicate wvolume to ghosts + + comm->forward_comm_fix(this); + // bond statistics int n = 0; @@ -251,6 +324,7 @@ double FixPeriNeigh::memory_usage() bytes += nmax*maxpartner * sizeof(int); bytes += nmax*maxpartner * sizeof(double); bytes += nmax * sizeof(double); + bytes += nmax * sizeof(double); return bytes; } @@ -267,6 +341,8 @@ void FixPeriNeigh::grow_arrays(int nmax) r0 = memory->grow_2d_double_array(r0,nmax,maxpartner,"peri_neigh:r0"); vinter = (double *) memory->srealloc(vinter,nmax*sizeof(double), "peri_neigh:vinter"); + wvolume = (double *) memory->srealloc(wvolume,nmax*sizeof(double), + "peri_neigh:wvolume"); } /* ---------------------------------------------------------------------- @@ -281,6 +357,7 @@ void FixPeriNeigh::copy_arrays(int i, int j) r0[j][m] = r0[i][m]; } vinter[j] = vinter[i]; + wvolume[j] = wvolume[i]; } /* ---------------------------------------------------------------------- @@ -301,6 +378,7 @@ int FixPeriNeigh::pack_exchange(int i, double *buf) buf[0] = m/2; buf[m++] = vinter[i]; + buf[m++] = wvolume[i]; return m; } @@ -317,9 +395,73 @@ int FixPeriNeigh::unpack_exchange(int nlocal, double *buf) r0[nlocal][n] = buf[m++]; } vinter[nlocal] = buf[m++]; + wvolume[nlocal] = buf[m++]; return m; } +/* ---------------------------------------------------------------------- */ + +int FixPeriNeigh::pack_comm(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,m; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = wvolume[j]; + } + return 1; +} + +/* ---------------------------------------------------------------------- */ + +void FixPeriNeigh::unpack_comm(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + wvolume[i] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- + pack entire state of Fix into one write +------------------------------------------------------------------------- */ + +void FixPeriNeigh::write_restart(FILE *fp) +{ + int n = 0; + double list[2]; + list[n++] = first; + list[n++] = maxpartner; + + if (comm->me == 0) { + int size = n * sizeof(double); + fwrite(&size,sizeof(int),1,fp); + fwrite(list,sizeof(double),n,fp); + } +} + +/* ---------------------------------------------------------------------- + use state info from restart file to restart the Fix +------------------------------------------------------------------------- */ + +void FixPeriNeigh::restart(char *buf) +{ + int n = 0; + double *list = (double *) buf; + + first = static_cast (list[n++]); + maxpartner = static_cast (list[n++]); + + // grow 2D arrays now, cannot change size of 2nd array index later + + grow_arrays(atom->nmax); +} + /* ---------------------------------------------------------------------- pack values in local atom-based arrays for restart file ------------------------------------------------------------------------- */ @@ -327,13 +469,14 @@ int FixPeriNeigh::unpack_exchange(int nlocal, double *buf) int FixPeriNeigh::pack_restart(int i, double *buf) { int m = 0; - buf[m++] = 2*npartner[i] + 2; + buf[m++] = 2*npartner[i] + 4; buf[m++] = npartner[i]; for (int n = 0; n < npartner[i]; n++) { buf[m++] = partner[i][n]; buf[m++] = r0[i][n]; } buf[m++] = vinter[i]; + buf[m++] = wvolume[i]; return m; } @@ -357,6 +500,7 @@ void FixPeriNeigh::unpack_restart(int nlocal, int nth) r0[nlocal][n] = extra[nlocal][m++]; } vinter[nlocal] = extra[nlocal][m++]; + wvolume[nlocal] = extra[nlocal][m++]; } /* ---------------------------------------------------------------------- @@ -365,7 +509,7 @@ void FixPeriNeigh::unpack_restart(int nlocal, int nth) int FixPeriNeigh::maxsize_restart() { - return 2*maxpartner + 3; + return 2*maxpartner + 4; } /* ---------------------------------------------------------------------- @@ -374,5 +518,5 @@ int FixPeriNeigh::maxsize_restart() int FixPeriNeigh::size_restart(int nlocal) { - return 2*npartner[nlocal] + 3; + return 2*npartner[nlocal] + 4; } diff --git a/src/PERI/fix_peri_neigh.h b/src/PERI/fix_peri_neigh.h index 9bd974ef85..239d18931f 100644 --- a/src/PERI/fix_peri_neigh.h +++ b/src/PERI/fix_peri_neigh.h @@ -26,6 +26,7 @@ namespace LAMMPS_NS { class FixPeriNeigh : public Fix { friend class PairPeriPMB; + friend class PairPeriLPS; friend class ComputeDamageAtom; public: @@ -41,10 +42,15 @@ class FixPeriNeigh : public Fix { void copy_arrays(int, int); int pack_exchange(int, double *); int unpack_exchange(int, double *); + void write_restart(FILE *); + void restart(char *); int pack_restart(int, double *); void unpack_restart(int, int); int size_restart(int); int maxsize_restart(); + int pack_comm(int, int *, double *, int, int *); + void unpack_comm(int, int, double *); + private: int first; // flag for first time initialization @@ -53,6 +59,7 @@ class FixPeriNeigh : public Fix { int **partner; // neighs for each atom, stored as global IDs double **r0; // initial distance to partners double *vinter; // sum of vfrac for bonded neighbors + double *wvolume; // weighted volume of particle class NeighList *list; }; diff --git a/src/PERI/pair_peri_lps.cpp b/src/PERI/pair_peri_lps.cpp new file mode 100644 index 0000000000..504e78f859 --- /dev/null +++ b/src/PERI/pair_peri_lps.cpp @@ -0,0 +1,745 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Mike Parks (SNL) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "values.h" +#include "stdlib.h" +#include "string.h" +#include "pair_peri_lps.h" +#include "atom.h" +#include "domain.h" +#include "lattice.h" +#include "force.h" +#include "update.h" +#include "modify.h" +#include "fix.h" +#include "fix_peri_neigh.h" +#include "comm.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "memory.h" +#include "error.h" +#include "update.h" + +using namespace LAMMPS_NS; + +#define MIN(a,b) ((a) < (b) ? (a) : (b)) +#define MAX(a,b) ((a) > (b) ? (a) : (b)) + +/* ---------------------------------------------------------------------- */ + +PairPeriLPS::PairPeriLPS(LAMMPS *lmp) : Pair(lmp) +{ + for (int i = 0; i < 6; i++) virial[i] = 0.0; + + ifix_peri = -1; + + nmax = 0; + s0_new = NULL; + theta = NULL; + + bulkmodulus = NULL; + shearmodulus = NULL; + s00 = alpha = NULL; + cut = NULL; + + // set comm size needed by this Pair + // comm_reverse not needed + comm_forward = 1; // For passing dilatation (theta) +} + +/* ---------------------------------------------------------------------- */ + +PairPeriLPS::~PairPeriLPS() +{ + if (ifix_peri >= 0) modify->delete_fix("PERI_NEIGH"); + + if (allocated) { + memory->destroy_2d_int_array(setflag); + memory->destroy_2d_double_array(cutsq); + memory->destroy_2d_double_array(bulkmodulus); + memory->destroy_2d_double_array(shearmodulus); + memory->destroy_2d_double_array(s00); + memory->destroy_2d_double_array(alpha); + memory->destroy_2d_double_array(cut); + memory->sfree(theta); + memory->sfree(s0_new); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairPeriLPS::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz; + double xtmp0,ytmp0,ztmp0,delx0,dely0,delz0,rsq0; + double rsq,r,dr,rk,evdwl,fpair,fbond; + int *ilist,*jlist,*numneigh,**firstneigh; + double d_ij,delta,stretch; + + evdwl = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + double **f = atom->f; + double **x = atom->x; + int *type = atom->type; + int nlocal = atom->nlocal; + + double *vfrac = atom->vfrac; + double *s0 = atom->s0; + double **x0 = atom->x0; + double **r0 = ((FixPeriNeigh *) modify->fix[ifix_peri])->r0; + int **partner = ((FixPeriNeigh *) modify->fix[ifix_peri])->partner; + int *npartner = ((FixPeriNeigh *) modify->fix[ifix_peri])->npartner; + double *wvolume = ((FixPeriNeigh *) modify->fix[ifix_peri])->wvolume; + + // lc = lattice constant + // init_style guarantees it's the same in x, y, and z + + double lc = domain->lattice->xlattice; + double half_lc = 0.5*lc; + double vfrac_scale = 1.0; + + // short-range forces + + int nall = atom->nlocal + atom->nghost; + int newton_pair = force->newton_pair; + int nonperiodic = domain->nonperiodic; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + double dt = update->dt; + + // loop over neighbors of my atoms + // need minimg() for x0 difference since not ghosted + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + xtmp0 = x0[i][0]; + ytmp0 = x0[i][1]; + ztmp0 = x0[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j %= nall; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + + rsq = delx*delx + dely*dely + delz*delz; + delx0 = xtmp0 - x0[j][0]; + dely0 = ytmp0 - x0[j][1]; + delz0 = ztmp0 - x0[j][2]; + if (nonperiodic == 0) domain->minimum_image(delx0,dely0,delz0); + rsq0 = delx0*delx0 + dely0*dely0 + delz0*delz0; + jtype = type[j]; + + r = sqrt(rsq); + + // short-range interaction distance based on initial particle position + // 0.9 and 1.35 are constants + + d_ij = MIN(0.9*sqrt(rsq0),1.35*lc); + + // short-range contact forces + // 15 is constant taken from the EMU Theory Manual + // Silling, 12 May 2005, p 18 + + if (r < d_ij) { + dr = r - d_ij; + + // kshort based upon short-range force constant of the bond-based theory used in PMB model + double kshort = ( 15.0 * 18.0 * bulkmodulus[itype][itype] ) / + ( 3.141592653589793 * cutsq[itype][jtype] * cutsq[itype][jtype] ); + rk = ( kshort * vfrac[j]) * (dr / sqrt(cutsq[itype][jtype])); + + if (r > 0.0) fpair = -(rk/r); + else fpair = 0.0; + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + } + + if (eflag) evdwl = 0.5*rk*dr; + if (evflag) ev_tally(i,j,nlocal,newton_pair, + evdwl,0.0,fpair,delx,dely,delz); + } + } + } + + // grow bond forces array if necessary + + if (atom->nmax > nmax) { + memory->sfree(s0_new); + memory->sfree(theta); + nmax = atom->nmax; + s0_new = (double *) memory->smalloc(nmax*sizeof(double),"pair:s0_new"); + theta = (double *) memory->smalloc(nmax*sizeof(double),"pair:theta"); + } + + // Compute the dilatation on each particle + compute_dilatation(); + + // communicate dilatation (theta) of each particle + comm->forward_comm_pair(this); + // communicate wighted volume (wvolume) upon every reneighbor + if (neighbor->ago == 0) + comm->forward_comm_fix(modify->fix[ifix_peri]); + + // Volume-dependent part of the energy + for (i = 0; i < nlocal; i++) { + itype = type[i]; + if (eflag) { + if (eflag_global) eng_vdwl += 0.5 * bulkmodulus[itype][itype] * (theta[i] * theta[i]); + if (eflag_atom) eatom[i] += 0.5 * bulkmodulus[itype][itype] * (theta[i] * theta[i]); + } + } + + // loop over my particles and their partners + // partner list contains all bond partners, so I-J appears twice + // if bond already broken, skip this partner + // first = true if this is first neighbor of particle i + + bool first; + double omega_minus, omega_plus; + + for (i = 0; i < nlocal; i++) { + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + xtmp0 = x0[i][0]; + ytmp0 = x0[i][1]; + ztmp0 = x0[i][2]; + itype = type[i]; + jnum = npartner[i]; + first = true; + + for (jj = 0; jj < jnum; jj++) { + if (partner[i][jj] == 0) continue; + j = atom->map(partner[i][jj]); + + // check if lost a partner without first breaking bond + + if (j < 0) { + partner[i][jj] = 0; + continue; + } + + // compute force density, add to PD equation of motion + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + if (nonperiodic == 0) domain->minimum_image(delx,dely,delz); + rsq = delx*delx + dely*dely + delz*delz; + delx0 = xtmp0 - x0[j][0]; + dely0 = ytmp0 - x0[j][1]; + delz0 = ztmp0 - x0[j][2]; + if (nonperiodic == 0) domain->minimum_image(delx0,dely0,delz0); + jtype = type[j]; + delta = sqrt(cutsq[itype][jtype]); + r = sqrt(rsq); + dr = r - r0[i][jj]; + + // avoid roundoff errors + + if (fabs(dr) < 2.2204e-016) dr = 0.0; + + // scale vfrac[j] if particle j near the horizon + + if ((fabs(r0[i][jj] - delta)) <= half_lc) + vfrac_scale = (-1.0/(2*half_lc))*(r0[i][jj]) + + (1.0 + ((delta - half_lc)/(2*half_lc) ) ); + else vfrac_scale = 1.0; + + omega_plus = influence_function(-1.0*delx0,-1.0*dely0,-1.0*delz0); + omega_minus = influence_function(delx0,dely0,delz0); + rk = ( (3.0 * bulkmodulus[itype][itype]) - (5.0 * shearmodulus[itype][itype]) ) * vfrac[j] * vfrac_scale * ( (omega_plus * theta[i] / wvolume[i]) + ( omega_minus * theta[j] / wvolume[j] ) ) * r0[i][jj] ; + rk += 15.0 * ( shearmodulus[itype][itype] * vfrac[j] * vfrac_scale ) * ( (omega_plus / wvolume[i]) + (omega_minus / wvolume[j]) ) * dr; + + if (r > 0.0) fbond = -(rk/r); + else fbond = 0.0; + + f[i][0] += delx*fbond; + f[i][1] += dely*fbond; + f[i][2] += delz*fbond; + + // since I-J is double counted, set newton off & use 1/2 factor and I,I + + double deviatoric_extension = dr - (theta[i]* r0[i][jj] / 3.0); + if (eflag) evdwl = 0.5 * 15 * (shearmodulus[itype][itype]/wvolume[i]) * omega_plus * ( deviatoric_extension * deviatoric_extension ) * vfrac[j] * vfrac_scale; + if (evflag) ev_tally(i,i,nlocal,0, + 0.5*evdwl,0.0,0.5*fbond,delx,dely,delz); + + // find stretch in bond I-J and break if necessary + // use s0 from previous timestep + + stretch = dr / r0[i][jj]; + if (stretch > MIN(s0[i],s0[j])) partner[i][jj] = 0; + + // update s0 for next timestep + + if (first) + s0_new[i] = s00[itype][jtype] - (alpha[itype][jtype] * stretch); + else + s0_new[i] = MAX(s0_new[i], + s00[itype][jtype] - (alpha[itype][jtype] * stretch)); + + first = false; + } + } + + if (vflag_fdotr) virial_compute(); + + // store new s0 + for (i = 0; i < nlocal; i++) s0[i] = s0_new[i]; + +} + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairPeriLPS::allocate() +{ + allocated = 1; + int n = atom->ntypes; + + setflag = memory->create_2d_int_array(n+1,n+1,"pair:setflag"); + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + setflag[i][j] = 0; + + cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq"); + bulkmodulus = memory->create_2d_double_array(n+1,n+1,"pair:bulkmodulus"); + shearmodulus = memory->create_2d_double_array(n+1,n+1,"pair:shearmodulus"); + s00 = memory->create_2d_double_array(n+1,n+1,"pair:s00"); + alpha = memory->create_2d_double_array(n+1,n+1,"pair:alpha"); + cut = memory->create_2d_double_array(n+1,n+1,"pair:cut"); + +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairPeriLPS::settings(int narg, char **arg) +{ + if (narg) error->all("Illegal pair_style command"); +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairPeriLPS::coeff(int narg, char **arg) +{ + if (narg != 7) error->all("Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo,ihi,jlo,jhi; + force->bounds(arg[0],atom->ntypes,ilo,ihi); + force->bounds(arg[1],atom->ntypes,jlo,jhi); + + double bulkmodulus_one = atof(arg[2]); + double shearmodulus_one = atof(arg[3]); + double cut_one = atof(arg[4]); + double s00_one = atof(arg[5]); + double alpha_one = atof(arg[6]); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + bulkmodulus[i][j] = bulkmodulus_one; + shearmodulus[i][j] = shearmodulus_one; + cut[i][j] = cut_one; + s00[i][j] = s00_one; + alpha[i][j] = alpha_one; + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) error->all("Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairPeriLPS::init_one(int i, int j) +{ + if (setflag[i][j] == 0) error->all("All pair coeffs are not set"); + + cutsq[i][j] = cut[i][j] * cut[i][j]; + cutsq[j][i] = cutsq[i][j]; + + // set other j,i parameters + + bulkmodulus[j][i] = bulkmodulus[i][j]; + shearmodulus[j][i] = shearmodulus[i][j]; + s00[j][i] = s00[i][j]; + alpha[j][i] = alpha[i][j]; + + return cut[i][j]; +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairPeriLPS::init_style() +{ + // error checks + + if (atom->map_style == 0) + error->all("Pair peri requires an atom map, see atom_modify"); + + if (atom->style_match("peri") == 0) + error->all("Pair style peri_lps requires atom style peri"); + + if (domain->lattice == NULL) + error->all("Pair peri requires a lattice be defined"); + if (domain->lattice->xlattice != domain->lattice->ylattice || + domain->lattice->xlattice != domain->lattice->zlattice || + domain->lattice->ylattice != domain->lattice->zlattice) + error->all("Pair peri lattice is not identical in x, y, and z"); + + // if first init, create Fix needed for storing fixed neighbors + + if (ifix_peri == -1) { + char **fixarg = new char*[3]; + fixarg[0] = (char *) "PERI_NEIGH"; + fixarg[1] = (char *) "all"; + fixarg[2] = (char *) "PERI_NEIGH"; + modify->add_fix(3,fixarg); + delete [] fixarg; + } + + // find associated PERI_NEIGH fix that must exist + // could have changed locations in fix list since created + + for (int i = 0; i < modify->nfix; i++) + if (strcmp(modify->fix[i]->style,"PERI_NEIGH") == 0) ifix_peri = i; + if (ifix_peri == -1) error->all("Fix peri neigh does not exist"); + + int irequest = neighbor->request(this); +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairPeriLPS::write_restart(FILE *fp) +{ + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + fwrite(&setflag[i][j],sizeof(int),1,fp); + if (setflag[i][j]) { + fwrite(&bulkmodulus[i][j],sizeof(double),1,fp); + fwrite(&shearmodulus[i][j],sizeof(double),1,fp); + fwrite(&s00[i][j],sizeof(double),1,fp); + fwrite(&alpha[i][j],sizeof(double),1,fp); + fwrite(&cut[i][j],sizeof(double),1,fp); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairPeriLPS::read_restart(FILE *fp) +{ + allocate(); + + int i,j; + int me = comm->me; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp); + MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); + if (setflag[i][j]) { + if (me == 0) { + fread(&bulkmodulus[i][j],sizeof(double),1,fp); + fread(&shearmodulus[i][j],sizeof(double),1,fp); + fread(&s00[i][j],sizeof(double),1,fp); + fread(&alpha[i][j],sizeof(double),1,fp); + fread(&cut[i][j],sizeof(double),1,fp); + } + MPI_Bcast(&bulkmodulus[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&shearmodulus[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&s00[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&alpha[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); + } + } +} + +/* ---------------------------------------------------------------------- */ + +double PairPeriLPS::single(int i, int j, int itype, int jtype, double rsq, + double factor_coul, double factor_lj, + double &fforce) +{ + double delx0,dely0,delz0,rsq0; + double d_ij,r,dr,rk,vfrac_scale; + + double *vfrac = atom->vfrac; + double **x0 = atom->x0; + double **r0 = ((FixPeriNeigh *) modify->fix[ifix_peri])->r0; + int **partner = ((FixPeriNeigh *) modify->fix[ifix_peri])->partner; + int *npartner = ((FixPeriNeigh *) modify->fix[ifix_peri])->npartner; + double *wvolume = ((FixPeriNeigh *) modify->fix[ifix_peri])->wvolume; + + double lc = domain->lattice->xlattice; + double half_lc = 0.5*lc; + + double kshort; + + delx0 = x0[i][0] - x0[j][0]; + dely0 = x0[i][1] - x0[j][1]; + delz0 = x0[i][2] - x0[j][2]; + if (domain->nonperiodic == 0) domain->minimum_image(delx0,dely0,delz0); + rsq0 = delx0*delx0 + dely0*dely0 + delz0*delz0; + + d_ij = MIN(0.9*sqrt(rsq0),1.35*lc); + r = sqrt(rsq); + + double energy = 0.0; + fforce = 0.0; + + if (r < d_ij) { + dr = r - d_ij; + // kshort resembles the short-range force constant of the bond-based theory in 3-D ! + kshort = (15.0 * 18.0 * bulkmodulus[itype][itype]) / + ( 3.141592653589793 * cutsq[itype][jtype] * cutsq[itype][jtype] ); + rk = ( kshort * vfrac[j]) * (dr / sqrt(cutsq[itype][jtype])); + if (r > 0.0) fforce += -(rk/r); + energy += 0.5*rk*dr; + } + + if (atom->nmax > nmax) { + memory->sfree(theta); + nmax = atom->nmax; + theta = (double *) memory->smalloc(nmax*sizeof(double),"pair:theta"); + } + + // Compute the dilatation on each particle + compute_dilatation(); + + // communicate dilatation (theta) of each particle + comm->forward_comm_pair(this); + // communicate wighted volume (wvolume) upon every reneighbor + if (neighbor->ago == 0) + comm->forward_comm_fix(modify->fix[ifix_peri]); + + double omega_plus, omega_minus; + + int jnum = npartner[i]; + for (int jj = 0; jj < jnum; jj++) { + if (partner[i][jj] == 0) continue; + if (j < 0) continue; + if (j == atom->map(partner[i][jj])) { + dr = r - r0[i][jj]; + if (fabs(dr) < 2.2204e-016) dr = 0.0; + + // scale vfrac[j] if particle j near the horizon + + if ( (fabs(r0[i][jj] - sqrt(cutsq[itype][jtype]))) <= half_lc) + vfrac_scale = (-1.0/(2*half_lc))*(r0[i][jj]) + + (1.0 + ((sqrt(cutsq[itype][jtype]) - half_lc)/(2*half_lc))); + else vfrac_scale = 1.0; + + omega_plus = influence_function(-1.0*delx0,-1.0*dely0,-1.0*delz0); + omega_minus = influence_function(delx0,dely0,delz0); + rk = (3.0* bulkmodulus[itype][itype] -5.0 * shearmodulus[itype][itype] )* vfrac[j] * vfrac_scale * ( (omega_plus * theta[i] / wvolume[i]) + ( omega_minus * theta[j] / wvolume[j] ) ) * r0[i][jj] ; + rk += 15.0 * ( shearmodulus[itype][itype] * vfrac[j] * vfrac_scale ) * ( (omega_plus / wvolume[i]) + (omega_minus / wvolume[j]) ) * dr; + + if (r > 0.0) fforce += -(rk/r); + energy += 0.5 * 15 * (shearmodulus[itype][itype]/wvolume[i]) * omega_plus * ( dr - theta[i]* r0[i][jj] / 3.0 ) * ( dr - theta[i]* r0[i][jj] / 3.0 ) * vfrac[j] * vfrac_scale; + + } + } + + return energy; +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based arrays +------------------------------------------------------------------------- */ + +double PairPeriLPS::memory_usage() +{ + double bytes = 2 * nmax * sizeof(double); + return bytes; +} + +/* ---------------------------------------------------------------------- + Influence function definition +------------------------------------------------------------------------- */ + +double PairPeriLPS::influence_function(double xi_x, double xi_y, double xi_z) +{ + double r = sqrt(xi_x*xi_x + xi_y*xi_y + xi_z*xi_z); + double omega; + + if (fabs(r) < 2.2204e-016) error->one("pair_peri_lps - influence_function: Divide by 0"); + + omega = 1.0/r; + + return omega; +} + +/* ---------------------------------------------------------------------- */ + +void PairPeriLPS::compute_dilatation() +{ + int i,j,jj,jnum,itype,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz; + double xtmp0,ytmp0,ztmp0,delx0,dely0,delz0; + double rsq,r,dr; + double delta; + + double **x = atom->x; + int *type = atom->type; + double **x0 = atom->x0; + int nlocal = atom->nlocal; + double *vfrac = atom->vfrac; + double vfrac_scale = 1.0; + + double lc = domain->lattice->xlattice; + double half_lc = 0.5*lc; + + double **r0 = ((FixPeriNeigh *) modify->fix[ifix_peri])->r0; + int **partner = ((FixPeriNeigh *) modify->fix[ifix_peri])->partner; + int *npartner = ((FixPeriNeigh *) modify->fix[ifix_peri])->npartner; + double *wvolume = ((FixPeriNeigh *) modify->fix[ifix_peri])->wvolume; + + int nonperiodic = domain->nonperiodic; + + // compute the dilatation theta + // loop over my particles + for (i = 0; i < nlocal; i++) { + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + xtmp0 = x0[i][0]; + ytmp0 = x0[i][1]; + ztmp0 = x0[i][2]; + jnum = npartner[i]; + theta[i] = 0.0; + itype = type[i]; + + // loop over partners of particle i + for (jj = 0; jj < jnum; jj++) { + + // if bond already broken, skip this partner + if (partner[i][jj] == 0) continue; + + // Look up local index of this partner particle + j = atom->map(partner[i][jj]); + + // Skip if particle is "lost" + if (j < 0) continue; + + // Compute force density and add to PD equation of motion + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + if (nonperiodic == 0) domain->minimum_image(delx,dely,delz); + rsq = delx*delx + dely*dely + delz*delz; + delx0 = xtmp0 - x0[j][0]; + dely0 = ytmp0 - x0[j][1]; + delz0 = ztmp0 - x0[j][2]; + if (nonperiodic == 0) domain->minimum_image(delx0,dely0,delz0); + + r = sqrt(rsq); + dr = r - r0[i][jj]; + if (fabs(dr) < 2.2204e-016) dr = 0.0; + + jtype = type[j]; + delta = sqrt(cutsq[itype][jtype]); + + // scale vfrac[j] if particle j near the horizon + + if ((fabs(r0[i][jj] - delta)) <= half_lc) + vfrac_scale = (-1.0/(2*half_lc))*(r0[i][jj]) + + (1.0 + ((delta - half_lc)/(2*half_lc) ) ); + else vfrac_scale = 1.0; + + theta[i] += influence_function(delx0, dely0, delz0) * r0[i][jj] * dr * vfrac[j] * vfrac_scale; + + } // end loop over all neighbors jj + + // If wvolume[i] is zero, then particle i has no bonds. Therefore, the dilatation is set to 0. + if (wvolume[i] != 0.0) theta[i] = (3.0/wvolume[i]) * theta[i]; + else theta[i] = 0; + + } // end loop over all particles i + +} + +/* ---------------------------------------------------------------------- + communication routines + ---------------------------------------------------------------------- */ + +int PairPeriLPS::pack_comm(int n, int *list, double *buf, int pbc_flag, int *pbc) +{ + + int i,j,m; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = theta[j]; + } + return 1; +} + +/* ---------------------------------------------------------------------- */ + +void PairPeriLPS::unpack_comm(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + theta[i] = buf[m++]; + } +} + + diff --git a/src/PERI/pair_peri_lps.h b/src/PERI/pair_peri_lps.h new file mode 100644 index 0000000000..47ec71901e --- /dev/null +++ b/src/PERI/pair_peri_lps.h @@ -0,0 +1,65 @@ +/* ---------------------------------------------------------------------- + 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 PAIR_CLASS + +PairStyle(peri/lps,PairPeriLPS) + +#else + +#ifndef LMP_PAIR_PERI_LPS_H +#define LMP_PAIR_PERI_LPS_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairPeriLPS : public Pair { + public: + PairPeriLPS(class LAMMPS *); + ~PairPeriLPS(); + int pack_comm(int, int *, double *, int, int *); + void unpack_comm(int, int, double *); + + void compute(int, int); + void settings(int, char **); + void coeff(int, char **); + double init_one(int, int); + void init_style(); + void write_restart(FILE *); + void read_restart(FILE *); + void write_restart_settings(FILE *) {} + void read_restart_settings(FILE *) {} + double single(int, int, int, int, double, double, double, double &); + double memory_usage(); + double influence_function(double, double, double); + void compute_dilatation(); + + protected: + int ifix_peri; + double **bulkmodulus; + double **shearmodulus; + double **s00, **alpha; + double **cut; + + double *s0_new; + double *theta; + int nmax; + + void allocate(); +}; + +} + +#endif +#endif diff --git a/src/PERI/pair_peri_pmb.cpp b/src/PERI/pair_peri_pmb.cpp index ffc0fe2812..1f3163f9ef 100644 --- a/src/PERI/pair_peri_pmb.cpp +++ b/src/PERI/pair_peri_pmb.cpp @@ -31,6 +31,7 @@ #include "comm.h" #include "neighbor.h" #include "neigh_list.h" +#include "neigh_request.h" #include "memory.h" #include "error.h" @@ -399,6 +400,7 @@ void PairPeriPMB::init_style() for (int i = 0; i < modify->nfix; i++) if (strcmp(modify->fix[i]->style,"PERI_NEIGH") == 0) ifix_peri = i; + if (ifix_peri == -1) error->all("Fix peri neigh does not exist"); int irequest = neighbor->request(this); }