Commit 60908eea authored by sjplimp's avatar sjplimp Committed by GitHub
Browse files

Merge pull request #246 from akohlmey/manybody-short-neighbor-list 

Manybody short neighbor list
parents 8214555b 02510ec3

src/MANYBODY/pair_sw.cpp

+49 −29
Original line number Diff line number Diff line
@@ -51,6 +51,9 @@ PairSW::PairSW(LAMMPS *lmp) : Pair(lmp) 
  params = NULL;

  elem2param = NULL;

  map = NULL;



  maxshort = 10;

  neighshort = NULL;

}



/* ----------------------------------------------------------------------
@@ -70,6 +73,7 @@ PairSW::~PairSW() 
  if (allocated) {

    memory->destroy(setflag);

    memory->destroy(cutsq);

    memory->destroy(neighshort);

    delete [] map;

  }

}
@@ -102,6 +106,8 @@ void PairSW::compute(int eflag, int vflag) 
  numneigh = list->numneigh;

  firstneigh = list->firstneigh;



  double fxtmp,fytmp,fztmp;



  // loop over full neighbor list of my atoms



  for (ii = 0; ii < inum; ii++) {
@@ -111,17 +117,36 @@ void PairSW::compute(int eflag, int vflag) 
    xtmp = x[i][0];

    ytmp = x[i][1];

    ztmp = x[i][2];

    fxtmp = fytmp = fztmp = 0.0;



    // two-body interactions, skip half of them



    jlist = firstneigh[i];

    jnum = numneigh[i];

    int numshort = 0;



    for (jj = 0; jj < jnum; jj++) {

      j = jlist[jj];

      j &= NEIGHMASK;

      jtag = tag[j];



      delx = xtmp - x[j][0];

      dely = ytmp - x[j][1];

      delz = ztmp - x[j][2];

      rsq = delx*delx + dely*dely + delz*delz;



      jtype = map[type[j]];

      ijparam = elem2param[itype][jtype][jtype];

      if (rsq >= params[ijparam].cutsq) {

        continue;

      } else {

        neighshort[numshort++] = j;

        if (numshort >= maxshort) {

          maxshort += maxshort/2;

          memory->grow(neighshort,maxshort,"pair:neighshort");

        }

      }



      jtag = tag[j];

      if (itag > jtag) {

        if ((itag+jtag) % 2 == 0) continue;

      } else if (itag < jtag) {
@@ -132,21 +157,11 @@ void PairSW::compute(int eflag, int vflag) 
        if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue;

      }



      jtype = map[type[j]];



      delx = xtmp - x[j][0];

      dely = ytmp - x[j][1];

      delz = ztmp - x[j][2];

      rsq = delx*delx + dely*dely + delz*delz;



      ijparam = elem2param[itype][jtype][jtype];

      if (rsq >= params[ijparam].cutsq) continue;



      twobody(&params[ijparam],rsq,fpair,eflag,evdwl);



      f[i][0] += delx*fpair;

      f[i][1] += dely*fpair;

      f[i][2] += delz*fpair;

      fxtmp += delx*fpair;

      fytmp += dely*fpair;

      fztmp += delz*fpair;

      f[j][0] -= delx*fpair;

      f[j][1] -= dely*fpair;

      f[j][2] -= delz*fpair;
@@ -155,22 +170,22 @@ void PairSW::compute(int eflag, int vflag) 
                           evdwl,0.0,fpair,delx,dely,delz);

    }



    jnumm1 = jnum - 1;

    jnumm1 = numshort - 1;



    for (jj = 0; jj < jnumm1; jj++) {

      j = jlist[jj];

      j &= NEIGHMASK;

      j = neighshort[jj];

      jtype = map[type[j]];

      ijparam = elem2param[itype][jtype][jtype];

      delr1[0] = x[j][0] - xtmp;

      delr1[1] = x[j][1] - ytmp;

      delr1[2] = x[j][2] - ztmp;

      rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2];

      if (rsq1 >= params[ijparam].cutsq) continue;



      for (kk = jj+1; kk < jnum; kk++) {

        k = jlist[kk];

        k &= NEIGHMASK;

      double fjxtmp,fjytmp,fjztmp;

      fjxtmp = fjytmp = fjztmp = 0.0;



      for (kk = jj+1; kk < numshort; kk++) {

        k = neighshort[kk];

        ktype = map[type[k]];

        ikparam = elem2param[itype][ktype][ktype];

        ijkparam = elem2param[itype][jtype][ktype];
@@ -179,24 +194,29 @@ void PairSW::compute(int eflag, int vflag) 
        delr2[1] = x[k][1] - ytmp;

        delr2[2] = x[k][2] - ztmp;

        rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];

        if (rsq2 >= params[ikparam].cutsq) continue;



        threebody(&params[ijparam],&params[ikparam],&params[ijkparam],

                  rsq1,rsq2,delr1,delr2,fj,fk,eflag,evdwl);



        f[i][0] -= fj[0] + fk[0];

        f[i][1] -= fj[1] + fk[1];

        f[i][2] -= fj[2] + fk[2];

        f[j][0] += fj[0];

        f[j][1] += fj[1];

        f[j][2] += fj[2];

        fxtmp -= fj[0] + fk[0];

        fytmp -= fj[1] + fk[1];

        fztmp -= fj[2] + fk[2];

        fjxtmp += fj[0];

        fjytmp += fj[1];

        fjztmp += fj[2];

        f[k][0] += fk[0];

        f[k][1] += fk[1];

        f[k][2] += fk[2];



        if (evflag) ev_tally3(i,j,k,evdwl,0.0,fj,fk,delr1,delr2);

      }

      f[j][0] += fjxtmp;

      f[j][1] += fjytmp;

      f[j][2] += fjztmp;

    }

    f[i][0] += fxtmp;

    f[i][1] += fytmp;

    f[i][2] += fztmp;

  }



  if (vflag_fdotr) virial_fdotr_compute();
@@ -211,7 +231,7 @@ void PairSW::allocate() 


  memory->create(setflag,n+1,n+1,"pair:setflag");

  memory->create(cutsq,n+1,n+1,"pair:cutsq");



  memory->create(neighshort,maxshort,"pair:neighshort");

  map = new int[n+1];

}

src/MANYBODY/pair_sw.h

+2 −0
Original line number Diff line number Diff line
@@ -55,6 +55,8 @@ class PairSW : public Pair { 
  int nparams;                  // # of stored parameter sets

  int maxparam;                 // max # of parameter sets

  Param *params;                // parameter set for an I-J-K interaction

  int maxshort;                 // size of short neighbor list array

  int *neighshort;              // short neighbor list array



  virtual void allocate();

  void read_file(char *);

src/MANYBODY/pair_tersoff.cpp

+57 −36
Original line number Diff line number Diff line
@@ -52,6 +52,9 @@ PairTersoff::PairTersoff(LAMMPS *lmp) : Pair(lmp) 
  params = NULL;

  elem2param = NULL;

  map = NULL;



  maxshort = 10;

  neighshort = NULL;

}



/* ----------------------------------------------------------------------
@@ -71,6 +74,7 @@ PairTersoff::~PairTersoff() 
  if (allocated) {

    memory->destroy(setflag);

    memory->destroy(cutsq);

    memory->destroy(neighshort);

    delete [] map;

  }

}
@@ -98,12 +102,15 @@ void PairTersoff::compute(int eflag, int vflag) 
  int *type = atom->type;

  int nlocal = atom->nlocal;

  int newton_pair = force->newton_pair;

  const double cutshortsq = cutmax*cutmax;



  inum = list->inum;

  ilist = list->ilist;

  numneigh = list->numneigh;

  firstneigh = list->firstneigh;



  double fxtmp,fytmp,fztmp;



  // loop over full neighbor list of my atoms



  for (ii = 0; ii < inum; ii++) {
@@ -113,17 +120,32 @@ void PairTersoff::compute(int eflag, int vflag) 
    xtmp = x[i][0];

    ytmp = x[i][1];

    ztmp = x[i][2];

    fxtmp = fytmp = fztmp = 0.0;



    // two-body interactions, skip half of them



    jlist = firstneigh[i];

    jnum = numneigh[i];

    int numshort = 0;



    for (jj = 0; jj < jnum; jj++) {

      j = jlist[jj];

      j &= NEIGHMASK;

      jtag = tag[j];



      delx = xtmp - x[j][0];

      dely = ytmp - x[j][1];

      delz = ztmp - x[j][2];

      rsq = delx*delx + dely*dely + delz*delz;



      if (rsq < cutshortsq) {

        neighshort[numshort++] = j;

        if (numshort >= maxshort) {

          maxshort += maxshort/2;

          memory->grow(neighshort,maxshort,"pair:neighshort");

        }

      }



      jtag = tag[j];

      if (itag > jtag) {

        if ((itag+jtag) % 2 == 0) continue;

      } else if (itag < jtag) {
@@ -135,20 +157,14 @@ void PairTersoff::compute(int eflag, int vflag) 
      }



      jtype = map[type[j]];



      delx = xtmp - x[j][0];

      dely = ytmp - x[j][1];

      delz = ztmp - x[j][2];

      rsq = delx*delx + dely*dely + delz*delz;



      iparam_ij = elem2param[itype][jtype][jtype];

      if (rsq > params[iparam_ij].cutsq) continue;

      if (rsq >= params[iparam_ij].cutsq) continue;



      repulsive(&params[iparam_ij],rsq,fpair,eflag,evdwl);



      f[i][0] += delx*fpair;

      f[i][1] += dely*fpair;

      f[i][2] += delz*fpair;

      fxtmp += delx*fpair;

      fytmp += dely*fpair;

      fztmp += delz*fpair;

      f[j][0] -= delx*fpair;

      f[j][1] -= dely*fpair;

      f[j][2] -= delz*fpair;
@@ -159,10 +175,10 @@ void PairTersoff::compute(int eflag, int vflag) 


    // three-body interactions

    // skip immediately if I-J is not within cutoff

    double fjxtmp,fjytmp,fjztmp;



    for (jj = 0; jj < jnum; jj++) {

      j = jlist[jj];

      j &= NEIGHMASK;

    for (jj = 0; jj < numshort; jj++) {

      j = neighshort[jj];

      jtype = map[type[j]];

      iparam_ij = elem2param[itype][jtype][jtype];
@@ -170,16 +186,16 @@ void PairTersoff::compute(int eflag, int vflag) 
      delr1[1] = x[j][1] - ytmp;

      delr1[2] = x[j][2] - ztmp;

      rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2];

      if (rsq1 > params[iparam_ij].cutsq) continue;

      if (rsq1 >= params[iparam_ij].cutsq) continue;



      // accumulate bondorder zeta for each i-j interaction via loop over k



      fjxtmp = fjytmp = fjztmp = 0.0;

      zeta_ij = 0.0;



      for (kk = 0; kk < jnum; kk++) {

      for (kk = 0; kk < numshort; kk++) {

        if (jj == kk) continue;

        k = jlist[kk];

        k &= NEIGHMASK;

        k = neighshort[kk];

        ktype = map[type[k]];

        iparam_ijk = elem2param[itype][jtype][ktype];
@@ -187,7 +203,7 @@ void PairTersoff::compute(int eflag, int vflag) 
        delr2[1] = x[k][1] - ytmp;

        delr2[2] = x[k][2] - ztmp;

        rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];

        if (rsq2 > params[iparam_ijk].cutsq) continue;

        if (rsq2 >= params[iparam_ijk].cutsq) continue;



        zeta_ij += zeta(&params[iparam_ijk],rsq1,rsq2,delr1,delr2);

      }
@@ -196,22 +212,21 @@ void PairTersoff::compute(int eflag, int vflag) 


      force_zeta(&params[iparam_ij],rsq1,zeta_ij,fpair,prefactor,eflag,evdwl);



      f[i][0] += delr1[0]*fpair;

      f[i][1] += delr1[1]*fpair;

      f[i][2] += delr1[2]*fpair;

      f[j][0] -= delr1[0]*fpair;

      f[j][1] -= delr1[1]*fpair;

      f[j][2] -= delr1[2]*fpair;

      fxtmp += delr1[0]*fpair;

      fytmp += delr1[1]*fpair;

      fztmp += delr1[2]*fpair;

      fjxtmp -= delr1[0]*fpair;

      fjytmp -= delr1[1]*fpair;

      fjztmp -= delr1[2]*fpair;



      if (evflag) ev_tally(i,j,nlocal,newton_pair,

                           evdwl,0.0,-fpair,-delr1[0],-delr1[1],-delr1[2]);



      // attractive term via loop over k



      for (kk = 0; kk < jnum; kk++) {

      for (kk = 0; kk < numshort; kk++) {

        if (jj == kk) continue;

        k = jlist[kk];

        k &= NEIGHMASK;

        k = neighshort[kk];

        ktype = map[type[k]];

        iparam_ijk = elem2param[itype][jtype][ktype];
@@ -219,24 +234,30 @@ void PairTersoff::compute(int eflag, int vflag) 
        delr2[1] = x[k][1] - ytmp;

        delr2[2] = x[k][2] - ztmp;

        rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];

        if (rsq2 > params[iparam_ijk].cutsq) continue;

        if (rsq2 >= params[iparam_ijk].cutsq) continue;



        attractive(&params[iparam_ijk],prefactor,

                   rsq1,rsq2,delr1,delr2,fi,fj,fk);



        f[i][0] += fi[0];

        f[i][1] += fi[1];

        f[i][2] += fi[2];

        f[j][0] += fj[0];

        f[j][1] += fj[1];

        f[j][2] += fj[2];

        fxtmp += fi[0];

        fytmp += fi[1];

        fztmp += fi[2];

        fjxtmp += fj[0];

        fjytmp += fj[1];

        fjztmp += fj[2];

        f[k][0] += fk[0];

        f[k][1] += fk[1];

        f[k][2] += fk[2];



        if (vflag_atom) v_tally3(i,j,k,fj,fk,delr1,delr2);

      }

      f[j][0] += fjxtmp;

      f[j][1] += fjytmp;

      f[j][2] += fjztmp;

    }

    f[i][0] += fxtmp;

    f[i][1] += fytmp;

    f[i][2] += fztmp;

  }



  if (vflag_fdotr) virial_fdotr_compute();
@@ -251,7 +272,7 @@ void PairTersoff::allocate() 


  memory->create(setflag,n+1,n+1,"pair:setflag");

  memory->create(cutsq,n+1,n+1,"pair:cutsq");



  memory->create(neighshort,maxshort,"pair:neighshort");

  map = new int[n+1];

}

src/MANYBODY/pair_tersoff.h

+2 −0
Original line number Diff line number Diff line
@@ -59,6 +59,8 @@ class PairTersoff : public Pair { 
  int nelements;                // # of unique elements

  int nparams;                  // # of stored parameter sets

  int maxparam;                 // max # of parameter sets

  int maxshort;                 // size of short neighbor list array

  int *neighshort;              // short neighbor list array



  virtual void allocate();

  virtual void read_file(char *);

src/MANYBODY/pair_vashishta.cpp

+53 −26
Original line number Diff line number Diff line
@@ -52,6 +52,10 @@ PairVashishta::PairVashishta(LAMMPS *lmp) : Pair(lmp) 
  params = NULL;

  elem2param = NULL;

  map = NULL;



  r0max = 0.0;

  maxshort = 10;

  neighshort = NULL;

}



/* ----------------------------------------------------------------------
@@ -69,6 +73,7 @@ PairVashishta::~PairVashishta() 
  if (allocated) {

    memory->destroy(setflag);

    memory->destroy(cutsq);

    memory->destroy(neighshort);

    delete [] map;

  }

}
@@ -95,12 +100,15 @@ void PairVashishta::compute(int eflag, int vflag) 
  int *type = atom->type;

  int nlocal = atom->nlocal;

  int newton_pair = force->newton_pair;

  const double cutshortsq = r0max*r0max;



  inum = list->inum;

  ilist = list->ilist;

  numneigh = list->numneigh;

  firstneigh = list->firstneigh;



  double fxtmp,fytmp,fztmp;



  // loop over full neighbor list of my atoms



  for (ii = 0; ii < inum; ii++) {
@@ -110,17 +118,32 @@ void PairVashishta::compute(int eflag, int vflag) 
    xtmp = x[i][0];

    ytmp = x[i][1];

    ztmp = x[i][2];

    fxtmp = fytmp = fztmp = 0.0;



    // two-body interactions, skip half of them



    jlist = firstneigh[i];

    jnum = numneigh[i];

    int numshort = 0;



    for (jj = 0; jj < jnum; jj++) {

      j = jlist[jj];

      j &= NEIGHMASK;

      jtag = tag[j];



      delx = xtmp - x[j][0];

      dely = ytmp - x[j][1];

      delz = ztmp - x[j][2];

      rsq = delx*delx + dely*dely + delz*delz;



      if (rsq < cutshortsq) {

        neighshort[numshort++] = j;

        if (numshort >= maxshort) {

          maxshort += maxshort/2;

          memory->grow(neighshort,maxshort,"pair:neighshort");

        }

      }



      jtag = tag[j];

      if (itag > jtag) {

        if ((itag+jtag) % 2 == 0) continue;

      } else if (itag < jtag) {
@@ -132,20 +155,14 @@ void PairVashishta::compute(int eflag, int vflag) 
      }



      jtype = map[type[j]];



      delx = xtmp - x[j][0];

      dely = ytmp - x[j][1];

      delz = ztmp - x[j][2];

      rsq = delx*delx + dely*dely + delz*delz;



      ijparam = elem2param[itype][jtype][jtype];

      if (rsq > params[ijparam].cutsq) continue;

      if (rsq >= params[ijparam].cutsq) continue;



      twobody(&params[ijparam],rsq,fpair,eflag,evdwl);



      f[i][0] += delx*fpair;

      f[i][1] += dely*fpair;

      f[i][2] += delz*fpair;

      fxtmp += delx*fpair;

      fytmp += dely*fpair;

      fztmp += delz*fpair;

      f[j][0] -= delx*fpair;

      f[j][1] -= dely*fpair;

      f[j][2] -= delz*fpair;
@@ -154,11 +171,10 @@ void PairVashishta::compute(int eflag, int vflag) 
      			   evdwl,0.0,fpair,delx,dely,delz);

    }



    jnumm1 = jnum - 1;

    jnumm1 = numshort - 1;



    for (jj = 0; jj < jnumm1; jj++) {

      j = jlist[jj];

      j &= NEIGHMASK;

      j = neighshort[jj];

      jtype = map[type[j]];

      ijparam = elem2param[itype][jtype][jtype];

      delr1[0] = x[j][0] - xtmp;
@@ -167,9 +183,11 @@ void PairVashishta::compute(int eflag, int vflag) 
      rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2];

      if (rsq1 >= params[ijparam].cutsq2) continue;



      for (kk = jj+1; kk < jnum; kk++) {

        k = jlist[kk];

        k &= NEIGHMASK;

      double fjxtmp,fjytmp,fjztmp;

      fjxtmp = fjytmp = fjztmp = 0.0;



      for (kk = jj+1; kk < numshort; kk++) {

        k = neighshort[kk];

        ktype = map[type[k]];

        ikparam = elem2param[itype][ktype][ktype];

        ijkparam = elem2param[itype][jtype][ktype];
@@ -183,19 +201,25 @@ void PairVashishta::compute(int eflag, int vflag) 
        threebody(&params[ijparam],&params[ikparam],&params[ijkparam],

                  rsq1,rsq2,delr1,delr2,fj,fk,eflag,evdwl);



        f[i][0] -= fj[0] + fk[0];

        f[i][1] -= fj[1] + fk[1];

        f[i][2] -= fj[2] + fk[2];

        f[j][0] += fj[0];

        f[j][1] += fj[1];

        f[j][2] += fj[2];

        fxtmp -= fj[0] + fk[0];

        fytmp -= fj[1] + fk[1];

        fztmp -= fj[2] + fk[2];

        fjxtmp += fj[0];

        fjytmp += fj[1];

        fjztmp += fj[2];

        f[k][0] += fk[0];

        f[k][1] += fk[1];

        f[k][2] += fk[2];



        if (evflag) ev_tally3(i,j,k,evdwl,0.0,fj,fk,delr1,delr2);

      }

      f[j][0] += fjxtmp;

      f[j][1] += fjytmp;

      f[j][2] += fjztmp;

    }

    f[i][0] += fxtmp;

    f[i][1] += fytmp;

    f[i][2] += fztmp;

  }



  if (vflag_fdotr) virial_fdotr_compute();
@@ -210,6 +234,7 @@ void PairVashishta::allocate() 


  memory->create(setflag,n+1,n+1,"pair:setflag");

  memory->create(cutsq,n+1,n+1,"pair:cutsq");

  memory->create(neighshort,maxshort,"pair:neighshort");



  map = new int[n+1];

}
@@ -529,13 +554,15 @@ void PairVashishta::setup_params() 
    params[m].c0 = params[m].cut*params[m].dvrc - params[m].vrc;

  }



  // set cutmax to max of all params

  // set cutmax to max of all cutoff params. r0max only for r0



  cutmax = 0.0;

  r0max = 0.0;

  for (m = 0; m < nparams; m++) {

    if (params[m].cut > cutmax) cutmax = params[m].cut;

    if (params[m].r0 > cutmax) cutmax = params[m].r0;

    if (params[m].r0 > r0max) r0max = params[m].r0;

  }

  if (r0max > cutmax) cutmax = r0max;

}



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

CRA Git | Maintained and supported by SUSTech CRA and CCSE