whitespace and formatting update

using namespace LAMMPS_NS;

// This is for debugging purposes. The ASSERT() macro is used in the code to check

// if everything runs as expected. Change this to #if 0 if you don't need the checking.

#if 0

        #define ASSERT(cond) ((!(cond)) ? my_failure(error,__FILE__,__LINE__) : my_noop())

        inline void my_noop() {}

        inline void my_failure(Error* error, const char* file, int line) {

                char str[1024];

                sprintf(str,"Assertion failure: File %s, line %i", file, line);

                error->one(FLERR,str);

        }

#else

#define ASSERT(cond)

#endif

#define MAXLINE 1024        // This sets the maximum line length in EAM input files.

PairEAMCD::PairEAMCD(LAMMPS *lmp, int _cdeamVersion) : PairEAM(lmp), PairEAMAlloy(lmp), cdeamVersion(_cdeamVersion)

PairEAMCD::PairEAMCD(LAMMPS *lmp, int _cdeamVersion)

  : PairEAM(lmp), PairEAMAlloy(lmp), cdeamVersion(_cdeamVersion)

{

  single_enable = 0;

  restartinfo = 0;

  hcoeff = NULL;

  // Set communication buffer sizes needed by this pair style.

  if (cdeamVersion == 1) {

    comm_forward = 4;

    comm_reverse = 3;

        }

        else if(cdeamVersion == 2) {

  } else if (cdeamVersion == 2) {

    comm_forward = 3;

    comm_reverse = 2;

        }

        else {

                error->all(FLERR,"Invalid CD-EAM potential version.");

  } else {

    error->all(FLERR,"Invalid eam/cd potential version.");

  }

}

  else evflag = vflag_fdotr = eflag_global = eflag_atom = 0;

  // Grow per-atom arrays if necessary

  if (atom->nmax > nmax) {

    memory->destroy(rho);

    memory->destroy(fp);

  firstneigh = list->firstneigh;

  // Zero out per-atom arrays.

  int m = nlocal + atom->nghost;

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

    rho[i] = 0.0;

  // Stage I

  // Compute rho and rhoB at each local atom site.

        // Additionally calculate the D_i values here if we are using the one-site formulation.

        // For the two-site formulation we have to calculate the D values in an extra loop (Stage II).

  // Additionally calculate the D_i values here if we are using the

  // one-site formulation.  For the two-site formulation we have to

  // calculate the D values in an extra loop (Stage II).

  for (ii = 0; ii < inum; ii++) {

    i = ilist[ii];

    xtmp = x[i][0];

        }

        if (cdeamVersion == 1 && itype != jtype) {

          // Note: if the i-j interaction is not concentration dependent (because either

          // i or j are not species A or B) then its contribution to D_i and D_j should

          // be ignored.

          // This if-clause is only required for a ternary.

                                        if((itype == speciesA && jtype == speciesB) || (jtype == speciesA && itype == speciesB)) {

          if ((itype == speciesA && jtype == speciesB)

              || (jtype == speciesA && itype == speciesB)) {

            double Phi_AB = PhiOfR(index, itype, jtype, 1.0 / r);

            D_values[i] += Phi_AB;

            if (newton_pair || j < nlocal)

  }

  // Communicate and sum densities.

  if (newton_pair) {

    communicationStage = 1;

    comm->reverse_comm_pair(this);

  // fp = derivative of embedding energy at each atom

  // phi = embedding energy at each atom

  for (ii = 0; ii < inum; ii++) {

    i = ilist[ii];

    EAMTableIndex index = rhoToTableIndex(rho[i]);

  // Communicate derivative of embedding function and densities

  // and D_values (this for one-site formulation only).

  communicationStage = 2;

  comm->forward_comm_pair(this);

        // The electron densities may not drop to zero because then the concentration would no longer be defined.

        // But the concentration is not needed anyway if there is no interaction with another atom, which is the case

        // if the electron density is exactly zero. That's why the following lines have been commented out.

  // The electron densities may not drop to zero because then the

  // concentration would no longer be defined.  But the concentration

  // is not needed anyway if there is no interaction with another atom,

  // which is the case if the electron density is exactly zero.

  // That's why the following lines have been commented out.

  //

  //for (i = 0; i < nlocal + atom->nghost; i++) {

  //        if (rho[i] == 0 && (type[i] == speciesA || type[i] == speciesB))

  // This is only required for the original two-site formulation of the CD-EAM potential.

  if (cdeamVersion == 2) {

    // Compute intermediate value D_i for each atom.

    for (ii = 0; ii < inum; ii++) {

      i = ilist[ii];

      xtmp = x[i][0];

      jnum = numneigh[i];

      // This code line is required for ternary alloys.

      if (itype != speciesA && itype != speciesB) continue;

      double x_i = rhoB[i] / rho[i];        // Concentration at atom i.

        if (itype == jtype) continue;

        // This code line is required for ternary alloys.

        if (jtype != speciesA && jtype != speciesB) continue;

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

          const EAMTableIndex index = radiusToTableIndex(r);

          // The concentration independent part of the cross pair potential.

          double Phi_AB = PhiOfR(index, itype, jtype, 1.0 / r);

          // Average concentration of two sites

          double x_ij = 0.5 * (x_i + rhoB[j]/rho[j]);

          // Calculate derivative of h(x_ij) polynomial function.

          double h_prime = evalHprime(x_ij);

          D_values[i] += h_prime * Phi_AB / (2.0 * rho[i] * rho[i]);

    }

    // Communicate and sum D values.

    if (newton_pair) {

      communicationStage = 3;

      comm->reverse_comm_pair(this);

  // Stage III

  // Compute force acting on each atom.

  for (ii = 0; ii < inum; ii++) {

    i = ilist[ii];

    xtmp = x[i][0];

    jnum = numneigh[i];

    // Concentration at site i

                double x_i = -1.0;                // The value -1 indicates: no concentration dependence for all interactions of atom i.

    // The value -1 indicates: no concentration dependence for all interactions of atom i.

    // It will be replaced by the concentration at site i if atom i is either A or B.

    double x_i = -1.0;

    double D_i, h_prime_i;

    // This if-clause is only required for ternary alloys.

    if ((itype == speciesA || itype == speciesB) && rho[i] != 0.0) {

      // Compute local concentration at site i.

      x_i = rhoB[i]/rho[i];

      ASSERT(x_i >= 0 && x_i<=1.0);

      if (cdeamVersion == 1) {

        // Calculate derivative of h(x_i) polynomial function.

        h_prime_i = evalHprime(x_i);

        D_i = D_values[i] * h_prime_i / (2.0 * rho[i] * rho[i]);

      } else if (cdeamVersion == 2) {

        // psip needs both fp[i] and fp[j] terms since r_ij appears in two

        //   terms of embed eng: Fi(sum rho_ij) and Fj(sum rho_ji)

        //   hence embed' = Fi(sum rho_ij) rhojp + Fj(sum rho_ji) rhoip

        rhoip = RhoPrimeOfR(index, itype, jtype);

        rhojp = RhoPrimeOfR(index, jtype, itype);

        fpair = fp[i]*rhojp + fp[j]*rhoip;

        recip = 1.0/r;

                                double x_j = -1;  // The value -1 indicates: no concentration dependence for this i-j pair

                                                  // because atom j is not of species A nor B.

        // The value -1 indicates: no concentration dependence for this

        // i-j pair because atom j is not of species A nor B.

        double x_j = -1;

        // This code line is required for ternary alloy.

        if (jtype == speciesA || jtype == speciesB) {

          ASSERT(rho[i] != 0.0);

          ASSERT(rho[j] != 0.0);

          // Compute local concentration at site j.

          x_j = rhoB[j]/rho[j];

          ASSERT(x_j >= 0 && x_j<=1.0);

          double D_j=0.0;

          if (cdeamVersion == 1) {

            // Calculate derivative of h(x_j) polynomial function.

            double h_prime_j = evalHprime(x_j);

            D_j = D_values[j] * h_prime_j / (2.0 * rho[j] * rho[j]);

          } else if (cdeamVersion == 2) {

        }

        // This if-clause is only required for a ternary alloy.

                                // Actually we don't need it at all because D_i should be zero anyway if

                                // atom i has no concentration dependent interactions (because it is not species A or B).

        // Actually we don't need it at all because D_i should be zero

        // anyway if atom i has no concentration dependent interactions

        // (because it is not species A or B).

        if (x_i != -1.0) {

          double t1 = -rhoB[i];

          if (jtype == speciesB) t1 += rho[i];

        double phip;

        double phi = PhiOfR(index, itype, jtype, recip, phip);

        if (itype == jtype || x_i == -1.0 || x_j == -1.0) {

          // Case of no concentration dependence.

          fpair += phip;

        } else {

          // We have a concentration dependence for the i-j interaction.

          double h=0.0;

          if (cdeamVersion == 1) {

            // Calculate h(x_i) polynomial function.

            double h_i = evalH(x_i);

            // Calculate h(x_j) polynomial function.

            double h_j = evalH(x_j);

            h = 0.5 * (h_i + h_j);

          } else if (cdeamVersion == 2) {

            // Average concentration.

            double x_ij = 0.5 * (x_i + x_j);

            // Calculate h(x_ij) polynomial function.

            h = evalH(x_ij);

          } else {

            ASSERT(false);

        }

        // Divide by r_ij and negate to get forces from gradient.

        fpair /= -r;

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

  PairEAMAlloy::coeff(narg, arg);

  // Make sure the EAM file is a CD-EAM binary alloy.

  if (setfl->nelements < 2)

    error->all(FLERR,"The EAM file must contain at least 2 elements to be used with the eam/cd pair style.");

  // Read in the coefficients of the h polynomial from the end of the EAM file.

  read_h_coeff(arg[2]);

        // Determine which atom type is the A species and which is the B species in the alloy.

        // By default take the first element (index 0) in the EAM file as the A species

        // and the second element (index 1) in the EAM file as the B species.

  // Determine which atom type is the A species and which is the B

  // species in the alloy.  By default take the first element (index 0)

  // in the EAM file as the A species and the second element (index 1)

  // in the EAM file as the B species.

  speciesA = -1;

  speciesB = -1;

  for (int i = 1; i <= atom->ntypes; i++) {

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

   Reads in the h(x) polynomial coefficients

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

void PairEAMCD::read_h_coeff(char *filename)

{

  if (comm->me == 0) {

    // Open potential file

    FILE *fptr;

    char line[MAXLINE];

    char nextline[MAXLINE];

    // h coefficients are stored at the end of the file.

    // Skip to last line of file.

    while(fgets(nextline, MAXLINE, fptr) != NULL) {

      strcpy(line, nextline);

    }

      error->one(FLERR,"Failed to read h(x) function coefficients from EAM file.");

    // Close the potential file.

    fclose(fptr);

  }

        buf[m++] = D_values[j];

      }

      return m;

                }

                else if(cdeamVersion == 2) {

    } else if (cdeamVersion == 2) {

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

        j = list[i];

        buf[m++] = fp[j];

        buf[m++] = rhoB[j];

      }

      return m;

                }

                else { ASSERT(false); return 0; }

        }

        else if(communicationStage == 4) {

    } else { ASSERT(false); return 0; }

  } else if (communicationStage == 4) {

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

      j = list[i];

      buf[m++] = D_values[j];

    }

    return m;

        }

        else return 0;

  } else return 0;

}

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

        rhoB[i] = buf[m++];

        D_values[i] = buf[m++];

      }

                }

                else if(cdeamVersion == 2) {

    } else if (cdeamVersion == 2) {

      for (i = first; i < last; i++) {

        fp[i] = buf[m++];

        rho[i] = buf[m++];

    } else {

      ASSERT(false);

    }

        }

        else if(communicationStage == 4) {

  } else if (communicationStage == 4) {

    for (i = first; i < last; i++) {

      D_values[i] = buf[m++];

    }

        buf[m++] = D_values[i];

      }

      return m;

                }

                else if(cdeamVersion == 2) {

    } else if (cdeamVersion == 2) {

      for (i = first; i < last; i++) {

        buf[m++] = rho[i];

        buf[m++] = rhoB[i];

      }

      return m;

                }

                else { ASSERT(false); return 0; }

        }

        else if(communicationStage == 3) {

    } else { ASSERT(false); return 0; }

  } else if (communicationStage == 3) {

    for (i = first; i < last; i++) {

      buf[m++] = D_values[i];

    }

    return m;

        }

        else return 0;

  } else return 0;

}

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

    } else {

      ASSERT(false);

    }

        }

        else if(communicationStage == 3) {

  } else if (communicationStage == 3) {

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

      j = list[i];

      D_values[j] += buf[m++];