diff --git a/src/MC/fix_gcmc.cpp b/src/MC/fix_gcmc.cpp index 80d2b3d48f..7da19241e9 100644 --- a/src/MC/fix_gcmc.cpp +++ b/src/MC/fix_gcmc.cpp @@ -234,6 +234,7 @@ void FixGCMC::options(int narg, char **arg) grouptypestrings = NULL; grouptypes = NULL; grouptypebits = NULL; + energy_intra = 0.0; int iarg = 0; while (iarg < narg) { @@ -320,6 +321,10 @@ void FixGCMC::options(int narg, char **arg) strcpy(grouptypestrings[ngrouptypes],arg[iarg+2]); ngrouptypes++; iarg += 3; + } else if (strcmp(arg[iarg],"intra_energy") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix gcmc command"); + energy_intra = force->numeric(FLERR,arg[iarg+1]); + iarg += 2; } else error->all(FLERR,"Illegal fix gcmc command"); } } @@ -1674,8 +1679,11 @@ void FixGCMC::attempt_molecule_deletion_full() if (force->kspace) force->kspace->qsum_qsq(); double energy_after = energy_full(); - if (random_equal->uniform() < - ngas*exp(beta*(energy_before - energy_after))/(zz*volume*natoms_per_molecule)) { + // energy_before corrected by energy_intra + + double deltaphi = ngas*exp(beta*((energy_before - energy_intra) - energy_after))/(zz*volume*natoms_per_molecule); + + if (random_equal->uniform() < deltaphi) { int i = 0; while (i < atom->nlocal) { if (atom->molecule[i] == deletion_molecule) { @@ -1827,8 +1835,12 @@ void FixGCMC::attempt_molecule_insertion_full() if (force->kspace) force->kspace->qsum_qsq(); double energy_after = energy_full(); - if (random_equal->uniform() < zz*volume*natoms_per_molecule* - exp(beta*(energy_before - energy_after))/(ngas + natoms_per_molecule)) { + // energy_after corrected by energy_intra + + double deltaphi = zz*volume*natoms_per_molecule* + exp(beta*(energy_before - (energy_after - energy_intra)))/(ngas + natoms_per_molecule); + + if (random_equal->uniform() < deltaphi) { ninsertion_successes += 1.0; energy_stored = energy_after; diff --git a/src/MC/fix_gcmc.h b/src/MC/fix_gcmc.h index 3b1e8c72b1..88d05e022c 100644 --- a/src/MC/fix_gcmc.h +++ b/src/MC/fix_gcmc.h @@ -112,6 +112,8 @@ class FixGCMC : public Fix { double **atom_coord; imageint imagetmp; + double energy_intra; + class Pair *pair; class RanPark *random_equal;