From 87023e98bb27d7b5ae8a8a64aba763214c1d1187 Mon Sep 17 00:00:00 2001
From: pscrozi <pscrozi@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Wed, 27 Mar 2013 22:35:13 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9717
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/MC/fix_gcmc.cpp | 13 +++++++++++++
 src/MC/fix_gcmc.h   |  2 ++
 2 files changed, 15 insertions(+)

diff --git a/src/MC/fix_gcmc.cpp b/src/MC/fix_gcmc.cpp
index 0eb8b5fe9e..7f4fa776d4 100644
--- a/src/MC/fix_gcmc.cpp
+++ b/src/MC/fix_gcmc.cpp
@@ -85,6 +85,9 @@ FixGCMC::FixGCMC(LAMMPS *lmp, int narg, char **arg) :
   rotation_group = 0;
   rotation_groupbit = 0;
   rotation_inversegroupbit = 0;
+  pressure_flag = false;
+  pressure = 0.0;
+  fugacity_coeff = 1.0;
 
   // read options from end of input line
 
@@ -199,6 +202,15 @@ void FixGCMC::options(int narg, char **arg)
       max_rotation_angle = atof(arg[iarg+1]);
       max_rotation_angle *= MY_PI/180;
       iarg += 2;
+    } else if (strcmp(arg[iarg],"pressure") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal fix gcmc command");
+      pressure = atof(arg[iarg+1]);
+      pressure_flag = true;
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"fugacity_coeff") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal fix gcmc command");
+      fugacity_coeff = atof(arg[iarg+1]);
+      iarg += 2;
     } else error->all(FLERR,"Illegal fix gcmc command");
   }
 }
@@ -335,6 +347,7 @@ void FixGCMC::init()
                         force->boltz*reservoir_temperature));
   sigma = sqrt(force->boltz*reservoir_temperature/gas_mass/force->mvv2e);
   zz = exp(beta*chemical_potential)/(pow(lambda,3.0));
+  if (pressure_flag) zz = pressure*fugacity_coeff*beta/force->nktv2p;
   
   imagetmp = ((tagint) IMGMAX << IMG2BITS) | 
              ((tagint) IMGMAX << IMGBITS) | IMGMAX;
diff --git a/src/MC/fix_gcmc.h b/src/MC/fix_gcmc.h
index dc7254f392..551e36e8d2 100644
--- a/src/MC/fix_gcmc.h
+++ b/src/MC/fix_gcmc.h
@@ -63,6 +63,7 @@ class FixGCMC : public Fix {
   int regionflag;           // 0 = anywhere in box, 1 = specific region
   int iregion;              // GCMC region
   char *idregion;           // GCMC region id
+  bool pressure_flag;       // true if user specified reservoir pressure, false otherwise
 
   int maxmol;               // largest molecule tag across all existing atoms
   int natoms_per_molecule;  // number of atoms in each gas molecule
@@ -84,6 +85,7 @@ class FixGCMC : public Fix {
   double displace;
   double max_rotation_angle;
   double beta,zz,sigma,volume;
+  double pressure,fugacity_coeff;
   double xlo,xhi,ylo,yhi,zlo,zhi;
   double region_xlo,region_xhi,region_ylo,region_yhi,region_zlo,region_zhi;
   double region_volume;