diff --git a/src/Makefile b/src/Makefile index 0fad2a5712..0fde23a842 100755 --- a/src/Makefile +++ b/src/Makefile @@ -16,7 +16,7 @@ OBJ = $(SRC:.cpp=.o) PACKAGE = asphere class2 colloid dipole dpd gpu granular \ kspace manybody meam molecule opt peri poems prd reax xtc -PACKUSER = user-ackland user-atc user-cg-cmm user-ewaldn user-smd +PACKUSER = user-ackland user-atc user-cd-eam user-cg-cmm user-ewaldn user-smd PACKALL = $(PACKAGE) $(PACKUSER) diff --git a/src/USER-CD-EAM/Install.csh b/src/USER-CD-EAM/Install.csh new file mode 100755 index 0000000000..60f847809a --- /dev/null +++ b/src/USER-CD-EAM/Install.csh @@ -0,0 +1,18 @@ +# Install/Uninstall package classes in LAMMPS + +if ($1 == 1) then + + cp style_user_cd_eam.h .. + + cp pair_cdeam.cpp .. + cp pair_cdeam.h .. + +else if ($1 == 0) then + + rm ../style_user_cd_eam.h + touch ../style_user_cd_eam.h + + rm ../pair_cdeam.cpp + rm ../pair_cdeam.h + +endif diff --git a/src/USER-CD-EAM/pair_cdeam.cpp b/src/USER-CD-EAM/pair_cdeam.cpp new file mode 100644 index 0000000000..b5c30b3c63 --- /dev/null +++ b/src/USER-CD-EAM/pair_cdeam.cpp @@ -0,0 +1,663 @@ +/* ---------------------------------------------------------------------- + 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 (stukowski at mm.tu-darmstadt.de) +------------------------------------------------------------------------- */ + +/* + Concentration-dependent EAM (CD-EAM) potential for multi-component + systems. + + This potential class implements an improved version of the original + CD-EAM formalism. The new version (a.k.a. one-site model; + cdeamVersion==1) has been published in + + A. Stukowski, B. Sadigh, P. Erhart and A. Caro + Efficient implementation of the concentration-dependent embedded + atom method for molecular-dynamics and Monte-Carlo simulations + Model. Simul. Mater. Sci. Eng., 2009, 075005 + + This new formulation is more efficient for MD and Monte-Carlo + simulations and is the default. + + The original formulation (a.k.a. two-site model; cdeamVersion==2) is + also implemented and has been published in + + A. Caro, D. A. Crowson and M. Caro + Classical Many-Body Potential for Concentrated Alloys and the + Inversion of Order in Iron-Chromium Alloys + Phys. Rev. Lett., APS, 2005, 95, 075702 +*/ + +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#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; +//using namespace std; + +// 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(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; + 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("Invalid CD-EAM potential version."); + } +} + +PairCDEAM::~PairCDEAM() +{ + memory->sfree(rhoB); + memory->sfree(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 = 0; + + // Grow per-atom arrays if necessary + if(atom->nmax > nmax) { + memory->sfree(rho); + memory->sfree(fp); + memory->sfree(rhoB); + memory->sfree(D_values); + nmax = atom->nmax; + rho = (double *)memory->smalloc(nmax*sizeof(double),"pair:rho"); + rhoB = (double *)memory->smalloc(nmax*sizeof(double),"pair:rhoB"); + fp = (double *)memory->smalloc(nmax*sizeof(double),"pair:fp"); + D_values = (double *)memory->smalloc(nmax*sizeof(double),"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]; + + 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->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("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]; + 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->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]; + + 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; + 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; + 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_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("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("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("The second element from the EAM file may only be mapped to a single atom type."); + speciesB = i; + } + } + if(speciesA < 0) + error->all("The first element from the EAM file must be mapped to exactly one atom type."); + if(speciesB < 0) + error->all("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 *fp; + char line[MAXLINE]; + char nextline[MAXLINE]; + fp = fopen(filename,"r"); + if (fp == NULL) { + char str[128]; + sprintf(str,"Cannot open EAM potential file %s", filename); + error->one(str); + } + + // h coefficients are stored at the end of the file. + // Skip to last line of file. + while(fgets(nextline, MAXLINE, fp) != 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("Failed to read h(x) function coefficients from EAM file."); + + // Close the potential file. + fclose(fp); + } + + 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_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 4; + } + 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 3; + } + else { ASSERT(false); return 0; } + } + else if(communicationStage == 4) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = D_values[j]; + } + return 1; + } + else return 0; +} + +/* ---------------------------------------------------------------------- */ + +void PairCDEAM::unpack_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 3; + } + else if(cdeamVersion == 2) { + for(i = first; i < last; i++) { + buf[m++] = rho[i]; + buf[m++] = rhoB[i]; + } + return 2; + } + else { ASSERT(false); return 0; } + } + else if(communicationStage == 3) { + for(i = first; i < last; i++) { + buf[m++] = D_values[i]; + } + return 1; + } + 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; +} diff --git a/src/USER-CD-EAM/pair_cdeam.h b/src/USER-CD-EAM/pair_cdeam.h new file mode 100644 index 0000000000..3f121d461b --- /dev/null +++ b/src/USER-CD-EAM/pair_cdeam.h @@ -0,0 +1,222 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifndef PAIR_CDEAM_H +#define PAIR_CDEAM_H + +#include "pair_eam_alloy.h" + +namespace LAMMPS_NS { + +class PairCDEAM : public PairEAMAlloy +{ +public: + /// Constructor. + PairCDEAM(class LAMMPS*, int cdeamVersion); + + /// Destructor. + virtual ~PairCDEAM(); + + /// Calculates the energies and forces for all atoms in the system. + void compute(int, int); + + /// Parses the pair_coeff command parameters for this pair style. + void coeff(int, char **); + + /// This is for MPI communication with neighbor nodes. + int pack_comm(int, int *, double *, int, int *); + void unpack_comm(int, int, double *); + int pack_reverse_comm(int, int, double *); + void unpack_reverse_comm(int, int *, double *); + + /// Reports the memory usage of this pair style to LAMMPS. + double memory_usage(); + + /// Parses the coefficients of the h polynomial from the end of the EAM file. + void read_h_coeff(char* filename); + + public: + // The public interface exposed by this potential class. + + // Evaluates the h(x) polynomial for a given local concentration x. + inline double evalH(double x) const { + double v = 0.0; + for(int i = nhcoeff-1; i >= 1; i--) { + v = (v + hcoeff[i]) * x; + } + return v + hcoeff[0]; + } + + // Calculates the derivative of the h(x) polynomial. + inline double evalHprime(double x) const { + double v = 0.0; + for(int i = nhcoeff-1; i >= 2; i--) { + v = (v + (double)i * hcoeff[i]) * x; + } + return v + hcoeff[1]; + } + + // We have two versions of the CD-EAM potential. The user can choose which one he wants to use: + // + // Version 1 - One-site concentration: The h(x_i) function depends only on the concentration at the atomic site i. + // This is a new version with a slight modification to the formula. It happens to be computationally more efficient. + // It has been published in the MSMSE 2009 paper. + // + // Version 2 - Two-site concentration: The h(x_ij) function depends on the concentrations at two atomic sites i and j. + // This is the original version from the 2005 PRL paper. + int cdeamVersion; + + // Per-atom arrays + + // The partial density of B atoms at each atom site. + double *rhoB; + + // The intermediate value D_i for each atom. + // The meaning of these values depend on the version of the CD-EAM potential used: + // + // For the one-site concentration version these are the v_i values defined in equation (21) + // of the MSMSE paper. + // + // For the old two-site concentration version these are the M_i values defined in equation (12) + // of the MSMSE paper. + double *D_values; + + // The atom type index that is considered to be the A species in the alloy. + int speciesA; + + // The atom type index that is considered to be the B species in the alloy. + int speciesB; + + protected: + + // Evaluation functions: + + // This structure specifies an entry in one of the EAM spline tables + // and the corresponding floating point part. + typedef struct { + int m; + double p; + } EAMTableIndex; + + // Converts a radius value to an index value to be used in a spline table lookup. + inline EAMTableIndex radiusToTableIndex(double r) const { + EAMTableIndex index; + index.p = r*rdr + 1.0; + index.m = static_cast(index.p); + index.m = index.m <= (nr-1) ? index.m : (nr-1); + index.p -= index.m; + index.p = index.p <= 1.0 ? index.p : 1.0; + return index; + } + + // Converts a density value to an index value to be used in a spline table lookup. + inline EAMTableIndex rhoToTableIndex(double rho) const { + EAMTableIndex index; + index.p = rho*rdrho + 1.0; + index.m = static_cast(index.p); + index.m = index.m <= (nrho-1) ? index.m : (nrho-1); + index.p -= index.m; + index.p = index.p <= 1.0 ? index.p : 1.0; + return index; + } + + // Computes the derivative of rho(r) + inline double RhoPrimeOfR(const EAMTableIndex& index, int itype, int jtype) const { + const double* coeff = rhor_spline[type2rhor[itype][jtype]][index.m]; + return (coeff[0]*index.p + coeff[1])*index.p + coeff[2]; + } + + // Computes rho(r) + inline double RhoOfR(const EAMTableIndex& index, int itype, int jtype) const { + const double* coeff = rhor_spline[type2rhor[itype][jtype]][index.m]; + return ((coeff[3]*index.p + coeff[4])*index.p + coeff[5])*index.p + coeff[6]; + } + + // Computes the derivative of F(rho) + inline double FPrimeOfRho(const EAMTableIndex& index, int itype) const { + const double* coeff = frho_spline[type2frho[itype]][index.m]; + return (coeff[0]*index.p + coeff[1])*index.p + coeff[2]; + } + + // Computes F(rho) + inline double FofRho(const EAMTableIndex& index, int itype) const { + const double* coeff = frho_spline[type2frho[itype]][index.m]; + return ((coeff[3]*index.p + coeff[4])*index.p + coeff[5])*index.p + coeff[6]; + } + + // Computes the derivative of z2(r) + inline double Z2PrimeOfR(const EAMTableIndex& index, int itype, int jtype) const { + const double* coeff = z2r_spline[type2z2r[itype][jtype]][index.m]; + return (coeff[0]*index.p + coeff[1])*index.p + coeff[2]; + } + + // Computes z2(r) + inline double Z2OfR(const EAMTableIndex& index, int itype, int jtype) const { + const double* coeff = z2r_spline[type2z2r[itype][jtype]][index.m]; + return ((coeff[3]*index.p + coeff[4])*index.p + coeff[5])*index.p + coeff[6]; + } + + // Computes pair potential V_ij(r). + inline double PhiOfR(const EAMTableIndex& index, int itype, int jtype, const double oneOverR) const { + // phi = pair potential energy + // z2 = phi * r + const double* coeff = z2r_spline[type2z2r[itype][jtype]][index.m]; + const double z2 = ((coeff[3]*index.p + coeff[4])*index.p + coeff[5])*index.p + coeff[6]; + return z2 * oneOverR; + } + + // Computes pair potential V_ij(r) and its derivative. + inline double PhiOfR(const EAMTableIndex& index, int itype, int jtype, const double oneOverR, double& phid) const { + // phi = pair potential energy + // phip = phi' + // z2 = phi * r + // z2p = (phi * r)' = (phi' r) + phi + const double* coeff = z2r_spline[type2z2r[itype][jtype]][index.m]; + const double z2p = (coeff[0]*index.p + coeff[1])*index.p + coeff[2]; + const double z2 = ((coeff[3]*index.p + coeff[4])*index.p + coeff[5])*index.p + coeff[6]; + const double phi = z2 * oneOverR; + phid = z2p * oneOverR - phi * oneOverR; + return phi; + } + + // Parameters + + // h() polynomial function coefficients + double* hcoeff; + // The number of coefficients in the polynomial. + int nhcoeff; + + // This specifies the calculation stage to let the pack/unpack communication routines know + // which data to exchange. + int communicationStage; +}; + +/// The one-site concentration formulation of CD-EAM. + class PairCDEAM_OneSite : public PairCDEAM + { + public: + /// Constructor. + PairCDEAM_OneSite(class LAMMPS* lmp) : PairEAM(lmp), PairCDEAM(lmp, 1) {} + }; + + /// The two-site concentration formulation of CD-EAM. + class PairCDEAM_TwoSite : public PairCDEAM + { + public: + /// Constructor. + PairCDEAM_TwoSite(class LAMMPS* lmp) : PairEAM(lmp), PairCDEAM(lmp, 2) {} + }; + +}; + +#endif diff --git a/src/USER-CD-EAM/style_user_cd_eam.h b/src/USER-CD-EAM/style_user_cd_eam.h new file mode 100755 index 0000000000..1c3bebd873 --- /dev/null +++ b/src/USER-CD-EAM/style_user_cd_eam.h @@ -0,0 +1,21 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifdef PairInclude +#include "pair_cdeam.h" +#endif + +#ifdef PairClass +PairStyle(eam/cd,PairCDEAM_OneSite) +PairStyle(eam/cd/old,PairCDEAM_TwoSite) +#endif diff --git a/src/USER-CG-CMM/Install.csh b/src/USER-CG-CMM/Install.csh index d17aa83c41..f47836d48f 100644 --- a/src/USER-CG-CMM/Install.csh +++ b/src/USER-CG-CMM/Install.csh @@ -19,7 +19,6 @@ if ($1 == 1) then cp pair_cg_cmm_coul_long.cpp .. cp pair_cg_cmm_coul_long.h .. - else if ($1 == 0) then rm ../style_user_cg_cmm.h diff --git a/src/USER-CG-CMM/style_user_cg_cmm.h b/src/USER-CG-CMM/style_user_cg_cmm.h index d623bb9429..525eee3dd8 100644 --- a/src/USER-CG-CMM/style_user_cg_cmm.h +++ b/src/USER-CG-CMM/style_user_cg_cmm.h @@ -11,10 +11,6 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -// add new include files in appropriate Include ifdef -// add new style keywords and class names in appropriate Class ifdef -// see style.h for examples - #ifdef AngleInclude #include "angle_cg_cmm.h" #endif diff --git a/src/style_user_packages.h b/src/style_user_packages.h index e11a715c56..70c68719e0 100644 --- a/src/style_user_packages.h +++ b/src/style_user_packages.h @@ -16,6 +16,7 @@ #include "style_user_ackland.h" #include "style_user_atc.h" +#include "style_user_cd_eam.h" #include "style_user_cg_cmm.h" #include "style_user_ewaldn.h" #include "style_user_smd.h"