Commit 82d9f5f5 authored by Anders Hafreager's avatar Anders Hafreager Committed by Axel Kohlmeyer
Browse files

Added 3-body neighbor list building for faster short range 3 body forces.

parent 944ebdcf
Loading
Loading
Loading
Loading
+35 −12
Original line number Diff line number Diff line
@@ -54,6 +54,10 @@ PairVashishta::PairVashishta(LAMMPS *lmp) : Pair(lmp)
  elem2param = NULL;
  map = NULL;

  neigh3BodyMax = 0;
  neigh3BodyCount = NULL; 
  neigh3Body = NULL;
  
  useTable = true;
  tableSize = 65536;
  forceTable = NULL;
@@ -78,8 +82,7 @@ PairVashishta::~PairVashishta()
    delete [] map;
  }
}
#include <iostream>
using namespace std;

void PairVashishta::modify_params(int narg, char **arg)
{
  if(narg < 2 || narg > 3) // We accept table yes/no [tableSize]
@@ -102,6 +105,16 @@ void PairVashishta::modify_params(int narg, char **arg)
  }
}

void PairVashishta::validateNeigh3Body() {
  if (atom->nlocal > neigh3BodyMax) {
        neigh3BodyMax = atom->nmax;
        memory->destroy(neigh3BodyCount);
        memory->destroy(neigh3Body);
        memory->create(neigh3BodyCount,neigh3BodyMax,  "pair:vashishta:neigh3BodyCount");
        memory->create(neigh3Body,neigh3BodyMax, 1000, "pair:vashishta:neigh3Body"); // TODO: pick this number more wisely? 
    }
}

void PairVashishta::createTable()
{
  int ntypes = atom->ntypes+1;
@@ -140,6 +153,7 @@ void PairVashishta::createTable()

void PairVashishta::compute(int eflag, int vflag)
{
  validateNeigh3Body();
  int i,j,k,ii,jj,kk,inum,jnum,jnumm1;
  int itype,jtype,ktype,ijparam,ikparam,ijkparam;
  tagint itag,jtag;
@@ -174,6 +188,8 @@ void PairVashishta::compute(int eflag, int vflag)
    ytmp = x[i][1];
    ztmp = x[i][2];

    neigh3BodyCount[i] = 0; // Reset the 3-body neighbor list

    // two-body interactions, skip half of them

    jlist = firstneigh[i];
@@ -184,6 +200,21 @@ void PairVashishta::compute(int eflag, int vflag)
      j &= NEIGHMASK;
      jtag = tag[j];

      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].cutsq2) {
          neigh3Body[i][neigh3BodyCount[i]] = j;
          neigh3BodyCount[i]++;
      }

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

      if (itag > jtag) {
        if ((itag+jtag) % 2 == 0) continue;
      } else if (itag < jtag) {
@@ -194,16 +225,6 @@ void PairVashishta::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, useTable);

      f[i][0] += delx*fpair;
@@ -217,6 +238,8 @@ void PairVashishta::compute(int eflag, int vflag)
      			   evdwl,0.0,fpair,delx,dely,delz);
    }

    jlist = neigh3Body[i];
    jnum = neigh3BodyCount[i];
    jnumm1 = jnum - 1;

    for (jj = 0; jj < jnumm1; jj++) {
+5 −0
Original line number Diff line number Diff line
@@ -53,6 +53,10 @@ class PairVashishta : public Pair {
  double ***forceTable;          // Table containing the forces
  double ***potentialTable;      // Table containing the potential energies

  int neigh3BodyMax;            // Max size of short neighborlist
  int *neigh3BodyCount;         // Number of neighbors in short range 3 particle forces neighbor list
  int **neigh3Body;             // Neighborlist for short range 3 particle forces

  double cutmax;                // max cutoff for all elements
  int nelements;                // # of unique elements
  char **elements;              // names of unique elements
@@ -66,6 +70,7 @@ class PairVashishta : public Pair {
  void read_file(char *);
  void setup_params();
  void createTable();
  void validateNeigh3Body();
  void twobody(Param *, double, double &, int, double &, bool);
  void threebody(Param *, Param *, Param *, double, double, double *, double *,
                 double *, double *, int, double &);
+21 −10
Original line number Diff line number Diff line
@@ -36,6 +36,8 @@ PairVashishtaOMP::PairVashishtaOMP(LAMMPS *lmp) :

void PairVashishtaOMP::compute(int eflag, int vflag)
{
  validateNeigh3Body();

  if (eflag || vflag) {
    ev_setup(eflag,vflag);
  } else evflag = vflag_fdotr = 0;
@@ -105,6 +107,8 @@ void PairVashishtaOMP::eval(int iifrom, int iito, ThrData * const thr)
    ztmp = x[i].z;
    fxtmp = fytmp = fztmp = 0.0;

    neigh3BodyCount[i] = 0;
    
    // two-body interactions, skip half of them

    jlist = firstneigh[i];
@@ -114,6 +118,21 @@ void PairVashishtaOMP::eval(int iifrom, int iito, ThrData * const thr)
      j = jlist[jj];
      j &= NEIGHMASK;
      jtag = tag[j];
      jtype = map[type[j]];

      delx = xtmp - x[j].x;
      dely = ytmp - x[j].y;
      delz = ztmp - x[j].z;
      rsq = delx*delx + dely*dely + delz*delz;

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

      if (rsq <= params[ijparam].cutsq2) {
          neigh3Body[i][neigh3BodyCount[i]] = j;
          neigh3BodyCount[i]++;
      }

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

      if (itag > jtag) {
        if ((itag+jtag) % 2 == 0) continue;
@@ -125,16 +144,6 @@ void PairVashishtaOMP::eval(int iifrom, int iito, ThrData * const thr)
        if (x[j].z == ztmp && x[j].y == ytmp && x[j].x < xtmp) continue;
      }

      jtype = map[type[j]];

      delx = xtmp - x[j].x;
      dely = ytmp - x[j].y;
      delz = ztmp - x[j].z;
      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, useTable);

      fxtmp += delx*fpair;
@@ -148,6 +157,8 @@ void PairVashishtaOMP::eval(int iifrom, int iito, ThrData * const thr)
                               evdwl,0.0,fpair,delx,dely,delz,thr);
    }

    jlist = neigh3Body[i];
    jnum = neigh3BodyCount[i];
    jnumm1 = jnum - 1;

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