From fc0a41fb71f49d115d9c18a3a75bf5b2db50b9d1 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Sat, 19 Apr 2025 18:15:43 -0600 Subject: [PATCH 1/7] Added atomic/energy keyword, EAM support, correct for examples/MC/in.sgcmc.eam in serial --- doc/src/fix_sgcmc.rst | 42 ++++++++--- src/MANYBODY/pair_eam.cpp | 65 +++++++++++++++++ src/MANYBODY/pair_eam.h | 1 + src/MC/fix_sgcmc.cpp | 145 ++++++++++++++++++++++++++++++++------ src/MC/fix_sgcmc.h | 11 +++ src/pair.cpp | 2 + src/pair.h | 3 + 7 files changed, 235 insertions(+), 34 deletions(-) diff --git a/doc/src/fix_sgcmc.rst b/doc/src/fix_sgcmc.rst index b9d933cec3..628c1e3760 100644 --- a/doc/src/fix_sgcmc.rst +++ b/doc/src/fix_sgcmc.rst @@ -30,7 +30,8 @@ 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 Examples """""""" @@ -127,6 +128,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 for this has only been implemented for EAM type potentials. + ------------ Restart, fix_modify, output, run start/stop, minimize info @@ -159,16 +168,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 +199,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..2b3ec0fcad 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/MC/fix_sgcmc.cpp b/src/MC/fix_sgcmc.cpp index c405421b7f..4f2c48f45b 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) { @@ -162,6 +163,7 @@ FixSemiGrandCanonicalMC::FixSemiGrandCanonicalMC(LAMMPS *_lmp, int narg, char ** if (iarg + 1 + atom->ntypes > narg) utils::missing_cmd_args(FLERR, "fix sgcmc variance", error); iarg++; + // printf("variance iarg = %d narg = %d\n", iarg, narg); kappa = utils::numeric(FLERR, arg[iarg], false, lmp); if (kappa < 0) error->all(FLERR, "Variance constraint parameter must not be negative."); @@ -175,7 +177,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 +185,7 @@ FixSemiGrandCanonicalMC::FixSemiGrandCanonicalMC(LAMMPS *_lmp, int narg, char ** utils::logmesg(lmp, " SGC - Target concentration of species {}: {}\n", i, targetConcentration[i]); } - + // printf("variance iarg = %d narg = %d\n", iarg, narg); } else if (strcmp(arg[iarg], "serial") == 0) { // Switch off second rejection. serialMode = true; @@ -193,6 +195,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 +251,41 @@ 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,13 +403,20 @@ 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. double dm = 0.0; if (serialMode && kappa != 0.0) { @@ -440,10 +471,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 +1024,70 @@ 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; + + // printf("In computeEnergyChangeEatom()\n"); + // 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; + + // printf("In computeEnergyChangeEatom() deltaE = %g\n", deltaE); + + return deltaE; +} + +/********************************************************************* + * + * + *********************************************************************/ + +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..837f532209 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 routne 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..16cb5fb08f 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) { From 75d620c0a77987901fbadf80dc1566270e8b1564 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Sat, 19 Apr 2025 18:21:39 -0600 Subject: [PATCH 2/7] Fixed doc page --- doc/src/fix_sgcmc.rst | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/doc/src/fix_sgcmc.rst b/doc/src/fix_sgcmc.rst index 628c1e3760..0a4156a886 100644 --- a/doc/src/fix_sgcmc.rst +++ b/doc/src/fix_sgcmc.rst @@ -32,6 +32,7 @@ Syntax 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 """""""" @@ -134,7 +135,7 @@ 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 for this has only been implemented for EAM type potentials. +So far this has only been implemented for EAM type potentials. ------------ From d5c4f9c158fde526fa436c98aa873b5b6a2e8570 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Sun, 20 Apr 2025 11:01:38 -0600 Subject: [PATCH 3/7] Fixed whitespace --- src/MANYBODY/pair_eam.cpp | 6 +++--- src/MC/fix_sgcmc.cpp | 25 ++++++++++++------------- src/pair.cpp | 2 +- tools/coding_standard/README | 5 +++++ 4 files changed, 21 insertions(+), 17 deletions(-) diff --git a/src/MANYBODY/pair_eam.cpp b/src/MANYBODY/pair_eam.cpp index 2b3ec0fcad..fcef67b985 100644 --- a/src/MANYBODY/pair_eam.cpp +++ b/src/MANYBODY/pair_eam.cpp @@ -351,7 +351,7 @@ double PairEAM::compute_atomic_energy(int i, NeighList *neighborList) 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]; @@ -382,7 +382,7 @@ double PairEAM::compute_atomic_energy(int i, NeighList *neighborList) // 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; @@ -393,7 +393,7 @@ double PairEAM::compute_atomic_energy(int i, NeighList *neighborList) } // 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)); diff --git a/src/MC/fix_sgcmc.cpp b/src/MC/fix_sgcmc.cpp index 4f2c48f45b..cc990f375a 100644 --- a/src/MC/fix_sgcmc.cpp +++ b/src/MC/fix_sgcmc.cpp @@ -195,12 +195,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]); } @@ -274,7 +274,7 @@ void FixSemiGrandCanonicalMC::init() 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); @@ -285,7 +285,6 @@ void FixSemiGrandCanonicalMC::init() 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; @@ -403,7 +402,7 @@ void FixSemiGrandCanonicalMC::doMC() deltaN[newSpecies] = +1; // Compute the energy difference that swapping this atom would cost or gain. - + // Atomic energy method: if(atomicenergyflag) { deltaE = computeEnergyChangeEatom(selectedAtom, oldSpecies, newSpecies); @@ -416,7 +415,7 @@ void FixSemiGrandCanonicalMC::doMC() deltaE = computeEnergyChangeGeneric(selectedAtom, oldSpecies, newSpecies); } } - + // Perform inner MC acceptance test. double dm = 0.0; if (serialMode && kappa != 0.0) { @@ -1037,9 +1036,9 @@ double FixSemiGrandCanonicalMC::computeEnergyChangeEatom(int flipAtom, int oldSp // 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]; @@ -1053,14 +1052,14 @@ double FixSemiGrandCanonicalMC::computeEnergyChangeEatom(int flipAtom, int oldSp 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; @@ -1071,8 +1070,8 @@ double FixSemiGrandCanonicalMC::computeEnergyChangeEatom(int flipAtom, int oldSp } /********************************************************************* - * - * + * Flips the type of one atom. + * This routine is for the per-atom energy case. *********************************************************************/ void FixSemiGrandCanonicalMC::flipAtomEatom(int flipAtom, int oldSpecies, int newSpecies) diff --git a/src/pair.cpp b/src/pair.cpp index 16cb5fb08f..78b2e13ee3 100644 --- a/src/pair.cpp +++ b/src/pair.cpp @@ -91,7 +91,7 @@ Pair::Pair(LAMMPS *lmp) : centroidstressflag = CENTROID_SAME; atomic_energy_enable = 0; - + // pair_modify settings compute_flag = 1; 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. + From 7bfc15ff2c1a847e96fa460ac7f9634178cfb6d9 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Mon, 21 Apr 2025 10:06:50 -0600 Subject: [PATCH 4/7] Address akohlmey comments --- doc/src/Modify_pair.rst | 4 ++++ src/MC/fix_sgcmc.cpp | 5 ----- 2 files changed, 4 insertions(+), 5 deletions(-) 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/src/MC/fix_sgcmc.cpp b/src/MC/fix_sgcmc.cpp index cc990f375a..c081d84f2d 100644 --- a/src/MC/fix_sgcmc.cpp +++ b/src/MC/fix_sgcmc.cpp @@ -163,7 +163,6 @@ FixSemiGrandCanonicalMC::FixSemiGrandCanonicalMC(LAMMPS *_lmp, int narg, char ** if (iarg + 1 + atom->ntypes > narg) utils::missing_cmd_args(FLERR, "fix sgcmc variance", error); iarg++; - // printf("variance iarg = %d narg = %d\n", iarg, narg); kappa = utils::numeric(FLERR, arg[iarg], false, lmp); if (kappa < 0) error->all(FLERR, "Variance constraint parameter must not be negative."); @@ -185,7 +184,6 @@ FixSemiGrandCanonicalMC::FixSemiGrandCanonicalMC(LAMMPS *_lmp, int narg, char ** utils::logmesg(lmp, " SGC - Target concentration of species {}: {}\n", i, targetConcentration[i]); } - // printf("variance iarg = %d narg = %d\n", iarg, narg); } else if (strcmp(arg[iarg], "serial") == 0) { // Switch off second rejection. serialMode = true; @@ -1032,7 +1030,6 @@ double FixSemiGrandCanonicalMC::computeEnergyChangeEatom(int flipAtom, int oldSp { double Eold, Enew, deltaE; - // printf("In computeEnergyChangeEatom()\n"); // Calculate old atomic energy of selected atom Eold = force->pair->compute_atomic_energy(flipAtom, neighborList); @@ -1064,8 +1061,6 @@ double FixSemiGrandCanonicalMC::computeEnergyChangeEatom(int flipAtom, int oldSp deltaE = Enew - Eold; - // printf("In computeEnergyChangeEatom() deltaE = %g\n", deltaE); - return deltaE; } From 04732e2efdf26ae20c57395305e7ce41c9100ae4 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Mon, 21 Apr 2025 10:08:27 -0600 Subject: [PATCH 5/7] Removed debug code from unrelated files in MANYBODY and MC --- src/MANYBODY/pair_extep.cpp | 1 - src/MC/fix_bond_break.cpp | 5 ----- 2 files changed, 6 deletions(-) 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/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]; From 385f35091810069c2be0c33d632c3390508c5a8d Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Mon, 21 Apr 2025 15:50:14 -0600 Subject: [PATCH 6/7] Removed debug code from unrelated files in MANYBODY and MC --- src/MANYBODY/pair_gw.h | 2 -- 1 file changed, 2 deletions(-) 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)); } From 2c27ea3706fbc1c838044f0fa6bbee21ae3a7b5a Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 26 Apr 2025 13:32:21 -0400 Subject: [PATCH 7/7] Fix spelling in comment Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- src/MC/fix_sgcmc.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/MC/fix_sgcmc.h b/src/MC/fix_sgcmc.h index 837f532209..16ed23391f 100644 --- a/src/MC/fix_sgcmc.h +++ b/src/MC/fix_sgcmc.h @@ -91,7 +91,7 @@ class FixSemiGrandCanonicalMC : public Fix { void flipAtomGeneric(int flipAtom, int oldSpecies, int newSpecies); // Flips the type of one atom. - // This routne is for the atomic energy method + // 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.