From d31eb1a08005d0963ab8538cd74617ad207e957a Mon Sep 17 00:00:00 2001 From: sjplimp Date: Wed, 3 Oct 2007 16:15:56 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@932 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/MANYBODY/pair_airebo.cpp | 112 +++++++++++++++++++++------------- src/MANYBODY/pair_airebo.h | 2 +- src/MANYBODY/pair_eam.cpp | 67 +++++++++++--------- src/MANYBODY/pair_eam.h | 2 +- src/MANYBODY/pair_sw.cpp | 65 ++++++++++++-------- src/MANYBODY/pair_tersoff.cpp | 78 +++++++++++++---------- 6 files changed, 193 insertions(+), 133 deletions(-) diff --git a/src/MANYBODY/pair_airebo.cpp b/src/MANYBODY/pair_airebo.cpp index 1edc70d603..7e37c5f3e1 100644 --- a/src/MANYBODY/pair_airebo.cpp +++ b/src/MANYBODY/pair_airebo.cpp @@ -22,10 +22,13 @@ #include "mpi.h" #include "pair_airebo.h" #include "atom.h" +#include "neighbor.h" +#include "neigh_request.h" #include "force.h" #include "comm.h" -#include "update.h" #include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" #include "memory.h" #include "error.h" @@ -43,8 +46,6 @@ using namespace LAMMPS_NS; PairAIREBO::PairAIREBO(LAMMPS *lmp) : Pair(lmp) { - neigh_half_every = 0; - neigh_full_every = 1; single_enable = 0; one_coeff = 1; @@ -89,9 +90,7 @@ void PairAIREBO::compute(int eflag, int vflag) eng_vdwl = 0.0; if (vflag) for (int i = 0; i < 6; i++) virial[i] = 0.0; - double **f; - if (vflag == 2) f = update->f_pair; - else f = atom->f; + double **f = atom->f; REBO_neigh(); FREBO(eflag,f); @@ -200,6 +199,30 @@ void PairAIREBO::coeff(int narg, char **arg) if (count == 0) error->all("Incorrect args for pair coefficients"); } +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairAIREBO::init_style() +{ + if (atom->tag_enable == 0) + error->all("Pair style AIREBO requires atom IDs"); + if (force->newton_pair == 0) + error->all("Pair style AIREBO requires newton pair on"); + + // need a full neighbor list + + int irequest = neighbor->request(this); + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->full = 1; + + // local REBO neighbor list memory + + pgsize = neighbor->pgsize; + oneatom = neighbor->oneatom; + if (maxlocal == 0) add_pages(0); +} + /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ @@ -262,31 +285,17 @@ double PairAIREBO::init_one(int i, int j) return cutmax; } -/* ---------------------------------------------------------------------- - init specific to this pair style -------------------------------------------------------------------------- */ - -void PairAIREBO::init_style() -{ - if (atom->tag_enable == 0) - error->all("Pair style AIREBO requires atom IDs"); - if (force->newton_pair == 0) - error->all("Pair style AIREBO requires newton pair on"); - - pgsize = neighbor->pgsize; - oneatom = neighbor->oneatom; - if (maxlocal == 0) add_pages(0); -} - /* ---------------------------------------------------------------------- create REBO neighbor list from main neighbor list + REBO neighbor list stores neighbors of ghost atoms ------------------------------------------------------------------------- */ void PairAIREBO::REBO_neigh() { - int i,j,k,n,itype,jtype,m,numneigh; + int i,j,ii,jj,m,n,inum,jnum,itype,jtype; double xtmp,ytmp,ztmp,delx,dely,delz,rsq,dS; - int *neighptr,*neighs; + int *ilist,*jlist,*numneigh,**firstneigh; + int *neighptr; double **x = atom->x; int *type = atom->type; @@ -307,6 +316,11 @@ void PairAIREBO::REBO_neigh() nH = new double[maxlocal]; } + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + // initialize ghost atom references to -1 for (i = nlocal; i < nall; i++) REBO_numneigh[i] = -1; @@ -320,7 +334,9 @@ void PairAIREBO::REBO_neigh() npage = 0; int npnt = 0; - for (i = 0; i < nlocal; i++) { + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + if (pgsize - npnt < oneatom) { npnt = 0; npage++; @@ -334,11 +350,11 @@ void PairAIREBO::REBO_neigh() ztmp = x[i][2]; itype = map[type[i]]; nC[i] = nH[i] = 0.0; - neighs = neighbor->firstneigh_full[i]; - numneigh = neighbor->numneigh_full[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; - for (k = 0; k < numneigh; k++) { - j = neighs[k]; + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; jtype = map[type[j]]; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; @@ -399,11 +415,11 @@ void PairAIREBO::REBO_neigh() else nH[i] += Sp(sqrt(rsq),rcmin[itype][jtype],rcmax[itype][jtype],dS); - neighs = neighbor->firstneigh_full[m]; - numneigh = neighbor->numneigh_full[m]; + jlist = firstneigh[i]; + jnum = numneigh[i]; - for (k = 0; k < numneigh; k++) { - j = neighs[k]; + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; if (j == i) continue; jtype = map[type[j]]; delx = xtmp - x[j][0]; @@ -503,14 +519,15 @@ void PairAIREBO::FREBO(int eflag, double **f) void PairAIREBO::FLJ(int eflag, double **f) { - int i,j,k,m,jj,kk,mm,itype,jtype,ktype,mtype; - int numneigh,atomi,atomj,atomk,atomm; + int i,j,k,m,ii,jj,kk,mm,inum,jnum,itype,jtype,ktype,mtype; + int atomi,atomj,atomk,atomm; int testpath,npath,done; double rsq,best,wik,wkm,cij,rij,dwij,dwik,dwkj,dwkm,dwmj; double delij[3],rijsq,delik[3],rik,delkj[3]; double rkj,wkj,dC,VLJ,dVLJ,fforce,VA,Str,dStr,Stb; double delkm[3],rkm,delmj[3],rmj,wmj,r2inv,r6inv,scale,delscale[3]; - int *neighs,*REBO_neighs_i,*REBO_neighs_k; + int *ilist,*jlist,*numneigh,**firstneigh; + int *REBO_neighs_i,*REBO_neighs_k; double delikS[3],delkjS[3],delkmS[3],delmjS[3]; double rikS,rkjS,rkmS,rmjS,wikS,dwikS; double wkjS,dwkjS,wkmS,dwkmS,wmjS,dwmjS; @@ -523,14 +540,23 @@ void PairAIREBO::FLJ(int eflag, double **f) int *type = atom->type; int nlocal = atom->nlocal; - for (i = 0; i < nlocal; i++) { + 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]; itype = map[type[i]]; - neighs = neighbor->firstneigh_full[i]; - numneigh = neighbor->numneigh_full[i]; - atomi=i; - for (jj = 0; jj < numneigh; jj++) { - j = neighs[jj]; - atomj=j; + atomi = i; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + atomj = j; if (tag[i] > tag[j]) continue; jtype = map[type[j]]; diff --git a/src/MANYBODY/pair_airebo.h b/src/MANYBODY/pair_airebo.h index 64cb4a50bf..c3a7accc48 100644 --- a/src/MANYBODY/pair_airebo.h +++ b/src/MANYBODY/pair_airebo.h @@ -25,8 +25,8 @@ class PairAIREBO : public Pair { void compute(int, int); void settings(int, char **); void coeff(int, char **); - double init_one(int, int); void init_style(); + double init_one(int, int); void write_restart(FILE *) {} void read_restart(FILE *) {} void write_restart_settings(FILE *) {} diff --git a/src/MANYBODY/pair_eam.cpp b/src/MANYBODY/pair_eam.cpp index 760167f8fb..8177be2d25 100644 --- a/src/MANYBODY/pair_eam.cpp +++ b/src/MANYBODY/pair_eam.cpp @@ -22,9 +22,9 @@ #include "pair_eam.h" #include "atom.h" #include "force.h" -#include "update.h" #include "comm.h" #include "neighbor.h" +#include "neigh_list.h" #include "memory.h" #include "error.h" @@ -124,12 +124,11 @@ PairEAM::~PairEAM() void PairEAM::compute(int eflag, int vflag) { - int i,j,k,m,numneigh,itype,jtype; + int i,j,ii,jj,m,inum,jnum,itype,jtype; double xtmp,ytmp,ztmp,delx,dely,delz; double rsq,r,p,fforce,rhoip,rhojp,z2,z2p,recip,phi,phip,psip; double *coeff; - int *neighs; - double **f; + int *ilist,*jlist,*numneigh,**firstneigh; // grow energy array if necessary @@ -144,13 +143,17 @@ void PairEAM::compute(int eflag, int vflag) eng_vdwl = 0.0; if (vflag) for (i = 0; i < 6; i++) virial[i] = 0.0; - if (vflag == 2) f = update->f_pair; - else f = atom->f; double **x = atom->x; + double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; int newton_pair = force->newton_pair; + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + // zero out density if (newton_pair) { @@ -161,16 +164,17 @@ void PairEAM::compute(int eflag, int vflag) // rho = density at each atom // loop over neighbors of my atoms - for (i = 0; i < nlocal; i++) { + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; itype = type[i]; - neighs = neighbor->firstneigh[i]; - numneigh = neighbor->numneigh[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; - for (k = 0; k < numneigh; k++) { - j = neighs[k]; + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; @@ -201,7 +205,8 @@ void PairEAM::compute(int eflag, int vflag) // fp = derivative of embedding energy at each atom // phi = embedding energy at each atom - for (i = 0; i < nlocal; i++) { + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; p = rho[i]*rdrho + 1.0; m = static_cast (p); m = MAX(1,MIN(m,nrho-1)); @@ -219,16 +224,18 @@ void PairEAM::compute(int eflag, int vflag) // compute forces on each atom // loop over neighbors of my atoms - for (i = 0; i < nlocal; i++) { + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; itype = type[i]; - neighs = neighbor->firstneigh[i]; - numneigh = neighbor->numneigh[i]; - for (k = 0; k < numneigh; k++) { - j = neighs[k]; + 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]; @@ -382,6 +389,20 @@ void PairEAM::coeff(int narg, char **arg) if (count == 0) error->all("Incorrect args for pair coefficients"); } +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairEAM::init_style() +{ + // convert read-in file(s) to arrays and spline them + + file2array(); + array2spline(); + + int irequest = neighbor->request(this); +} + /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ @@ -404,18 +425,6 @@ double PairEAM::init_one(int i, int j) return cutmax; } -/* ---------------------------------------------------------------------- - init specific to this pair style -------------------------------------------------------------------------- */ - -void PairEAM::init_style() -{ - // convert read-in file(s) to arrays and spline them - - file2array(); - array2spline(); -} - /* ---------------------------------------------------------------------- read potential values from a DYNAMO single element funcfl file ------------------------------------------------------------------------- */ diff --git a/src/MANYBODY/pair_eam.h b/src/MANYBODY/pair_eam.h index 29cebe433a..bcf187a646 100644 --- a/src/MANYBODY/pair_eam.h +++ b/src/MANYBODY/pair_eam.h @@ -25,8 +25,8 @@ class PairEAM : public Pair { void compute(int, int); void settings(int, char **); virtual void coeff(int, char **); - double init_one(int, int); void init_style(); + double init_one(int, int); void single(int, int, int, int, double, double, double, int, One &); void single_embed(int, int, double &); diff --git a/src/MANYBODY/pair_sw.cpp b/src/MANYBODY/pair_sw.cpp index 857064231d..0ee132407c 100755 --- a/src/MANYBODY/pair_sw.cpp +++ b/src/MANYBODY/pair_sw.cpp @@ -21,11 +21,13 @@ #include "string.h" #include "pair_sw.h" #include "atom.h" +#include "neighbor.h" +#include "neigh_request.h" #include "force.h" #include "comm.h" -#include "update.h" #include "memory.h" #include "neighbor.h" +#include "neigh_list.h" #include "memory.h" #include "error.h" @@ -38,8 +40,6 @@ using namespace LAMMPS_NS; PairSW::PairSW(LAMMPS *lmp) : Pair(lmp) { - neigh_half_every = 0; - neigh_full_every = 1; single_enable = 0; one_coeff = 1; @@ -74,23 +74,26 @@ PairSW::~PairSW() void PairSW::compute(int eflag, int vflag) { - int i,j,k,m,n,itag,jtag,itype,jtype,ktype,iparam,numneigh; + int i,j,k,ii,jj,kk,inum,jnum,jnumm1,itag,jtag,itype,jtype,ktype,iparam; double xtmp,ytmp,ztmp,delx,dely,delz; double rsq,rsq1,rsq2,eng,fforce; double delr1[3],delr2[3],fj[3],fk[3]; - int *neighs; - double **f; + int *ilist,*jlist,*numneigh,**firstneigh; eng_vdwl = 0.0; if (vflag) for (i = 0; i < 6; i++) virial[i] = 0.0; - if (vflag == 2) f = update->f_pair; - else f = atom->f; double **x = atom->x; + double **f = atom->f; int *tag = atom->tag; int *type = atom->type; int nlocal = atom->nlocal; + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + // loop over full neighbor list of my atoms for (i = 0; i < nlocal; i++) { @@ -102,11 +105,11 @@ void PairSW::compute(int eflag, int vflag) // two-body interactions, skip half of them - neighs = neighbor->firstneigh_full[i]; - numneigh = neighbor->numneigh_full[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; - for (m = 0; m < numneigh; m++) { - j = neighs[m]; + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; jtag = tag[j]; if (itag > jtag) { @@ -145,12 +148,14 @@ void PairSW::compute(int eflag, int vflag) // cannot test I-J distance against cutoff outside of 2nd loop // b/c must use I-J-K cutoff for both rij and rik - for (m = 0; m < numneigh-1; m++) { - j = neighs[m]; + jnumm1 = jnum - 1; + + for (jj = 0; jj < jnumm1; jj++) { + j = jlist[jj]; jtype = map[type[j]]; - for (n = m+1; n < numneigh; n++) { - k = neighs[n]; + for (kk = jj+1; kk < jnum; kk++) { + k = jlist[kk]; ktype = map[type[k]]; iparam = elem2param[itype][jtype][ktype]; @@ -278,6 +283,24 @@ void PairSW::coeff(int narg, char **arg) if (count == 0) error->all("Incorrect args for pair coefficients"); } +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairSW::init_style() +{ + if (atom->tag_enable == 0) + error->all("Pair style Stillinger-Weber requires atom IDs"); + if (force->newton_pair == 0) + error->all("Pair style Stillinger-Weber requires newton pair on"); + + // need a full neighbor list + + int irequest = neighbor->request(this); + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->full = 1; +} + /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ @@ -291,16 +314,6 @@ double PairSW::init_one(int i, int j) /* ---------------------------------------------------------------------- */ -void PairSW::init_style() -{ - if (atom->tag_enable == 0) - error->all("Pair style Stillinger-Weber requires atom IDs"); - if (force->newton_pair == 0) - error->all("Pair style Stillinger-Weber requires newton pair on"); -} - -/* ---------------------------------------------------------------------- */ - void PairSW::read_file(char *file) { int params_per_line = 13; diff --git a/src/MANYBODY/pair_tersoff.cpp b/src/MANYBODY/pair_tersoff.cpp index 453928b7c4..3ac5061855 100755 --- a/src/MANYBODY/pair_tersoff.cpp +++ b/src/MANYBODY/pair_tersoff.cpp @@ -21,11 +21,11 @@ #include "string.h" #include "pair_tersoff.h" #include "atom.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" #include "force.h" #include "comm.h" -#include "update.h" -#include "memory.h" -#include "neighbor.h" #include "memory.h" #include "error.h" @@ -38,8 +38,6 @@ using namespace LAMMPS_NS; PairTersoff::PairTersoff(LAMMPS *lmp) : Pair(lmp) { - neigh_half_every = 0; - neigh_full_every = 1; single_enable = 0; one_coeff = 1; @@ -77,27 +75,32 @@ PairTersoff::~PairTersoff() void PairTersoff::compute(int eflag, int vflag) { - int i,j,k,m,n,itag,jtag,itype,jtype,ktype,iparam_ij,iparam_ijk,numneigh; + int i,j,k,ii,jj,kk,inum,jnum; + int itag,jtag,itype,jtype,ktype,iparam_ij,iparam_ijk; double xtmp,ytmp,ztmp,delx,dely,delz; double rsq,rsq1,rsq2,eng,fforce; double delr1[3],delr2[3],fi[3],fj[3],fk[3]; double zeta_ij,prefactor; - int *neighs; - double **f; + int *ilist,*jlist,*numneigh,**firstneigh; eng_vdwl = 0.0; if (vflag) for (i = 0; i < 6; i++) virial[i] = 0.0; - if (vflag == 2) f = update->f_pair; - else f = atom->f; double **x = atom->x; + double **f = atom->f; int *tag = atom->tag; int *type = atom->type; int nlocal = atom->nlocal; + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + // loop over full neighbor list of my atoms - for (i = 0; i < nlocal; i++) { + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; itag = tag[i]; itype = map[type[i]]; xtmp = x[i][0]; @@ -106,11 +109,11 @@ void PairTersoff::compute(int eflag, int vflag) // two-body interactions, skip half of them - neighs = neighbor->firstneigh_full[i]; - numneigh = neighbor->numneigh_full[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; - for (m = 0; m < numneigh; m++) { - j = neighs[m]; + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; jtag = tag[j]; if (itag > jtag) { @@ -148,8 +151,8 @@ void PairTersoff::compute(int eflag, int vflag) // three-body interactions // skip immediately if I-J is not within cutoff - for (m = 0; m < numneigh; m++) { - j = neighs[m]; + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; jtype = map[type[j]]; iparam_ij = elem2param[itype][jtype][jtype]; @@ -162,9 +165,10 @@ void PairTersoff::compute(int eflag, int vflag) // accumulate bondorder zeta for each i-j interaction via loop over k zeta_ij = 0.0; - for (n = 0; n < numneigh; n++) { - if (n == m) continue; - k = neighs[n]; + + for (kk = 0; kk < jnum; kk++) { + if (jj == kk) continue; + k = jlist[kk]; ktype = map[type[k]]; iparam_ijk = elem2param[itype][jtype][ktype]; @@ -191,9 +195,9 @@ void PairTersoff::compute(int eflag, int vflag) // attractive term via loop over k - for (n = 0; n < numneigh; n++) { - if (n == m) continue; - k = neighs[n]; + for (kk = 0; kk < jnum; kk++) { + if (jj == kk) continue; + k = jlist[kk]; ktype = map[type[k]]; iparam_ijk = elem2param[itype][jtype][ktype]; @@ -315,6 +319,24 @@ void PairTersoff::coeff(int narg, char **arg) if (count == 0) error->all("Incorrect args for pair coefficients"); } +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairTersoff::init_style() +{ + if (atom->tag_enable == 0) + error->all("Pair style Tersoff requires atom IDs"); + if (force->newton_pair == 0) + error->all("Pair style Tersoff requires newton pair on"); + + // need a full neighbor list + + int irequest = neighbor->request(this); + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->full = 1; +} + /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ @@ -328,16 +350,6 @@ double PairTersoff::init_one(int i, int j) /* ---------------------------------------------------------------------- */ -void PairTersoff::init_style() -{ - if (atom->tag_enable == 0) - error->all("Pair style Tersoff requires atom IDs"); - if (force->newton_pair == 0) - error->all("Pair style Tersoff requires newton pair on"); -} - -/* ---------------------------------------------------------------------- */ - void PairTersoff::read_file(char *file) { int params_per_line = 15;