diff --git a/LAMMPS_plugin/install.sh b/LAMMPS_plugin/install.sh index 7431a30..f070724 100755 --- a/LAMMPS_plugin/install.sh +++ b/LAMMPS_plugin/install.sh @@ -22,3 +22,8 @@ ln -sf $(pwd)/src/pair_hybrid_overlay_mlml.h ${lammps}/src/pair_hybrid_overlay_m ln -sf $(pwd)/src/pair_hybrid_overlay_mlml.cpp ${lammps}/src/pair_hybrid_overlay_mlml.cpp ln -sf $(pwd)/src/KOKKOS/pair_hybrid_overlay_mlml_kokkos.h ${lammps}/src/KOKKOS/pair_hybrid_overlay_mlml_kokkos.h ln -sf $(pwd)/src/KOKKOS/pair_hybrid_overlay_mlml_kokkos.cpp ${lammps}/src/KOKKOS/pair_hybrid_overlay_mlml_kokkos.cpp + +ln -sf $(pwd)/src/pair_eam_mlml.h ${lammps}/src/pair_eam_mlml.h +ln -sf $(pwd)/src/pair_eam_mlml.cpp ${lammps}/src/pair_eam_mlml.cpp +ln -sf $(pwd)/src/pair_eam_fs_mlml.h ${lammps}/src/pair_eam_fs_mlml.h +ln -sf $(pwd)/src/pair_eam_fs_mlml.cpp ${lammps}/src/pair_eam_fs_mlml.cpp diff --git a/LAMMPS_plugin/src/pair_eam_fs_mlml.cpp b/LAMMPS_plugin/src/pair_eam_fs_mlml.cpp new file mode 100644 index 0000000..cddcd74 --- /dev/null +++ b/LAMMPS_plugin/src/pair_eam_fs_mlml.cpp @@ -0,0 +1,372 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, 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 authors: + originally + Tim Lau (MIT) + removing communication within compute + David Immel (d.immel@fz-juelich.de, FZJ, Germany) +------------------------------------------------------------------------- */ + +#include "pair_eam_fs_mlml.h" + +#include "atom.h" +#include "comm.h" +#include "error.h" +#include "memory.h" +#include "potential_file_reader.h" + +#include + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +PairEAMFSMLML::PairEAMFSMLML(LAMMPS *lmp) : PairEAMMLML(lmp) +{ + one_coeff = 1; + manybody_flag = 1; + he_flag = 0; +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs + read EAM Finnis-Sinclair file +------------------------------------------------------------------------- */ + +void PairEAMFSMLML::coeff(int narg, char **arg) +{ + int i,j; + + if (!allocated) allocate(); + + if (narg != 3 + atom->ntypes) + error->all(FLERR,"Incorrect args for pair coefficients"); + + // insure I,J args are * * + + if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0) + error->all(FLERR,"Incorrect args for pair coefficients"); + + // read EAM Finnis-Sinclair file + + if (fs) { + for (i = 0; i < fs->nelements; i++) delete [] fs->elements[i]; + delete [] fs->elements; + memory->destroy(fs->mass); + memory->destroy(fs->frho); + memory->destroy(fs->rhor); + memory->destroy(fs->z2r); + delete fs; + } + fs = new Fs(); + read_file(arg[2]); + + // read args that map atom types to elements in potential file + // map[i] = which element the Ith atom type is, -1 if "NULL" + + for (i = 3; i < narg; i++) { + if (strcmp(arg[i],"NULL") == 0) { + map[i-2] = -1; + continue; + } + for (j = 0; j < fs->nelements; j++) + if (strcmp(arg[i],fs->elements[j]) == 0) break; + if (j < fs->nelements) map[i-2] = j; + else error->all(FLERR,"No matching element in EAM potential file"); + } + + // clear setflag since coeff() called once with I,J = * * + + int n = atom->ntypes; + for (i = 1; i <= n; i++) + for (j = i; j <= n; j++) + setflag[i][j] = 0; + + // set setflag i,j for type pairs where both are mapped to elements + // set mass of atom type if i = j + + int count = 0; + for (i = 1; i <= n; i++) { + for (j = i; j <= n; j++) { + if (map[i] >= 0 && map[j] >= 0) { + setflag[i][j] = 1; + if (i == j) atom->set_mass(FLERR,i,fs->mass[map[i]]); + count++; + } + scale[i][j] = 1.0; + } + } + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + read a multi-element DYNAMO setfl file +------------------------------------------------------------------------- */ + +void PairEAMFSMLML::read_file(char *filename) +{ + Fs *file = fs; + + // read potential file + if (comm->me == 0) { + PotentialFileReader reader(lmp, filename, he_flag ? "eam/he" : "eam/fs", + unit_convert_flag); + + // transparently convert units for supported conversions + + int unit_convert = reader.get_unit_convert(); + double conversion_factor = utils::get_conversion_factor(utils::ENERGY, + unit_convert); + try { + reader.skip_line(); + reader.skip_line(); + reader.skip_line(); + + // extract element names from nelements line + ValueTokenizer values = reader.next_values(1); + file->nelements = values.next_int(); + + if ((int)values.count() != file->nelements + 1) + error->one(FLERR,"Incorrect element names in EAM potential file"); + + file->elements = new char*[file->nelements]; + for (int i = 0; i < file->nelements; i++) { + const std::string word = values.next_string(); + const int n = word.length() + 1; + file->elements[i] = new char[n]; + strcpy(file->elements[i], word.c_str()); + } + + // + + if (he_flag) values = reader.next_values(6); + else values = reader.next_values(5); + file->nrho = values.next_int(); + file->drho = values.next_double(); + file->nr = values.next_int(); + file->dr = values.next_double(); + file->cut = values.next_double(); + if (he_flag) rhomax = values.next_double(); + + if ((file->nrho <= 0) || (file->nr <= 0) || (file->dr <= 0.0)) + error->one(FLERR,"Invalid EAM potential file"); + + memory->create(file->mass, file->nelements, "pair:mass"); + memory->create(file->frho, file->nelements, file->nrho + 1, "pair:frho"); + memory->create(file->rhor, file->nelements, file->nelements, file->nr + 1, "pair:rhor"); + memory->create(file->z2r, file->nelements, file->nelements, file->nr + 1, "pair:z2r"); + + for (int i = 0; i < file->nelements; i++) { + values = reader.next_values(2); + values.next_int(); // ignore + file->mass[i] = values.next_double(); + + reader.next_dvector(&file->frho[i][1], file->nrho); + if (unit_convert) { + for (int j = 1; j <= file->nrho; ++j) + file->frho[i][j] *= conversion_factor; + } + + for (int j = 0; j < file->nelements; j++) { + reader.next_dvector(&file->rhor[i][j][1], file->nr); + } + } + + for (int i = 0; i < file->nelements; i++) { + for (int j = 0; j <= i; j++) { + reader.next_dvector(&file->z2r[i][j][1], file->nr); + if (unit_convert) { + for (int k = 1; k <= file->nr; ++k) + file->z2r[i][j][k] *= conversion_factor; + } + } + } + } catch (TokenizerException &e) { + error->one(FLERR, e.what()); + } + } + + // broadcast potential information + MPI_Bcast(&file->nelements, 1, MPI_INT, 0, world); + + MPI_Bcast(&file->nrho, 1, MPI_INT, 0, world); + MPI_Bcast(&file->drho, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&file->nr, 1, MPI_INT, 0, world); + MPI_Bcast(&file->dr, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&file->cut, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&rhomax, 1, MPI_DOUBLE, 0, world); + + // allocate memory on other procs + if (comm->me != 0) { + file->elements = new char*[file->nelements]; + for (int i = 0; i < file->nelements; i++) file->elements[i] = nullptr; + memory->create(file->mass, file->nelements, "pair:mass"); + memory->create(file->frho, file->nelements, file->nrho + 1, "pair:frho"); + memory->create(file->rhor, file->nelements, file->nelements, file->nr + 1, "pair:rhor"); + memory->create(file->z2r, file->nelements, file->nelements, file->nr + 1, "pair:z2r"); + } + + // broadcast file->elements string array + for (int i = 0; i < file->nelements; i++) { + int n; + if (comm->me == 0) n = strlen(file->elements[i]) + 1; + MPI_Bcast(&n, 1, MPI_INT, 0, world); + if (comm->me != 0) file->elements[i] = new char[n]; + MPI_Bcast(file->elements[i], n, MPI_CHAR, 0, world); + } + + // broadcast file->mass, frho, rhor + for (int i = 0; i < file->nelements; i++) { + MPI_Bcast(&file->mass[i], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&file->frho[i][1], file->nrho, MPI_DOUBLE, 0, world); + + for (int j = 0; j < file->nelements; j++) { + MPI_Bcast(&file->rhor[i][j][1], file->nr, MPI_DOUBLE, 0, world); + } + } + + // broadcast file->z2r + for (int i = 0; i < file->nelements; i++) { + for (int j = 0; j <= i; j++) { + MPI_Bcast(&file->z2r[i][j][1], file->nr, MPI_DOUBLE, 0, world); + } + } +} + +/* ---------------------------------------------------------------------- + copy read-in setfl potential to standard array format +------------------------------------------------------------------------- */ + +void PairEAMFSMLML::file2array() +{ + int i,j,m,n; + int ntypes = atom->ntypes; + + // set function params directly from fs file + + nrho = fs->nrho; + nr = fs->nr; + drho = fs->drho; + dr = fs->dr; + if (he_flag) rhomin = rhomax - (nrho-1) * drho; + else rhomax = (nrho-1) * drho; + + // ------------------------------------------------------------------ + // setup frho arrays + // ------------------------------------------------------------------ + + // allocate frho arrays + // nfrho = # of fs elements + 1 for zero array + + nfrho = fs->nelements + 1; + memory->destroy(frho); + memory->create(frho,nfrho,nrho+1,"pair:frho"); + + // copy each element's frho to global frho + + for (i = 0; i < fs->nelements; i++) + for (m = 1; m <= nrho; m++) frho[i][m] = fs->frho[i][m]; + + // add extra frho of zeroes for non-EAM types to point to (pair hybrid) + // this is necessary b/c fp is still computed for non-EAM atoms + + for (m = 1; m <= nrho; m++) frho[nfrho-1][m] = 0.0; + + // type2frho[i] = which frho array (0 to nfrho-1) each atom type maps to + // if atom type doesn't point to element (non-EAM atom in pair hybrid) + // then map it to last frho array of zeroes + + for (i = 1; i <= ntypes; i++) + if (map[i] >= 0) type2frho[i] = map[i]; + else type2frho[i] = nfrho-1; + + // ------------------------------------------------------------------ + // setup rhor arrays + // ------------------------------------------------------------------ + + // allocate rhor arrays + // nrhor = square of # of fs elements + + nrhor = fs->nelements * fs->nelements; + memory->destroy(rhor); + memory->create(rhor,nrhor,nr+1,"pair:rhor"); + + // copy each element pair rhor to global rhor + + n = 0; + for (i = 0; i < fs->nelements; i++) + for (j = 0; j < fs->nelements; j++) { + for (m = 1; m <= nr; m++) rhor[n][m] = fs->rhor[i][j][m]; + n++; + } + + // type2rhor[i][j] = which rhor array (0 to nrhor-1) each type pair maps to + // for fs files, there is a full NxN set of rhor arrays + // OK if map = -1 (non-EAM atom in pair hybrid) b/c type2rhor not used + + for (i = 1; i <= ntypes; i++) + for (j = 1; j <= ntypes; j++) + type2rhor[i][j] = map[i] * fs->nelements + map[j]; + + // ------------------------------------------------------------------ + // setup z2r arrays + // ------------------------------------------------------------------ + + // allocate z2r arrays + // nz2r = N*(N+1)/2 where N = # of fs elements + + nz2r = fs->nelements * (fs->nelements+1) / 2; + memory->destroy(z2r); + memory->create(z2r,nz2r,nr+1,"pair:z2r"); + + // copy each element pair z2r to global z2r, only for I >= J + + n = 0; + for (i = 0; i < fs->nelements; i++) + for (j = 0; j <= i; j++) { + for (m = 1; m <= nr; m++) z2r[n][m] = fs->z2r[i][j][m]; + n++; + } + + // type2z2r[i][j] = which z2r array (0 to nz2r-1) each type pair maps to + // set of z2r arrays only fill lower triangular Nelement matrix + // value = n = sum over rows of lower-triangular matrix until reach irow,icol + // swap indices when irow < icol to stay lower triangular + // if map = -1 (non-EAM atom in pair hybrid): + // type2z2r is not used by non-opt + // but set type2z2r to 0 since accessed by opt + + int irow,icol; + for (i = 1; i <= ntypes; i++) { + for (j = 1; j <= ntypes; j++) { + irow = map[i]; + icol = map[j]; + if (irow == -1 || icol == -1) { + type2z2r[i][j] = 0; + continue; + } + if (irow < icol) { + irow = map[j]; + icol = map[i]; + } + n = 0; + for (m = 0; m < irow; m++) n += m + 1; + n += icol; + type2z2r[i][j] = n; + } + } +} diff --git a/LAMMPS_plugin/src/pair_eam_fs_mlml.h b/LAMMPS_plugin/src/pair_eam_fs_mlml.h new file mode 100644 index 0000000..8d61430 --- /dev/null +++ b/LAMMPS_plugin/src/pair_eam_fs_mlml.h @@ -0,0 +1,44 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, 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 PAIR_CLASS +// clang-format off +PairStyle(eam/fs/mlml,PairEAMFSMLML); +// clang-format on +#else + +#ifndef LMP_PAIR_EAM_FS_MLML_H +#define LMP_PAIR_EAM_FS_MLML_H + +#include "pair_eam_mlml.h" + +namespace LAMMPS_NS { + +// need virtual public b/c of how eam/fs/opt inherits from it + +class PairEAMFSMLML : virtual public PairEAMMLML { + public: + PairEAMFSMLML(class LAMMPS *); + + void coeff(int, char **) override; + + protected: + void read_file(char *) override; + void file2array() override; + int he_flag; +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/LAMMPS_plugin/src/pair_eam_mlml.cpp b/LAMMPS_plugin/src/pair_eam_mlml.cpp new file mode 100644 index 0000000..ab3384e --- /dev/null +++ b/LAMMPS_plugin/src/pair_eam_mlml.cpp @@ -0,0 +1,838 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, 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 authors: + originally + Stephen Foiles (SNL), Murray Daw (SNL) + removing communication within compute + David Immel (d.immel@fz-juelich.de, FZJ, Germany) +------------------------------------------------------------------------- */ + +#include "pair_eam_mlml.h" + +#include "atom.h" +#include "comm.h" +#include "error.h" +#include "force.h" +#include "memory.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "potential_file_reader.h" +#include "update.h" + +#include +#include + +using namespace LAMMPS_NS; + +#define MAXLINE 1024 + +/* ---------------------------------------------------------------------- */ + +PairEAMMLML::PairEAMMLML(LAMMPS *lmp) : Pair(lmp) +{ + restartinfo = 0; + manybody_flag = 1; + unit_convert_flag = utils::get_supported_conversions(utils::ENERGY); + + nmax = 0; + rho = nullptr; + fp = nullptr; + numforce = nullptr; + type2frho = nullptr; + + nfuncfl = 0; + funcfl = nullptr; + + setfl = nullptr; + fs = nullptr; + + frho = nullptr; + rhor = nullptr; + z2r = nullptr; + scale = nullptr; + + rhomax = rhomin = 0.0; + + frho_spline = nullptr; + rhor_spline = nullptr; + z2r_spline = nullptr; +} + +/* ---------------------------------------------------------------------- + check if allocated, since class can be destructed when incomplete +------------------------------------------------------------------------- */ + +PairEAMMLML::~PairEAMMLML() +{ + if (copymode) return; + + memory->destroy(rho); + memory->destroy(fp); + memory->destroy(numforce); + + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + delete [] type2frho; + type2frho = nullptr; + memory->destroy(type2rhor); + memory->destroy(type2z2r); + memory->destroy(scale); + } + + if (funcfl) { + for (int i = 0; i < nfuncfl; i++) { + delete [] funcfl[i].file; + memory->destroy(funcfl[i].frho); + memory->destroy(funcfl[i].rhor); + memory->destroy(funcfl[i].zr); + } + memory->sfree(funcfl); + funcfl = nullptr; + } + + if (setfl) { + for (int i = 0; i < setfl->nelements; i++) delete [] setfl->elements[i]; + delete [] setfl->elements; + memory->destroy(setfl->mass); + memory->destroy(setfl->frho); + memory->destroy(setfl->rhor); + memory->destroy(setfl->z2r); + delete setfl; + setfl = nullptr; + } + + if (fs) { + for (int i = 0; i < fs->nelements; i++) delete [] fs->elements[i]; + delete [] fs->elements; + memory->destroy(fs->mass); + memory->destroy(fs->frho); + memory->destroy(fs->rhor); + memory->destroy(fs->z2r); + delete fs; + fs = nullptr; + } + + memory->destroy(frho); + memory->destroy(rhor); + memory->destroy(z2r); + + memory->destroy(frho_spline); + memory->destroy(rhor_spline); + memory->destroy(z2r_spline); +} + +/* ---------------------------------------------------------------------- */ + +void PairEAMMLML::compute(int eflag, int vflag) +{ + // start timers + int i,j,ii,jj,m,inum,jnum,itype,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair,fpair_cl; + double rsq,r,p,rhoip,rhojp,z2,z2p,recip,phip,psip,phi,psip_cl; + double *coeff; + int *ilist,*jlist,*numneigh,**firstneigh; + + evdwl = 0.0; + ev_init(eflag,vflag); + + // grow energy and fp arrays if necessary + // need to be atom->nmax in length + + if (atom->nmax > nmax) { + memory->destroy(rho); + memory->destroy(fp); + memory->destroy(numforce); + nmax = atom->nmax; + memory->create(rho,nmax,"pair:rho"); + memory->create(fp,nmax,"pair:fp"); + memory->create(numforce,nmax,"pair:numforce"); + } + + double **x = atom->x; + double **f = atom->f; + + int *type = atom->type; + int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; + int newton_pair = force->newton_pair; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // zero out density + + for (i = 0; i < nlocal; i++) rho[i] = 0.0; + + // rho = density at each atom + // loop over neighbors of my atoms + + 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; + + // avoid double counting due to full neighbour list for two local particles + // j < i implies j < nlocal + if (j < i && i < nlocal) { 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) { + jtype = type[j]; + p = sqrt(rsq)*rdr + 1.0; + m = static_cast (p); + m = MIN(m,nr-1); + p -= m; + p = MIN(p,1.0); + coeff = rhor_spline[type2rhor[jtype][itype]][m]; + rho[i] += ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]; + + // do not calculate rho for ghost atoms + if (j < nlocal) { + coeff = rhor_spline[type2rhor[itype][jtype]][m]; + rho[j] += ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]; + } + } + } + } + + // fp = derivative of embedding energy at each atom + // phi = embedding energy at each atom + // if rho > rhomax (e.g. due to close approach of two atoms), + // will exceed table, so add linear term to conserve energy + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + p = rho[i]*rdrho + 1.0; + m = static_cast (p); + m = MAX(1,MIN(m,nrho-1)); + p -= m; + p = MIN(p,1.0); + coeff = frho_spline[type2frho[type[i]]][m]; + fp[i] = (coeff[0]*p + coeff[1])*p + coeff[2]; + if (eflag) { + phi = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]; + if (rho[i] > rhomax) phi += fp[i] * (rho[i]-rhomax); + phi *= scale[type[i]][type[i]]; + ev_tally_full(i, 2.0 * phi, 0.0, 0.0, 0.0, 0.0, 0.0); + } + } + + // compute forces on each atom + // loop over neighbors of my atoms + + 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]; + numforce[i] = 0; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + // avoid double counting due to full neighbour list for two local particles + // j < i implies j < nlocal + if (j < i && i < nlocal) { 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) { + ++numforce[i]; + jtype = type[j]; + r = sqrt(rsq); + p = r*rdr + 1.0; + m = static_cast (p); + m = MIN(m,nr-1); + p -= m; + p = MIN(p,1.0); + + // rhoip = derivative of (density at atom j due to atom i) + // rhojp = derivative of (density at atom i due to atom j) + // phi = pair potential energy + // phip = phi' + // z2 = phi * r + // z2p = (phi * r)' = (phi' r) + phi + // 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 + // scale factor can be applied by thermodynamic integration + + coeff = rhor_spline[type2rhor[jtype][itype]][m]; + rhojp = (coeff[0]*p + coeff[1])*p + coeff[2]; + coeff = z2r_spline[type2z2r[itype][jtype]][m]; + z2p = (coeff[0]*p + coeff[1])*p + coeff[2]; + z2 = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]; + + recip = 1.0/r; + phi = z2*recip; + phip = z2p*recip - phi*recip; + if (j < nlocal) { + coeff = rhor_spline[type2rhor[itype][jtype]][m]; + rhoip = (coeff[0]*p + coeff[1])*p + coeff[2]; + psip = fp[i]*rhojp + fp[j]*rhoip + phip; + } else { + // processor of j calculates the remaining terms (compared to psip in if case) + psip = fp[i]*rhojp + phip*0.5; + } + + fpair = -scale[itype][jtype]*psip*recip; + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + + if (eflag) { + evdwl = scale[itype][jtype]*phi; + ev_tally_full(i, evdwl, 0.0, 0.0, 0.0, 0.0, 0.0); + if (j < nlocal) + ev_tally_full(j, evdwl, 0.0, 0.0, 0.0, 0.0, 0.0); + } + if(vflag) ev_tally(i, j, nlocal, newton_pair, 0.0, 0.0, fpair, delx, dely, delz); + } + } + } + + if (vflag_fdotr) virial_fdotr_compute(); +} + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairEAMMLML::allocate() +{ + allocated = 1; + int n = atom->ntypes; + + memory->create(setflag,n+1,n+1,"pair:setflag"); + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + setflag[i][j] = 0; + + memory->create(cutsq,n+1,n+1,"pair:cutsq"); + + delete[] map; + map = new int[n+1]; + for (int i = 1; i <= n; i++) map[i] = -1; + + type2frho = new int[n+1]; + memory->create(type2rhor,n+1,n+1,"pair:type2rhor"); + memory->create(type2z2r,n+1,n+1,"pair:type2z2r"); + memory->create(scale,n+1,n+1,"pair:scale"); +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairEAMMLML::settings(int narg, char ** arg) +{ + if (narg > 0) error->all(FLERR, "Illegal pair_style command"); +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs + read DYNAMO funcfl file +------------------------------------------------------------------------- */ + +void PairEAMMLML::coeff(int narg, char **arg) +{ + if (!allocated) allocate(); + + if (narg != 3) error->all(FLERR,"Incorrect args for pair coefficients"); + + // parse pair of atom types + + int ilo,ihi,jlo,jhi; + utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error); + utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error); + + // read funcfl file if hasn't already been read + // store filename in Funcfl data struct + + int ifuncfl; + for (ifuncfl = 0; ifuncfl < nfuncfl; ifuncfl++) + if (strcmp(arg[2],funcfl[ifuncfl].file) == 0) break; + + if (ifuncfl == nfuncfl) { + nfuncfl++; + funcfl = (Funcfl *) + memory->srealloc(funcfl,nfuncfl*sizeof(Funcfl),"pair:funcfl"); + read_file(arg[2]); + funcfl[ifuncfl].file = utils::strdup(arg[2]); + } + + // set setflag and map only for i,i type pairs + // set mass of atom type if i = j + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + if (i == j) { + setflag[i][i] = 1; + map[i] = ifuncfl; + atom->set_mass(FLERR,i,funcfl[ifuncfl].mass); + count++; + } + scale[i][j] = 1.0; + } + } + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairEAMMLML::init_style() +{ + // convert read-in file(s) to arrays and spline them + if (force->newton_pair == 0) error->all(FLERR, "Pair style eam/mlml requires newton pair on"); + + file2array(); + array2spline(); + + // communication during computation should be avoided + // -> do not exchange the derivative of the embedding function by default + neighbor->add_request(this, NeighConst::REQ_FULL); +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairEAMMLML::init_one(int i, int j) +{ + // single global cutoff = max of cut from all files read in + // for funcfl could be multiple files + // for setfl or fs, just one file + + if (setflag[i][j] == 0) scale[i][j] = 1.0; + scale[j][i] = scale[i][j]; + + if (funcfl) { + cutmax = 0.0; + for (int m = 0; m < nfuncfl; m++) + cutmax = MAX(cutmax,funcfl[m].cut); + } else if (setfl) cutmax = setfl->cut; + else if (fs) cutmax = fs->cut; + + cutforcesq = cutmax*cutmax; + + return cutmax; +} + +/* ---------------------------------------------------------------------- + read potential values from a DYNAMO single element funcfl file +------------------------------------------------------------------------- */ + +void PairEAMMLML::read_file(char *filename) +{ + Funcfl *file = &funcfl[nfuncfl-1]; + + // read potential file + if (comm->me == 0) { + PotentialFileReader reader(lmp, filename, "eam", unit_convert_flag); + + // transparently convert units for supported conversions + + int unit_convert = reader.get_unit_convert(); + double conversion_factor = utils::get_conversion_factor(utils::ENERGY, + unit_convert); + try { + reader.skip_line(); + + ValueTokenizer values = reader.next_values(2); + values.next_int(); // ignore + file->mass = values.next_double(); + + values = reader.next_values(5); + file->nrho = values.next_int(); + file->drho = values.next_double(); + file->nr = values.next_int(); + file->dr = values.next_double(); + file->cut = values.next_double(); + + if ((file->nrho <= 0) || (file->nr <= 0) || (file->dr <= 0.0)) + error->one(FLERR,"Invalid EAM potential file"); + + memory->create(file->frho, (file->nrho+1), "pair:frho"); + memory->create(file->rhor, (file->nr+1), "pair:rhor"); + memory->create(file->zr, (file->nr+1), "pair:zr"); + + reader.next_dvector(&file->frho[1], file->nrho); + reader.next_dvector(&file->zr[1], file->nr); + reader.next_dvector(&file->rhor[1], file->nr); + + if (unit_convert) { + const double sqrt_conv = sqrt(conversion_factor); + for (int i = 1; i <= file->nrho; ++i) + file->frho[i] *= conversion_factor; + for (int j = 1; j <= file->nr; ++j) + file->zr[j] *= sqrt_conv; + } + } catch (TokenizerException &e) { + error->one(FLERR, e.what()); + } + } + + MPI_Bcast(&file->mass, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&file->nrho, 1, MPI_INT, 0, world); + MPI_Bcast(&file->drho, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&file->nr, 1, MPI_INT, 0, world); + MPI_Bcast(&file->dr, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&file->cut, 1, MPI_DOUBLE, 0, world); + + if (comm->me != 0) { + memory->create(file->frho, (file->nrho+1), "pair:frho"); + memory->create(file->rhor, (file->nr+1), "pair:rhor"); + memory->create(file->zr, (file->nr+1), "pair:zr"); + } + + MPI_Bcast(&file->frho[1], file->nrho, MPI_DOUBLE, 0, world); + MPI_Bcast(&file->zr[1], file->nr, MPI_DOUBLE, 0, world); + MPI_Bcast(&file->rhor[1], file->nr, MPI_DOUBLE, 0, world); +} + +/* ---------------------------------------------------------------------- + convert read-in funcfl potential(s) to standard array format + interpolate all file values to a single grid and cutoff +------------------------------------------------------------------------- */ + +void PairEAMMLML::file2array() +{ + int i,j,k,m,n; + int ntypes = atom->ntypes; + double sixth = 1.0/6.0; + + // determine max function params from all active funcfl files + // active means some element is pointing at it via map + + int active; + double rmax; + dr = drho = rmax = rhomax = 0.0; + + for (int i = 0; i < nfuncfl; i++) { + active = 0; + for (j = 1; j <= ntypes; j++) + if (map[j] == i) active = 1; + if (active == 0) continue; + Funcfl *file = &funcfl[i]; + dr = MAX(dr,file->dr); + drho = MAX(drho,file->drho); + rmax = MAX(rmax,(file->nr-1) * file->dr); + rhomax = MAX(rhomax,(file->nrho-1) * file->drho); + } + + // set nr,nrho from cutoff and spacings + // 0.5 is for round-off in divide + + nr = static_cast (rmax/dr + 0.5); + nrho = static_cast (rhomax/drho + 0.5); + + // ------------------------------------------------------------------ + // setup frho arrays + // ------------------------------------------------------------------ + + // allocate frho arrays + // nfrho = # of funcfl files + 1 for zero array + + nfrho = nfuncfl + 1; + memory->destroy(frho); + memory->create(frho,nfrho,nrho+1,"pair:frho"); + + // interpolate each file's frho to a single grid and cutoff + + double r,p,cof1,cof2,cof3,cof4; + + n = 0; + for (i = 0; i < nfuncfl; i++) { + Funcfl *file = &funcfl[i]; + for (m = 1; m <= nrho; m++) { + r = (m-1)*drho; + p = r/file->drho + 1.0; + k = static_cast (p); + k = MIN(k,file->nrho-2); + k = MAX(k,2); + p -= k; + p = MIN(p,2.0); + cof1 = -sixth*p*(p-1.0)*(p-2.0); + cof2 = 0.5*(p*p-1.0)*(p-2.0); + cof3 = -0.5*p*(p+1.0)*(p-2.0); + cof4 = sixth*p*(p*p-1.0); + frho[n][m] = cof1*file->frho[k-1] + cof2*file->frho[k] + + cof3*file->frho[k+1] + cof4*file->frho[k+2]; + } + n++; + } + + // add extra frho of zeroes for non-EAM types to point to (pair hybrid) + // this is necessary b/c fp is still computed for non-EAM atoms + + for (m = 1; m <= nrho; m++) frho[nfrho-1][m] = 0.0; + + // type2frho[i] = which frho array (0 to nfrho-1) each atom type maps to + // if atom type doesn't point to file (non-EAM atom in pair hybrid) + // then map it to last frho array of zeroes + + for (i = 1; i <= ntypes; i++) + if (map[i] >= 0) type2frho[i] = map[i]; + else type2frho[i] = nfrho-1; + + // ------------------------------------------------------------------ + // setup rhor arrays + // ------------------------------------------------------------------ + + // allocate rhor arrays + // nrhor = # of funcfl files + + nrhor = nfuncfl; + memory->destroy(rhor); + memory->create(rhor,nrhor,nr+1,"pair:rhor"); + + // interpolate each file's rhor to a single grid and cutoff + + n = 0; + for (i = 0; i < nfuncfl; i++) { + Funcfl *file = &funcfl[i]; + for (m = 1; m <= nr; m++) { + r = (m-1)*dr; + p = r/file->dr + 1.0; + k = static_cast (p); + k = MIN(k,file->nr-2); + k = MAX(k,2); + p -= k; + p = MIN(p,2.0); + cof1 = -sixth*p*(p-1.0)*(p-2.0); + cof2 = 0.5*(p*p-1.0)*(p-2.0); + cof3 = -0.5*p*(p+1.0)*(p-2.0); + cof4 = sixth*p*(p*p-1.0); + rhor[n][m] = cof1*file->rhor[k-1] + cof2*file->rhor[k] + + cof3*file->rhor[k+1] + cof4*file->rhor[k+2]; + } + n++; + } + + // type2rhor[i][j] = which rhor array (0 to nrhor-1) each type pair maps to + // for funcfl files, I,J mapping only depends on I + // OK if map = -1 (non-EAM atom in pair hybrid) b/c type2rhor not used + + for (i = 1; i <= ntypes; i++) + for (j = 1; j <= ntypes; j++) + type2rhor[i][j] = map[i]; + + // ------------------------------------------------------------------ + // setup z2r arrays + // ------------------------------------------------------------------ + + // allocate z2r arrays + // nz2r = N*(N+1)/2 where N = # of funcfl files + + nz2r = nfuncfl*(nfuncfl+1)/2; + memory->destroy(z2r); + memory->create(z2r,nz2r,nr+1,"pair:z2r"); + + // create a z2r array for each file against other files, only for I >= J + // interpolate zri and zrj to a single grid and cutoff + // final z2r includes unit conversion of 27.2 eV/Hartree and 0.529 Ang/Bohr + + double zri,zrj; + + n = 0; + for (i = 0; i < nfuncfl; i++) { + Funcfl *ifile = &funcfl[i]; + for (j = 0; j <= i; j++) { + Funcfl *jfile = &funcfl[j]; + + for (m = 1; m <= nr; m++) { + r = (m-1)*dr; + + p = r/ifile->dr + 1.0; + k = static_cast (p); + k = MIN(k,ifile->nr-2); + k = MAX(k,2); + p -= k; + p = MIN(p,2.0); + cof1 = -sixth*p*(p-1.0)*(p-2.0); + cof2 = 0.5*(p*p-1.0)*(p-2.0); + cof3 = -0.5*p*(p+1.0)*(p-2.0); + cof4 = sixth*p*(p*p-1.0); + zri = cof1*ifile->zr[k-1] + cof2*ifile->zr[k] + + cof3*ifile->zr[k+1] + cof4*ifile->zr[k+2]; + + p = r/jfile->dr + 1.0; + k = static_cast (p); + k = MIN(k,jfile->nr-2); + k = MAX(k,2); + p -= k; + p = MIN(p,2.0); + cof1 = -sixth*p*(p-1.0)*(p-2.0); + cof2 = 0.5*(p*p-1.0)*(p-2.0); + cof3 = -0.5*p*(p+1.0)*(p-2.0); + cof4 = sixth*p*(p*p-1.0); + zrj = cof1*jfile->zr[k-1] + cof2*jfile->zr[k] + + cof3*jfile->zr[k+1] + cof4*jfile->zr[k+2]; + + z2r[n][m] = 27.2*0.529 * zri*zrj; + } + n++; + } + } + + // type2z2r[i][j] = which z2r array (0 to nz2r-1) each type pair maps to + // set of z2r arrays only fill lower triangular Nelement matrix + // value = n = sum over rows of lower-triangular matrix until reach irow,icol + // swap indices when irow < icol to stay lower triangular + // if map = -1 (non-EAM atom in pair hybrid): + // type2z2r is not used by non-opt + // but set type2z2r to 0 since accessed by opt + + int irow,icol; + for (i = 1; i <= ntypes; i++) { + for (j = 1; j <= ntypes; j++) { + irow = map[i]; + icol = map[j]; + if (irow == -1 || icol == -1) { + type2z2r[i][j] = 0; + continue; + } + if (irow < icol) { + irow = map[j]; + icol = map[i]; + } + n = 0; + for (m = 0; m < irow; m++) n += m + 1; + n += icol; + type2z2r[i][j] = n; + } + } +} + +/* ---------------------------------------------------------------------- */ + +void PairEAMMLML::array2spline() +{ + rdr = 1.0/dr; + rdrho = 1.0/drho; + + memory->destroy(frho_spline); + memory->destroy(rhor_spline); + memory->destroy(z2r_spline); + + memory->create(frho_spline,nfrho,nrho+1,7,"pair:frho"); + memory->create(rhor_spline,nrhor,nr+1,7,"pair:rhor"); + memory->create(z2r_spline,nz2r,nr+1,7,"pair:z2r"); + + for (int i = 0; i < nfrho; i++) + interpolate(nrho,drho,frho[i],frho_spline[i]); + + for (int i = 0; i < nrhor; i++) + interpolate(nr,dr,rhor[i],rhor_spline[i]); + + for (int i = 0; i < nz2r; i++) + interpolate(nr,dr,z2r[i],z2r_spline[i]); +} + +/* ---------------------------------------------------------------------- */ + +void PairEAMMLML::interpolate(int n, double delta, double *f, double **spline) +{ + for (int m = 1; m <= n; m++) spline[m][6] = f[m]; + + spline[1][5] = spline[2][6] - spline[1][6]; + spline[2][5] = 0.5 * (spline[3][6]-spline[1][6]); + spline[n-1][5] = 0.5 * (spline[n][6]-spline[n-2][6]); + spline[n][5] = spline[n][6] - spline[n-1][6]; + + for (int m = 3; m <= n-2; m++) + spline[m][5] = ((spline[m-2][6]-spline[m+2][6]) + + 8.0*(spline[m+1][6]-spline[m-1][6])) / 12.0; + + for (int m = 1; m <= n-1; m++) { + spline[m][4] = 3.0*(spline[m+1][6]-spline[m][6]) - + 2.0*spline[m][5] - spline[m+1][5]; + spline[m][3] = spline[m][5] + spline[m+1][5] - + 2.0*(spline[m+1][6]-spline[m][6]); + } + + spline[n][4] = 0.0; + spline[n][3] = 0.0; + + for (int m = 1; m <= n; m++) { + spline[m][2] = spline[m][5]/delta; + spline[m][1] = 2.0*spline[m][4]/delta; + spline[m][0] = 3.0*spline[m][3]/delta; + } +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based arrays +------------------------------------------------------------------------- */ + +double PairEAMMLML::memory_usage() +{ + double bytes = (double)maxeatom * sizeof(double); + bytes += (double)maxvatom*6 * sizeof(double); + bytes += (double)2 * nmax * sizeof(double); + return bytes; +} + +/* ---------------------------------------------------------------------- + swap fp array with one passed in by caller +------------------------------------------------------------------------- */ + +void PairEAMMLML::swap_eam(double *fp_caller, double **fp_caller_hold) +{ + double *tmp = fp; + fp = fp_caller; + *fp_caller_hold = tmp; +} + +/* ---------------------------------------------------------------------- */ + +void *PairEAMMLML::extract(const char *str, int &dim) +{ + dim = 2; + if (strcmp(str,"scale") == 0) return (void *) scale; + return nullptr; +} diff --git a/LAMMPS_plugin/src/pair_eam_mlml.h b/LAMMPS_plugin/src/pair_eam_mlml.h new file mode 100644 index 0000000..befd439 --- /dev/null +++ b/LAMMPS_plugin/src/pair_eam_mlml.h @@ -0,0 +1,109 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, 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 PAIR_CLASS +// clang-format off +PairStyle(eam/mlml,PairEAMMLML); +// clang-format on +#else + +#ifndef LMP_PAIR_EAM_MLML_H +#define LMP_PAIR_EAM_MLML_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairEAMMLML : public Pair { + public: + friend class FixSemiGrandCanonicalMC; // Alex Stukowski option + + // public variables so ATC package can access them + + double cutmax; + + // potentials as array data + + int nrho, nr; + int nfrho, nrhor, nz2r; + double **frho, **rhor, **z2r; + int *type2frho, **type2rhor, **type2z2r; + + // potentials in spline form used for force computation + + double dr, rdr, drho, rdrho, rhomax, rhomin; + double ***rhor_spline, ***frho_spline, ***z2r_spline; + + PairEAMMLML(class LAMMPS *); + ~PairEAMMLML() override; + void compute(int, int) override; + void settings(int, char **) override; + void coeff(int, char **) override; + void init_style() override; + double init_one(int, int) override; + void *extract(const char *, int &) override; + + double memory_usage() override; + void swap_eam(double *, double **) override; + + protected: + int nmax; // allocated size of per-atom arrays + double cutforcesq; + double **scale; + + // per-atom arrays + + double *rho, *fp; + int *numforce; + + // potentials as file data + + struct Funcfl { + char *file; + int nrho, nr; + double drho, dr, cut, mass; + double *frho, *rhor, *zr; + }; + Funcfl *funcfl; + int nfuncfl; + + struct Setfl { + char **elements; + int nelements, nrho, nr; + double drho, dr, cut; + double *mass; + double **frho, **rhor, ***z2r; + }; + Setfl *setfl; + + struct Fs { + char **elements; + int nelements, nrho, nr; + double drho, dr, cut; + double *mass; + double **frho, ***rhor, ***z2r; + }; + Fs *fs; + + virtual void allocate(); + virtual void array2spline(); + void interpolate(int, double, double *, double **); + + virtual void read_file(char *); + virtual void file2array(); +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/docs/source/commands/pair_eam_mlml.rst b/docs/source/commands/pair_eam_mlml.rst new file mode 100644 index 0000000..a0dc995 --- /dev/null +++ b/docs/source/commands/pair_eam_mlml.rst @@ -0,0 +1,125 @@ +.. index:: pair_style eam/mlml +.. index:: pair_style eam/fs/mlml + +pair_style eam/mlml command +============================= + +For constant precision optimised variant: *eam* + +pair_style eam/fs/mlml command +================================ + +For constant precision optimised variant: *eam/fs* + +Syntax +"""""" + +.. code-block:: LAMMPS + + pair_style eam/mlml + pair_style eam/fs/mlml + +Examples +"""""""" + +.. code-block:: LAMMPS + + pair_style hybrid/overlay/mlml pace eam/fs/mlml + pair_coeff * * pace 1 W.yace W + pair_coeff * * eam/fs/full 2 W.eam.fs W + + +Description +""""""""""" +Style *eam* computes pairwise interactions for metals and metal alloys +using embedded-atom method (EAM) potentials :ref:`(Daw) `. The total +energy :math:`E_i` of an atom :math:`i` is given by + +.. math:: + + E_i^\text{EAM} = F_\alpha \left(\sum_{j \neq i}\ \rho_\beta (r_{ij})\right) + + \frac{1}{2} \sum_{j \neq i} \phi_{\alpha\beta} (r_{ij}) + +where :math:`F` is the embedding energy which is a function of the atomic +electron density :math:`\rho`, :math:`\phi` is a pair potential interaction, +and :math:`\alpha` and :math:`\beta` are the element types of atoms +:math:`i` and :math:`j`. The multi-body nature of the EAM potential is a +result of the embedding energy term. Both summations in the formula are over +all neighbors :math:`j` of atom :math:`i` within the cutoff distance. +EAM is documented in detail in :doc:`pair_style eam `. + +This pair style is optimized for the usage with ML-MIX: Communication within +the force computation routine is prevented at the cost of some double +computations. Thereby, one can efficiently balance the computational load +of MLML-EAM simulations via :doc:`fix balance `. + +.. warning:: + + The double computations computed by this pair style are disadvantageous + when force-mixing via :doc:`pair hybrid/overlay/mlml ` + is not used. + +Force-mixing provides the force :math:`\pmb{F}_i` on an atom :math:`i` that +is interpolated between two models. + +.. math:: + + \pmb{F}_i^\text{MIX} = \lambda_i \pmb{F}_i^\text{(fast)} + (1-\lambda_i) + \pmb{F}_i^\text{(precise)}\,, + +where the switching parameter :math:`\lambda_i` is computed +dynamically during a simulation by :doc:`fix mlml `. + +The pair style *eam/fs/mlml* computes the force +:math:`\pmb{F}_i^\text{EAM}` and should be combined +with a precise method like +:doc:`pair_style pace ` via +:doc:`pair_style hybrid/overlay/mlml `. + +Mixing, shift, table, tail correction, restart, rRESPA info +""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" + +For atom type pairs I,J and I != J, where types I and J correspond to +two different element types, mixing is performed by LAMMPS as +described above with the individual styles. You never need to specify +a pair_coeff command with I != J arguments for the eam/mlml styles. + +This pair style does not support the :doc:`pair_modify ` +shift, table, and tail options. + +The eam/mlml pair styles do not write their information to :doc:`binary +restart files `, since it is stored in tabulated potential +files. Thus, you need to re-specify the pair_style and pair_coeff +commands in an input script that reads a restart file. + +The eam/mlml pair styles can only be used via the *pair* keyword of the +:doc:`run_style respa ` command. They do not support the +*inner*, *middle*, *outer* keywords. + +---------- + +Restrictions +"""""""""""" + +This pair styles are part of the MLML package. It is only enabled if +LAMMPS was built with that package. See the :doc:`Build package +` page for more info. + +Related commands +"""""""""""""""" + +:doc:`pair_style eam `, +:doc:`pair_style hybrid/overlay/mlml `, +:doc:`fix mlml ` + +Default +""""""" + +none + +---------- + +.. _Daw3: + +**(Daw)** Daw, Baskes, Phys Rev Lett, 50, 1285 (1983). +Daw, Baskes, Phys Rev B, 29, 6443 (1984).