Commit d31eb1a0 authored by sjplimp's avatar sjplimp
Browse files

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@932 f3b2605a-c512-4ea7-a41b-209d697bcdaa
parent bc528eda
Loading
Loading
Loading
Loading
+69 −43
Original line number Diff line number Diff line
@@ -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,13 +540,22 @@ 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];
    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]];
+1 −1
Original line number Diff line number Diff line
@@ -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 *) {}
+38 −29
Original line number Diff line number Diff line
@@ -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<int> (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
------------------------------------------------------------------------- */
+1 −1
Original line number Diff line number Diff line
@@ -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 &);
+39 −26
Original line number Diff line number Diff line
@@ -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];

@@ -279,24 +284,32 @@ void PairSW::coeff(int narg, char **arg)
}

/* ----------------------------------------------------------------------
   init for one type pair i,j and corresponding j,i
   init specific to this pair style
------------------------------------------------------------------------- */

double PairSW::init_one(int i, int j)
{
  if (setflag[i][j] == 0) error->all("All pair coeffs are not set");

  return cutmax;
}

/* ---------------------------------------------------------------------- */

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
------------------------------------------------------------------------- */

double PairSW::init_one(int i, int j)
{
  if (setflag[i][j] == 0) error->all("All pair coeffs are not set");

  return cutmax;
}

/* ---------------------------------------------------------------------- */
Loading