Commit 2c656594 authored by Sebastian Hütter's avatar Sebastian Hütter
Browse files

MEAM/C: implement scaling factor for reversible scaling calculations

parent 204529bc

src/USER-MEAMC/meam.h

+2 −2
Original line number Diff line number Diff line
@@ -247,9 +247,9 @@ public: 
  void meam_dens_init(int i, int ntype, int* type, int* fmap, double** x, int numneigh, int* firstneigh,

                      int numneigh_full, int* firstneigh_full, int fnoffset);

  void meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_atom, double* eng_vdwl,

                       double* eatom, int ntype, int* type, int* fmap, int& errorflag);

                       double* eatom, int ntype, int* type, int* fmap, double** scale, int& errorflag);

  void meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int vflag_atom, double* eng_vdwl,

                  double* eatom, int ntype, int* type, int* fmap, double** x, int numneigh, int* firstneigh,

                  double* eatom, int ntype, int* type, int* fmap, double** scale, double** x, int numneigh, int* firstneigh,

                  int numneigh_full, int* firstneigh_full, int fnoffset, double** f, double** vatom);

};

src/USER-MEAMC/meam_dens_final.cpp

+4 −1
Original line number Diff line number Diff line
@@ -4,18 +4,20 @@ using namespace LAMMPS_NS; 


void

MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_atom, double* eng_vdwl,

                      double* eatom, int /*ntype*/, int* type, int* fmap, int& errorflag)

                      double* eatom, int /*ntype*/, int* type, int* fmap, double** scale, int& errorflag)

