Description of the program Gibbs: Calculation of liquid phase equilibria with the cavity-biased Gibbs ensemble method.

Molec. Simulation, Vol 9, pp 257-267 (1992) Click HERE to download a copy

Mihaly Mezei

Department of Pharmacological Sciences,

Icahn School of Medicine at Mount Sinai

New York, NY 10029

Mihaly.Mezei@mssm.edu

Nov. 15, 2004.

I. General description

II. Installation

III. Input syntax

IV. Sample input

V. File formats

VI. References

I. General Description

The program calculates phase equilibria of a one-component molecular fluid by simulating the system in the Gibbs ensemble of Panagiotpoulus REF. In the Gibbs ensemble, two copies of the system are simulated in equilibrium with each other. The equilibrium is obtained by periodic attempt to transfer either some volume or a molecule from one system to an other. Such exchanges are accepted usiong rules similar to those used in (T,P,N) or (T,V,mu) ensemble simulations, respectively.

The transfer of molecules into a liquid is facilitated by the cavity-biasig of Mezei REF, originally introduced for simulations in the grand-canonical ensemble REF.

The program can also be used for the simulation of two (T,P,N) or (T,V,mu) ensemble fluids.

The molecules are assumed to be rigid. They can interact with a combination of Lennard-Jones and electrostatic terms. The long-range contributions (that are of greater significance here than for constant N simulations) can be estimated by the reaction field method for the electrostatic part and using the formula of Mezei & Bencsath for the Lennard-Jones part REF.

II. Installation

The program is written in standard Fortran-77. Before compilation, certain symbols have to be changed in the source file gibbs.for. These symbols determine the maximum dimensions of arrays used.

change '#IJ' to '<(maxmol-1)*maxmol/2>'
change '#MO' to '<maxmol>' (limits nmmg and nmmr)
change '#NA' to '<maxslv*maxmol>'
change '#RG' to '<maximumn number of rdf grids>' (limits nrgrid)
change '#SV' to '<maxslv>' (limits nslv)
change '#CV' to '<maximum number of cavities>'
change '#GX' to '<maximum number of (limits ngrid) cavity grids in each direction>'

Furthermore, certain platform-dependent lines may have to be uncommented as follows:

For IBM VM/MVS systems, change 'C@BM' to ''
For VAX/VMS systems, change 'C@VX' to ''
For Cray, change 'C@RY' to ''

There is a c-shell script pre.csh that performs the dimension-related substitution. It can be edited to accomodated the system to be simulated.

II. Input syntax

The endings G and R in the variable names refer to the grid or random insertion. Grid insertion is done in the liquid and random insertion in the vapor phase. For some quantities, zero input will result in using default values, shown below enclosed between braces.

