From 192008e7c7fc86655371522e51aa0814f5f23ea3 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 24 Jul 2008 20:46:29 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@1965 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/Makefile | 2 +- src/PERI/Install.csh | 32 ++ src/PERI/atom_vec_peri.cpp | 613 +++++++++++++++++++++++++++++++ src/PERI/atom_vec_peri.h | 54 +++ src/PERI/compute_damage_atom.cpp | 132 +++++++ src/PERI/compute_damage_atom.h | 37 ++ src/PERI/fix_peri_neigh.cpp | 365 ++++++++++++++++++ src/PERI/fix_peri_neigh.h | 55 +++ src/PERI/pair_peri_pmb.cpp | 520 ++++++++++++++++++++++++++ src/PERI/pair_peri_pmb.h | 51 +++ src/PERI/style_peri.h | 45 +++ src/atom.cpp | 8 +- src/atom.h | 4 +- src/set.cpp | 55 ++- src/set.h | 1 + src/style.h | 1 + src/update.cpp | 21 ++ 17 files changed, 1987 insertions(+), 9 deletions(-) create mode 100644 src/PERI/Install.csh create mode 100644 src/PERI/atom_vec_peri.cpp create mode 100755 src/PERI/atom_vec_peri.h create mode 100644 src/PERI/compute_damage_atom.cpp create mode 100644 src/PERI/compute_damage_atom.h create mode 100644 src/PERI/fix_peri_neigh.cpp create mode 100644 src/PERI/fix_peri_neigh.h create mode 100644 src/PERI/pair_peri_pmb.cpp create mode 100644 src/PERI/pair_peri_pmb.h create mode 100644 src/PERI/style_peri.h diff --git a/src/Makefile b/src/Makefile index f026d4adb5..360eba0a75 100755 --- a/src/Makefile +++ b/src/Makefile @@ -14,7 +14,7 @@ OBJ = $(SRC:.cpp=.o) # Package variables PACKAGE = asphere class2 colloid dipole dpd granular \ - kspace manybody meam molecule opt poems xtc + kspace manybody meam molecule opt peri poems xtc PACKUSER = user-ackland user-cg-cmm user-ewaldn user-smd diff --git a/src/PERI/Install.csh b/src/PERI/Install.csh new file mode 100644 index 0000000000..4321be0177 --- /dev/null +++ b/src/PERI/Install.csh @@ -0,0 +1,32 @@ +# Install/unInstall package classes in LAMMPS + +if ($1 == 1) then + + cp style_peri.h .. + + cp atom_vec_peri.cpp .. + cp pair_peri_pmb.cpp .. + cp fix_peri_neigh.cpp .. + cp compute_damage_atom.cpp .. + + cp atom_vec_peri.h .. + cp pair_peri_pmb.h .. + cp fix_peri_neigh.h .. + cp compute_damage_atom.h .. + +else if ($1 == 0) then + + rm -f ../style_peri.h + touch ../style_peri.h + + rm -f ../atom_vec_peri.cpp + rm -f ../pair_peri_pmb.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 ../fix_peri_neigh.h + rm -f ../compute_damage_atom.h + +endif diff --git a/src/PERI/atom_vec_peri.cpp b/src/PERI/atom_vec_peri.cpp new file mode 100644 index 0000000000..dbc87887e1 --- /dev/null +++ b/src/PERI/atom_vec_peri.cpp @@ -0,0 +1,613 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + www.cs.sandia.gov/~sjplimp/lammps.html + Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories + + 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 "values.h" +#include "stdlib.h" +#include "atom_vec_peri.h" +#include "atom.h" +#include "domain.h" +#include "modify.h" +#include "fix.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +#define DELTA 10000 + +/* ---------------------------------------------------------------------- */ + +AtomVecPeri::AtomVecPeri(LAMMPS *lmp, int narg, char **arg) : +AtomVec(lmp, narg, arg) +{ + comm_x_only = 0; + size_comm = 4; + size_reverse = 3; + size_border = 13; + size_data_atom = 7; + size_data_vel = 4; + xcol_data = 3; + + atom->rmass_flag = atom->vfrac_flag = atom->s0_flag = 1; + atom->vinter_flag = 1; +} + +/* ---------------------------------------------------------------------- + grow atom arrays + n = 0 grows arrays by DELTA + n > 0 allocates arrays to size n +------------------------------------------------------------------------- */ + +void AtomVecPeri::grow(int n) +{ + if (n == 0) nmax += DELTA; + else nmax = n; + atom->nmax = nmax; + + tag = atom->tag = (int *) + memory->srealloc(atom->tag,nmax*sizeof(int),"atom:tag"); + type = atom->type = (int *) + memory->srealloc(atom->type,nmax*sizeof(int),"atom:type"); + mask = atom->mask = (int *) + memory->srealloc(atom->mask,nmax*sizeof(int),"atom:mask"); + image = atom->image = (int *) + memory->srealloc(atom->image,nmax*sizeof(int),"atom:image"); + x = atom->x = memory->grow_2d_double_array(atom->x,nmax,3,"atom:x"); + v = atom->v = memory->grow_2d_double_array(atom->v,nmax,3,"atom:v"); + f = atom->f = memory->grow_2d_double_array(atom->f,nmax,3,"atom:f"); + + vfrac = atom->vfrac = (double *) + memory->srealloc(atom->vfrac,nmax*sizeof(double),"atom:vfrac"); + rmass = atom->rmass = (double *) + memory->srealloc(atom->rmass,nmax*sizeof(double),"atom:rmass"); + s0 = atom->s0 = (double *) + memory->srealloc(atom->s0,nmax*sizeof(double),"atom:s0"); + x0 = atom->x0 = memory->grow_2d_double_array(atom->x0,nmax,3,"atom:x0"); + vinter = atom->vinter = (double *) + memory->srealloc(atom->vinter,nmax*sizeof(double),"atom:vinter"); + + if (atom->nextra_grow) + for (int iextra = 0; iextra < atom->nextra_grow; iextra++) + modify->fix[atom->extra_grow[iextra]]->grow_arrays(nmax); +} + +/* ---------------------------------------------------------------------- */ + +void AtomVecPeri::copy(int i, int j) +{ + tag[j] = tag[i]; + type[j] = type[i]; + mask[j] = mask[i]; + image[j] = image[i]; + x[j][0] = x[i][0]; + x[j][1] = x[i][1]; + x[j][2] = x[i][2]; + v[j][0] = v[i][0]; + v[j][1] = v[i][1]; + v[j][2] = v[i][2]; + + vfrac[j] = vfrac[i]; + rmass[j] = rmass[i]; + s0[j] = s0[i]; + x0[j][0] = x0[i][0]; + x0[j][1] = x0[i][1]; + x0[j][2] = x0[i][2]; + vinter[j] = vinter[i]; + + if (atom->nextra_grow) + for (int iextra = 0; iextra < atom->nextra_grow; iextra++) + modify->fix[atom->extra_grow[iextra]]->copy_arrays(i,j); +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecPeri::pack_comm(int n, int *list, double *buf, + int pbc_flag, int *pbc) + +{ + int i,j,m; + double dx,dy,dz; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; + buf[m++] = s0[j]; + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0]*domain->xprd; + dy = pbc[1]*domain->yprd; + dz = pbc[2]*domain->zprd; + } else { + dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz; + dy = pbc[1]*domain->yprd + pbc[3]*domain->yz; + dz = pbc[2]*domain->zprd; + } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = s0[j]; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecPeri::pack_comm_one(int i, double *buf) +{ + buf[0] = s0[i]; + return 1; +} + +/* ---------------------------------------------------------------------- */ + +void AtomVecPeri::unpack_comm(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + s0[i] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecPeri::unpack_comm_one(int i, double *buf) +{ + s0[i] = buf[0]; + return 1; +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecPeri::pack_reverse(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + buf[m++] = f[i][0]; + buf[m++] = f[i][1]; + buf[m++] = f[i][2]; + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +void AtomVecPeri::unpack_reverse(int n, int *list, double *buf) +{ + int i,j,m; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + f[j][0] += buf[m++]; + f[j][1] += buf[m++]; + f[j][2] += buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecPeri::pack_border(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,m; + double dx,dy,dz; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; + buf[m++] = vfrac[j]; + buf[m++] = rmass[j]; + buf[m++] = s0[j]; + buf[m++] = x0[j][0]; + buf[m++] = x0[j][1]; + buf[m++] = x0[j][2]; + buf[m++] = vinter[j]; + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0]*domain->xprd; + dy = pbc[1]*domain->yprd; + dz = pbc[2]*domain->zprd; + } else { + dx = pbc[0]; + dy = pbc[1]; + dz = pbc[2]; + } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; + buf[m++] = vfrac[j]; + buf[m++] = rmass[j]; + buf[m++] = s0[j]; + buf[m++] = x0[j][0]; + buf[m++] = x0[j][1]; + buf[m++] = x0[j][2]; + buf[m++] = vinter[j]; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecPeri::pack_border_one(int i, double *buf) +{ + buf[0] = vfrac[i]; + buf[1] = rmass[i]; + buf[2] = s0[i]; + buf[3] = x0[i][0]; + buf[4] = x0[i][1]; + buf[5] = x0[i][2]; + buf[6] = vinter[i]; + return 7; +} + +/* ---------------------------------------------------------------------- */ + +void AtomVecPeri::unpack_border(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + if (i == nmax) grow(0); + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + tag[i] = static_cast (buf[m++]); + type[i] = static_cast (buf[m++]); + mask[i] = static_cast (buf[m++]); + vfrac[i] = buf[m++]; + rmass[i] = buf[m++]; + s0[i] = buf[m++]; + x0[i][0] = buf[m++]; + x0[i][1] = buf[m++]; + x0[i][2] = buf[m++]; + vinter[i] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecPeri::unpack_border_one(int i, double *buf) +{ + vfrac[i] = buf[0]; + rmass[i] = buf[1]; + s0[i] = buf[2]; + x0[i][0] = buf[3]; + x0[i][1] = buf[4]; + x0[i][2] = buf[5]; + vinter[i]= buf[6]; + return 7; +} + +/* ---------------------------------------------------------------------- + pack data for atom I for sending to another proc + xyz must be 1st 3 values, so comm::exchange() can test on them +------------------------------------------------------------------------- */ + +int AtomVecPeri::pack_exchange(int i, double *buf) +{ + int m = 1; + buf[m++] = x[i][0]; + buf[m++] = x[i][1]; + buf[m++] = x[i][2]; + buf[m++] = v[i][0]; + buf[m++] = v[i][1]; + buf[m++] = v[i][2]; + buf[m++] = tag[i]; + buf[m++] = type[i]; + buf[m++] = mask[i]; + buf[m++] = image[i]; + + buf[m++] = vfrac[i]; + buf[m++] = rmass[i]; + buf[m++] = s0[i]; + buf[m++] = x0[i][0]; + buf[m++] = x0[i][1]; + buf[m++] = x0[i][2]; + buf[m++] = vinter[i]; + + if (atom->nextra_grow) + for (int iextra = 0; iextra < atom->nextra_grow; iextra++) + m += modify->fix[atom->extra_grow[iextra]]->pack_exchange(i,&buf[m]); + + buf[0] = m; + return m; +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecPeri::unpack_exchange(double *buf) +{ + int nlocal = atom->nlocal; + if (nlocal == nmax) grow(0); + + int m = 1; + x[nlocal][0] = buf[m++]; + x[nlocal][1] = buf[m++]; + x[nlocal][2] = buf[m++]; + v[nlocal][0] = buf[m++]; + v[nlocal][1] = buf[m++]; + v[nlocal][2] = buf[m++]; + tag[nlocal] = static_cast (buf[m++]); + type[nlocal] = static_cast (buf[m++]); + mask[nlocal] = static_cast (buf[m++]); + image[nlocal] = static_cast (buf[m++]); + + vfrac[nlocal] = buf[m++]; + rmass[nlocal] = buf[m++]; + s0[nlocal] = buf[m++]; + x0[nlocal][0] = buf[m++]; + x0[nlocal][1] = buf[m++]; + x0[nlocal][2] = buf[m++]; + vinter[nlocal] = buf[m++]; + + if (atom->nextra_grow) + for (int iextra = 0; iextra < atom->nextra_grow; iextra++) + m += modify->fix[atom->extra_grow[iextra]]-> + unpack_exchange(nlocal,&buf[m]); + + atom->nlocal++; + return m; +} + + +/* ---------------------------------------------------------------------- + size of restart data for all atoms owned by this proc + include extra data stored by fixes +------------------------------------------------------------------------- */ + +int AtomVecPeri::size_restart() +{ + int i; + + int nlocal = atom->nlocal; + int n = 18 * nlocal; + + if (atom->nextra_restart) + for (int iextra = 0; iextra < atom->nextra_restart; iextra++) + for (i = 0; i < nlocal; i++) + n += modify->fix[atom->extra_restart[iextra]]->size_restart(i); + + return n; +} + +/* ---------------------------------------------------------------------- + pack atom I's data for restart file including extra quantities + xyz must be 1st 3 values, so that read_restart can test on them + molecular types may be negative, but write as positive +------------------------------------------------------------------------- */ + +int AtomVecPeri::pack_restart(int i, double *buf) +{ + int m = 1; + buf[m++] = x[i][0]; + buf[m++] = x[i][1]; + buf[m++] = x[i][2]; + buf[m++] = tag[i]; + buf[m++] = type[i]; + buf[m++] = mask[i]; + buf[m++] = image[i]; + buf[m++] = v[i][0]; + buf[m++] = v[i][1]; + buf[m++] = v[i][2]; + + buf[m++] = rmass[i]; + buf[m++] = vfrac[i]; + buf[m++] = s0[i]; + buf[m++] = x0[i][0]; + buf[m++] = x0[i][1]; + buf[m++] = x0[i][2]; + buf[m++] = vinter[i]; + + if (atom->nextra_restart) + for (int iextra = 0; iextra < atom->nextra_restart; iextra++) + m += modify->fix[atom->extra_restart[iextra]]->pack_restart(i,&buf[m]); + + buf[0] = m; + return m; +} + +/* ---------------------------------------------------------------------- + unpack data for one atom from restart file including extra quantities +------------------------------------------------------------------------- */ + +int AtomVecPeri::unpack_restart(double *buf) +{ + int nlocal = atom->nlocal; + if (nlocal == nmax) { + grow(0); + if (atom->nextra_store) + atom->extra = memory->grow_2d_double_array(atom->extra,nmax, + atom->nextra_store, + "atom:extra"); + } + + int m = 1; + x[nlocal][0] = buf[m++]; + x[nlocal][1] = buf[m++]; + x[nlocal][2] = buf[m++]; + tag[nlocal] = static_cast (buf[m++]); + type[nlocal] = static_cast (buf[m++]); + mask[nlocal] = static_cast (buf[m++]); + image[nlocal] = static_cast (buf[m++]); + v[nlocal][0] = buf[m++]; + v[nlocal][1] = buf[m++]; + v[nlocal][2] = buf[m++]; + + rmass[nlocal] = buf[m++]; + vfrac[nlocal] = buf[m++]; + s0[nlocal] = buf[m++]; + x0[nlocal][0] = buf[m++]; + x0[nlocal][1] = buf[m++]; + x0[nlocal][2] = buf[m++]; + vinter[nlocal] = buf[m++]; + + double **extra = atom->extra; + if (atom->nextra_store) { + int size = static_cast (buf[0]) - m; + for (int i = 0; i < size; i++) extra[nlocal][i] = buf[m++]; + } + + atom->nlocal++; + return m; +} + +/* ---------------------------------------------------------------------- + create one atom of itype at coord + set other values to defaults +------------------------------------------------------------------------- */ + +void AtomVecPeri::create_atom(int itype, double *coord) +{ + int nlocal = atom->nlocal; + if (nlocal == nmax) grow(0); + + tag[nlocal] = 0; + type[nlocal] = itype; + x[nlocal][0] = coord[0]; + x[nlocal][1] = coord[1]; + x[nlocal][2] = coord[2]; + mask[nlocal] = 1; + image[nlocal] = (512 << 20) | (512 << 10) | 512; + v[nlocal][0] = 0.0; + v[nlocal][1] = 0.0; + v[nlocal][2] = 0.0; + + vfrac[nlocal] = 1.0; + rmass[nlocal] = 1.0; + s0[nlocal] = MAXDOUBLE; + x0[nlocal][0] = x[nlocal][0]; + x0[nlocal][1] = x[nlocal][1]; + x0[nlocal][2] = x[nlocal][2]; + vinter[nlocal] = 0.0; + + atom->nlocal++; +} + +/* ---------------------------------------------------------------------- + unpack one line from Atoms section of data file + initialize other atom quantities +------------------------------------------------------------------------- */ + +void AtomVecPeri::data_atom(double *coord, int imagetmp, char **values) +{ + int nlocal = atom->nlocal; + if (nlocal == nmax) grow(0); + + tag[nlocal] = atoi(values[0]); + if (tag[nlocal] <= 0) + error->one("Invalid atom ID in Atoms section of data file"); + + type[nlocal] = atoi(values[1]); + if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes) + error->one("Invalid atom type in Atoms section of data file"); + + vfrac[nlocal] = atof(values[5]); + rmass[nlocal] = atof(values[6]); + s0[nlocal] = MAXDOUBLE; + vinter[nlocal] = 0.0; + + x[nlocal][0] = coord[0]; + x[nlocal][1] = coord[1]; + x[nlocal][2] = coord[2]; + + image[nlocal] = imagetmp; + + mask[nlocal] = 1; + v[nlocal][0] = 0.0; + v[nlocal][1] = 0.0; + v[nlocal][2] = 0.0; + + atom->nlocal++; +} + + +/* ---------------------------------------------------------------------- + unpack hybrid quantities from one line in Atoms section of data file + initialize other atom quantities for this sub-style +------------------------------------------------------------------------- */ + +int AtomVecPeri::data_atom_hybrid(int nlocal, char **values) +{ + vfrac[nlocal] = atof(values[0]); + rmass[nlocal] = atof(values[1]); + + s0[nlocal] = MAXDOUBLE; + vinter[nlocal] = 0.0; + v[nlocal][0] = 0.0; + v[nlocal][1] = 0.0; + v[nlocal][2] = 0.0; + + return 2; +} + +/* ---------------------------------------------------------------------- + return # of bytes of allocated memory +------------------------------------------------------------------------- */ + +double AtomVecPeri::memory_usage() +{ + double bytes = 0.0; + + if (atom->memcheck("tag")) bytes += nmax * sizeof(int); + if (atom->memcheck("type")) bytes += nmax * sizeof(int); + if (atom->memcheck("mask")) bytes += nmax * sizeof(int); + if (atom->memcheck("image")) bytes += nmax * sizeof(int); + if (atom->memcheck("x")) bytes += nmax*3 * sizeof(double); + if (atom->memcheck("v")) bytes += nmax*3 * sizeof(double); + if (atom->memcheck("f")) bytes += nmax*3 * sizeof(double); + + if (atom->memcheck("vfrac")) bytes += nmax * sizeof(double); + if (atom->memcheck("rmass")) bytes += nmax * sizeof(double); + if (atom->memcheck("s0")) bytes += nmax * sizeof(double); + if (atom->memcheck("x0")) bytes += nmax*3 * sizeof(double); + if (atom->memcheck("vinter")) bytes += nmax * sizeof(double); + + return bytes; +} diff --git a/src/PERI/atom_vec_peri.h b/src/PERI/atom_vec_peri.h new file mode 100755 index 0000000000..0c785c2400 --- /dev/null +++ b/src/PERI/atom_vec_peri.h @@ -0,0 +1,54 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + www.cs.sandia.gov/~sjplimp/lammps.html + Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories + + 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. +------------------------------------------------------------------------- */ + +#ifndef ATOM_VEC_PERI_H +#define ATOM_VEC_PERI_H + +#include "atom_vec.h" + +namespace LAMMPS_NS { + +class AtomVecPeri : public AtomVec { + public: + AtomVecPeri(class LAMMPS *, int, char **); + void grow(int); + void copy(int, int); + int pack_comm(int, int *, double *, int, int *); + int pack_comm_one(int, double *); + void unpack_comm(int, int, double *); + int unpack_comm_one(int, double *); + int pack_reverse(int, int, double *); + void unpack_reverse(int, int *, double *); + int pack_border(int, int *, double *, int, int *); + int pack_border_one(int, double *); + void unpack_border(int, int, double *); + int unpack_border_one(int, double *); + int pack_exchange(int, double *); + int unpack_exchange(double *); + int size_restart(); + int pack_restart(int, double *); + int unpack_restart(double *); + void create_atom(int, double *); + void data_atom(double *, int, char **); + int data_atom_hybrid(int, char **); + double memory_usage(); + + private: + int *tag,*type,*mask,*image; + double **x,**v,**f; + double *rmass,*vfrac,*s0,**x0,*vinter; +}; + +} + +#endif diff --git a/src/PERI/compute_damage_atom.cpp b/src/PERI/compute_damage_atom.cpp new file mode 100644 index 0000000000..b3549031dc --- /dev/null +++ b/src/PERI/compute_damage_atom.cpp @@ -0,0 +1,132 @@ +/* ---------------------------------------------------------------------- + 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 "string.h" +#include "compute_damage_atom.h" +#include "atom.h" +#include "modify.h" +#include "comm.h" +#include "force.h" +#include "pair_peri_pmb.h" +#include "fix_peri_neigh.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +ComputeDamageAtom::ComputeDamageAtom(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if (narg != 3) error->all("Illegal compute damage/atom command"); + + peratom_flag = 1; + size_peratom = 0; + + nmax = 0; + damage = NULL; +} + +/* ---------------------------------------------------------------------- */ + +ComputeDamageAtom::~ComputeDamageAtom() +{ + memory->sfree(damage); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeDamageAtom::init() +{ + int count = 0; + for (int i = 0; i < modify->ncompute; i++) + if (strcmp(modify->compute[i]->style,"damage/peri") == 0) count++; + if (count > 1 && comm->me == 0) + error->warning("More than one compute damage/atom"); + + // find associated PERI_NEIGH fix that must exist + + ifix_peri == -1; + 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("Compute damage/atom requires peridynamic potential"); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeDamageAtom::compute_peratom() +{ + // grow damage array if necessary + + if (atom->nlocal > nmax) { + memory->sfree(damage); + nmax = atom->nmax; + damage = (double *) memory->smalloc(nmax*sizeof(double), + "compute/damage/atom:damage"); + scalar_atom = damage; + } + + // compute damage for each atom in group + + int nlocal = atom->nlocal; + int *mask = atom->mask; + double **x = atom->x; + double *vinter = atom->vinter; + double *vfrac = atom->vfrac; + int **partner = ((FixPeriNeigh *) modify->fix[ifix_peri])->partner; + int *npartner = ((FixPeriNeigh *) modify->fix[ifix_peri])->npartner; + int i,j,jj,jnum; + + Pair *anypair = force->pair_match("peri_pmb"); + double **cutsq = anypair->cutsq; + double damage_temp; + + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + jnum = npartner[i]; + damage_temp = 0.0; + + for (jj = 0; jj < jnum; jj++) { + if (partner[i][jj] == 0) continue; + + // look up local index of this partner particle + // skip if particle is "lost" + + j = atom->map(partner[i][jj]); + if (j < 0) continue; + + damage_temp += vfrac[j]; + } + } + else damage_temp = vinter[i]; + + if (vinter[i] != 0.0) damage[i] = 1.0 - damage_temp/vinter[i]; + else damage[i] = 0.0; + } +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based array +------------------------------------------------------------------------- */ + +double ComputeDamageAtom::memory_usage() +{ + double bytes = nmax * sizeof(double); + return bytes; +} diff --git a/src/PERI/compute_damage_atom.h b/src/PERI/compute_damage_atom.h new file mode 100644 index 0000000000..47346549e4 --- /dev/null +++ b/src/PERI/compute_damage_atom.h @@ -0,0 +1,37 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifndef COMPUTE_DAMAGE_ATOM_H +#define COMPUTE_DAMAGE_ATOM_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeDamageAtom : public Compute { + public: + ComputeDamageAtom(class LAMMPS *, int, char **); + ~ComputeDamageAtom(); + void init(); + void compute_peratom(); + double memory_usage(); + + private: + int nmax; + double *damage; + int ifix_peri; +}; + +} + +#endif diff --git a/src/PERI/fix_peri_neigh.cpp b/src/PERI/fix_peri_neigh.cpp new file mode 100644 index 0000000000..509dfab3eb --- /dev/null +++ b/src/PERI/fix_peri_neigh.cpp @@ -0,0 +1,365 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + www.cs.sandia.gov/~sjplimp/lammps.html + Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories + + 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 "fix_peri_neigh.h" +#include "atom.h" +#include "force.h" +#include "pair.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.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) +#define MAX(A,B) ((A) > (B)) ? (A) : (B) + +/* ---------------------------------------------------------------------- */ + +FixPeriNeigh::FixPeriNeigh(LAMMPS *lmp,int narg, char **arg) : + Fix(lmp, narg, arg) +{ + restart_peratom = 1; + first = 1; + + // perform initial allocation of atom-based arrays + // register with atom class + // set maxpartner = 1 as placeholder + + maxpartner = 1; + npartner = NULL; + partner = NULL; + r0 = NULL; + grow_arrays(atom->nmax); + atom->add_callback(0); + atom->add_callback(1); + + // initialize npartner to 0 so atom migration is OK the 1st time + + int nlocal = atom->nlocal; + for (int i = 0; i < nlocal; i++) npartner[i] = 0; +} + +/* ---------------------------------------------------------------------- */ + +FixPeriNeigh::~FixPeriNeigh() +{ + // if atom class still exists: + // unregister this fix so atom class doesn't invoke it any more + + if (atom) atom->delete_callback(id,0); + if (atom) atom->delete_callback(id,1); + + // delete locally stored arrays + + memory->sfree(npartner); + memory->destroy_2d_int_array(partner); + memory->destroy_2d_double_array(r0); +} + +/* ---------------------------------------------------------------------- */ + +int FixPeriNeigh::setmask() +{ + int mask = 0; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixPeriNeigh::init() +{ + // need a full neighbor list once + + int irequest = neighbor->request((void *) this); + neighbor->requests[irequest]->pair = 0; + neighbor->requests[irequest]->fix = 1; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->full = 1; + neighbor->requests[irequest]->occasional = 1; +} + +/* ---------------------------------------------------------------------- */ + +void FixPeriNeigh::init_list(int id, NeighList *ptr) +{ + list = ptr; +} + +/* ---------------------------------------------------------------------- + create initial list of neighbor partners via call to neighbor->build() + must be done in setup (not init) since fix init comes before neigh init +------------------------------------------------------------------------- */ + +void FixPeriNeigh::setup(int vflag) +{ + int i,j,ii,jj,itype,jtype,inum,jnum; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int *ilist,*jlist,*numneigh; + int **firstneigh; + + double **x = atom->x; + double *vfrac = atom->vfrac; + double *vinter = atom->vinter; + int *type = atom->type; + int *tag = atom->tag; + int nlocal = atom->nlocal; + + // only build list of bonds on very first run + + if (!first) return; + first = 0; + + // invoke full neighbor list (will copy or build if necessary) + + neighbor->build_one(list->index); + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // scan neighbor list to set maxpartner + + Pair *anypair = force->pair_match("peri"); + double **cutsq = anypair->cutsq; + + double *s0 = atom->s0; + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + + 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 (rsq < cutsq[itype][jtype]) npartner[i]++; + } + } + + maxpartner = 0; + for (i = 0; i < nlocal; i++) maxpartner = MAX(maxpartner,npartner[i]); + int maxall; + MPI_Allreduce(&maxpartner,&maxall,1,MPI_INT,MPI_MAX,world); + + // realloc arrays stored by fix with correct value for maxpartner + + memory->sfree(npartner); + memory->destroy_2d_int_array(partner); + memory->destroy_2d_double_array(r0); + npartner = NULL; + partner = NULL; + r0 = NULL; + grow_arrays(atom->nmax); + + // create partner list and r0 values from neighbor list + + for (i = 0; i < nlocal; i++) npartner[i] = 0; + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + vinter[i] = 0.0; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + + 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 (rsq <= cutsq[itype][jtype]) { + partner[i][npartner[i]] = tag[j]; + r0[i][npartner[i]] = sqrt(rsq); + npartner[i]++; + vinter[i] += vfrac[j]; + } + } + } + + // bond statistics + + int n = 0; + for (i = 0; i < nlocal; i++) n += npartner[i]; + int nall; + MPI_Allreduce(&n,&nall,1,MPI_INT,MPI_SUM,world); + + if (comm->me == 0) { + if (screen) { + fprintf(screen,"Peridynamic bonds:\n"); + fprintf(screen," total # of bonds = %d\n",nall); + fprintf(screen," bonds/atom = %g\n",nall/atom->natoms); + } + if (logfile) { + fprintf(logfile,"Peridynamic bonds:\n"); + fprintf(logfile," total # of bonds = %d\n",nall); + fprintf(logfile," bonds/atom = %g\n",nall/atom->natoms); + } + } +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based arrays +------------------------------------------------------------------------- */ + +double FixPeriNeigh::memory_usage() +{ + int nmax = atom->nmax; + int bytes = nmax * sizeof(int); + bytes += nmax*maxpartner * sizeof(int); + bytes += nmax*maxpartner * sizeof(double); + return bytes; +} + +/* ---------------------------------------------------------------------- + allocate local atom-based arrays +------------------------------------------------------------------------- */ + +void FixPeriNeigh::grow_arrays(int nmax) +{ + npartner = (int *) memory->srealloc(npartner,nmax*sizeof(int), + "peri_neigh:npartner"); + partner = memory->grow_2d_int_array(partner,nmax,maxpartner, + "peri_neigh:partner"); + r0 = memory->grow_2d_double_array(r0,nmax,maxpartner,"peri_neigh:r0"); +} + +/* ---------------------------------------------------------------------- + copy values within local atom-based arrays +------------------------------------------------------------------------- */ + +void FixPeriNeigh::copy_arrays(int i, int j) +{ + npartner[j] = npartner[i]; + for (int m = 0; m < npartner[j]; m++) { + partner[j][m] = partner[i][m]; + r0[j][m] = r0[i][m]; + } +} + +/* ---------------------------------------------------------------------- + pack values in local atom-based arrays for exchange with another proc +------------------------------------------------------------------------- */ + +int FixPeriNeigh::pack_exchange(int i, double *buf) +{ + // compact list by eliminating partner = 0 entries + // set buf[0] after compaction + + int m = 1; + for (int n = 0; n < npartner[i]; n++) { + if (partner[i][n] == 0) continue; + buf[m++] = partner[i][n]; + buf[m++] = r0[i][n]; + } + buf[0] = m/2; + return m; +} + +/* ---------------------------------------------------------------------- + unpack values in local atom-based arrays from exchange with another proc +------------------------------------------------------------------------- */ + +int FixPeriNeigh::unpack_exchange(int nlocal, double *buf) +{ + int m = 0; + npartner[nlocal] = static_cast (buf[m++]); + for (int n = 0; n < npartner[nlocal]; n++) { + partner[nlocal][n] = static_cast (buf[m++]); + r0[nlocal][n] = buf[m++]; + } + return m; +} + +/* ---------------------------------------------------------------------- + pack values in local atom-based arrays for restart file +------------------------------------------------------------------------- */ + +int FixPeriNeigh::pack_restart(int i, double *buf) +{ + int m = 0; + buf[m++] = 2*npartner[i] + 2; + buf[m++] = npartner[i]; + for (int n = 0; n < npartner[i]; n++) { + buf[m++] = partner[i][n]; + buf[m++] = r0[i][n]; + } + return m; +} + +/* ---------------------------------------------------------------------- + unpack values from atom->extra array to restart the fix +------------------------------------------------------------------------- */ + +void FixPeriNeigh::unpack_restart(int nlocal, int nth) +{ + double **extra = atom->extra; + + // skip to Nth set of extra values + + int m = 0; + for (int i = 0; i < nth; i++) m += static_cast (extra[nlocal][m]); + m++; + + npartner[nlocal] = static_cast (extra[nlocal][m++]); + for (int n = 0; n < npartner[nlocal]; n++) { + partner[nlocal][n] = static_cast (extra[nlocal][m++]); + r0[nlocal][n] = extra[nlocal][m++]; + } +} + +/* ---------------------------------------------------------------------- + maxsize of any atom's restart data +------------------------------------------------------------------------- */ + +int FixPeriNeigh::maxsize_restart() +{ + return 2*maxpartner + 2; +} + +/* ---------------------------------------------------------------------- + size of atom nlocal's restart data +------------------------------------------------------------------------- */ + +int FixPeriNeigh::size_restart(int nlocal) +{ + return 2*npartner[nlocal] + 2; +} diff --git a/src/PERI/fix_peri_neigh.h b/src/PERI/fix_peri_neigh.h new file mode 100644 index 0000000000..85f50e27b2 --- /dev/null +++ b/src/PERI/fix_peri_neigh.h @@ -0,0 +1,55 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + www.cs.sandia.gov/~sjplimp/lammps.html + Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories + + 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. +------------------------------------------------------------------------- */ + +#ifndef FIX_PERI_NEIGH_H +#define FIX_PERI_NEIGH_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixPeriNeigh : public Fix { + friend class PairPeriPMB; + friend class ComputeDamageAtom; + + public: + FixPeriNeigh(class LAMMPS *,int, char **); + ~FixPeriNeigh(); + int setmask(); + void init(); + void init_list(int, class NeighList *); + void setup(int); + + double memory_usage(); + void grow_arrays(int); + void copy_arrays(int, int); + int pack_exchange(int, double *); + int unpack_exchange(int, double *); + int pack_restart(int, double *); + void unpack_restart(int, int); + int size_restart(int); + int maxsize_restart(); + + private: + int first; + int maxpartner; // max # of peridynamic neighs for any atom + int *npartner; // # of neighbors for each atom + int **partner; // neighs for each atom, stored as global IDs + double **r0; // initial distance to partners + + class NeighList *list; +}; + +} + +#endif diff --git a/src/PERI/pair_peri_pmb.cpp b/src/PERI/pair_peri_pmb.cpp new file mode 100644 index 0000000000..3d0b38f049 --- /dev/null +++ b/src/PERI/pair_peri_pmb.cpp @@ -0,0 +1,520 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + www.cs.sandia.gov/~sjplimp/lammps.html + Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories + + 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_pmb.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" + +using namespace LAMMPS_NS; + +#define MIN(a,b) ((a) < (b) ? (a) : (b)) +#define MAX(a,b) ((a) > (b) ? (a) : (b)) + +/* ---------------------------------------------------------------------- */ + +PairPeriPMB::PairPeriPMB(LAMMPS *lmp) : Pair(lmp) +{ + for (int i = 0; i < 6; i++) virial[i] = 0.0; + + ifix_peri = -1; + + nmax = 0; + s0_new = NULL; + + kspring = NULL; + s00 = NULL; + alpha = NULL; + cut = NULL; +} + +/* ---------------------------------------------------------------------- */ + +PairPeriPMB::~PairPeriPMB() +{ + 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(kspring); + memory->destroy_2d_double_array(s00); + memory->destroy_2d_double_array(alpha); + memory->destroy_2d_double_array(cut); + memory->sfree(s0_new); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairPeriPMB::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; + + // 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; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + + 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]; + 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]; + domain->minimum_image(delx0,dely0,delz0); + rsq0 = delx0*delx0 + dely0*dely0 + delz0*delz0; + jtype = type[j]; + + r = sqrt(rsq); + + // shortrange interaction distance based on initial particle position + // 0.9 and 1.35 are constants + + d_ij = MIN(0.9*sqrt(rsq0),1.35*lc); + + // apply short-range contact forces + // 15 is a constant taken from the EMU Theory Manual + // Silling, 12 May 2005, p 18 + + if (r < d_ij) { + dr = r - d_ij; + + rk = (15.0 * kspring[itype][jtype] * 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 = rk*dr; + if (evflag) ev_tally(i,j,nlocal,newton_pair, + evdwl,0.0,fpair,delx,dely,delz); + } + } + } + + // bond forces + + if (atom->nmax > nmax) { + memory->sfree(s0_new); + nmax = atom->nmax; + s0_new = (double *) memory->smalloc(nmax*sizeof(double),"pair:s0_new"); + } + + // first = flag indicating if this is first neighbor of particle i + + bool first; + + // loop over my particles and their partners + // if bond already broken, skip this partner + + for (i = 0; i < nlocal; i++) { + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jnum = npartner[i]; + s0_new[i] = MAXDOUBLE; + first = true; + + for (jj = 0; jj < jnum; jj++) { + if (partner[i][jj] == 0) continue; + j = atom->map(partner[i][jj]); + + // check if we've lost a partner without first breaking bond + + if (j < 0) { + partner[i][jj] = 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]; + domain->minimum_image(delx,dely,delz); + rsq = delx*delx + dely*dely + delz*delz; + jtype = type[j]; + delta = sqrt(cutsq[itype][jtype]); // the horizon + r = sqrt(rsq); + dr = r - r0[i][jj]; + + // avoid roundoff errors + + if (fabs(dr) < 2.2204e-016) dr = 0.0; + + // apply scaling for 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; + + stretch = dr / r0[i][jj]; + rk = (kspring[itype][jtype] * vfrac[j]) * vfrac_scale * stretch; + 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; + + if (eflag) evdwl = rk*dr; + if (evflag) ev_tally(i,j,nlocal,newton_pair, + evdwl,0.0,fbond,delx,dely,delz); + + // find stretch in bond i-j and break if necessary + // use s0 from previous timestep + + 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 neighbor of particle i has now been examined + + first = false; + } + } + + // update with newly computed s0 + + for (i = 0; i < nlocal; i++) s0[i] = s0_new[i]; +} + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairPeriPMB::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"); + kspring = memory->create_2d_double_array(n+1,n+1,"pair:kspring"); + 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 PairPeriPMB::settings(int narg, char **arg) +{ + if (narg) error->all("Illegal pair_style command"); +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairPeriPMB::coeff(int narg, char **arg) +{ + if (narg != 6) 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 kspring_one = atof(arg[2]); + double cut_one = atof(arg[3]); + double s00_one = atof(arg[4]); + double alpha_one = atof(arg[5]); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + kspring[i][j] = kspring_one; + s00[i][j] = s00_one; + alpha[i][j] = alpha_one; + cut[i][j] = cut_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 PairPeriPMB::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 + kspring[j][i] = kspring[i][j]; + alpha[j][i] = alpha[i][j]; + s00[j][i] = s00[i][j]; + + return cut[i][j]; +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairPeriPMB::init_style() +{ + // error checks + + if (atom->map_style == 0) + error->all("Must use atom_modify to define a map"); + + if (atom->style_match("peri") == 0) + error->all("Pair style peri_pmb 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; + + int irequest = neighbor->request(this); +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairPeriPMB::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(&kspring[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 PairPeriPMB::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(&kspring[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(&kspring[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 PairPeriPMB::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 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; + + + delx0 = x0[i][0] - x0[j][0]; + dely0 = x0[i][1] - x0[j][1]; + delz0 = x0[i][2] - x0[j][2]; + 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; + rk = (15.0 * kspring[itype][jtype] * vfrac[j]) * + (dr / sqrt(cutsq[itype][jtype])); + if (r > 0.0) fforce += -(rk/r); + energy += rk*dr; + } + + 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; + 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; + rk = (kspring[itype][jtype] * vfrac[j] * vfrac_scale) * + (dr / r0[i][jj]); + if (r > 0.0) fforce += -(rk/r); + energy += rk*dr; + } + } + + return energy; +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based arrays +------------------------------------------------------------------------- */ + +double PairPeriPMB::memory_usage() +{ + double bytes = nmax * sizeof(double); + return bytes; +} diff --git a/src/PERI/pair_peri_pmb.h b/src/PERI/pair_peri_pmb.h new file mode 100644 index 0000000000..40cbad2667 --- /dev/null +++ b/src/PERI/pair_peri_pmb.h @@ -0,0 +1,51 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + www.cs.sandia.gov/~sjplimp/lammps.html + Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories + + 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. +------------------------------------------------------------------------- */ + +#ifndef PAIR_PERI_PMB_H +#define PAIR_PERI_PMB_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairPeriPMB : public Pair { + public: + PairPeriPMB(class LAMMPS *); + ~PairPeriPMB(); + 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(); + + protected: + int ifix_peri; + double **kspring; + double **s00, **alpha; + double **cut; + + double *s0_new; + int nmax; + + void allocate(); +}; + +} + +#endif diff --git a/src/PERI/style_peri.h b/src/PERI/style_peri.h new file mode 100644 index 0000000000..1d840255a0 --- /dev/null +++ b/src/PERI/style_peri.h @@ -0,0 +1,45 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + www.cs.sandia.gov/~sjplimp/lammps.html + Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories + + 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 AtomInclude +#include "atom_vec_peri.h" +#endif + +#ifdef AtomClass +AtomStyle(peri,AtomVecPeri) +# endif + +#ifdef FixInclude +#include "fix_peri_neigh.h" +#endif + +#ifdef FixClass +FixStyle(PERI_NEIGH,FixPeriNeigh) +#endif + +#ifdef PairInclude +#include "pair_peri_pmb.h" +#endif + +#ifdef PairClass +PairStyle(peri/pmb,PairPeriPMB) +#endif + +#ifdef ComputeInclude +#include "compute_damage_atom.h" +#endif + +#ifdef ComputeClass +ComputeStyle(damage/atom,ComputeDamageAtom) +#endif + diff --git a/src/atom.cpp b/src/atom.cpp index 0c6487357f..f13b2cbc0e 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -64,7 +64,8 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp) q = NULL; mu = NULL; quat = omega = angmom = torque = NULL; - radius = density = rmass = vfrac = NULL; + radius = density = rmass = vfrac = s0 = vinter = NULL; + x0 = NULL; maxspecial = 1; nspecial = NULL; @@ -90,7 +91,7 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp) molecule_flag = 0; q_flag = mu_flag = 0; quat_flag = omega_flag = angmom_flag = torque_flag = 0; - radius_flag = density_flag = rmass_flag = vfrac_flag = 0; + radius_flag = density_flag = rmass_flag = vfrac_flag = s0_flag = vinter_flag = 0; // ntype-length arrays @@ -162,6 +163,9 @@ Atom::~Atom() memory->sfree(density); memory->sfree(rmass); memory->sfree(vfrac); + memory->sfree(s0); + memory->sfree(vinter); + memory->destroy_2d_double_array(x0); memory->sfree(molecule); diff --git a/src/atom.h b/src/atom.h index 658dc35f28..4b9905b933 100644 --- a/src/atom.h +++ b/src/atom.h @@ -48,7 +48,7 @@ class Atom : protected Pointers { int *molecule; double *q,**mu; double **quat,**omega,**angmom,**torque; - double *radius,*density,*rmass,*vfrac; + double *radius,*density,*rmass,*vfrac,*s0,**x0,*vinter; int maxspecial; int **nspecial,**special; @@ -75,7 +75,7 @@ class Atom : protected Pointers { int molecule_flag; int q_flag,mu_flag; int quat_flag,omega_flag,angmom_flag,torque_flag; - int radius_flag,density_flag,rmass_flag,vfrac_flag; + int radius_flag,density_flag,rmass_flag,vfrac_flag,s0_flag,vinter_flag; // extra peratom info in restart file destined for fix & diag diff --git a/src/set.cpp b/src/set.cpp index b9a3872368..12f9557c10 100644 --- a/src/set.cpp +++ b/src/set.cpp @@ -33,13 +33,17 @@ using namespace LAMMPS_NS; enum{ATOM,GROUP,REGION}; enum{TYPE,TYPE_FRACTION,MOLECULE, - X,Y,Z,VX,VY,VZ,CHARGE, - DIPOLE,DIPOLE_RANDOM,QUAT,QUAT_RANDOM, - BOND,ANGLE,DIHEDRAL,IMPROPER}; + X,Y,Z,VX,VY,VZ,CHARGE, + DIPOLE,DIPOLE_RANDOM,QUAT,QUAT_RANDOM, + DIAMETER,DENSITY,VOLUME, + BOND,ANGLE,DIHEDRAL,IMPROPER}; /* ---------------------------------------------------------------------- */ -Set::Set(LAMMPS *lmp) : Pointers(lmp) {} +Set::Set(LAMMPS *lmp) : Pointers(lmp) +{ + PI = 4.0*atan(1.0); +} /* ---------------------------------------------------------------------- */ @@ -171,6 +175,29 @@ void Set::command(int narg, char **arg) if (ivalue <= 0) error->all("Invalid random number seed in set command"); setrandom(QUAT_RANDOM); iarg += 2; + } else if (strcmp(arg[iarg],"diameter") == 0) { + if (iarg+2 > narg) error->all("Illegal set command"); + dvalue = atof(arg[iarg+1]); + if (!atom->radius_flag) + error->all("Cannot set this attribute for this atom style"); + set(DIAMETER); + iarg += 2; + } else if (strcmp(arg[iarg],"density") == 0) { + if (iarg+2 > narg) error->all("Illegal set command"); + dvalue = atof(arg[iarg+1]); + /* + if (!atom->density_flag) + error->all("Cannot set this attribute for this atom style"); + */ + set(DENSITY); + iarg += 2; + } else if (strcmp(arg[iarg],"volume") == 0) { + if (iarg+2 > narg) error->all("Illegal set command"); + dvalue = atof(arg[iarg+1]); + if (!atom->vfrac_flag) + error->all("Cannot set this attribute for this atom style"); + set(VOLUME); + iarg += 2; } else if (strcmp(arg[iarg],"bond") == 0) { if (iarg+2 > narg) error->all("Illegal set command"); ivalue = atoi(arg[iarg+1]); @@ -291,6 +318,26 @@ void Set::set(int keyword) else if (keyword == VY) atom->v[i][1] = dvalue; else if (keyword == VZ) atom->v[i][2] = dvalue; else if (keyword == CHARGE) atom->q[i] = dvalue; + else if (keyword == DIAMETER) { + atom->radius[i] = 0.5 * dvalue; + if (domain->dimension == 3) + atom->rmass[i] = 4.0*PI/3.0 * + atom->radius[i]*atom->radius[i]*atom->radius[i] * atom->density[i]; + else + atom->rmass[i] = PI * + atom->radius[i]*atom->radius[i] * atom->density[i]; + } else if (keyword == DENSITY) { + atom->rmass[i] = dvalue; + /* + atom->density[i] = dvalue; + if (domain->dimension == 3) + atom->rmass[i] = 4.0*PI/3.0 * + atom->radius[i]*atom->radius[i]*atom->radius[i] * atom->density[i]; + else + atom->rmass[i] = PI * + atom->radius[i]*atom->radius[i] * atom->density[i]; + */ + } else if (keyword == VOLUME) atom->vfrac[i] = dvalue; else if (keyword == DIPOLE) { if (atom->dipole[atom->type[i]] > 0.0) { diff --git a/src/set.h b/src/set.h index 1136b8b0b7..d39dfe1428 100644 --- a/src/set.h +++ b/src/set.h @@ -28,6 +28,7 @@ class Set : protected Pointers { int *select; int style,ivalue,newtype,count; double dvalue,xvalue,yvalue,zvalue,wvalue,fraction; + double PI; void selection(int); void set(int); diff --git a/src/style.h b/src/style.h index fe84d72531..245e1c5c81 100644 --- a/src/style.h +++ b/src/style.h @@ -360,6 +360,7 @@ RegionStyle(union,RegUnion) #include "style_meam.h" #include "style_molecule.h" #include "style_opt.h" +#include "style_peri.h" #include "style_poems.h" #include "style_xtc.h" diff --git a/src/update.cpp b/src/update.cpp index 9cbdcc3b94..5abc8acb61 100644 --- a/src/update.cpp +++ b/src/update.cpp @@ -127,6 +127,27 @@ void Update::set_units(const char *style) force->vxmu2f = 0.6241509647; dt = 0.001; neighbor->skin = 2.0; + + } else if (strcmp(style,"si") == 0) { + force->boltz = 1.3806504e-23; + force->mvv2e = 1.0; + force->ftm2v = 1.0; + force->nktv2p = 1.0; + force->qqr2e = 8.9876e9; + force->qe2f = 1.0; + dt = 1.0e-8; + neighbor->skin = 0.001; + + } else if (strcmp(style,"cgs") == 0) { + force->boltz = 1.3806504e-16; + force->mvv2e = 1.0; + force->ftm2v = 1.0; + force->nktv2p = 1.0; + force->qqr2e = 1.0; + force->qe2f = 1.0; + dt = 1.0e-8; + neighbor->skin = 0.1; + } else error->all("Illegal units command"); delete [] unit_style;