jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@1965 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2008-07-24 20:46:29 +00:00
parent 0c081ffa4b
commit 192008e7c7
17 changed files with 1987 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/Makefile
View File

@ -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

32
src/PERI/Install.csh Normal file
View File

@ -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

613
src/PERI/atom_vec_peri.cpp Normal file
View File

@ -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<int> (buf[m++]);
type[i] = static_cast<int> (buf[m++]);
mask[i] = static_cast<int> (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<int> (buf[m++]);
type[nlocal] = static_cast<int> (buf[m++]);
mask[nlocal] = static_cast<int> (buf[m++]);
image[nlocal] = static_cast<int> (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<int> (buf[m++]);
type[nlocal] = static_cast<int> (buf[m++]);
mask[nlocal] = static_cast<int> (buf[m++]);
image[nlocal] = static_cast<int> (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<int> (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;
}

54
src/PERI/atom_vec_peri.h Executable file
View File

@ -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

132
src/PERI/compute_damage_atom.cpp Normal file
View File

@ -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;
}

37
src/PERI/compute_damage_atom.h Normal file
View File

@ -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

365
src/PERI/fix_peri_neigh.cpp Normal file
View File

@ -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<int> (buf[m++]);
for (int n = 0; n < npartner[nlocal]; n++) {
partner[nlocal][n] = static_cast<int> (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<int> (extra[nlocal][m]);
m++;
npartner[nlocal] = static_cast<int> (extra[nlocal][m++]);
for (int n = 0; n < npartner[nlocal]; n++) {
partner[nlocal][n] = static_cast<int> (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;
}

55
src/PERI/fix_peri_neigh.h Normal file
View File

@ -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

520
src/PERI/pair_peri_pmb.cpp Normal file
View File

@ -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;
}

51
src/PERI/pair_peri_pmb.h Normal file
View File

@ -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

45
src/PERI/style_peri.h Normal file
View File

@ -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

8
src/atom.cpp
View File

@ -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);

4
src/atom.h
View File

@ -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

55
src/set.cpp
View File

@ -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) {

1
src/set.h
View File

@ -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);

1
src/style.h
View File

@ -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"

21
src/update.cpp
View File

@ -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;