diff --git a/doc/src/Modify_pair.rst b/doc/src/Modify_pair.rst index 64831e726f..84f06a3acd 100644 --- a/doc/src/Modify_pair.rst +++ b/doc/src/Modify_pair.rst @@ -46,6 +46,8 @@ Here is a brief list of some the class methods in the Pair class that +---------------------------------+------------------------------------------------------------------------+ | compute_inner/middle/outer | versions of compute used by rRESPA | +---------------------------------+------------------------------------------------------------------------+ +| compute_atomic_energy | energy of one atom, equivalent to per-atom energy | ++---------------------------------+------------------------------------------------------------------------+ | memory_usage | return estimated amount of memory used by the pair style | +---------------------------------+------------------------------------------------------------------------+ | modify_params | process arguments to pair_modify command | @@ -122,3 +124,5 @@ setting. +---------------------------------+-------------------------------------------------------------+---------+ | spinflag | 1 if compatible with spin kspace_style | 0 | +---------------------------------+-------------------------------------------------------------+---------+ +| atomic_energy_enable | 1 if compute_atomic_energy() routine exists | 0 | ++---------------------------------+-------------------------------------------------------------+---------+ diff --git a/doc/src/fix_sgcmc.rst b/doc/src/fix_sgcmc.rst index b9d933cec3..0a4156a886 100644 --- a/doc/src/fix_sgcmc.rst +++ b/doc/src/fix_sgcmc.rst @@ -30,7 +30,9 @@ Syntax N = number of times sampling window is moved during one MC cycle *window_size* frac frac = size of sampling window (must be between 0.5 and 1.0) - + *atomic/energy* yes/no + yes = use the atomic energy method to calculate energy changes + no = use the default method to calculate energy changes Examples """""""" @@ -127,6 +129,14 @@ The number of times the window is moved during a MC cycle is set using the parameter *window_moves* (see Sect. III.B in :ref:`Sadigh1 ` for details). +The *atomic/energy* keyword controls which method is used for calculating +the energy change when atom types are swapped. A value of *no* +uses the default method, see discussion below in Restrictions section. +A value of *yes* uses the atomic energy method, +if the method has been implemented for the LAMMPS energy model, +otherwise LAMMPS will exit with an error message. +So far this has only been implemented for EAM type potentials. + ------------ Restart, fix_modify, output, run start/stop, minimize info @@ -159,16 +169,26 @@ page for more info. This fix style requires an :doc:`atom style ` with per atom type masses. -At present the fix provides optimized subroutines for EAM type -potentials (see above) that calculate potential energy changes due to -*local* atom type swaps very efficiently. Other potentials are -supported by using the generic potential functions. This, however, will -lead to exceedingly slow simulations since it implies that the -energy of the *entire* system is recomputed at each MC trial step. If -other potentials are to be used it is strongly recommended to modify and -optimize the existing generic potential functions for this purpose. -Also, the generic energy calculation can not be used for parallel -execution i.e. it only works with a single MPI process. +The fix provides three methods for calculating the potential energy +change due to atom type swaps. For EAM type potentials, the default +method is a carefully optimized local energy change calculation that +is part of the source code for this fix. It takes advantage of the +specific computational and communication requirements of EAM. Customizing +the local method to handle other energy models such as Tersoff has been done, +but these cases are not supported in the public LAMMPS code. +For all other LAMMPS energy models, the default method calculates +the *total* potential energy of the system before and after each +atom type swap. This method does not depend on the details of the +energy model and so is guaranteed to be correct. It is also +orders of magnitude slower than the custom EAM calculation. +In addition, it can not be used with parallel execution i.e. only +a single MPI process is allowed. +The third method uses the *atomic/energy* keyword described above. +This allows parallel execution and it is also a local calculation, +making it only a bit slower than a fully-optimized local calculation. +So far, this has been implemented for EAM type potentials. +It is straightforward to extend this to other potentials, +requiring adding an atomic energy method to the pair style. ------------ @@ -180,6 +200,7 @@ The optional parameters default to the following values: * *randseed* = 324234 * *window_moves* = 8 * *window_size* = automatic +* *atomic/energy* = no ------------ diff --git a/src/MANYBODY/pair_eam.cpp b/src/MANYBODY/pair_eam.cpp index b5a0ca0f77..fcef67b985 100644 --- a/src/MANYBODY/pair_eam.cpp +++ b/src/MANYBODY/pair_eam.cpp @@ -39,6 +39,7 @@ PairEAM::PairEAM(LAMMPS *lmp) : Pair(lmp) { restartinfo = 0; manybody_flag = 1; + atomic_energy_enable = 1; embedstep = -1; unit_convert_flag = utils::get_supported_conversions(utils::ENERGY); @@ -340,6 +341,70 @@ void PairEAM::compute(int eflag, int vflag) if (vflag_fdotr) virial_fdotr_compute(); } +/********************************************************************* + * Calculates the atomic energy of atom i + *********************************************************************/ +double PairEAM::compute_atomic_energy(int i, NeighList *neighborList) +{ + double p; + int m; + double* coeff; + double Ei = 0.0; + double rhoi = 0.0; + + double xi = atom->x[i][0]; + double yi = atom->x[i][1]; + double zi = atom->x[i][2]; + int itype = atom->type[i]; + + // loop over all neighbors of the selected atom. + + int* jlist = neighborList->firstneigh[i]; + int jnum = neighborList->numneigh[i]; + + for(int jj = 0; jj < jnum; jj++) { + int j = jlist[jj]; + + double delx = xi - atom->x[j][0]; + double dely = yi - atom->x[j][1]; + double delz = zi - atom->x[j][2]; + double rsq = delx*delx + dely*dely + delz*delz; + if(rsq >= cutforcesq) continue; + + int jtype = atom->type[j]; + double r = sqrt(rsq); + + p = r * rdr + 1.0; + m = static_cast(p); + m = MIN(m, nr - 1); + p -= m; + p = MIN(p, 1.0); + + // sum pair energy ij + // divide by 2 to avoid double counting energy + + coeff = z2r_spline[type2z2r[jtype][itype]][m]; + double z2 = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]; + Ei += 0.5*z2 / r; + + // sum rho_ij to rho_i + coeff = rhor_spline[type2rhor[jtype][itype]][m]; + rhoi += ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]; + } + + // compute the change in embedding energy of atom i. + + p = rhoi * 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[itype]][m]; + Ei += ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]; + + return Ei; +} + /* ---------------------------------------------------------------------- allocate all arrays ------------------------------------------------------------------------- */ diff --git a/src/MANYBODY/pair_eam.h b/src/MANYBODY/pair_eam.h index 24221a07ce..6a112a3221 100644 --- a/src/MANYBODY/pair_eam.h +++ b/src/MANYBODY/pair_eam.h @@ -47,6 +47,7 @@ class PairEAM : public Pair { PairEAM(class LAMMPS *); ~PairEAM() override; void compute(int, int) override; + double compute_atomic_energy(int, NeighList *) override; void settings(int, char **) override; void coeff(int, char **) override; void init_style() override; diff --git a/src/MANYBODY/pair_extep.cpp b/src/MANYBODY/pair_extep.cpp index 623dc510c4..1d8242127d 100644 --- a/src/MANYBODY/pair_extep.cpp +++ b/src/MANYBODY/pair_extep.cpp @@ -161,7 +161,6 @@ void PairExTeP::SR_neigh() } } } - //printf("SR_neigh : N[%d] = %f\n",i,N[i]); ipage->vgot(n); if (ipage->status()) diff --git a/src/MANYBODY/pair_gw.h b/src/MANYBODY/pair_gw.h index 15ddaa21a8..f098bfa082 100644 --- a/src/MANYBODY/pair_gw.h +++ b/src/MANYBODY/pair_gw.h @@ -94,8 +94,6 @@ class PairGW : public Pair { const double gw_d = param->d * param->d; const double hcth = param->h - costheta; - //printf("gw_gijk: gw_c=%f gw_d=%f hcth=%f=%f-%f\n", gw_c, gw_d, hcth, param->h, costheta); - return param->gamma * (1.0 + gw_c / gw_d - gw_c / (gw_d + hcth * hcth)); } diff --git a/src/MC/fix_bond_break.cpp b/src/MC/fix_bond_break.cpp index 94ec5a89bb..f5b34ac356 100644 --- a/src/MC/fix_bond_break.cpp +++ b/src/MC/fix_bond_break.cpp @@ -439,11 +439,6 @@ void FixBondBreak::update_topology() ndihedrals = 0; nimpropers = 0; - //printf("NBREAK %d: ",nbreak); - //for (i = 0; i < nbreak; i++) - // printf(" %d %d,",broken[i][0],broken[i][1]); - //printf("\n"); - for (i = 0; i < nlocal; i++) { influenced = 0; slist = special[i]; diff --git a/src/MC/fix_sgcmc.cpp b/src/MC/fix_sgcmc.cpp index c405421b7f..c081d84f2d 100644 --- a/src/MC/fix_sgcmc.cpp +++ b/src/MC/fix_sgcmc.cpp @@ -124,6 +124,7 @@ FixSemiGrandCanonicalMC::FixSemiGrandCanonicalMC(LAMMPS *_lmp, int narg, char ** // Default values for optional parameters (where applicable). numSamplingWindowMoves = 8; seed = 324234; + atomicenergyflag = 0; // Parse extra/optional parameters while (iarg < narg) { @@ -175,7 +176,7 @@ FixSemiGrandCanonicalMC::FixSemiGrandCanonicalMC(LAMMPS *_lmp, int narg, char ** targetConcentration[i] = utils::numeric(FLERR, arg[iarg], false, lmp); targetConcentration[1] -= targetConcentration[i]; } - for (int i = 1; i <= atom->ntypes; i++, iarg++) { + for (int i = 1; i <= atom->ntypes; i++) { if ((targetConcentration[i] < 0.0) || (targetConcentration[i] > 1.0)) error->all(FLERR, "Target concentration {} for species {} is out of range", targetConcentration[i], i); @@ -183,7 +184,6 @@ FixSemiGrandCanonicalMC::FixSemiGrandCanonicalMC(LAMMPS *_lmp, int narg, char ** utils::logmesg(lmp, " SGC - Target concentration of species {}: {}\n", i, targetConcentration[i]); } - } else if (strcmp(arg[iarg], "serial") == 0) { // Switch off second rejection. serialMode = true; @@ -193,6 +193,12 @@ FixSemiGrandCanonicalMC::FixSemiGrandCanonicalMC(LAMMPS *_lmp, int narg, char ** if (comm->nprocs != 1) error->all(FLERR, "Cannot use serial mode Monte Carlo in a parallel simulation."); + + } else if (strcmp(arg[iarg],"atomic/energy") == 0) { + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, fmt::format("fix {} atomic/energy", style), error); + atomicenergyflag = utils::logical(FLERR,arg[iarg+1],false,lmp); + iarg += 2; + } else { error->all(FLERR, "Unknown fix sgcmc keyword: {}", arg[iarg]); } @@ -243,25 +249,40 @@ void FixSemiGrandCanonicalMC::init() if (modify->get_fix_by_style("sgcmc").size() > 1) error->all(FLERR, "More than one fix sgcmc defined."); - // Save a pointer to the EAM potential. - pairEAM = dynamic_cast(force->pair); - if (!pairEAM) { - if (comm->me == 0) - utils::logmesg(lmp, " SGC - Using naive total energy calculation for MC -> SLOW!\n"); + if (atomicenergyflag) { + // Save a pointer to the EAM potential. + if (comm->me == 0) utils::logmesg(lmp, " SGC - Using atomic energy method for SGCMC\n"); + if (!force->pair->atomic_energy_enable) { + error->all(FLERR, "SGC - Pair style does not support atomic energy method"); + } + } else { + // Save a pointer to the EAM potential. + pairEAM = dynamic_cast(force->pair); + if (!pairEAM) { - if (comm->nprocs > 1) - error->all(FLERR, "Can not run fix sgcmc with naive total energy calculation " - "and more than one MPI process."); + if (comm->me == 0) + utils::logmesg(lmp, " SGC - Using naive total energy calculation for MC -> SLOW!\n"); - // Get reference to a compute that will provide the total energy of the system. - // This is needed by computeTotalEnergy(). - compute_pe = modify->get_compute_by_id("thermo_pe"); + if (comm->nprocs > 1) + error->all(FLERR, "Can not run fix sgcmc with naive total energy calculation " + "and more than one MPI process."); + + // Get reference to a compute that will provide the total energy of the system. + // This is needed by computeTotalEnergy(). + compute_pe = modify->get_compute_by_id("thermo_pe"); + } } + interactionRadius = force->pair->cutforce; if (comm->me == 0) utils::logmesg(lmp, " SGC - Interaction radius: {}\n", interactionRadius); + // This fix needs a full neighbor list. - neighbor->add_request(this, NeighConst::REQ_FULL); + if (atomicenergyflag) + // for atomic energy method, need ghost neighbors + neighbor->add_request(this, NeighConst::REQ_FULL | NeighConst::REQ_GHOST); + else + neighbor->add_request(this, NeighConst::REQ_FULL); // Count local number of atoms from each species. const int *type = atom->type; @@ -379,11 +400,18 @@ void FixSemiGrandCanonicalMC::doMC() deltaN[newSpecies] = +1; // Compute the energy difference that swapping this atom would cost or gain. - if (pairEAM) { - deltaE = computeEnergyChangeEAM(selectedAtom, selectedAtomNL, oldSpecies, newSpecies); + + // Atomic energy method: + if(atomicenergyflag) { + deltaE = computeEnergyChangeEatom(selectedAtom, oldSpecies, newSpecies); } else { - // Generic case: - deltaE = computeEnergyChangeGeneric(selectedAtom, oldSpecies, newSpecies); + // EAM method: + if (pairEAM) { + deltaE = computeEnergyChangeEAM(selectedAtom, selectedAtomNL, oldSpecies, newSpecies); + } else { + // Generic method + deltaE = computeEnergyChangeGeneric(selectedAtom, oldSpecies, newSpecies); + } } // Perform inner MC acceptance test. @@ -440,10 +468,14 @@ void FixSemiGrandCanonicalMC::doMC() // Make accepted atom swap permanent. if (selectedAtom >= 0) { - if (pairEAM) - flipAtomEAM(selectedAtom, selectedAtomNL, oldSpecies, newSpecies); - else - flipAtomGeneric(selectedAtom, oldSpecies, newSpecies); + if(atomicenergyflag) { + flipAtomEatom(selectedAtom, oldSpecies, newSpecies); + } else { + if (pairEAM) + flipAtomEAM(selectedAtom, selectedAtomNL, oldSpecies, newSpecies); + else + flipAtomGeneric(selectedAtom, oldSpecies, newSpecies); + } nAcceptedSwapsLocal++; } else { nRejectedSwapsLocal++; @@ -989,6 +1021,67 @@ void FixSemiGrandCanonicalMC::flipAtomGeneric(int flipAtom, int oldSpecies, int changedAtoms[flipAtom] = true; } +/********************************************************************* + * Calculates the change in energy that swapping the given + * atom would produce. This routine uses a per-atom energy calculation + *********************************************************************/ + +double FixSemiGrandCanonicalMC::computeEnergyChangeEatom(int flipAtom, int oldSpecies, int newSpecies) +{ + double Eold, Enew, deltaE; + + // Calculate old atomic energy of selected atom + + Eold = force->pair->compute_atomic_energy(flipAtom, neighborList); + + // calculate the old per-atom energy of neighbors + + int* jlist = neighborList->firstneigh[flipAtom]; + int jnum = neighborList->numneigh[flipAtom]; + + for(int jj = 0; jj < jnum; jj++) { + int j = jlist[jj]; + Eold += force->pair->compute_atomic_energy(j, neighborList); + } + + // Calculate new per-atom energy of selected atom + + atom->type[flipAtom] = newSpecies; + + Enew = force->pair->compute_atomic_energy(flipAtom, neighborList); + + // calculate the new per-atom energy of neighbors + + for(int jj = 0; jj < jnum; jj++) { + int j = jlist[jj]; + Enew += force->pair->compute_atomic_energy(j, neighborList); + } + + atom->type[flipAtom] = oldSpecies; + + deltaE = Enew - Eold; + + return deltaE; +} + +/********************************************************************* + * Flips the type of one atom. + * This routine is for the per-atom energy case. + *********************************************************************/ + +void FixSemiGrandCanonicalMC::flipAtomEatom(int flipAtom, int oldSpecies, int newSpecies) +{ + atom->type[flipAtom] = newSpecies; + + // Rescale particle velocity vector to conserve kinetic energy. + double vScaleFactor = sqrt(atom->mass[oldSpecies] / atom->mass[newSpecies]); + atom->v[flipAtom][0] *= vScaleFactor; + atom->v[flipAtom][1] *= vScaleFactor; + atom->v[flipAtom][2] *= vScaleFactor; + + changedAtoms[flipAtom] = true; +} + /********************************************************************* * Lets the fix report one of its internal state variables to LAMMPS. *********************************************************************/ diff --git a/src/MC/fix_sgcmc.h b/src/MC/fix_sgcmc.h index 9ad3f65214..16ed23391f 100644 --- a/src/MC/fix_sgcmc.h +++ b/src/MC/fix_sgcmc.h @@ -75,6 +75,10 @@ class FixSemiGrandCanonicalMC : public Fix { // This routine should only be used for debugging purposes. double computeEnergyChangeGeneric(int flipAtom, int oldSpecies, int newSpecies); + // Calculates the change in energy that swapping the given atom would produce. + // This uses the atomic energy method + double computeEnergyChangeEatom(int flipAtom, int oldSpecies, int newSpecies); + // Lets LAMMPS calculate the total potential energy of the system. double computeTotalEnergy(); @@ -86,6 +90,10 @@ class FixSemiGrandCanonicalMC : public Fix { // This routine is for the generic case. void flipAtomGeneric(int flipAtom, int oldSpecies, int newSpecies); + // Flips the type of one atom. + // This routine is for the atomic energy method + void flipAtomEatom(int flipAtom, int oldSpecies, int newSpecies); + // Transfers the locally changed electron densities and atom types to the neighbors. void communicateRhoAndTypes(); @@ -181,6 +189,9 @@ class FixSemiGrandCanonicalMC : public Fix { // A compute used to compute the total potential energy of the system. class Compute *compute_pe; + + // Indicate whether or not atomic energy method is used + int atomicenergyflag; }; } // namespace LAMMPS_NS diff --git a/src/pair.cpp b/src/pair.cpp index 896957c087..78b2e13ee3 100644 --- a/src/pair.cpp +++ b/src/pair.cpp @@ -90,6 +90,8 @@ Pair::Pair(LAMMPS *lmp) : reinitflag = 1; centroidstressflag = CENTROID_SAME; + atomic_energy_enable = 0; + // pair_modify settings compute_flag = 1; diff --git a/src/pair.h b/src/pair.h index 461d0c121b..1a4cc1cbf3 100644 --- a/src/pair.h +++ b/src/pair.h @@ -54,6 +54,8 @@ class Pair : protected Pointers { int single_enable; // 1 if single() routine exists int born_matrix_enable; // 1 if born_matrix() routine exists int single_hessian_enable; // 1 if single_hessian() routine exists + int atomic_energy_enable; // 1 if compute_atomic_energy() routine exists + int restartinfo; // 1 if pair style writes restart info int respa_enable; // 1 if inner/middle/outer rRESPA routines int one_coeff; // 1 if allows only one coeff * * call @@ -156,6 +158,7 @@ class Pair : protected Pointers { virtual void compute_inner() {} virtual void compute_middle() {} virtual void compute_outer(int, int) {} + virtual double compute_atomic_energy(int, NeighList *) { return 0.0; } virtual double single(int, int, int, int, double, double, double, double &fforce) { diff --git a/tools/coding_standard/README b/tools/coding_standard/README index b19d7241f4..da4552e081 100644 --- a/tools/coding_standard/README +++ b/tools/coding_standard/README @@ -6,3 +6,8 @@ whitespace.py detects TAB characters and trailing whitespace homepage.py detects outdated LAMMPS homepage URLs (pointing to sandia.gov instead of lammps.org) errordocs.py detects deprecated error docs in header files versiontags.py detects .. versionadded:: and .. versionchanged:: with pending version date + +Example usage: python3 versiontags.py + +where is the path of a LAMMPS directory containing subdirs src/, python/, doc/src/, etc. +