diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
index 671b9341d524834e6e9ef8c0ee30a674818eaacb..d51108ba7aed1aae0a6040e7bda6f9d8dab20555 100644
--- a/cmake/CMakeLists.txt
+++ b/cmake/CMakeLists.txt
@@ -181,7 +181,6 @@ endmacro()
pkg_depends(MPIIO MPI)
pkg_depends(USER-ATC MANYBODY)
pkg_depends(USER-LB MPI)
-pkg_depends(USER-MISC MANYBODY)
pkg_depends(USER-PHONON KSPACE)
######################################################
diff --git a/doc/src/Section_commands.txt b/doc/src/Section_commands.txt
index 7948e0de32857547e1b9d01a6fa4bbf44194e342..daedb00a0a6f6b431356d04fece2dfaeddbb9f23 100644
--- a/doc/src/Section_commands.txt
+++ b/doc/src/Section_commands.txt
@@ -969,6 +969,7 @@ OPT.
"dsmc"_pair_dsmc.html,
"eam (gikot)"_pair_eam.html,
"eam/alloy (gikot)"_pair_eam.html,
+"eam/cd (o)"_pair_eam.html,
"eam/fs (gikot)"_pair_eam.html,
"eim (o)"_pair_eim.html,
"gauss (go)"_pair_gauss.html,
@@ -1069,7 +1070,6 @@ package"_Section_start.html#start_3.
"coul/shield"_pair_coul_shield.html,
"dpd/fdt"_pair_dpd_fdt.html,
"dpd/fdt/energy (k)"_pair_dpd_fdt.html,
-"eam/cd (o)"_pair_eam.html,
"edip (o)"_pair_edip.html,
"edip/multi"_pair_edip.html,
"edpd"_pair_meso.html,
diff --git a/doc/src/pair_eam.txt b/doc/src/pair_eam.txt
index 8d4d11341cd45dae75a2e48c7ec9d9e0e207055a..b3c770179fb948329d4b746b863041a101c83436 100644
--- a/doc/src/pair_eam.txt
+++ b/doc/src/pair_eam.txt
@@ -413,15 +413,10 @@ The eam pair styles can only be used via the {pair} keyword of the
[Restrictions:]
-All of these styles except the {eam/cd} style are part of the MANYBODY
-package. They are only enabled if LAMMPS was built with that package.
+All of these styles are part of the MANYBODY package. They are only
+enabled if LAMMPS was built with that package.
See the "Making LAMMPS"_Section_start.html#start_3 section for more info.
-The {eam/cd} style is part of the USER-MISC package and also requires
-the MANYBODY package. It is only enabled if LAMMPS was built with
-those packages. See the "Making LAMMPS"_Section_start.html#start_3
-section for more info.
-
[Related commands:]
"pair_coeff"_pair_coeff.html
diff --git a/src/MANYBODY/pair_eam_cd.cpp b/src/MANYBODY/pair_eam_cd.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..66ebad6244a9d16b939bfb61c63862742e709cc8
--- /dev/null
+++ b/src/MANYBODY/pair_eam_cd.cpp
@@ -0,0 +1,677 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ Contributing author: Alexander Stukowski
+ Technical University of Darmstadt,
+ Germany Department of Materials Science
+------------------------------------------------------------------------- */
+
+#include <cmath>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include "pair_eam_cd.h"
+#include "atom.h"
+#include "force.h"
+#include "comm.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+#define ASSERT(cond)
+#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)
+{
+ single_enable = 0;
+ restartinfo = 0;
+
+ rhoB = NULL;
+ D_values = NULL;
+ hcoeff = NULL;
+
+ // Set communication buffer sizes needed by this pair style.
+
+ if (cdeamVersion == 1) {
+ comm_forward = 4;
+ comm_reverse = 3;
+ } else if (cdeamVersion == 2) {
+ comm_forward = 3;
+ comm_reverse = 2;
+ } else {
+ error->all(FLERR,"Invalid eam/cd potential version.");
+ }
+}
+
+PairEAMCD::~PairEAMCD()
+{
+ memory->destroy(rhoB);
+ memory->destroy(D_values);
+ if (hcoeff) delete[] hcoeff;
+}
+
+void PairEAMCD::compute(int eflag, int vflag)
+{
+ int i,j,ii,jj,inum,jnum,itype,jtype;
+ double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
+ double rsq,rhoip,rhojp,recip,phi;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+
+ evdwl = 0.0;
+ if (eflag || vflag) ev_setup(eflag,vflag);
+ 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);
+ memory->destroy(rhoB);
+ memory->destroy(D_values);
+ nmax = atom->nmax;
+ memory->create(rho,nmax,"pair:rho");
+ memory->create(rhoB,nmax,"pair:rhoB");
+ memory->create(fp,nmax,"pair:fp");
+ memory->create(D_values,nmax,"pair:D_values");
+ }
+
+ double **x = atom->x;
+ double **f = atom->f;
+ int *type = atom->type;
+ int nlocal = atom->nlocal;
+ int newton_pair = force->newton_pair;
+
+ inum = list->inum;
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+
+ // Zero out per-atom arrays.
+
+ int m = nlocal + atom->nghost;
+ for (i = 0; i < m; i++) {
+ rho[i] = 0.0;
+ rhoB[i] = 0.0;
+ D_values[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).
+
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+ xtmp = x[i][0];
+ ytmp = x[i][1];
+ ztmp = x[i][2];
+ itype = type[i];
+ jlist = firstneigh[i];
+ jnum = numneigh[i];
+
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ j &= NEIGHMASK;
+
+ delx = xtmp - x[j][0];
+ dely = ytmp - x[j][1];
+ delz = ztmp - x[j][2];
+ rsq = delx*delx + dely*dely + delz*delz;
+
+ if (rsq < cutforcesq) {
+ jtype = type[j];
+ double r = sqrt(rsq);
+ const EAMTableIndex index = radiusToTableIndex(r);
+ double localrho = RhoOfR(index, jtype, itype);
+ rho[i] += localrho;
+ if (jtype == speciesB) rhoB[i] += localrho;
+ if (newton_pair || j < nlocal) {
+ localrho = RhoOfR(index, itype, jtype);
+ rho[j] += localrho;
+ if (itype == speciesB) rhoB[j] += localrho;
+ }
+
+ 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)) {
+ double Phi_AB = PhiOfR(index, itype, jtype, 1.0 / r);
+ D_values[i] += Phi_AB;
+ if (newton_pair || j < nlocal)
+ D_values[j] += Phi_AB;
+ }
+ }
+ }
+ }
+ }
+
+ // 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]);
+ fp[i] = FPrimeOfRho(index, type[i]);
+ if (eflag) {
+ phi = FofRho(index, type[i]);
+ if (eflag_global) eng_vdwl += phi;
+ if (eflag_atom) eatom[i] += phi;
+ }
+ }
+
+ // 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.
+ //
+ //for (i = 0; i < nlocal + atom->nghost; i++) {
+ // if (rho[i] == 0 && (type[i] == speciesA || type[i] == speciesB))
+ // error->one(FLERR,"CD-EAM potential routine: Detected atom with zero electron density.");
+ //}
+
+ // Stage II
+ // 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];
+ ytmp = x[i][1];
+ ztmp = x[i][2];
+ itype = type[i];
+ jlist = firstneigh[i];
+ 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.
+
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ j &= NEIGHMASK;
+ jtype = type[j];
+ if (itype == jtype) continue;
+
+ // This code line is required for ternary alloys.
+
+ if (jtype != speciesA && jtype != speciesB) continue;
+
+ delx = xtmp - x[j][0];
+ dely = ytmp - x[j][1];
+ delz = ztmp - x[j][2];
+ rsq = delx*delx + dely*dely + delz*delz;
+
+ if (rsq < cutforcesq) {
+ double r = sqrt(rsq);
+ 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]);
+ if (newton_pair || j < nlocal)
+ D_values[j] += h_prime * Phi_AB / (2.0 * rho[j] * rho[j]);
+ }
+ }
+ }
+
+ // Communicate and sum D values.
+
+ if (newton_pair) {
+ communicationStage = 3;
+ comm->reverse_comm_pair(this);
+ }
+ communicationStage = 4;
+ comm->forward_comm_pair(this);
+ }
+
+ // Stage III
+
+ // Compute force acting on each atom.
+
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+ xtmp = x[i][0];
+ ytmp = x[i][1];
+ ztmp = x[i][2];
+ itype = type[i];
+
+ jlist = firstneigh[i];
+ jnum = numneigh[i];
+
+ // Concentration at site 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) {
+ D_i = D_values[i];
+ } else {
+ ASSERT(false);
+ }
+ }
+
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ j &= NEIGHMASK;
+
+ delx = xtmp - x[j][0];
+ dely = ytmp - x[j][1];
+ delz = ztmp - x[j][2];
+ rsq = delx*delx + dely*dely + delz*delz;
+
+ if (rsq < cutforcesq) {
+ jtype = type[j];
+ double r = sqrt(rsq);
+ const EAMTableIndex index = radiusToTableIndex(r);
+
+ // rhoip = derivative of (density at atom j due to atom i)
+ // rhojp = derivative of (density at atom i due to atom j)
+ // 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;
+
+ // 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) {
+ D_j = D_values[j];
+ } else {
+ ASSERT(false);
+ }
+ double t2 = -rhoB[j];
+ if (itype == speciesB) t2 += rho[j];
+ fpair += D_j * rhoip * t2;
+ }
+
+ // 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).
+
+ if (x_i != -1.0) {
+ double t1 = -rhoB[i];
+ if (jtype == speciesB) t1 += rho[i];
+ fpair += D_i * rhojp * t1;
+ }
+
+ 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);
+ }
+ fpair += h * phip;
+ phi *= h;
+ }
+
+ // Divide by r_ij and negate to get forces from gradient.
+
+ fpair /= -r;
+
+ f[i][0] += delx*fpair;
+ f[i][1] += dely*fpair;
+ f[i][2] += delz*fpair;
+ if (newton_pair || j < nlocal) {
+ f[j][0] -= delx*fpair;
+ f[j][1] -= dely*fpair;
+ f[j][2] -= delz*fpair;
+ }
+
+ if (eflag) evdwl = phi;
+ if (evflag) ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,fpair,delx,dely,delz);
+ }
+ }
+ }
+
+ if (vflag_fdotr) virial_fdotr_compute();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairEAMCD::coeff(int narg, char **arg)
+{
+ 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.
+
+ speciesA = -1;
+ speciesB = -1;
+ for (int i = 1; i <= atom->ntypes; i++) {
+ if (map[i] == 0) {
+ if (speciesA >= 0)
+ error->all(FLERR,"The first element from the EAM file may only be mapped to a single atom type.");
+ speciesA = i;
+ }
+ if (map[i] == 1) {
+ if (speciesB >= 0)
+ error->all(FLERR,"The second element from the EAM file may only be mapped to a single atom type.");
+ speciesB = i;
+ }
+ }
+ if (speciesA < 0)
+ error->all(FLERR,"The first element from the EAM file must be mapped to exactly one atom type.");
+ if (speciesB < 0)
+ error->all(FLERR,"The second element from the EAM file must be mapped to exactly one atom type.");
+}
+
+/* ----------------------------------------------------------------------
+ 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];
+ fptr = force->open_potential(filename);
+ if (fptr == NULL) {
+ char str[128];
+ sprintf(str,"Cannot open EAM potential file %s", filename);
+ error->one(FLERR,str);
+ }
+
+ // 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);
+ }
+ char* ptr = strtok(line, " \t\n\r\f");
+ int degree = atoi(ptr);
+ nhcoeff = degree+1;
+ hcoeff = new double[nhcoeff];
+ int i = 0;
+ while((ptr = strtok(NULL," \t\n\r\f")) != NULL && i < nhcoeff) {
+ hcoeff[i++] = atof(ptr);
+ }
+ if (i != nhcoeff || nhcoeff < 1)
+ error->one(FLERR,"Failed to read h(x) function coefficients from EAM file.");
+
+ // Close the potential file.
+
+ fclose(fptr);
+ }
+
+ MPI_Bcast(&nhcoeff, 1, MPI_INT, 0, world);
+ if (comm->me != 0) hcoeff = new double[nhcoeff];
+ MPI_Bcast(hcoeff, nhcoeff, MPI_DOUBLE, 0, world);
+}
+
+
+/* ---------------------------------------------------------------------- */
+
+int PairEAMCD::pack_forward_comm(int n, int *list, double *buf,
+ int pbc_flag, int *pbc)
+{
+ int i,j,m;
+
+ m = 0;
+ if (communicationStage == 2) {
+ if (cdeamVersion == 1) {
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ buf[m++] = fp[j];
+ buf[m++] = rho[j];
+ buf[m++] = rhoB[j];
+ buf[m++] = D_values[j];
+ }
+ return m;
+ } else if (cdeamVersion == 2) {
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ buf[m++] = fp[j];
+ buf[m++] = rho[j];
+ buf[m++] = rhoB[j];
+ }
+ return m;
+ } 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;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairEAMCD::unpack_forward_comm(int n, int first, double *buf)
+{
+ int i,m,last;
+
+ m = 0;
+ last = first + n;
+ if (communicationStage == 2) {
+ if (cdeamVersion == 1) {
+ for (i = first; i < last; i++) {
+ fp[i] = buf[m++];
+ rho[i] = buf[m++];
+ rhoB[i] = buf[m++];
+ D_values[i] = buf[m++];
+ }
+ } else if (cdeamVersion == 2) {
+ for (i = first; i < last; i++) {
+ fp[i] = buf[m++];
+ rho[i] = buf[m++];
+ rhoB[i] = buf[m++];
+ }
+ } else {
+ ASSERT(false);
+ }
+ } else if (communicationStage == 4) {
+ for (i = first; i < last; i++) {
+ D_values[i] = buf[m++];
+ }
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+int PairEAMCD::pack_reverse_comm(int n, int first, double *buf)
+{
+ int i,m,last;
+
+ m = 0;
+ last = first + n;
+
+ if (communicationStage == 1) {
+ if (cdeamVersion == 1) {
+ for (i = first; i < last; i++) {
+ buf[m++] = rho[i];
+ buf[m++] = rhoB[i];
+ buf[m++] = D_values[i];
+ }
+ return m;
+ } 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) {
+ for (i = first; i < last; i++) {
+ buf[m++] = D_values[i];
+ }
+ return m;
+ } else return 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairEAMCD::unpack_reverse_comm(int n, int *list, double *buf)
+{
+ int i,j,m;
+
+ m = 0;
+ if (communicationStage == 1) {
+ if (cdeamVersion == 1) {
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ rho[j] += buf[m++];
+ rhoB[j] += buf[m++];
+ D_values[j] += buf[m++];
+ }
+ } else if (cdeamVersion == 2) {
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ rho[j] += buf[m++];
+ rhoB[j] += buf[m++];
+ }
+ } else {
+ ASSERT(false);
+ }
+ } else if (communicationStage == 3) {
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ D_values[j] += buf[m++];
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ memory usage of local atom-based arrays
+------------------------------------------------------------------------- */
+double PairEAMCD::memory_usage()
+{
+ double bytes = 2 * nmax * sizeof(double);
+ return PairEAMAlloy::memory_usage() + bytes;
+}
diff --git a/src/USER-MISC/pair_cdeam.h b/src/MANYBODY/pair_eam_cd.h
similarity index 94%
rename from src/USER-MISC/pair_cdeam.h
rename to src/MANYBODY/pair_eam_cd.h
index 934b7601a4bfac37bbad057f94ed8a373a6848fe..15486aed6d96fafd13014296638bcdb82af30962 100644
--- a/src/USER-MISC/pair_cdeam.h
+++ b/src/MANYBODY/pair_eam_cd.h
@@ -13,8 +13,8 @@
#ifdef PAIR_CLASS
-PairStyle(eam/cd,PairCDEAM_OneSite)
-PairStyle(eam/cd/old,PairCDEAM_TwoSite)
+PairStyle(eam/cd,PairEAMCD_OneSite)
+PairStyle(eam/cd/old,PairEAMCD_TwoSite)
#else
@@ -25,14 +25,14 @@ PairStyle(eam/cd/old,PairCDEAM_TwoSite)
namespace LAMMPS_NS {
-class PairCDEAM : public PairEAMAlloy
+class PairEAMCD : public PairEAMAlloy
{
public:
/// Constructor.
- PairCDEAM(class LAMMPS*, int cdeamVersion);
+ PairEAMCD(class LAMMPS*, int cdeamVersion);
/// Destructor.
- virtual ~PairCDEAM();
+ virtual ~PairEAMCD();
/// Calculates the energies and forces for all atoms in the system.
virtual void compute(int, int);
@@ -211,19 +211,19 @@ public:
};
/// The one-site concentration formulation of CD-EAM.
- class PairCDEAM_OneSite : public PairCDEAM
+ class PairEAMCD_OneSite : public PairEAMCD
{
public:
/// Constructor.
- PairCDEAM_OneSite(class LAMMPS* lmp) : PairEAM(lmp), PairCDEAM(lmp, 1) {}
+ PairEAMCD_OneSite(class LAMMPS* lmp) : PairEAM(lmp), PairEAMCD(lmp, 1) {}
};
/// The two-site concentration formulation of CD-EAM.
- class PairCDEAM_TwoSite : public PairCDEAM
+ class PairEAMCD_TwoSite : public PairEAMCD
{
public:
/// Constructor.
- PairCDEAM_TwoSite(class LAMMPS* lmp) : PairEAM(lmp), PairCDEAM(lmp, 2) {}
+ PairEAMCD_TwoSite(class LAMMPS* lmp) : PairEAM(lmp), PairEAMCD(lmp, 2) {}
};
}
diff --git a/src/USER-MISC/Install.sh b/src/USER-MISC/Install.sh
deleted file mode 100755
index 2d42125ec3a79cdb00cc88cb31d12b8be3959f4d..0000000000000000000000000000000000000000
--- a/src/USER-MISC/Install.sh
+++ /dev/null
@@ -1,40 +0,0 @@
-# Install/unInstall package files in LAMMPS
-# mode = 0/1/2 for uninstall/install/update
-
-mode=$1
-
-# enforce using portable C locale
-LC_ALL=C
-export LC_ALL
-
-# arg1 = file, arg2 = file it depends on
-
-action () {
- if (test $mode = 0) then
- rm -f ../$1
- elif (! cmp -s $1 ../$1) then
- if (test -z "$2" || test -e ../$2) then
- cp $1 ..
- if (test $mode = 2) then
- echo " updating src/$1"
- fi
- fi
- elif (test ! -n "$2") then
- if (test ! -e ../$2) then
- rm -f ../$1
- fi
- fi
-}
-
-# all package files
-# only a few files have dependencies
-
-for file in *.cpp *.h; do
- if (test $file = "pair_cdeam.cpp") then
- action pair_cdeam.cpp pair_eam_alloy.cpp
- elif (test $file = "pair_cdeam.h") then
- action pair_cdeam.h pair_eam_alloy.cpp
- else
- test -f ${file} && action $file
- fi
-done
diff --git a/src/USER-MISC/README b/src/USER-MISC/README
index 68a6252d8dc692285b5577c6281723c4515f8d54..0f9e7bf383bf3e005b828672ad24b360d6f9b440 100644
--- a/src/USER-MISC/README
+++ b/src/USER-MISC/README
@@ -65,7 +65,6 @@ pair_style buck/mdf, Paolo Raiteri, p.raiteri at curtin.edu.au, 2 Dec 15
pair_style coul/diel, Axel Kohlmeyer, akohlmey at gmail.com, 1 Dec 11
pair_style dipole/sf, Mario Orsi, orsimario at gmail.com, 8 Aug 11
pair_style edip, Luca Ferraro, luca.ferraro at caspur.it, 15 Sep 11
-pair_style eam/cd, Alexander Stukowski, stukowski at mm.tu-darmstadt.de, 7 Nov 09
pair_style extep, Jaap Kroes (Radboud U), jaapkroes at gmail dot com, 28 Nov 17
pair_style gauss/cut, Axel Kohlmeyer, akohlmey at gmail.com, 1 Dec 11
pair_style lennard/mdf, Paolo Raiteri, p.raiteri at curtin.edu.au, 2 Dec 15
diff --git a/src/USER-MISC/pair_cdeam.cpp b/src/USER-MISC/pair_cdeam.cpp
deleted file mode 100644
index 53d9036a61d44da205f921c3272ccf41b9075cde..0000000000000000000000000000000000000000
--- a/src/USER-MISC/pair_cdeam.cpp
+++ /dev/null
@@ -1,644 +0,0 @@
-/* ----------------------------------------------------------------------
- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
- http://lammps.sandia.gov, Sandia National Laboratories
- Steve Plimpton, sjplimp@sandia.gov
-
- Copyright (2003) Sandia Corporation. Under the terms of Contract
- DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
- certain rights in this software. This software is distributed under
- the GNU General Public License.
-
- See the README file in the top-level LAMMPS directory.
-------------------------------------------------------------------------- */
-
-/* ----------------------------------------------------------------------
- Contributing author: Alexander Stukowski
- Technical University of Darmstadt,
- Germany Department of Materials Science
-------------------------------------------------------------------------- */
-
-#include <cmath>
-#include <cstdio>
-#include <cstdlib>
-#include <cstring>
-#include "pair_cdeam.h"
-#include "atom.h"
-#include "force.h"
-#include "comm.h"
-#include "neighbor.h"
-#include "neigh_list.h"
-#include "memory.h"
-#include "error.h"
-
-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.
-
-PairCDEAM::PairCDEAM(LAMMPS *lmp, int _cdeamVersion) : PairEAM(lmp), PairEAMAlloy(lmp), cdeamVersion(_cdeamVersion)
-{
- single_enable = 0;
- restartinfo = 0;
-
- rhoB = NULL;
- D_values = NULL;
- hcoeff = NULL;
-
- // Set communication buffer sizes needed by this pair style.
- if(cdeamVersion == 1) {
- comm_forward = 4;
- comm_reverse = 3;
- }
- else if(cdeamVersion == 2) {
- comm_forward = 3;
- comm_reverse = 2;
- }
- else {
- error->all(FLERR,"Invalid CD-EAM potential version.");
- }
-}
-
-PairCDEAM::~PairCDEAM()
-{
- memory->destroy(rhoB);
- memory->destroy(D_values);
- if(hcoeff) delete[] hcoeff;
-}
-
-void PairCDEAM::compute(int eflag, int vflag)
-{
- int i,j,ii,jj,inum,jnum,itype,jtype;
- double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
- double rsq,rhoip,rhojp,recip,phi;
- int *ilist,*jlist,*numneigh,**firstneigh;
-
- evdwl = 0.0;
- if (eflag || vflag) ev_setup(eflag,vflag);
- 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);
- memory->destroy(rhoB);
- memory->destroy(D_values);
- nmax = atom->nmax;
- memory->create(rho,nmax,"pair:rho");
- memory->create(rhoB,nmax,"pair:rhoB");
- memory->create(fp,nmax,"pair:fp");
- memory->create(D_values,nmax,"pair:D_values");
- }
-
- double **x = atom->x;
- double **f = atom->f;
- int *type = atom->type;
- int nlocal = atom->nlocal;
- int newton_pair = force->newton_pair;
-
- inum = list->inum;
- ilist = list->ilist;
- numneigh = list->numneigh;
- firstneigh = list->firstneigh;
-
- // Zero out per-atom arrays.
- int m = nlocal + atom->nghost;
- for(i = 0; i < m; i++) {
- rho[i] = 0.0;
- rhoB[i] = 0.0;
- D_values[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).
- for(ii = 0; ii < inum; ii++) {
- i = ilist[ii];
- xtmp = x[i][0];
- ytmp = x[i][1];
- ztmp = x[i][2];
- itype = type[i];
- jlist = firstneigh[i];
- jnum = numneigh[i];
-
- for(jj = 0; jj < jnum; jj++) {
- j = jlist[jj];
- j &= NEIGHMASK;
-
- delx = xtmp - x[j][0];
- dely = ytmp - x[j][1];
- delz = ztmp - x[j][2];
- rsq = delx*delx + dely*dely + delz*delz;
-
- if(rsq < cutforcesq) {
- jtype = type[j];
- double r = sqrt(rsq);
- const EAMTableIndex index = radiusToTableIndex(r);
- double localrho = RhoOfR(index, jtype, itype);
- rho[i] += localrho;
- if(jtype == speciesB) rhoB[i] += localrho;
- if(newton_pair || j < nlocal) {
- localrho = RhoOfR(index, itype, jtype);
- rho[j] += localrho;
- if(itype == speciesB) rhoB[j] += localrho;
- }
-
- 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)) {
- double Phi_AB = PhiOfR(index, itype, jtype, 1.0 / r);
- D_values[i] += Phi_AB;
- if(newton_pair || j < nlocal)
- D_values[j] += Phi_AB;
- }
- }
- }
- }
- }
-
- // 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]);
- fp[i] = FPrimeOfRho(index, type[i]);
- if(eflag) {
- phi = FofRho(index, type[i]);
- if (eflag_global) eng_vdwl += phi;
- if (eflag_atom) eatom[i] += phi;
- }
- }
-
- // 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.
- //
- //for(i = 0; i < nlocal + atom->nghost; i++) {
- // if(rho[i] == 0 && (type[i] == speciesA || type[i] == speciesB))
- // error->one(FLERR,"CD-EAM potential routine: Detected atom with zero electron density.");
- //}
-
- // Stage II
- // 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];
- ytmp = x[i][1];
- ztmp = x[i][2];
- itype = type[i];
- jlist = firstneigh[i];
- 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.
-
- for(jj = 0; jj < jnum; jj++) {
- j = jlist[jj];
- j &= NEIGHMASK;
- jtype = type[j];
- if(itype == jtype) continue;
-
- // This code line is required for ternary alloys.
- if(jtype != speciesA && jtype != speciesB) continue;
-
- delx = xtmp - x[j][0];
- dely = ytmp - x[j][1];
- delz = ztmp - x[j][2];
- rsq = delx*delx + dely*dely + delz*delz;
-
- if(rsq < cutforcesq) {
- double r = sqrt(rsq);
- 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]);
- if(newton_pair || j < nlocal)
- D_values[j] += h_prime * Phi_AB / (2.0 * rho[j] * rho[j]);
- }
- }
- }
-
- // Communicate and sum D values.
- if(newton_pair) {
- communicationStage = 3;
- comm->reverse_comm_pair(this);
- }
- communicationStage = 4;
- comm->forward_comm_pair(this);
- }
-
- // Stage III
-
- // Compute force acting on each atom.
- for(ii = 0; ii < inum; ii++) {
- i = ilist[ii];
- xtmp = x[i][0];
- ytmp = x[i][1];
- ztmp = x[i][2];
- itype = type[i];
-
- jlist = firstneigh[i];
- 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.
- // It will be replaced by the concentration at site i if atom i is either A or B.
-
- 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) {
- D_i = D_values[i];
- } else {
- ASSERT(false);
- }
- }
-
- for(jj = 0; jj < jnum; jj++) {
- j = jlist[jj];
- j &= NEIGHMASK;
-
- delx = xtmp - x[j][0];
- dely = ytmp - x[j][1];
- delz = ztmp - x[j][2];
- rsq = delx*delx + dely*dely + delz*delz;
-
- if(rsq < cutforcesq) {
- jtype = type[j];
- double r = sqrt(rsq);
- const EAMTableIndex index = radiusToTableIndex(r);
-
- // rhoip = derivative of (density at atom j due to atom i)
- // rhojp = derivative of (density at atom i due to atom j)
- // 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.
-
- // 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) {
- D_j = D_values[j];
- } else {
- ASSERT(false);
- }
- double t2 = -rhoB[j];
- if(itype == speciesB) t2 += rho[j];
- fpair += D_j * rhoip * t2;
- }
-
- // 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).
- if(x_i != -1.0) {
- double t1 = -rhoB[i];
- if(jtype == speciesB) t1 += rho[i];
- fpair += D_i * rhojp * t1;
- }
-
- 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);
- }
- fpair += h * phip;
- phi *= h;
- }
-
- // Divide by r_ij and negate to get forces from gradient.
- fpair /= -r;
-
- f[i][0] += delx*fpair;
- f[i][1] += dely*fpair;
- f[i][2] += delz*fpair;
- if(newton_pair || j < nlocal) {
- f[j][0] -= delx*fpair;
- f[j][1] -= dely*fpair;
- f[j][2] -= delz*fpair;
- }
-
- if(eflag) evdwl = phi;
- if(evflag) ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,fpair,delx,dely,delz);
- }
- }
- }
-
- if(vflag_fdotr) virial_fdotr_compute();
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairCDEAM::coeff(int narg, char **arg)
-{
- 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.
- speciesA = -1;
- speciesB = -1;
- for(int i = 1; i <= atom->ntypes; i++) {
- if(map[i] == 0) {
- if(speciesA >= 0)
- error->all(FLERR,"The first element from the EAM file may only be mapped to a single atom type.");
- speciesA = i;
- }
- if(map[i] == 1) {
- if(speciesB >= 0)
- error->all(FLERR,"The second element from the EAM file may only be mapped to a single atom type.");
- speciesB = i;
- }
- }
- if(speciesA < 0)
- error->all(FLERR,"The first element from the EAM file must be mapped to exactly one atom type.");
- if(speciesB < 0)
- error->all(FLERR,"The second element from the EAM file must be mapped to exactly one atom type.");
-}
-
-/* ----------------------------------------------------------------------
- Reads in the h(x) polynomial coefficients
-------------------------------------------------------------------------- */
-void PairCDEAM::read_h_coeff(char *filename)
-{
- if(comm->me == 0) {
- // Open potential file
- FILE *fptr;
- char line[MAXLINE];
- char nextline[MAXLINE];
- fptr = force->open_potential(filename);
- if (fptr == NULL) {
- char str[128];
- sprintf(str,"Cannot open EAM potential file %s", filename);
- error->one(FLERR,str);
- }
-
- // 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);
- }
- char* ptr = strtok(line, " \t\n\r\f");
- int degree = atoi(ptr);
- nhcoeff = degree+1;
- hcoeff = new double[nhcoeff];
- int i = 0;
- while((ptr = strtok(NULL," \t\n\r\f")) != NULL && i < nhcoeff) {
- hcoeff[i++] = atof(ptr);
- }
- if(i != nhcoeff || nhcoeff < 1)
- error->one(FLERR,"Failed to read h(x) function coefficients from EAM file.");
-
- // Close the potential file.
- fclose(fptr);
- }
-
- MPI_Bcast(&nhcoeff, 1, MPI_INT, 0, world);
- if(comm->me != 0) hcoeff = new double[nhcoeff];
- MPI_Bcast(hcoeff, nhcoeff, MPI_DOUBLE, 0, world);
-}
-
-
-/* ---------------------------------------------------------------------- */
-
-int PairCDEAM::pack_forward_comm(int n, int *list, double *buf,
- int pbc_flag, int *pbc)
-{
- int i,j,m;
-
- m = 0;
- if(communicationStage == 2) {
- if(cdeamVersion == 1) {
- for (i = 0; i < n; i++) {
- j = list[i];
- buf[m++] = fp[j];
- buf[m++] = rho[j];
- buf[m++] = rhoB[j];
- buf[m++] = D_values[j];
- }
- return m;
- }
- else if(cdeamVersion == 2) {
- for (i = 0; i < n; i++) {
- j = list[i];
- buf[m++] = fp[j];
- buf[m++] = rho[j];
- buf[m++] = rhoB[j];
- }
- return m;
- }
- 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;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairCDEAM::unpack_forward_comm(int n, int first, double *buf)
-{
- int i,m,last;
-
- m = 0;
- last = first + n;
- if(communicationStage == 2) {
- if(cdeamVersion == 1) {
- for(i = first; i < last; i++) {
- fp[i] = buf[m++];
- rho[i] = buf[m++];
- rhoB[i] = buf[m++];
- D_values[i] = buf[m++];
- }
- }
- else if(cdeamVersion == 2) {
- for(i = first; i < last; i++) {
- fp[i] = buf[m++];
- rho[i] = buf[m++];
- rhoB[i] = buf[m++];
- }
- } else {
- ASSERT(false);
- }
- }
- else if(communicationStage == 4) {
- for(i = first; i < last; i++) {
- D_values[i] = buf[m++];
- }
- }
-}
-
-/* ---------------------------------------------------------------------- */
-int PairCDEAM::pack_reverse_comm(int n, int first, double *buf)
-{
- int i,m,last;
-
- m = 0;
- last = first + n;
-
- if(communicationStage == 1) {
- if(cdeamVersion == 1) {
- for(i = first; i < last; i++) {
- buf[m++] = rho[i];
- buf[m++] = rhoB[i];
- buf[m++] = D_values[i];
- }
- return m;
- }
- 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) {
- for(i = first; i < last; i++) {
- buf[m++] = D_values[i];
- }
- return m;
- }
- else return 0;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairCDEAM::unpack_reverse_comm(int n, int *list, double *buf)
-{
- int i,j,m;
-
- m = 0;
- if(communicationStage == 1) {
- if(cdeamVersion == 1) {
- for(i = 0; i < n; i++) {
- j = list[i];
- rho[j] += buf[m++];
- rhoB[j] += buf[m++];
- D_values[j] += buf[m++];
- }
- } else if(cdeamVersion == 2) {
- for(i = 0; i < n; i++) {
- j = list[i];
- rho[j] += buf[m++];
- rhoB[j] += buf[m++];
- }
- } else {
- ASSERT(false);
- }
- }
- else if(communicationStage == 3) {
- for(i = 0; i < n; i++) {
- j = list[i];
- D_values[j] += buf[m++];
- }
- }
-}
-
-/* ----------------------------------------------------------------------
- memory usage of local atom-based arrays
-------------------------------------------------------------------------- */
-double PairCDEAM::memory_usage()
-{
- double bytes = 2 * nmax * sizeof(double);
- return PairEAMAlloy::memory_usage() + bytes;
-}