Option array : iop (40i2)
- iop(1) = 1 : create inital configuration
- iop(1) = 2 : input from unit 13, but update number of molecules
- iop(2) = 0 : apply volume change to the R system
- iop(2) = 1 : apply volume change to the G system
- iop(3) = 0 : Gibbs ensemble simulation
- iop(3) = 1 : two (T,P,N) ensemble simulations
- iop(3) = 2 : two (T,P,mu) ensemble simulations
- iop(4) = 0 : binary input on unit 13
- iop(4) = 1 : alphanumeric input on unit 13
- iop(6) = 0 : no particle exchange
- iop(6) = 1 : use random insertion
- iop(6) = 2 : use cavity-biased insertion
- iop(8) = 0 : don't use the long-range correction
- iop(8) = 1 : use the long-range correction
- iop(9) = 0 : spherical cut-off for the electrostatic interactions
- iop(9) = 1 : reaction field correction with eps(rf)=infinity (applicable only if iop(11)=6)
- iop(10) > 0 : restart, input from unit 10
- iop(11) = 0 : MCY water
- iop(11) = 1 : ST2 water
- iop(11) = 2 : SPC/TIPS water
- iop(11) = 3 : TIPS2/TIPS4P water
- iop(11) = 4 : LJ fluid
- iop(11) = 5 : general 12-6 (LJ) molecule
- iop(11) = 6 : general 12-6-1 molecule
- iop(12) = 0 : com synchronization at 10⁵th MC step
- iop(12) > 0 : com synchronization at 10^(iop12-1)-th MC step
- iop(13) > 0 : echo inputted data
- iop(15) > 0 : save initial config on unit 13 and restart
- iop(15) = 2 : save initial config on unit 13 in alphanum and restart
- iop(16) > 0 : create two free-format coordinate file and stop
- iop(17) > 0 : calculate g(r) and angular correlations for the G system
- iop(17) > 1 : calculate g(r) and angular correlations for the R system
- iop(17) > 2 : calculate the absolute value of the correlations
If iop(10) > 0 (run continued from a checkpoint file), skip items 2-8.
id (20a4/20a4)
nrgrid, icolcr, rgrid, alj, blj (2i5,3f20.0)
- nrgrid {#RG}: number of rdf girdpoints (in each dimension)
- icolcr {1}: the colum of the orientation matrix to use for correlation calc
- rgrid {0.05}: radial gridsize in A
- alj, blj: LJ parameters for TIPS/SPC type water models.
nmmg, nmmr, temp, edgea, prsag, prsar (2i5,4f15.5)
- nmmg, nmmr: If nmmx<nmolx, the balance is deleted; if nmmx>nmolx, the balance of molecules is created randomly and the new configuration is saved on unit 13 and the calculation stops. Here nmolx is the number of molecules in the G and R systems as read from unit 13.
- temp: the temperature in Kelvin
- edgea: the cell size in A
- prsag, prsar: When iop(3) > 0 (i.e., two (T,P,N) simulations are run), the pressures in the G and R systems, respectively, in A.
cedgs, rotax, delvol (20x,3f20.0)
- cedgs: displacements along each coordinate axis will be limited to +/-cedgs/2
- rotax: rotations will be limited to +/-rotax/2
- delvol: volume changes will be limited to +/-delvol
nslv, ngrid, cutof, rsph (2i5,3f20.0)
- nslv: the number of solvent atoms
- ngrid {#GX}: number of cb grids
- cutof: potential cutoff (in A)
- rsph: the cavity diameter
nslv records (one for each atom in the simulated molecule):
- If iop(11) < 5 then islv, x, y, z, q (i5,3f15.10,3f10.0)
  - islv: atomic number
  - x, y, z: coordinates of he atom in a local frame
  - q: partial charge
- If iop(11) > 5 then islv, x, y, z, q, sigma, eps (i5,3f15.10,3f10.0)
  - sigma, eps: the LJ sigma (in A) and epsilon (in kcal/mol)
nmcmax, nmcrep, nrecd, nmcvch, nmccnt (7i10)
- nmcmax: number of MC steps to run
- nmcrep: one-line report frequency
- nrecd: check-poinf file saving frequency
- nmcvch: volume change attempt frequency
- nmccnt: control function saving frequency
When iop(10)=0, the input terminates here.
nmcmax, nmcrep (7i10)

IV. Sample input

The input for methylene fluoride REF at T= 225 K is given below, first properly formatted and next with line-by-line annotations:

 2 0 0 1 0 2 0 1 1 0 6 0 0 0 0 0 0
 GIBBS   methylene fluoride at 225K
 10 million MC steps
    0    0         0.0                0.0                0.0
  511    2         225.0        32.50000
             0.3                15.0                400.
    5  100   15.                 3.0
    6   0.00000         0.00000      0.00000         0.300     3.150   0.10844
    1   0.91213         0.00000     -0.59676         0.075     2.170   0.01985
    1  -0.91213         0.00000     -0.59676         0.075     2.170   0.01985
    9   0.00000         1.10464      0.79333        -0.225     2.975   0.07944
    9   0.00000        -1.10464      0.79333        -0.225     2.975   0.07944
  10000000     10000    100000       500    50000

Option array: iop
2 0 0 1 0 2 0 1 1 0 6 0 0 0 0 0 0
id
GIBBS methylene fluoride at 225K
10 million MC steps
nrgrid, icolcr, rgrid, alj, blj
0 0 0.0 0.0 0.0
nmmg, nmmr, temp, edgea, prsag, prsar
511 2 225.0 32.50000
cedgs, rotax, delvol
0.3 15.0 400.
nslv, ngrid, cutof, rsph
5 100 15. 3.0
nslv records (iop(11) > 5): islv, x, y, z, q, sigma, eps
6 0.00000 0.00000 0.00000 0.300 3.150 0.10844
1 0.91213 0.00000 -0.59676 0.075 2.170 0.01985
1 -0.91213 0.00000 -0.59676 0.075 2.170 0.01985
9 0.00000 1.10464 0.79333 -0.225 2.975 0.07944
9 0.00000 -1.10464 0.79333 -0.225 2.975 0.07944
nmcmax, nmcrep, nrecd, nmcvch, nmccnt
10000000 10000 100000 500 50000

V. File formats

Unit 10: checkpoint file
Unit 13: configuration file
The configuration file contains the description of both simulated systems. Depending on iop(4) it is either an ASCII file or a Fortran binary file.
- When iop(4)=0, the file is in binary format, written by the following Fortran statements:
```
   if (iop(4) .eq. 0) then
     write (13) nmolg,nmolr,ix,icopvg,volg,volr,
-      ((crmg(k,im),k=1,3),im=1,nmolg),
-      ((crmr(k,im),k=1,3),im=1,nmolr)

     write (13)
-      ((cg(k,ia,1),k=1,3),ia=1,natg,1),
-      ((cr(k,ia,1),k=1,3),ia=1,natr)
```
- When iop(4)=1, the file is ASCII format, containing the following records:
  - nmolg,nmolr,ix,icopvg,volg,volr (2i6,i12,i2,2e15.7)
  - ((crmg(k,im),k=1,3),im=1,nmolg) (3f20.10)
  - ((crmr(k,im),k=1,3),im=1,nmolr) (3f20.10)
  - ((cg(k,ia,icopvg),k=1,3),ia=1,natg) (3f20.10)
  - ((cr(k,ia,icopvr),k=1,3),ia=1,natr) (3f20.10)
Here
- nmolg,nmolr are the number of molecule in the G and R systems,
- ix is the current random number seed
- volg,volr are the volumes of the G and R systems
- crmg,crmr are the center of masses of the molecules in the G and R systems
- cg,cr are the atomic coordinates of the molecules in the G and R systems

VI. References

Original description of the Gibbs ensemble method: A.Z. Panagiotopoulos, Molec. Phys., 61, p 813 (1987)
Original description of the cavity-biased Gibbs ensemble method: M. Mezei, Molec. Simulation, 9, pp 257-267 (1992)
Click HERE
Cavity-biased insertion - grid-based formalism: M. Mezei, Mol. Phys., 61, 565-582 (1987); Erratum, 67,1207-1208 (1989).
Long-range correction for L-J fluids using a center of mass cutoff: M. Mezei, and K. Bencsath, J. Math. Chem., 16, 61-71 (1994).
Phase equilibria of methylene fluoride: P. Jedlovszky, and M. Mezei, J. Chem., Phys., in print (1999).

I.	General description
II.	Installation
III.	Input syntax
IV.	Sample input
V.	File formats
VI.	References