{

  int i, elti;

  int m;

  double rhob, G, dG, Gbar, dGbar, gam, shp[3], Z;

  double denom, rho_bkgd, Fl;

  double scaleii;



  //     Complete the calculation of density



  for (i = 0; i < nlocal; i++) {

    elti = fmap[type[i]];

    if (elti >= 0) {

      scaleii = scale[type[i]][type[i]];

      rho1[i] = 0.0;

      rho2[i] = -1.0 / 3.0 * arho2b[i] * arho2b[i];

      rho3[i] = 0.0;
@@ -113,6 +115,7 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_ 
      Fl = embedding(this->A_meam[elti], this->Ec_meam[elti][elti], rhob, frhop[i]);



      if (eflag_either != 0) {

        Fl *= scaleii;

        if (eflag_global != 0) {

          *eng_vdwl = *eng_vdwl + Fl;

        }

src/USER-MEAMC/meam_force.cpp

+14 −4
Original line number Diff line number Diff line
@@ -8,7 +8,7 @@ using namespace LAMMPS_NS; 


void

MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int vflag_atom, double* eng_vdwl,

                 double* eatom, int /*ntype*/, int* type, int* fmap, double** x, int numneigh, int* firstneigh,

                 double* eatom, int /*ntype*/, int* type, int* fmap, double** scale, double** x, int numneigh, int* firstneigh,

                 int numneigh_full, int* firstneigh_full, int fnoffset, double** f, double** vatom)

{

  int j, jn, k, kn, kk, m, n, p, q;
@@ -42,6 +42,7 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int 
  double arg1i1, arg1j1, arg1i2, arg1j2, arg1i3, arg1j3, arg3i3, arg3j3;

  double dsij1, dsij2, force1, force2;

  double t1i, t2i, t3i, t1j, t2j, t3j;

  double scaleij;



  third = 1.0 / 3.0;

  sixth = 1.0 / 6.0;
@@ -59,6 +60,7 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int 
  for (jn = 0; jn < numneigh; jn++) {

    j = firstneigh[jn];

    eltj = fmap[type[j]];

    scaleij = scale[type[i]][type[j]];



    if (!iszero(scrfcn[fnoffset + jn]) && eltj >= 0) {
@@ -84,11 +86,12 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int 
        recip = 1.0 / r;



        if (eflag_either != 0) {

          double phi_sc = phi * scaleij;

          if (eflag_global != 0)

            *eng_vdwl = *eng_vdwl + phi * sij;

            *eng_vdwl = *eng_vdwl + phi_sc * sij;

          if (eflag_atom != 0) {

            eatom[i] = eatom[i] + 0.5 * phi * sij;

            eatom[j] = eatom[j] + 0.5 * phi * sij;

            eatom[i] = eatom[i] + 0.5 * phi_sc * sij;

            eatom[j] = eatom[j] + 0.5 * phi_sc * sij;

          }

        }
@@ -379,6 +382,13 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int 
        for (m = 0; m < 3; m++) {

          dUdrijm[m] = frhop[i] * drhodrm1[m] + frhop[j] * drhodrm2[m];

        }

        if (!isone(scaleij)) {

          dUdrij *= scaleij;

          dUdsij *= scaleij;

          dUdrijm[0] *= scaleij;

          dUdrijm[1] *= scaleij;

          dUdrijm[2] *= scaleij;

        }



        //     Add the part of the force due to dUdrij and dUdsij

src/USER-MEAMC/pair_meamc.cpp

+22 −5
Original line number Diff line number Diff line
@@ -56,6 +56,7 @@ PairMEAMC::PairMEAMC(LAMMPS *lmp) : Pair(lmp) 
  elements = NULL;

  mass = NULL;

  meam_inst = new MEAM(memory);

  scale = NULL;



  // set comm size needed by this Pair
@@ -79,6 +80,7 @@ PairMEAMC::~PairMEAMC() 
  if (allocated) {

    memory->destroy(setflag);

    memory->destroy(cutsq);

    memory->destroy(scale);

    delete [] map;

  }

}
@@ -143,7 +145,7 @@ void PairMEAMC::compute(int eflag, int vflag) 
  comm->reverse_comm_pair(this);



  meam_inst->meam_dens_final(nlocal,eflag_either,eflag_global,eflag_atom,

                   &eng_vdwl,eatom,ntype,type,map,errorflag);

                   &eng_vdwl,eatom,ntype,type,map,scale,errorflag);

  if (errorflag) {

    char str[128];

    sprintf(str,"MEAM library error %d",errorflag);
@@ -164,7 +166,7 @@ void PairMEAMC::compute(int eflag, int vflag) 
  for (ii = 0; ii < inum_half; ii++) {

    i = ilist_half[ii];

    meam_inst->meam_force(i,eflag_either,eflag_global,eflag_atom,

                vflag_atom,&eng_vdwl,eatom,ntype,type,map,x,

                vflag_atom,&eng_vdwl,eatom,ntype,type,map,scale,x,

                numneigh_half[i],firstneigh_half[i],

                numneigh_full[i],firstneigh_full[i],

                offset,f,vptr);
@@ -183,6 +185,7 @@ void PairMEAMC::allocate() 


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

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

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



  map = new int[n+1];

}
@@ -267,13 +270,16 @@ void PairMEAMC::coeff(int narg, char **arg) 
  // set mass for i,i in atom class



  int count = 0;

  for (int i = 1; i <= n; i++)

    for (int j = i; j <= n; j++)

  for (int i = 1; i <= n; i++) {

    for (int j = i; j <= n; j++) {

      if (map[i] >= 0 && map[j] >= 0) {

        setflag[i][j] = 1;

        if (i == j) atom->set_mass(FLERR,i,mass[map[i]]);

        count++;

      }

      scale[i][j] = 1.0;

    }

  }



  if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");

}
@@ -312,8 +318,10 @@ void PairMEAMC::init_list(int id, NeighList *ptr) 
   init for one type pair i,j and corresponding j,i

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



double PairMEAMC::init_one(int /*i*/, int /*j*/)

double PairMEAMC::init_one(int i, int j)

{

  if (setflag[i][j] == 0) scale[i][j] = 1.0;

  scale[j][i] = scale[i][j];

  return cutmax;

}
@@ -773,3 +781,12 @@ void PairMEAMC::neigh_strip(int inum, int *ilist, 
    for (j = 0; j < jnum; j++) jlist[j] &= NEIGHMASK;

  }

}



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



void *PairMEAMC::extract(const char *str, int &dim)

{

  dim = 2;

  if (strcmp(str,"scale") == 0) return (void *) scale;

  return NULL;

}

src/USER-MEAMC/pair_meamc.h

+2 −0
Original line number Diff line number Diff line
@@ -35,6 +35,7 @@ class PairMEAMC : public Pair { 
  void init_style();

  void init_list(int, class NeighList *);

  double init_one(int, int);

  virtual void *extract(const char *, int &);



  int pack_forward_comm(int, int *, double *, int, int *);

  void unpack_forward_comm(int, int, double *);
@@ -50,6 +51,7 @@ class PairMEAMC : public Pair { 
  double *mass;                 // mass of each element



  int *map;                     // mapping from atom types (1-indexed) to elements (1-indexed)

  double **scale;               // scaling factor for adapt



  void allocate();

  void read_files(char *, char *);

CRA Git | Maintained and supported by SUSTech CRA and CCSE