Introduction
Described below are methods to grow molecules while using the gcmc or Gibbs program. These methods tell the computer how to form molecules with many interaction sites. Each interaction site is referred to as a bead. Each method takes as input integer parameters and real parameters. The parameters required for each method are conveyed in the following format:
Method ( I1, I2, ... , In; R1, R2, ... , Rn )
The method is listed followed by the required integer parameters, denoted
Im,
and the required real parameters, denoted Rm..
The units that the real parameters need to be inputted as follows along
with any restrictions. The method and parameters are then described in
the text that follows. A //
symbol indicates a line change in the input file. Please note that the
method names must be spelled exactly as seen here, including letter case,
in the input file.
Random
Random
Random
chooses a X, Y, and Z coordinate for the bead between zero and the box
length. This method is used to place the first bead for all molecules.
Sphere
Sphere ( I1; R1 )
R1 [=] Å
Sphere randomly chooses a X, Y, and Z coordinate on the surface of a sphere with radius R1 and a center located at the position of bead I1.
Example: Non-polar diatomic molecules are often modeled with two Lennard-Jones interaction sites separated by a distance L. The growth of this molecule would consist of two steps. First, one of the beads would randomly be placed in the simulation box using Random. Next, Sphere would be used to select a position for the second bead.
1. Random
2. Sphere ( 1; L )
Cone
Cone ( I1, I2, I3; R1_1 .. R1_I3, R2_1 .. R2_I3, R3_1 .. R3_I3 )
R1 [=] Degrees
R2 [=] Degrees
R3 [=] Å
R2_1 = 0.0
Cone randomly places I3 beads on a circle. Bead i is placed a distance R3_i from bead I2 at an angle R1_i formed by beads I1, I2, and i. Beads 2 to I3 are rotated R2_j(1>j>=I3) degrees from the first bead around an axis formed by beads I1 and I2. Figure 1 and Figure 2 display a case in which two beads, B1 and B2, are grown from the existing beads I1 and I2. In Figure 1 the I1-I2 axis exists in the plane of the paper. Figure 2 displays a view looking down the I1-I2 axis.
Figure 1: A plane formed by I1, I2, and B1.
Figure 2: Looking down the I1-I2
axis.
Match
Match ( I1 )
Match gives the bead being grown the same coordinates as bead I1. This method is useful when a site interacts with a Lennard-Jones and coulombic potential.
Example: The EMP2 model for carbon dioxide contains a carbon center and two oxygen centers. Each of the centers interact via a Lennard-Jones and coulombic potential. The C-O distance is fixed at 1.149 Å and the O-C-O angle is rigid. Let’s designate beads one to three as the Lennard-Jones beads for carbon, oxygen, and oxygen respectively. Similarly, the coulombic beads for carbon, oxygen, and oxygen will be numbered four through six respectively. Growth of the molecule will start by randomly placing the Lennard-Jones bead for carbon using Random. One of the Lennard-Jones oxygen beads will then be placed using Sphere. The second oxygen bead will then be placed using Cone. The molecule will be completed by matching the coulombic beads to the existing Lennard-Jones beads.
1. Random
2. Sphere ( 1; 1.149 )
3. Cone ( 2, 1, 1; -180.0, 0.0, 1.149 )
4. Match ( 1 )
5. Match ( 2 )
6. Match ( 3 )
Complete
Complete
Complete is used in combination with Cone. In an input file a method has to be entered for each bead. Using the Cone method more than one bead can be grown at a time. When the Cone method is used to grow more than one bead Complete should be listed as the method for the beads that have already been grown.
Example: Imagine a potential for ammonia that consists of four centers, one for the nitrogen core and three for the hydrogen centers. The nitrogen center will contain both a Lennard-Jones and Coulombic bead, however the hydrogen centers will contain only a coulombic bead. For simplicity the N-H distance and the H-N-H angle are fixed to 1.0124 Å and 106.7o respectively. Growth of this molecule consists of five steps. First the Lennard-Jones nitrogen bead is placed randomly in the simulation box using Random. The coulombic nitrogen bead is then matched to the existing nitrogen bead using Match. Next, one of the hydrogen beads is placed on the surface of a sphere using Sphere. Then, the last two hydrogens are placed using Cone in combination with Complete.
1. Random
2. Match ( 1 )
3. Sphere ( 1; 1.0124 )
4. Cone ( 3, 1, 2; 106.7, 106.7, 0.0, 113.72, 1.0124, 1.0124 )
5. Complete
FxBend
FxBend ( I1, I2 // I3_1 .. I3_I1; R1, R2_1, R3_1 .. R2_I1, R3_I1 )
R1 [=] Å
R2 [=] Kelvin rad-2
R3 [=] degrees
This method places a bead, B1, a distance R1 from bead I2. A position is selected by satisfying I1 bending potentials. The bending angles, qi, are defined by beads I3_i,I2, and B1. The bending energy for each angle is calculated from:
(1)
where kq is equal to R2_i and q0 is given by R3_i.
This is accomplished using an iterative process. First a X, Y, Z position is randomly chosen on the surface of a sphere with radius R1 centered at bead I2. The angles I3_i-I2-B1 are then calculated along with ubend. The position is then accepted with probability:
(2)
If rejected, the process is repeated until a position has been accepted.
NOTE: Using this method to grow a bead that must satisfy multiple
bending angles is not strictly correct. In this case alternative
methods should be employed.
FxBendTor
This method places a bead, B1, a distance R1 from bead I4.A position is selected by satisfying I1 bending potentials and I2 torsion potentials. The bending angles, qi, are defined by beads I5_i, I4, and B1. The torsion angles, fi, are defined by beads I6_i, I7_i, I4, and B1. The bending potential used is described in the FxBend method section. Two torsion potentials exist, the OPLS potential and the de Pablo potential. The OPLS model is selected by setting I3 equal to 1. For the de Pablo potential I3 is set to 2.
OPLS Torsion Potential
FxBendTor ( I1, I2, I3, I4 // I5_1 .. I5_I1, I6_1, I7_1 .. I6_I2, I7_I2; R1, R2_1, R3_1 .. R2_I1, R3_I1, R4_1, R5_1, R6_1, R7_1, R8_1, R9_1 .. R4_I2, R5_I2, R6_I2, R7_I2, R8_I2, R9_I2 )
I3 = 1
R1 [=] Å
R2 [=] Kelvin rad-2
R3 [=] degrees
R4 to R9 [=] Kelvin
The OPLS potential is:
(3)
De Pablo Torsion Potential
FxBendTor ( I1, I2, I3, I4 // I5_1 .. I5_I1, I6_1 .. I6_I2; R1, R2_1, R3_1 .. R2_I1, R3_I1, R4_1, R5_1, R6_1, R7_1 .. R4_I2, R5_I2, R6_I2, R7_I2 )
I3 = 2
R1 [=] Å
R2 [=] Kelvin rad-2
R3 [=] degrees
R4 to R7 [=] Kelvin
The de Pablo potential is:
(4)
The values of kq and q0 in equation (1) are R2_i and R3_i respectively. The ci values in equations (3) and (4) are given by R4 to R9 and R4 to R7 respectively.
Analogous to FxBend a X, Y, Z position is found using an iterative process. First, Sphere is used to randomly generate a X, Y, Z site a distance R3 from bead I4. The bending angles and the torsion angles are then calculated along with ubend and utor. The position is then accepted with probability:
(5)
If rejected, the process is repeated until a position has been accepted.
NOTE: Using this method to grow a bead that must satisfy multiple bending and/or torsion angles is not strictly correct. In this case alternative methods should be employed.
Example: The S-K-S model for n-alkanes consists of a series of beads representing the CH3 and CH2 groups. The bond lengths are kept fixed at 1.53 Å. The bending angles are flexible with a modulus of 62,500 Kelvin rad-2 and an equilibrium angle of 114o. The de Pablo potential is used with parameters c1 = 355 K, c2 = -68.19 K, and c3 = 791.3 K. To grow hexane a CH3 bead would be randomly placed within the simulation box using Random. Then a CH2 bead would randomly be placed on the surface of a sphere using Sphere. Next, FxBend would find a suitable position for the next CH2 bead. Finally, the remaining three beads would be placed using FxBendTor.
1. Random
2. Sphere ( 1; 1.53 )
3. FxBend ( 1, 2 // 1; 1.53, 62500.0, 114.0 )
4. FxBendTor ( 1, 1, 2, 3 // 2, 1, 2; 1.53, 62500.0, 114.0, 0.0, 355.0, -68.19, 791.3 )
5. FxBendTor ( 1, 1, 2, 4 // 3, 2, 3; 1.53, 62500.0, 114.0, 0.0, 355.0, -68.19, 791.3 )
6. FxBendTor ( 1,
1, 2, 5 // 4, 3, 4; 1.53, 62500.0, 114.0, 0.0, 355.0, -68.19, 791.3 )
Stretch
Stretch ( I1; R1, R2 )
R1 [=] Kelvin Å -2
R2 [=] Å
Stretch randomly chooses a X, Y, and Z coordinate on the surface of a sphere with a radius selected from the stretching potential:
(6)
where kr
is equal to R1
and rO
is equal to R2.
StBend
StBend ( I1, I2 // I3_1 .. I3_I1; R1, R2, R3_1, R4_1 .. R3_I1, R4_I1 )
R1 [=] Kelvin Å -2
R2 [=] Å
R3 [=] Kelvin rad-2
R4 [=] degrees
The method is analogous to FxBend. The difference is that the bond length is selected according to a stretching potential instead of having a fixed length. As a result one additional parameter is needed. The stretching modulus is added as the first real parameter. The rest of the parameters follow the same as in FxBend.
NOTE: Using this method to grow a bead that must satisfy multiple
bending angles is not strictly correct. In this case alternative
methods should be employed.
StBendTor
OPLS Torsion Potential
StBendTor ( I1, I2, I3, I4 // I5_1 .. I5_I1, I6_1, I7_1 .. I6_I2, I7_I2; R1, R2, R3_1, R4_1 .. R3_I1, R4_I1, R5_1, R6_1, R7_1, R8_1, R9_1, R10_1 .. R5_I2, R6_I2, R7_I2, R8_I2, R9_I2, R10_I2 )
I3 = 1
R1 [=] Kelvin Å -2
R2 [=] Å
R3 [=] Kelvin rad-2
R4 [=] degrees
R5 to R10 [=] Kelvin
De Pablo Torsion Potential
StBendTor ( I1, I2, I3, I4 // I5_1 .. I5_I1, I6_1, I7_1 .. I6_I2, I7_I2; R1, R2, R3_1, R4_1 .. R3_I1, R4_I1, R5_1, R6_1, R7_1, R8_1 .. R5_I2, R6_I2, R7_I2, R8_I2 )
I3 = 2
R1 [=] Kelvin Å -2
R2 [=] Å
R3 [=] Kelvin rad-2
R4 [=] degrees
R5 to R8 [=] Kelvin
The method is analogous to FxBendTor. The difference is that the bond length is selected according to a stretching potential instead of having a fixed length. As a result one additional parameter is needed. The stretching modulus is added as the first real parameter. The rest of the parameters follow the same as in FxBendTor.
NOTE: Using this method to grow a bead that must satisfy multiple
bending and/or torsion angles is not strictly correct. In this case
alternative methods should be employed.