ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'develop' into talinke/develop

Browse Source
This commit is contained in:
Axel Kohlmeyer
2025-04-29 17:20:15 -04:00
parent 1062a8fb3a 440e24c60e
commit 0ce5d68455
12 changed files with 238 additions and 41 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
doc/src/Modify_pair.rst
View File

@ -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 |
+---------------------------------+-------------------------------------------------------------+---------+

43
doc/src/fix_sgcmc.rst
View File

@ -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
<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 <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
------------

65
src/MANYBODY/pair_eam.cpp
View File

@ -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<int>(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<int>(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
------------------------------------------------------------------------- */

1
src/MANYBODY/pair_eam.h
View File

@ -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;

1
src/MANYBODY/pair_extep.cpp
View File

@ -161,7 +161,6 @@ void PairExTeP::SR_neigh()
}
}
}
//printf("SR_neigh : N[%d] = %f\n",i,N[i]);
ipage->vgot(n);
if (ipage->status())

2
src/MANYBODY/pair_gw.h
View File

@ -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));
}

5
src/MC/fix_bond_break.cpp
View File

@ -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];

137
src/MC/fix_sgcmc.cpp
View File

@ -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<PairEAM*>(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<PairEAM*>(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.
*********************************************************************/

11
src/MC/fix_sgcmc.h
View File

@ -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

2
src/pair.cpp
View File

@ -90,6 +90,8 @@ Pair::Pair(LAMMPS *lmp) :
reinitflag = 1;
centroidstressflag = CENTROID_SAME;
atomic_energy_enable = 0;
// pair_modify settings
compute_flag = 1;

3
src/pair.h
View File

@ -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)
{

5
tools/coding_standard/README
View File

@ -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 <DIRECTORY>
where <DIRECTORY> is the path of a LAMMPS directory containing subdirs src/, python/, doc/src/, etc.