Unverified Commit e4d4f3a7 authored by Axel Kohlmeyer's avatar Axel Kohlmeyer Committed by GitHub
Browse files

Merge pull request #1132 from martok/meam-init

Small fixes to MEAM/C
parents de010551 e5ddc909
Loading
Loading
Loading
Loading
+2 −3
Original line number Diff line number Diff line
@@ -101,7 +101,7 @@ MEAM::getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double**
  double xjtmp, yjtmp, zjtmp, delxik, delyik, delzik, rik2 /*,rik*/;
  double xktmp, yktmp, zktmp, delxjk, delyjk, delzjk, rjk2 /*,rjk*/;
  double xik, xjk, sij, fcij, sfcij, dfcij, sikj, dfikj, cikj;
  double Cmin, Cmax, delc, /*ebound,*/ rbound, a, coef1, coef2;
  double Cmin, Cmax, delc, /*ebound,*/ a, coef1, coef2;
  double dCikj;
  double rnorm, fc, dfc, drinv;

@@ -129,6 +129,7 @@ MEAM::getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double**
    rij2 = delxij * delxij + delyij * delyij + delzij * delzij;
    rij = sqrt(rij2);

    const double rbound = this->ebound_meam[elti][eltj] * rij2;
    if (rij > this->rc_meam) {
      fcij = 0.0;
      dfcij = 0.0;
@@ -138,7 +139,6 @@ MEAM::getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double**
      sij = 1.0;

      //     if rjk2 > ebound*rijsq, atom k is definitely outside the ellipse
      const double rbound = this->ebound_meam[elti][eltj] * rij2;
      for (kn = 0; kn < numneigh_full; kn++) {
        k = firstneigh_full[kn];
        eltk = fmap[type[k]];
@@ -193,7 +193,6 @@ MEAM::getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double**
    if (iszero(sfcij) || iszero(sfcij - 1.0))
      goto LABEL_100;

    rbound = this->ebound_meam[elti][eltj] * rij2;
    for (kn = 0; kn < numneigh_full; kn++) {
      k = firstneigh_full[kn];
      if (k == j) continue;
+1 −1
Original line number Diff line number Diff line
@@ -87,7 +87,7 @@ MEAM::dG_gam(const double gamma, const int ibar, double& dG) const
        //         e.g. gsmooth_factor is 99, {:
        //         gsmooth_switchpoint = -0.99
        //         G = 0.01*(-0.99/gamma)**99
        double G = 1 / (gsmooth_factor + 1) * pow((gsmooth_switchpoint / gamma), gsmooth_factor);
        G = 1 / (gsmooth_factor + 1) * pow((gsmooth_switchpoint / gamma), gsmooth_factor);
        G = sqrt(G);
        dG = -gsmooth_factor * G / (2.0 * gamma);
        return G;
+13 −0
Original line number Diff line number Diff line
@@ -34,6 +34,19 @@ MEAM::MEAM(Memory* mem)

  maxneigh = 0;
  scrfcn = dscrfcn = fcpair = NULL;
  
  neltypes = 0;
  for (int i = 0; i < maxelt; i++) {
    Omega_meam[i] = Z_meam[i] = A_meam[i] = rho0_meam[i] = beta0_meam[i] = 
      beta1_meam[i]= beta2_meam[i] = beta3_meam[i] = 
      t0_meam[i] = t1_meam[i] = t2_meam[i] = t3_meam[i] = 
      rho_ref_meam[i] = ibar_meam[i] = ielt_meam[i] = 0.0;
    for (int j = 0; j < maxelt; j++) {
      lattce_meam[i][j] = FCC;
      Ec_meam[i][j] = re_meam[i][j] = alpha_meam[i][j] = delta_meam[i][j] = Ec_meam[i][j] = ebound_meam[i][j] = attrac_meam[i][j] = repuls_meam[i][j] = 0.0;
      nn2_meam[i][j] = zbl_meam[i][j] = eltind[i][j] = 0;
    }
  }
}

MEAM::~MEAM()
+5 −4
Original line number Diff line number Diff line
@@ -206,7 +206,7 @@ void PairMEAMC::settings(int narg, char **/*arg*/)

void PairMEAMC::coeff(int narg, char **arg)
{
  int i,j,m,n;
  int m,n;

  if (!allocated) allocate();

@@ -222,7 +222,7 @@ void PairMEAMC::coeff(int narg, char **arg)
  // elements = list of unique element names

  if (nelements) {
    for (i = 0; i < nelements; i++) delete [] elements[i];
    for (int i = 0; i < nelements; i++) delete [] elements[i];
    delete [] elements;
    delete [] mass;
  }
@@ -231,7 +231,7 @@ void PairMEAMC::coeff(int narg, char **arg)
  elements = new char*[nelements];
  mass = new double[nelements];

  for (i = 0; i < nelements; i++) {
  for (int i = 0; i < nelements; i++) {
    n = strlen(arg[i+3]) + 1;
    elements[i] = new char[n];
    strcpy(elements[i],arg[i+3]);
@@ -247,8 +247,9 @@ void PairMEAMC::coeff(int narg, char **arg)
  // read args that map atom types to MEAM elements
  // map[i] = which element the Ith atom type is, -1 if not mapped

  for (i = 4 + nelements; i < narg; i++) {
  for (int i = 4 + nelements; i < narg; i++) {
    m = i - (4+nelements) + 1;
    int j;
    for (j = 0; j < nelements; j++)
      if (strcmp(arg[i],elements[j]) == 0) break;
    if (j < nelements) map[m] = j;