From d2f90c7e0c6dc48cc3293b8298f690f1180d0fc3 Mon Sep 17 00:00:00 2001
From: athomps <athomps@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Mon, 11 May 2015 23:16:44 +0000
Subject: [PATCH] Added intra_energy keyword

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@13459 f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
 src/MC/fix_gcmc.cpp | 20 ++++++++++++++++----
 src/MC/fix_gcmc.h   |  2 ++
 2 files changed, 18 insertions(+), 4 deletions(-)

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;