Input File Formats

In order to run simulation in GOMC, the following files need to be provided:

  • GOMC executable

  • PDB file(s)

  • PSF file(s)

  • Parameter file

  • Input file “NAME.conf” (proprietary control file)

In order to restart a simulation in GOMC from exactly where it left off, the following files also need to be provided:

  • XSC file(s)

  • COOR file(s)

  • CHK file

In order to run a hybrid MCMD simulation, the following files also need to be provided:

  • XSC file(s)

  • COOR file(s)

  • VEL files(s)

  • CHK file

PDB File

GOMC requires only one PDB file for NVT and NPT ensembles. However, GOMC requires two PDB files for GEMC and GCMC ensembles.

What is PDB file

The term PDB can refer to the Protein Data Bank (http://www.rcsb.org/pdb/), to a data file provided there, or to any file following the PDB format. Files in the PDB include various information such as the name of the compound, the ATOM and HETATM records containing the coordinates of the molecules, and etc. PDB widely used by NAMD, GROMACS, CHARMM, ACEMD, and Amber. GOMC ignore everything in a PDB file except for the REMARK, CRYST1, ATOM, and END records. An overview of the PDB standard can be found here:

http://www.wwpdb.org/documentation/file-format-content/format33/sect2.html#HEADER

http://www.wwpdb.org/documentation/file-format-content/format33/sect8.html#CRYST1

http://www.wwpdb.org/documentation/file-format-content/format33/sect9.html#ATOM

PDB contains four major parts; REMARK, CRYST1, ATOM, and END. Here is the definition of each field and how GOMC is using them to get the information it requires.

  • REMARK: This header records present experimental details, annotations, comments, and information not included in other records (for more information, click here).

    However, GOMC uses this header to print simulation informations.

    • Max Displacement (Å)

    • Max Rotation (Degree)

    • Max volume exchange (\(\AA^3\))

    • Monte Carlo Steps (MC)

  • CRYST1: This header records the unit cell dimension parameters.

    • Lattice constant: a,b,c (Å)

    • Lattice angles: \(\alpha, \beta, \gamma\) (Degree)

  • ATOM: The ATOM records present the atomic coordinates for standard amino acids and nucleotides. They also present the occupancy and temperature factor for each atom.

    • ATOM: Record name

    • serial: Atom serial number.

    • name: Atom name.

    • resName: Residue name.

    • chainID: Chain identifier.

    • resSeq: Residue sequence number.

    • x: Coordinates for X (Å).

    • y: Coordinates for Y (Å).

    • z: Coordinates for Z (Å).

    • occupancy: GOMC uses to define which atoms belong to which box.

    • beta: Beta or Temperature factor. GOMC uses this value to define the mobility of the atoms. element: Element symbol.

  • END: A frame in the PDB file is terminated with the keyword.

Here are the PDB output of GOMC for the first molecule of isobutane:

REMARK    GOMC   122.790    3.14159     3439.817     1000000
CRYST1  35.245    35.245    35.245    90.00   90.00   90.00
ATOM    1     C1    ISB     1     0.911    -0.313    0.000    0.00    0.00    C
ATOM    2     C1    ISB     1     1.424    -1.765    0.000    0.00    0.00    C
ATOM    3     C1    ISB     1    -0.629    -0.313    0.000    0.00    0.00    C
ATOM    4     C1    ISB     1     1.424     0.413   -1.257    0.00    0.00    C
END

The fields seen here in order from left to right are the record type, atom ID, atom name, residue name, residue ID, x, y, and z coordinates, occupancy, temperature factor (called beta), and segment name.

The atom name is “C1” and residue name is “ISB”. The PSF file (next section) contains a lookup table of atoms. These contain the atom name from the PDB and the name of the atom kind in the parameter file it corresponds to. As multiple different atom names will all correspond to the same parameter, these can be viewed “atom aliases” of sorts. The chain letter (in this case ‘A’) is sometimes used when packing a number of PDBs into a single PDB file.

Important

  • VMD requires a constant number of ATOMs in a multi-frame PDB (multiple records terminated by “END” in a single file). To compensate for this, all atoms from all boxes in the system are written to the output PDBs of this code.

  • For atoms not currently in a box, the coordinates are set to < 0.00, 0.00, 0.00 >. The occupancy is commonly just set to “1.00” and is left unused by many codes. We recycle this legacy parameter by using it to denote, in our output PDBs, the box a molecule is in (box 0 occupancy=0.00 ; box 1 occupancy=1.00)

  • The beta value in GOMC code is used to define the mobility of the molecule.

    • Beta = 0.00: molecule can move and transfer within and between boxes.

    • Beta = 1.00: molecule is fixed in its position.

    • Beta = 2.00: molecule can move within the box but cannot be transferred between boxes.

Generating PDB file

With that overview of the format in mind, the following steps describe how a PDB file is typically built.

  1. A single molecule PDB is obtained. In this example, the GaussView was used to draw the molecule, which was then edited by hand to adhere to the PDB spec properly. There are many open-source software that can build a molecule for you, such as Avagadro , molefacture in VMD and more. The end result is a PDB for a single molecule:

REMARK   1 File created by GaussView 5.0.8
ATOM     1  C1   ISB  1   0.911  -0.313    0.000  C
ATOM     2  C1   ISB  1   1.424  -1.765    0.000  C
ATOM     3  C1   ISB  1  -0.629  -0.313    0.000  C
ATOM     4  C1   ISB  1   1.424   0.413   -1.257  C
END
  1. Next, packings are calculated to place the simulation in a region of vapor-liquid coexistence. There are a couple of ways to do this in Gibbs ensemble:

  • Pack both boxes to a single middle density, which is an average of the liquid and vapor densities.

  • Same as previous method, but add a modest amount to axis of one box (e.g. 10-30 A). This technique can be handy in the constant pressure Gibbs ensemble.

  • Pack one box to the predicted liquid density and the other to the vapor density.

    A good reference for getting the information needed to estimate packing is the NIST Web Book database of pure compounds:

    http://webbook.nist.gov/chemistry/

  1. After packing is determined, a basic pack can be performed with a Packmol script. Here is the example of packing 1000 isobutane in 70 A cubic box:

tolerance   3.0
filetype    pdb
output      STEP2_ISB_packed_BOX 0.pdb
structure   isobutane.pdb
number      1000
inside cube 0.1   0.1   0.1   70.20
end     structure

Copy the above text into “pack_isobutane.inp” file, save it and run the script by typing the following line into the terminal:

$ ./packmol < pack_isobutane.inp

PSF File

GOMC requires only one PSF file for NVT and NPT ensembles. However, GOMC requires two PSF files for GEMC and GCMC ensembles.

What is PSF file

Protein structure file (PSF), contains all of the molecule-specific information needed to apply a particular force field to a molecular system. The CHARMM force field is divided into a topology file, which is needed to generate the PSF file, and a parameter file, which supplies specific numerical values for the generic CHARMM potential function. The topology file defines the atom types used in the force field; the atom names, types, bonds, and partial charges of each residue type; and any patches necessary to link or otherwise mutate these basic residues. The parameter file provides a mapping between bonded and nonbonded interactions involving the various combinations of atom types found in the topology file and specific spring constants and similar parameters for all of the bond, angle, dihedral, improper, and van der Waals terms in the CHARMM potential function. PSF file widely used by by NAMD, CHARMM, and X-PLOR.

The PSF file contains six main sections: remarks, atoms, bonds, angles, dihedrals, and impropers (dihedral force terms used to maintain planarity). Each section starts with a specific header described bellow:

  • NTITLE: remarks on the file. The following is taken from a PSF file for isobutane:

    PSF
          3  !NTITLE
    REMARKS  original generated structure x-plor psf file
    REMARKS  topology ./Top_Branched_Alkanes.inp
    REMARKS  segment ISB { first NONE; last NONE; auto angles dihedrals }
    
  • NATOM: Defines the atom names, types, and partial charges of each residue type.

    atom    ID
    segment name
    residue ID
    residue name
    atom    name
    atom    type
    atom    charge
    atom    mass
    

    The following is taken from a PSF file for isobutane:

    4000 !NATOM
    1    ISB  1  ISB    C1    CH1    0.000000   13.0190  0
    2    ISB  1  ISB    C2    CH3    0.000000   15.0350  0
    3    ISB  1  ISB    C3    CH3    0.000000   15.0350  0
    4    ISB  1  ISB    C4    CH3    0.000000   15.0350  0
    5    ISB  2  ISB    C1    CH1    0.000000   13.0190  0
    6    ISB  2  ISB    C2    CH3    0.000000   15.0350  0
    7    ISB  2  ISB    C3    CH3    0.000000   15.0350  0
    8    ISB  2  ISB    C4    CH3    0.000000   15.0350  0
    

    The fields in the atom section, from left to right are atom ID, segment name, residue ID, residue name, atom name, atom type, charge, mass, and an unused 0.

  • NBOND: The covalent bond section lists four pairs of atoms per line. The following is taken from a PSF file for isobutane:

    3000   !BOND:     bonds
    1   2   1   3   1   4   5   6
    5   7   5   8
    
  • NTHETA: The angle section lists three triples of atoms per line. The following is taken from a PSF file for isobutane:

    3000   !NTHETA:   angles
    2   1   4   2   1   3   3   1   4
    6   5   8   6   5   7   7   5   8
    
  • NPHI: The dihedral sections list two quadruples of atoms per line.

  • NIMPHI: The improper sections list two quadruples of atoms per line. GOMC currently does not support improper. For the molecules without dihedral or improper, PDF file look like the following:

    0   !NPHI: dihedrals
    0   !NIMPHI: impropers
    
  • (other sections such as cross terms)

Important

  • The PSF file format is a highly redundant file format. It repeats identical topology of thousands of molecules of a common kind in some cases. GOMC follows the same approach as NAMD, allowing this excess information externally and compiling it in the code.

  • Other sections (e.g. cross terms) contain unsupported or legacy parameters and are ignored.

  • Following the restriction of VMD, the order of the atoms in PSF file must match the order of the atoms in the PDB file.

  • Improper entries are read and stored, but are not currently used. Support will eventually be added for this.

Generating PSF file

The PSF file is typically generated using PSFGen. It is convenient to make a script, such as the example below, to do this:

package require psfgen
topology  ./Top_branched_Alaknes.inp
segment ISB {
  pdb   ./STEP2_ISB_packed_BOX 0.pdb
  first   none
  last  none
}

coordpdb ./STEP2_ISB_packed_BOX 0.pdb ISB

writepsf ./STEP3_START_ISB_sys_BOX_0.psf
writepdb ./STEP3_START_ISB_sys_BOX_0.pdb

Typically, one script is run per box to generate a finalized PDB/PSF for that box. The script requires one additional file, the NAMD-style topology file. While GOMC does not directly read or interact with this file, it’s typically used to generate the PSF and, hence, is considered one of the integral file types. It will be briefly discussed in the following section.

Topology File

A CHARMM forcefield topology file contains all of the information needed to convert a list of residue names into a complete PSF structure file. The topology is a whitespace separated file format, which contains a list of atoms and their corresponding masses, and a list of residue information (charges, composition, and topology). Essentially, it is a non-redundant lookup table equivalent to the PSF file.

This is followed by a series of residues, which tell PSFGen what atoms are bonded to a given atom. Each residue is comprised of four key elements:

  • A header beginning with the keyword RESI with the residue name and net charge

  • A body with multiple ATOM entries (not to be confused with the PDB-style entries of the same name), which list the partial charge on the particle and what kind of atom each named atom in a specific molecule/residue is.

  • A section of lines starting with the word BOND contains pairs of bonded atoms (typically 3 per line)

  • A closing section with instructions for PSFGen.

Here’s an example of topology file for isobutane:

* Custom top file -- branched alkanes *
11
!
MASS 1 CH3 15.035 C !
MASS 2 CH1 13.019 C !

AUTOGENERATE ANGLES DIHEDRALS

RESI ISB    0.00 !  isobutane - TraPPE
GROUP
ATOM  C1  CH1   0.00 !  C3\
ATOM  C2  CH3   0.00 !     C1-C2
ATOM  C3  CH3   0.00 !  C4/
ATOM  C4  CH3   0.00 !
BOND  C1  C2  C1  C3  C1  C4
PATCHING FIRS NONE LAST NONE

END

Note

The keyword END must be used to terminate this file and keywords related to the auto-generation process must be placed near the top of the file, after the MASS definitions.

Parameter File

Currently, GOMC uses a single parameter file and the user has the two kinds of parameter file choices:

  • CHARMM (Chemistry at Harvard Molecular Mechanics) compatible parameter file

  • EXOTIC or Mie parameter file

If the parameter file type is not specified or if the chosen file is missing, an error will result.

Both force field file options are whitespace separated files with sections preceded by a tag. When a known tag (representing a molecular interaction in the model) is encountered, reading of that section of the force field begins. Comments (anything after a * or !) and whitespace are ignored. Reading concludes when the end of the file is reached or another section tag is encountered.

CHARMM format parameter file

CHARMM contains a widely used model for describing energies in Monte Carlo and molecular dynamics simulations. It is intended to be compatible with other codes that use such a format, such as NAMD. See here for a general overview of the CHARMM force field.

Here’s the basic CHARMM contributions that are supported in GOMC:

\[\begin{split}U_{\texttt{bond}}&=\sum_{\texttt{bonds}} K_b(b-b_0)^2\\ U_{\texttt{angle}}&=\sum_{\texttt{angles}} K_{\theta}(\theta-\theta_0)^2\\ U_{\texttt{dihedral}}&=\sum_{\texttt{dihedrals}} K_{\phi} [1+\cos(n\phi - \delta)]\\ U_{\texttt{LJ}}&=\sum_{\texttt{nonbonded}} \epsilon_{ij}\left[\left(\frac{R_{min_{ij}}}{r_{ij}}\right)^{12}-2\left(\frac{R_{min_{ij}}}{r_{ij}}\right)^6\right]+ \frac{q_i q_j}{\epsilon r_{ij}}\end{split}\]

As seen above, the following are recognized, read and used:

  • BONDS - Quadratic expression describing bond stretching based on bond length (b) in Angstrom – Typically, it is ignored as bonds are rigid for Monte Carlo simulations.

    Note

    GOMC does not sample bond stretch. To ignore the relative bond energy, set the \(K_b\) to a large value i.e. “999999999999”.

    images/bonds.png

    Oscillations about the equilibrium bond length

  • ANGLES - Describe the conformational behavior of an angle (\(\delta\)) between three atoms, one of which is shared branch point to the other two.

    Note

    To fix any angle and ignore the related angle energy, set the \(K_\theta\) to a large value i.e. “999999999999”.

    images/angle.png

    Oscillations of 3 atoms about an equilibrium bond angle

  • DIHEDRALS - Describes crankshaft-like rotation behavior about a central bond in a series of three consecutive bonds (rotation is given as \(\phi\)).

    images/dihedrals.png

    Torsional rotation of 4 atoms about a central bond

  • NONBONDED - This tag name only should be used if CHARMM force files are being used. This section describes 12-6 (Lennard-Jones) non-bonded interactions. Non-bonded parameters are assigned by specifying atom type name followed by polarizabilities (which will be ignored), minimum energy, and (minimum radius)/2. In order to modify 1-4 interaction, a second polarizability (again, will be ignored), minimum energy, and (minimum radius)/2 need to be defined; otherwise, the same parameter will be considered for 1-4 interaction.

    images/nonbonded.png

    Non-bonded energy terms (electrostatics and Lennard-Jones)

  • NBFIX - This tag name only should be used if CHARMM force field is being used. This section allows in- teraction between two pairs of atoms to be modified, done by specifying two atom type names followed by minimum energy and minimum radius. In order to modify 1-4 interaction, a second minimum energy and minimum radius need to be defined; otherwise, the same parameter will be considered for 1-4 interaction.

    Note

    Please pay attention that in this section we define minimum radius, not (minimum radius)/2 as it is defined in the NONBONDED section.

    Note

    This does not modify the 1-4 electrostatic interactions.

    Currently, supported sections of the CHARMM compliant file include BONDS, ANGLES, DIHEDRALS, NONBONDED, NBFIX. Other sections such as CMAP are not currently read or supported.

BONDS

(“bond stretching”) is one key section of the CHARMM-compliant file. Units for the \(K_b\) variable in this section are in kcal/mol; the \(b_0\) section (which represents the equilibrium bond length for that kind of pair) is measured in Angstroms.

\[\begin{split}U_{\texttt{bond}}&=\sum_{\texttt{bonds}} K_b(b-b_0)^2\\\end{split}\]
BONDS
!V(bond) = Kb(b - b0)**2
!
!Kb:  kcal/mole/A**2
!b0:  A
!
!Kb (kcal/mol) = Kb (K) * Boltz.  const.;
!
!atom type      Kb          b0     description
CH3   CH1   9999999999    1.540 !  TraPPE 2

Note

The \(K_b\) value may appear odd, but this is because a larger value corresponds to a more rigid bond. As Monte Carlo force fields (e.g. TraPPE) typically treat molecules as rigid constructs, \(K_b\) is set to a large value - 9999999999. Sampling bond stretch is not supported in GOMC.

ANGLES

(“bond bending”), where \(\theta\) is the measured bond angle and \(\theta_0\) is the equilibrium bond angle for that kind of pair, are commonly measured in degrees and \(K_\theta\) is the force constant measured in kcal/mol/K. These values, in literature, are often expressed in Kelvin (K).

To convert Kelvin to kcal/mol/K, multiply by the Boltzmann constant – \(K_\theta\), 0.0019872041 kcal/mol. In order to fix the angle, it requires to set a large value for \(K_\theta\). By assigning a large value like 9999999999, specified angle will be fixed and energy of that angle will considered to be zero.

\[\begin{split}U_{\texttt{angle}}&=\sum_{\texttt{angles}} K_{\theta}(\theta-\theta_0)^2\\\end{split}\]

Here is an example of what is necessary for isobutane:

ANGLES
!
!V(angle) = Ktheta(Theta - Theta0)**2
!
!V(Urey-Bradley) = Kub(S - S0)**2
!
!Ktheta:  kcal/mole/rad**2
!Theta0:  degrees
!S0:  A
!
!Ktheta (kcal/mol) = Ktheta (K) * Boltz.  const.
!
!atom types         Ktheta        Theta0
CH3   CH1   CH3     62.100125     112.00 !  TraPPE 2

Some CHARMM ANGLES section entries include Urey-Bradley potentials (\(K_{ub}\), \(b_{ub}\)), in addition to the standard quadratic angle potential. The constants related to this potential function are currently read, but the logic has not been added to calculate this potential function. Support for this potential function will be added in later versions of the code.

DIHEDRALS

The final major bonded interactions section of the CHARMM compliant parameter file are the DIHEDRALS. Dihedral energies were represented by a cosine series where \(\phi\) is the dihedral angle, \(C_n\) are dihedral force constants, \(n\) is the multiplicity, and \(\delta_n\) is the phase shift. Often, there are 4 to 6 terms in a dihedral. Angles for the dihedrals’ deltas are given in degrees.

\[\begin{split}U_{\texttt{dihedral}}&= C_0 + \sum_{\texttt{n = 1}} C_n [1+\cos(n\phi_i - \delta_n)]\\\end{split}\]

Since isobutane has no dihedral, here are the parameters pertaining to 2,3-dimethylbutane:

DIHEDRALS
!
!V(dihedral) = Kchi(1 + cos(n(chi) - delta))
!
!Kchi:  kcal/mole
!n:  multiplicity
!delta:  degrees
!
!Kchi (kcal/mol) = Kchi (K) * Boltz.  const.
!
!atom types             Kchi    n     delta   description
X   CH1   CH1   X    -0.498907  0     0.0   !  TraPPE 2
X   CH1   CH1   X     0.851974  1     0.0   !  TraPPE 2
X   CH1   CH1   X    -0.222269  2   180.0   !  TraPPE 2
X   CH1   CH1   X     0.876894  3     0.0   !  TraPPE 2

Note

The code allows the use of ‘X’ to indicate ambiguous positions on the ends. This is useful because this kind is often determined solely by the two middle atoms in the middle of the dihedral, according to literature.

Note

If a dihedral parameter was defined with multiplicity value of zero (\(n\) = 0), GOMC will automatically assign the phase shift value to 90 (\(\delta_n\) = 90) to recover the above dihedral expresion.

IMPROPERS

Energy parameters used to describe out-of-plane rocking are currently read, but unused. The section is often blank. If it becomes necessary, algorithms to calculate the improper energy will need to be added.

NONBONDED

The next section of the CHARMM style parameter file is the NONBONDED. The nonbonded energy in CHARMM is presented as 12-6 potential where, \(r_{ij}\), \(\epsilon_{ij}\), \({R_{min}}_{ij}\) are the separation, minimum potential, and minimum potential distance, respectively. In order to use TraPPE this section of the CHARMM compliant file is critical.

\[\begin{split}U_{\texttt{LJ}}&=\sum_{\texttt{nonbonded}} \epsilon_{ij}\left[\left(\frac{R_{min_{ij}}}{r_{ij}}\right)^{12}-2\left(\frac{R_{min_{ij}}}{r_{ij}}\right)^6\right] \\\end{split}\]

Here’s an example with our isobutane potential model:

NONBONDED
!
!V(Lennard-Jones) = Eps,i,j[(Rmin,i,j/ri,j)**12 - 2(Rmin,i,j/ri,j)**6]
!
!atom ignored epsilon         Rmin/2        ignored   eps,1-4     Rmin/2,1-4
CH3   0.0     -0.194745992  2.10461634058     0.0       0.0       0.0 !  TraPPE 1
CH1   0.0     -0.019872040  2.62656119304     0.0       0.0       0.0 !  TraPPE 2
End

Note

The \(R_{min}\) is the potential well-depth, where the attraction is maximum. However, \(\sigma\) is the particle diameter, where the interaction energy is zero. To convert \(\sigma\) to \(R_{min}\), simply multiply \(\sigma\) by 0.56123102415.

Important

If no parameter was defined for 1-4 interaction e.g (\(\epsilon_{1-4}, Rmin_{1-4}/2\)), GOMC will use the \(\epsilon, Rmin/2\) for 1-4 interaction.

NBFIX

The last section of the CHARMM style parameter file is the NBFIX. In this section, individual pair interaction will be modified. First, pseudo non-bonded parameters have to be defined in NONBONDED and modified in NBFIX. Here iss an example if it is required to modify interaction between CH3 and CH1 atoms:

NBFIX
!V(Lennard-Jones) = Eps,i,j[(Rmin,i,j/ri,j)**12 - 2(Rmin,i,j/ri,j)**6]
!
!atom atom  epsilon         Rmin          eps,1-4   Rmin,1-4
CH3   CH1   -0.294745992    1.10461634058 !
End

Important

If no parameter was defined for 1-4 interaction e.g (\(\epsilon_{1-4}, Rmin_{1-4}\)), GOMC will use the \(\epsilon, Rmin\) for 1-4 interaction.

Exotic or Mie Parameter File

The Mie file is intended for use with nonstandard/specialty models of molecular interaction, which are not included in CHARMM standard.

Mie Potential

\[E_{ij} = C_{n_{ij}} \epsilon_{ij} \bigg[\bigg(\frac{\sigma_{ij}}{r_{ij}}\bigg)^{n_{ij}} - \bigg(\frac{\sigma_{ij}}{r_{ij}}\bigg)^6\bigg]\]

where \(r_{ij}\), \(\epsilon_{ij}\), and \(\sigma_{ij}\) are, respectively, the separation, minimum potential, and collision diameter for the pair of interaction sites \(i\) and \(j\). The constant \(C_n\) is a normalization factor such that the minimum of the potential remains at \(-\epsilon_{ij}\) for all \(n_{ij}\). In the 12-6 potential, \(C_n\) reduces to the familiar value of 4.

\[C_{n_{ij}} = \bigg(\frac{n_{ij}}{n_{ij} - 6} \bigg)\bigg(\frac{n_{ij}}{6} \bigg)^{6/(n_{ij} - 6)}\]

Buckingham Potential (Exp-6)

\[\begin{split}E_{ij} = \begin{cases} \frac{\alpha_{ij}\epsilon_{ij}}{\alpha_{ij}-6} \bigg[\frac{6}{\alpha_{ij}} exp\bigg(\alpha_{ij} \bigg[1-\frac{r_{ij}}{R_{min,ij}} \bigg]\bigg) - {\bigg(\frac{R_{min,ij}}{r_{ij}}\bigg)}^6 \bigg] & r_{ij} \geq R_{max,ij} \\ \infty & r_{ij} < R_{max,ij} \end{cases}\end{split}\]

where \(r_{ij}\), \(\epsilon_{ij}\), and \(R_{min,ij}\) are, respectively, the separation, minimum potential, and minimum potential distance for the pair of interaction sites \(i\) and \(j\). The constant \(\alpha_{ij}\) is an exponential-6 parameter. The cutoff distance \(R_{max,ij}\) is the smallest positive value for which \(\frac{dE_{ij}}{dr_{ij}}=0\).

Note

In order to use Mie or Exotice potential file format for Buckingham potential, instead of defining \(R_{min}\), we define \(\sigma\) (collision diameter or the distance, where potential is zero) and GOMC will calculate the \(R_{min}\) and \(R_{max}\) using Buckingham potential equation.

Currently, two custom interaction are included:

  • NONBODED_MIE This section describes n-6 (Lennard-Jones) or Exp-6 (Buckingham) non-bonded interactions. The Lennard-Jones potential (12-6) is a subset of Mie potential. Non-bonded parameters are assigned by specifying the following fields in order:

    1. Atom type name

    2. Minimum energy (\(\epsilon\))

    3. Atom diameter (\(\sigma\))

    4. Repulsion exponent (\(n\)) in Mie potential or \(\alpha\) in Buckingham potential.

    The 1-4 interaction can be modified by specifying the following fields in order:

    1. Minimum energy (\(\epsilon_{1-4}\))

    2. Atom diameter (\(\sigma_{1-4}\))

    3. Repulsion exponent (\(n_{1-4}\)) in Mie potential or \(\alpha_{1-4}\) in Buckingham potential.

    Note

    If no parameter is provided for 1-4 interaction, same parameters (item 2, 3, 4) would be considered for 1-4 interaction.

  • NBFIX_MIE This section allows n-6 (Lennard-Jones) or Exp-6 (Buckingham) interaction between two pairs of atoms to be modified. Interaction between two pairs of atoms can be modified by specifying the following fields in order:

    1. Atom type 1 name

    2. Atom type 2 name

    3. Minimum energy (\(\epsilon\))

    4. Atom diameter (\(\sigma\))

    5. Repulsion exponent (\(n\)) in Mie potential or \(\alpha\) in Buckingham potential.

    The 1-4 interaction between two pairs of atoms can be modified by specifying the following fields in order:

    1. Minimum energy (\(\epsilon_{1-4}\))

    2. Atom diameter (\(\sigma_{1-4}\))

    3. Repulsion exponent (\(n_{1-4}\)) in Mie potential or \(\alpha_{1-4}\) in Buckingham potential.

    Note

    If no parameter is provided for 1-4 interaction, same parameters (item 3, 4, 5) would be considered for 1-4 interaction.

Note

In Mie or Buckingham potential, the definition of atom diameter(\(\sigma\)) is same for both NONBONDED_MIE and NBFIX_MIE.

Important

If no parameter was defined for 1-4 interaction e.g (\(\epsilon_{1-4}, \sigma_{1-4}, n_{1-4}\)), GOMC will use the \(\epsilon, \sigma, n\) for 1-4 interaction.

Otherwise, the Mie file reuses the same geometry section headings - BONDS / ANGLES / DIHEDRALS / etc. The only difference in these sections versus in the CHARMM format force field file is that the energies are in Kelvin (‘K’), the unit most commonly found for parameters in Monte Carlo chemical simulation literature. This precludes the need to convert to kcal/mol, the energy unit used in CHARMM. The most frequently used section of the Mie files in the Mie potential section is NONBONDED_MIE.

Here is the example of Mie or Exotic parameters file format that are used to simulate alkanes with Mie potential:

NONBONDED_MIE
!
!V(Mie) = const*eps*((sig/r)^n-(sig/r)^6)
!
!atom eps       sig     n     eps,1-4   sig,1-4   n,1-4
CH4   161.00    3.740   14    0.0       0.0       0.0 ! Potoff, et al. '09
CH3   121.25    3.783   16    0.0       0.0       0.0 ! Potoff, et al. '09
CH2    61.00    3.990   16    0.0       0.0       0.0 ! Potoff, et al. '09

NBFIX_MIE
!V(Mie) = const*eps*((sig/r)^n-(sig/r)^6)
!
!atom atom  epsilon  sig     n     eps,1-4   sig,1-4   n,1-4
CH3   CH2   100.00   3.8     16    0.0       0.0       0.0 !
End

Here is the example of Mie or Exotic parameters file format that are used to simulate water with Buckingham potential:

NONBONDED_MIE
!
!V(exp-6) = ((eps-ij * alpha)/(alpha - 6)) * ((6 / alpha) * exp(alpha * [1 - (r / rmin)]) - (rmin / r)^6))
!
!atom eps       sig     alpha     eps,1-4   sig,1-4   n,1-4
OT    159.78    3.195   12        0.0       0.0       0.0 ! Errington, et al. 1998
HT      0.0     0.0      0        0.0       0.0       0.0 ! Errington, et al. 1998

NBFIX_MIE
!V(exp-6) = ((eps-ij * alpha)/(alpha - 6)) * ((6 / alpha) * exp(alpha * [1 - (r / rmin)]) - (rmin / r)^6))
!
!atom atom  epsilon  sig     alpha     eps,1-4   sig,1-4   n,1-4
HT   OT      0.00    0.0     0         0.0       0.0       0.0 !
End

Note

Although the units (Angstroms) are the same, the Mie file uses \(\sigma\), not the \(R_{min}\) used by CHARMM. The energy in the exotic file are expressed in Kelvin (K), as this is the standard convention in the literature.

Control File (*.conf)

The control file is GOMC’s proprietary input file. It contains key settings. The settings generally fall under three categories:

  • Input/Simulation Setup

  • System Settings for During Run

  • Output Settings

Note

The control file is designed to recognize logic values, such as “yes/true/on” or “no/false/off”. The keyword in control file is not case sensitive.

Input/Simulation Setup

In this section, input file names are listed. In addition, if you want to restart your simulation or use integer seed for running your simulation, you need to modify this section according to your purpose.

Restart

Determines whether to restart the simulation from restart file (*_restart.pdb) or not.

  • Value 1: Boolean - True if restart, false otherwise.

ExpertMode

Determines whether to perform error checking of move selection to ensure correct ensemble is sampled. This allows the user to run a simulation with no volume moves in NPT, NPT-GEMC; no molecule transfers in GCMC, GEMC.

  • Value 1: Boolean - True if enable expert mode; false otherwise.

Checkpoint

Determines whether to restart the simulation from checkpoint file or not. Restarting the simulation with would result in an identitcal outcome, as if previous simulation was continued. This is required for hybrid Monte-Carlo Molecular Dyanamics in open-ensembles (GCMC/GEMC) to concatenate trajectory files since the molecular transfers rearranges the order of the molecules. Checkpointing will ensure the molecules are loaded in the same order each cycle.

  • Value 1: Boolean - True if restart with checkpoint file, false otherwise.

  • Value 2: String - Sets the name of the checkpoint file.

    Checkpoint   true        AR_KR_continued.chk
    
PRNG

Dictates how to start the pseudo-random number generator (PRNG)

  • Value 1: String

    • RANDOM: Randomizes Mersenne Twister PRNG with random bits based on the system time.

    #################################
    # kind {RANDOM, INTSEED}
    #################################
    PRNG   RANDOM
    
    • INTSEED: This option “seeds” the Mersenne Twister PRNG with a standard integer. When the same integer is used, the generated PRNG stream should be the same every time, which is helpful in tracking down bugs.

Random_Seed

Defines the seed number. If “INTSEED” is chosen, seed number needs to be specified; otherwise, the program will terminate.

  • Value 1: ULONG - If “INTSEED” option is selected for PRNG (See bellow example)

#################################
# kind {RANDOM, INTSEED}
#################################
PRNG          INTSEED
Random_Seed    50
ParaTypeCHARMM

Sets force field type to CHARMM style.

  • Value 1: Boolean - True if it is CHARMM forcefield, false otherwise.

#################################
# FORCE FIELD TYPE
#################################
ParaTypeCHARMM    true
ParaTypeEXOTIC or ParaTypeMie

Sets force field type to Mie style.

  • Value 1: Boolean - True if it is Mie forcefield, false otherwise.

#################################
# FORCE FIELD TYPE
#################################
ParaTypeEXOTIC    true
ParaTypeMARTINI

Sets force field type to MARTINI style.

  • Value 1: Boolean - True if it is MARTINI forcefield, false otherwise.

#################################
# FORCE FIELD TYPE
#################################
ParaTypeMARTINI     true
Parameters

Provides the name and location of the parameter file to use for the simulation.

  • Value 1: String - Sets the name of the parameter file.

#################################
# FORCE FIELD TYPE
#################################
ParaTypeCHARMM    yes
Parameters        ../../common/Par_TraPPE_Alkanes.inp
Parameters        ../../common/Par_TIP4P.inp

Note

More than one parameter file can be provided.

Coordinates

Defines the PDB file names (coordinates) and location for each box in the system.

  • Value 1: Integer - Sets box number (starts from ‘0’).

  • Value 2: String - Sets the name of PDB file.

Note

NVT and NPT ensembles requires only one PDB file and GEMC/GCMC requires two PDB files. If the number of PDB files is not compatible with the simulation type, the program will terminate.

Example of NVT or NPT ensemble:

#############################################
# INPUT PDB FILES - NVT or NPT ensemble
#############################################
Coordinates   0   STEP3_START_ISB_sys.pdb

Example of Gibbs or GC ensemble:

#############################################
# INPUT PDB FILES - Gibbs or GCMC ensemble
#############################################
Coordinates   0   STEP3_START_ISB_sys_BOX_0.pdb
Coordinates   1   STEP3_START_ISB_sys_BOX_1.pdb

Note

In case of Restart true, the restart PDB output file from GOMC (OutputName_BOX_N_restart.pdb) can be used for each box.

Example of Gibbs ensemble when Restart mode is active:

#################################
# INPUT PDB FILES
#################################
Coordinates   0   ISB_T_270_k_BOX_0_restart.pdb
Coordinates   1   ISB_T_270_k_BOX_1_restart.pdb
Structures

Defines the PSF filenames (structures) for each box in the system.

  • Value 1: Integer - Sets box number (start from ‘0’)

  • Value 2: String - Sets the name of PSF file.

Note

NVT and NPT ensembles requires only one PSF file and GEMC/GCMC requires two PSF files. If the number of PSF files is not compatible with the simulation type, the program will terminate.

Example of NVT or NPT ensemble:

#################################
# INPUT PSF FILES
#################################
Structure   0   STEP3_START_ISB_sys.psf

Example of Gibbs or GC ensemble:

#################################
# INPUT PSF FILES
#################################
Structure   0   STEP3_START_ISB_sys_BOX_0.psf
Structure   1   STEP3_START_ISB_sys_BOX_1.psf

Note

In case of Restart true, the PSF output file from GOMC (OutputName_BOX_N_restart.psf) can be used for both boxes.

Example of Gibbs ensemble when Restart mode is active:

#################################
# INPUT PSF FILES
#################################
Structure   0   ISB_T_270_k_BOX_0_restart.psf
Structure   1   ISB_T_270_k_BOX_1_restart.psf
binCoordinates

Defines the DCD file names (coordinates) and location for each box in the system.

  • Value 1: Integer - Sets box number (starts from ‘0’).

  • Value 2: String - Sets the name of PDB file.

Note

NVT and NPT ensembles requires only one DCD file and GEMC/GCMC requires only one PDB files, although loading two is supported. This is different from PDB files, for which two are required in GEMC/GCMC. This allows the user to only load binary coordinates for one box.

Example of NVT or NPT ensemble:

#############################################
# INPUT PDB FILES - NVT or NPT ensemble
#############################################
binCoordinates   0   STEP3_START_ISB_sys.coor

Example of Gibbs or GC ensemble:

#############################################
# INPUT PDB FILES - Gibbs or GCMC ensemble
#############################################
binCoordinates   0   STEP3_START_ISB_sys_BOX_0.coor
binCoordinates   1   STEP3_START_ISB_sys_BOX_1.coor

Note

In case of Restart, the restart DCD output file from GOMC (OutputName_BOX_N_restart.coor) can be used for each box.

Example of Gibbs ensemble when Restart mode is active:

#################################
# INPUT PDB FILES
#################################
binCoordinates   0   ISB_T_270_k_BOX_0_restart.coor
binCoordinates   1   ISB_T_270_k_BOX_1_restart.coor
binVelocities

Defines the VEL file names (velocities) and location for each box in the system.

  • Value 1: Integer - Sets box number (starts from ‘0’).

  • Value 2: String - Sets the name of VEL file.

Note

Originate from a Molecular Dynamics softwrae such as NAMD. GOMC will only output a velocity restart file if it is provided one using this keyword.

Note

In hybrid Monte-Carlo Molecular Dynamics, the velocities of the atoms should be preserved across cycles to increase accuracy. These files are not used internally by GOMC, only maintained. If a molecular transfer occurs, a new velocity is generated by Langevin dynamics.

Example of NVT or NPT ensemble:

#############################################
# INPUT PDB FILES - NVT or NPT ensemble
#############################################
binVelocities   0   STEP3_START_ISB_sys.vel

Example of Gibbs or GC ensemble:

#############################################
# INPUT PDB FILES - Gibbs or GCMC ensemble
#############################################
binVelocities   0   STEP3_START_ISB_sys_BOX_0.vel
binVelocities   1   STEP3_START_ISB_sys_BOX_1.vel

Note

In case of Restart, the restart VEL output file from GOMC (OutputName_BOX_N_restart.vel) can be used for each box.

Example of Gibbs ensemble when Restart mode is active:

#################################
# INPUT PDB FILES
#################################
binVelocities   0   ISB_T_270_k_BOX_0_restart.vel
binVelocities   1   ISB_T_270_k_BOX_1_restart.vel
extendedSystem

Defines the XSC file names (box dimensions and origin) and location for each box in the system.

  • Value 1: Integer - Sets box number (starts from ‘0’).

  • Value 2: String - Sets the name of XSC file.

Note

Previously, this information was stored in plain-text format at the top of restart PDB files. This will be deprecated in favor of binary XSC.

Example of NVT or NPT ensemble:

#############################################
# INPUT PDB FILES - NVT or NPT ensemble
#############################################
extendedSystem   0   STEP3_START_ISB_sys.xsc

Example of Gibbs or GC ensemble:

#############################################
# INPUT PDB FILES - Gibbs or GCMC ensemble
#############################################
extendedSystem   0   STEP3_START_ISB_sys_BOX_0.xsc
extendedSystem   1   STEP3_START_ISB_sys_BOX_1.xsc

Note

In case of Restart, the restart XSC output file from GOMC (OutputName_BOX_N_restart.xsc) can be used for each box.

Example of Gibbs ensemble when Restart mode is active:

#################################
# INPUT PDB FILES
#################################
extendedSystem   0   ISB_T_270_k_BOX_0_restart.xsc
extendedSystem   1   ISB_T_270_k_BOX_1_restart.xsc
MultiSimFolderName

The name of the folder to be created which contains output from the multisim.

  • Value 1: String - Name of the folder to contain output

MultiSimFolderName  outputFolderName

System Settings for During Run Setup

This section contains all the variables not involved in the output of data during the simulation, or in the reading of input files at the start of the simulation. In other words, it contains settings related to the moves, the thermodynamic constants (based on choice of ensemble), and the length of the simulation. Note that some tags, or entries for tags, are only used in certain ensembles (e.g. Gibbs ensemble). These cases are denoted with colored text.

GEMC

(For Gibbs Ensemble runs only) Defines the type of Gibbs Ensemble simulation you want to run. If neglected in Gibbs Ensemble, it simply defaults to const volume (NVT) Gibbs Ensemble.

  • Value 1: String - Allows you to pick between isovolumetric (“NVT”) and isobaric (“NPT”) Gibbs ensemble simulations.

Note

The default value for GEMC is NVT.

#################################
# GEMC TYPE (DEFAULT IS NVT GEMC)
#################################
GEMC    NVT
Pressure

For NPT or NPT-GEMC simulation, imposed pressure (in bar) needs to be specified; otherwise, the program will terminate.

  • Value 1: Double - Constant pressure in bar.

#################################
# GEMC TYPE (DEFAULT IS NVT GEMC)
#################################
GEMC        NPT
Pressure    5.76
Temperature

Sets the temperature at which the system will run.

  • Value 1: Double - Constant temperature of simulation in degrees Kelvin.

#################################
# SIMULATION CONDITION
#################################
Temperature   270.00

(MPI-GOMC Only)

  • Value 1: List of Doubles - A list of constant temperatures for simulations in degrees Kelvin.

#################################
# SIMULATION CONDITION
#################################
Temperature   270.00    280.00    290.00    300.00

Note

To use more than one temperature, GOMC must be compiled in MPI mode. Also, if GOMC is compiled in MPI mode, more than one temperature is required. To use only one temperature, use standard GOMC.

Rcut

Sets a specific radius that non-bonded interaction energy and force will be considered and calculated using defined potential function.

  • Value 1: Double - The distance to truncate the Lennard-Jones potential at.

RcutLow

Sets a specific minimum possible in angstrom that reject any move that places any atom closer than specified distance.

  • Value 1: Double - The minimum possible distance between any atoms.

RcutCoulomb

Sets a specific radius for each box in the system that short range electrostatic energy will be calculated.

  • Value 1: Integer - Sets box number (start from ‘0’)

  • Value 2: Double - The distance to truncate the short rage electrostatic energy at.

Note

The default value for RcutCoulomb is the value of Rcut

Important

  • In Ewald Summation method, at constant Tolerance and box volume, increasing RcutCoulomb would result is decreasing reciprocal vector [Fincham 1993]. Decreasing the reciprocal vector decreases the computation time in long range electrostatic calculation.

  • Increasing the RcutCoulomb results in increasing the computation time in short range electrostatic calculation.

  • Parallelization of Ewald summation method is done on reciprocal vector loop, rather than molecule loop. So, in case of running on multiple CPU threads or GPU, it is better to use the lower value for RcutCoulomb, to maximize the parallelization efficiency.

  • There is an optimum value for RcutCoulomb, where result in maximum effeciency of the method. We encourage to run a short simulation with various RcutCoulomb to find the optimum value.

LRC

Defines whether or not long range corrections are used.

  • Value 1: Boolean - True to consider long range correction.

Note

In case of using SHIFT or SWITCH potential functions, LRC will be ignored.

Exclude

Defines which pairs of bonded atoms should be excluded from non-bonded interactions.

  • Value 1: String - Allows you to choose between “1-2”, “1-3”, and “1-4”.

    • 1-2: All interactions pairs of bonded atoms, except the ones that separated with one bond, will be considered and modified using 1-4 parameters defined in parameter file.

    • 1-3: All interaction pairs of bonded atoms, except the ones that separated with one or two bonds, will be considered and modified using 1-4 parameters defined in parameter file.

    • 1-4: All interaction pairs of bonded atoms, except the ones that separated with one, two or three bonds, will be considered using non-bonded parameters defined in parameter file.

    Note

    The default value for Exclude is “1-4”.

    Note

    In CHARMM force field, the 1-4 interaction needs to be considered. Choosing “Exclude 1-3” will modify 1-4 interaction based on 1-4 parameters in parameter file. If a kind of force field is used, where 1-4 interaction needs to be ignored, such as TraPPE, either “Exclude 1-4” needs to be chosen or 1-4 parameter needs to be assigned to zero in the parameter file.

Potential

Defines the potential function type to calculate non-bonded interaction energy and force between atoms.

  • Value 1: String - Allows you to pick between “VDW”, “EXP6”, “SHIFT” and “SWITCH”.

    • VDW: Nonbonded interaction energy and force calculated based on n-6 (Lennard-Johns) equation. This function will be discussed further in the Intermolecular energy and Virial calculation section.

      #################################
      # SIMULATION CONDITION
      #################################
      Temperature   270.00
      Potential     VDW
      LRC           true
      Rcut          10
      Exclude       1-4
      
    • EXP6: Nonbonded interaction energy and force calculated based on exp-6 (Buckingham potential) equation. This function will be discussed further in the Intermolecular energy and Virial calculation section.

      #################################
      # SIMULATION CONDITION
      #################################
      Temperature   270.00
      Potential     EXP6
      LRC           true
      Rcut          10
      Exclude       1-4
      
    • SHIFT: This option forces the potential energy to be zero at Rcut distance. This function will be discussed further in the Intermolecular energy and Virial calculation section.

      #################################
      # SIMULATION CONDITION
      #################################
      Temperature     270.00
      Potential       SHIFT
      LRC             false
      Rcut            10
      Exclude         1-4
      RcutCoulomb  0  12.0
      RcutCoulomb  1  20.0
      
    • SWITCH: This option smoothly forces the potential energy to be zero at Rcut distance and starts modifying the potential at Rswitch distance. Depending on force field type, specific potential function will be applied. These functions will be discussed further in the Intermolecular energy and Virial calculation section.

    Rswitch

    In the case of choosing “SWITCH” as potential function, a distance is set in which non-bonded interaction energy is truncated smoothly at Rcut distance.

    • Value 1: Double - Define switch distance in angstrom. If the “SWITCH” function is chosen, Rswitch needs to be defined; otherwise, the program will be terminated.

    #################################
    # SIMULATION CONDITION
    #################################
    Temperature   270.00
    Potential     SWITCH
    LRC           false
    Rcut          12
    Rswitch       9
    Exclude       1-4
    
VDWGeometricSigma

Use geometric mean, as required by OPLS force field, to combining Lennard-Jones sigma parameters for different atom types.

  • Value 1: Boolean - True, uses geometric mean to combine L-J sigmas

    Note

    The default setting of VDWGeometricSigma is false to use arithmetic mean when combining Lennard-Jones sigma parameters for different atom types.

ElectroStatic

Considers coulomb interaction or not. This function will be discussed further in the Inter- molecular energy and Virial calculation section.

  • Value 1: Boolean - True if coulomb interaction needs to be considered and false if not.

    Note

    To simulate the polar molecule in MARTINI force field, ElectroStatic needs to be turn on. MARTINI force field uses short range coulomb interaction with constant Dielectric 15.0.

Ewald

Implements the standard Ewald summation method for electrostatic calculation. This function will be discussed further in the Intermolecular energy and Virial calculation section.

  • Value 1: Double - True if Ewald summation calculation needs to be considered and false if not.

    Note

    By default, ElectroStatic will be set to true if the Ewald summation method was used to calculate coulomb interaction.

CachedFourier

Stores the reciprocal terms for Ewald summation calculation to improve code performance. This option increases code performance at the cost of more memory usage. The effects of this option depend on the system size and the value of RcutCoulomb, as the number of k-vectors increases with the value of RcutCoulomb. Small systems may see no or minimal change in performance and memory usage.

  • Value 1: Boolean - True to store reciprocal terms of Ewald summation calculation for reuse and false to not store them.

    Note

    The default, the value of CachedFourier is set to true.

    Note

    The CachedFourier option is ignored when running on the GPU.

    Warning

    Monte Carlo moves, such as MEMC-1, MEMC-2, MEMC-3, IntraMEMC-1, IntraMEMC-2, IntraMEMC-3 do not support CachedFourier.

Tolerance

Specifies the accuracy of the Ewald summation calculation. Ewald separation parameter and number of reciprocal vectors for the Ewald summation are determined based on the accuracy parameter.

  • Value 1: Double - Sets the accuracy in Ewald summation calculation.

    Note

    • A reasonable value for the accuracy is 0.00001.

    • If “Ewald” was chosen and no value was set for Tolerance, the program will terminate.

Dielectric

Defines dielectric constant for coulomb interaction in MARTINI force field.

  • Value 1: Double - Sets dielectric value used in coulomb interaction.

    Note

    • In MARTINI force field, Dielectric needs to be set to 15.0.

    • If MARTINI force field was chosen and Dielectric was not specified, a default value of 15.0 will be assigned.

PressureCalc

Considers to calculate the pressure or not. If it is set to true, the frequency of pressure calculation need to be set.

  • Value 1: Boolean - True enabling pressure calculation during the simulation, false disabling pressure calculation.

  • Value 2: Ulong - The frequency of calculating the pressure.

1-4scaling

Defines constant factor to modify intra-molecule coulomb interaction.

  • Value 1: Double - A fraction number between 0.0 and 1.0.

    Note

    CHARMM force field uses a value between 0.0 and 1.0. In MARTINI force field, it needs to be set to 1.0 because 1-4 interaction will not be modified in this force field.

    #################################
    # SIMULATION CONDITION
    #################################
    ElectroStatic   true
    Ewald           true
    Tolerance       0.00001
    CachedFourier   false
    1-4scaling      0.0
    
RunSteps

Sets the total number of steps to run (one move is performed for each step) (cycles = this value / number of molecules in the system)

  • Value 1: Ulong - Total run steps

Note

RunSteps is a delta.

Important

Setting the RunSteps to zero, and activating Restart simulation, will recalculate the energy of stored simulation’s snapshots.

EqSteps

Sets the number of steps necessary to equilibrate the system; averaging will begin at this step.

  • Value 1: Ulong - Equilibration steps

Note

EqSteps is not a delta. If restarting a simulation with a start step greater than EqSteps, no equilibration is performed.

Note

In GCMC simulation, the Histogram files will be outputed at EqSteps.

AdjSteps

Sets the number of steps per adjustment of the parameter associated with each move (e.g. maximum translate distance, maximum rotation, maximum volume exchange, etc.)

  • Value 1: Ulong - Number of steps per move adjustment

InitStep

Sets the first step of the simulation.

  • Value 1: Ulong - Number of first step of simulation.

Note

Hybrid Monte-Carlo Molecular Dynamics (py-MCMD) requires resetting start step to 0 for combination of NAMD and GOMC data.

#################################
# STEPS
#################################
RunSteps    25000000
EqSteps     5000000
AdjSteps    1000
InitStep    0
ChemPot

For Grand Canonical (GC) ensemble runs only: Chemical potential at which simulation is run.

  • Value 1: String - The residue name to apply this chemical potential.

  • Value 2: Double - The chemical potential value in degrees Kelvin (should be negative).

Note

  • For binary systems, include multiple copies of the tag (one per residue kind).

  • If there is a molecule kind that cannot be transfer between boxes (in PDB file the beta value is set to 1.00 or 2.00), an arbitrary value (e.g. 0.00) can be assigned to the residue name.

#################################
# Mol.  Name Chem.  Pot.  (K)
#################################
ChemPot   ISB     -968
Fugacity

For Grand Canonical (GC) ensemble runs only: Fugacity at which simulation is run.

  • Value 1: String - The residue to apply this fugacity.

  • Value 2: Double - The fugacity value in bar.

Note

  • For binary systems, include multiple copies of the tag (one per residue kind).

  • If there is a molecule kind that cannot be transfer between boxes (in PDB file the beta value is set to 1.00 or 2.00) an arbitrary value e.g. 0.00 can be assigned to the residue name.

#################################
# Mol.  Name Fugacity (bar)
#################################
Fugacity  ISB   10.0
Fugacity  Si     0.0
Fugacity  O      0.0
DisFreq

Fractional percentage at which displacement move will occur.

  • Value 1: Double - % Displacement

RotFreq

Fractional percentage at which rigid rotation move will occur.

  • Value 1: Double - % Rotatation

IntraSwapFreq

Fractional percentage at which molecule will be removed from a box and inserted into the same box using coupled-decoupled configurational-bias algorithm.

  • Value 1: Double - % Intra molecule swap

Note

The default value for IntraSwapFreq is 0.000

IntraTargetedSwapFreq

Fractional percentage at which molecule will be removed from the box and inserted into a subvolume within the same box, or deleted from the subvolume and inserted into the same box using coupled-decoupled configurational-bias algorithm.

  • Value 1: Double - % Intra molecule swap

Note

The default value for IntraTargetedSwapFreq is 0.000

RegrowthFreq

Fractional percentage at which part of the molecule will be deleted and then regrown using coupled-decoupled configurational-bias algorithm.

  • Value 1: Double - % Molecular growth

Note

The default value for RegrowthFreq is 0.000

CrankShaftFreq

Fractional percentage at which crankshaft move will occur. In this move, two atoms that are forming angle or dihedral are selected randomely and form a shaft. Then any atoms or group that are within these two selected atoms, will rotate around the shaft to sample intramolecular degree of freedom.

  • Value 1: Double - % Crankshaft

Note

The default value for CrankShaftFreq is 0.000

MultiParticleFreq

Fractional percentage at which multi-particle move will occur. In this move, all molecules in the selected simulation box will be rigidly rotated or displaced simultaneously, along the calculated torque or force, respectively.

  • Value 1: Double - % Multiparticle

Note

The default value for MultiParticleFreq is 0.000

MultiParticleBrownianFreq

Fractional percentage at which multi-particle brownian move will occur. In this move, all molecules in the selected simulation box will be rigidly rotated or displaced simultaneously, along the calculated torque or force, respectively.

  • Value 1: Double - % Multiparticle

Note

The default value for MultiParticleBrownianFreq is 0.000

IntraMEMC-1Freq

Fractional percentage at which specified number of small molecule kind will be exchanged with a specified large molecule kind in defined sub-volume within same simulation box.

  • Value 1: Double - % Molecular exchange

Note

  • The default value for IntraMEMC-1Freq is 0.000

  • This move need additional information such as ExchangeVolumeDim, ExchangeRatio, ExchangeSmallKind, and ExchangeLargeKind, which will be explained later.

  • For more information about this move, please refere to MEMC-GCMC and MEMC-GEMC papers.

IntraMEMC-2Freq

Fractional percentage at which specified number of small molecule kind will be exchanged with a specified large molecule kind in defined sub-volume within same simulation box. Backbone of small and large molecule kind will be used to insert the large molecule more efficiently.

  • Value 1: Double - % Molecular exchange

Note

  • The default value for IntraMEMC-2Freq is 0.000

  • This move need additional information such as ExchangeVolumeDim, ExchangeRatio, ExchangeSmallKind, ExchangeLargeKind, SmallKindBackBone, and LargeKindBackBone, which will be explained later.

  • For more information about this move, please refere to MEMC-GCMC and MEMC-GEMC papers.

IntraMEMC-3Freq

Fractional percentage at which specified number of small molecule kind will be exchanged with a specified large molecule kind in defined sub-volume within same simulation box. Specified atom of the large molecule kind will be used to insert the large molecule using coupled-decoupled configurational-bias.

  • Value 1: Double - % Molecular exchange

Note

  • The default value for IntraMEMC-3Freq is 0.000

  • This move need additional information such as ExchangeVolumeDim, ExchangeRatio, ExchangeSmallKind, ExchangeLargeKind, and LargeKindBackBone, which will be explained later.

  • For more information about this move, please refere to MEMC-GCMC and MEMC-GEMC papers.

MEMC-1Freq

For Gibbs and Grand Canonical (GC) ensemble runs only: Fractional percentage at which specified number of small molecule kind will be exchanged with a specified large molecule kind in defined sub-volume in dense simulation box.

  • Value 1: Double - % Molecular exchange

Note

  • The default value for MEMC-1Freq is 0.000

  • This move need additional information such as ExchangeVolumeDim, ExchangeRatio, ExchangeSmallKind, and ExchangeLargeKind, which will be explained later.

  • For more information about this move, please refere to MEMC-GCMC and MEMC-GEMC papers.

MEMC-2Freq

For Gibbs and Grand Canonical (GC) ensemble runs only: Fractional percentage at which specified number of small molecule kind will be exchanged with a specified large molecule kind in defined sub-volume in dense simulation box. Backbone of small and large molecule kind will be used to insert the large molecule more efficiently.

  • Value 1: Double - % Molecular exchange

Note

  • The default value for MEMC-2Freq is 0.000

  • This move need additional information such as ExchangeVolumeDim, ExchangeRatio, ExchangeSmallKind, ExchangeLargeKind, SmallKindBackBone, and LargeKindBackBone, which will be explained later.

  • For more information about this move, please refere to MEMC-GCMC and MEMC-GEMC papers.

MEMC-3Freq

For Gibbs and Grand Canonical (GC) ensemble runs only: Fractional percentage at which specified number of small molecule kind will be exchanged with a specified large molecule kind in defined sub-volume in dense simulation box. Specified atom of the large molecule kind will be used to insert the large molecule using coupled-decoupled configurational-bias.

  • Value 1: Double - % Molecular exchange

Note

  • The default value for MEMC-3Freq is 0.000

  • This move need additional information such as ExchangeVolumeDim, ExchangeRatio, ExchangeSmallKind, ExchangeLargeKind, and LargeKindBackBone, which will be explained later.

  • For more information about this move, please refere to MEMC-GCMC and MEMC-GEMC papers.

SwapFreq

For Gibbs and Grand Canonical (GC) ensemble runs only: Fractional percentage at which molecule swap move will occur using coupled-decoupled configurational-bias.

  • Value 1: Double - % Molecule swaps

TargetedSwapFreq

For Gibbs and Grand Canonical (GC) ensemble runs only: Fractional percentage at which targeted molecule swap move will occur using coupled-decoupled configurational-bias in the sub-volumes specified.

  • Value 1: Double - % Molecule targeted swaps

VolFreq

For isobaric-isothermal ensemble and Gibbs ensemble runs only: Fractional percentage at which molecule will be removed from one box and inserted into the other box using configurational bias algorithm.

  • Value 1: Double - % Volume swaps

#################################
# MOVE FREQEUNCY
#################################
DisFreq         0.39
RotFreq         0.10
IntraSwapFreq   0.10
RegrowthFreq    0.10
CrankShaftFreq  0.10
SwapFreq        0.20
VolFreq         0.01

Warning

All move percentages should add up to 1.0; otherwise, the program will terminate.

ExchangeVolumeDim

To use all variation of MEMC and IntraMEMC Monte Carlo moves, the exchange sub-volume must be defined. The exchange sub-volume is defined as an orthogonal box with x-, y-, and z-dimensions, where small molecule/molecules kind will be selected from to be exchanged with a large molecule kind.

  • Value 1: Double - X dimension in \(Å\)

  • Value 2: Double - Y dimension in \(Å\)

  • Value 3: Double - Z dimension in \(Å\)

Note

  • Currently, the X and Y dimension cannot be set independently (X = Y = max(X, Y))

  • A heuristic for setting good values of the x-, y-, and z-dimensions is to use the geometric size of the large molecule plus 1-2 Å in each dimension.

  • In case of exchanging 1 small molecule kind with 1 large molecule kind in IntraMEMC-2, IntraMEMC-3, MEMC-2, MEMC-3 Monte Carlo moves, the sub-volume dimension has no effect on acceptance rate.

ExchangeSmallKind

To use all variation of MEMC and IntraMEMC Monte Carlo moves, the small molecule kind to be exchanged with a large molecule kind must be defined. Multiple small molecule kind can be specified.

  • Value 1: String - Small molecule kind (resname) to be exchanged.

ExchangeLargeKind

To use all variation of MEMC and IntraMEMC Monte Carlo moves, the large molecule kind to be exchanged with small molecule kind must be defined. Multiple large molecule kind can be specified.

  • Value 1: String - Large molecule kind (resname) to be exchanged.

ExchangeRatio

To use all variation of MEMC and IntraMEMC Monte Carlo moves, the exchange ratio must be defined. The exchange ratio defines how many small molecule will be exchanged with 1 large molecule. For each large-small molecule pairs, one exchange ratio must be defined.

  • Value 1: Integer - Ratio of exchanging small molecule/molecules with 1 large molecule.

LargeKindBackBone

To use MEMC-2, MEMC-3, IntraMEMC-2, and IntraMEMC-3 Monte Carlo moves, the large molecule backbone must be defined. The backbone of the molecule is defined as a vector that connects two atoms belong to the large molecule. The large molecule backbone will be used to align the sub-volume in MEMC-2 and IntraMEMC-2 moves, while in MEMC-3 and IntraMEMC-3 moves, it uses the atom name to start growing the large molecule using coupled-decoupled configurational-bias. For each large-small molecule pairs, two atom names must be defined.

  • Value 1: String - Atom name 1 belong to the large molecule’s backbone

  • Value 2: String - Atom name 2 belong to the large molecule’s backbone

Important

  • In MEMC-3 and IntraMEMC-3 Monte Carlo moves, both atom names must be same, otherwise program will be terminated.

  • If large molecule has only one atom (mono atomic molecules), same atom name must be used for Value 1 and Value 2 of the LargeKindBackBone.

SmallKindBackBone

To use MEMC-2, and IntraMEMC-2 Monte Carlo moves, the small molecule backbone must be defined. The backbone of the molecule is defined as a vector that connects two atoms belong to the small molecule and will be used to align the sub-volume. For each large-small molecule pairs, two atom names must be defined.

  • Value 1: String - Atom name 1 belong to the small molecule’s backbone

  • Value 2: String - Atom name 2 belong to the small molecule’s backbone

Important

  • If small molecule has only one atom (mono atomic molecules), same atom name must be used for Value 1 and Value 2 of the SmallKindBackBone.

Here is the example of MEMC-2 Monte Carlo moves, where 7 large-small molecule pairs are defined with an exchange ratio of 1:1: (ethane, methane), (propane, ethane), (n-butane, propane), (n-pentane, nbutane), (n-hexane, n-pentane), (n-heptane, n-hexane), and (noctane, n-heptane).

######################################################################
# MEMC PARAMETER
######################################################################
ExchangeVolumeDim   1.0   1.0   1.0
ExchangeRatio       1       1       1      1      1      1      1
ExchangeLargeKind   C8P    C7P    C6P    C5P    C4P    C3P    C2P
ExchangeSmallKind   C7P    C6P    C5P    C4P    C3P    C2P    C1P
LargeKindBackBone   C1 C8  C1 C7  C1 C6  C1 C5  C1 C4  C1 C3  C1 C2
SmallKindBackBone   C1 C7  C1 C6  C1 C5  C1 C4  C1 C3  C1 C2  C1 C1

Here is the example of MEMC-3 Monte Carlo moves, where 7 large-small molecule pairs are defined with an exchange ratio of 1:1: (ethane, methane), (propane, ethane), (n-butane, propane), (n-pentane, nbutane), (n-hexane, n-pentane), (n-heptane, n-hexane), and (noctane, n-heptane).

######################################################################
# MEMC PARAMETER
######################################################################
ExchangeVolumeDim   1.0   1.0   1.0
ExchangeRatio       1       1       1      1      1      1      1
ExchangeLargeKind   C8P    C7P    C6P    C5P    C4P    C3P    C2P
ExchangeSmallKind   C7P    C6P    C5P    C4P    C3P    C2P    C1P
LargeKindBackBone   C4 C4  C4 C4  C3 C3  C3 C3  C2 C2  C2 C2  C1 C1
SmallKindBackBone   C1 C7  C1 C6  C1 C5  C1 C4  C1 C3  C1 C2  C1 C1

Here is the example of MEMC-2 Monte Carlo moves, where 1 large-small molecule pair is defined with an exchange ratio of 1:2: (xenon, methane).

######################################################################
# MEMC PARAMETER
######################################################################
ExchangeVolumeDim   5.0   5.0   5.0
ExchangeRatio       2
ExchangeLargeKind   XE
ExchangeSmallKind   C1P
LargeKindBackBone   Xe Xe
SmallKindBackBone   C1 C1
SubVolumeBox

Define which box the subvolume occupies. - Value 1: Integer - Sub-volume id. - Value 2: Integer - Sets box number (first box is box ‘0’).

SubVolumeCenter

Define the center of the static subvolume. - Value 1: Integer - Sub-volume id. - Value 2: Double - x value of SubVolumeCenter \(Å\). - Value 3: Double - y value of SubVolumeCenter \(Å\). - Value 4: Double - z value of SubVolumeCenter \(Å\).

SubVolumePBC

Define which dimensions periodic box wrapping is applied in the subvolume. - Value 1: Integer - Sub-volume id. - Value 2: String - X, Y, Z, XY, XZ, YZ, XYZ (Axes should have no spaced between them)

SubVolumeCenterList

Define the center of the dynamic subvolume by defining the atoms to use for the geometric mean calculation. - Value 1: Integer - Sub-volume id. - Value 2: Integer Range - Atom indices used to calculate geometric center of subvolume.

SubVolumeDim

Define the dimensions of the subvolume. - Value 1: Integer - Sub-volume id. - Value 2: Double - x value of SubVolumeDim \(Å\). - Value 3: Double - y value of SubVolumeDim \(Å\). - Value 4: Double - z value of SubVolumeDim \(Å\).

SubVolumeResidueKind

Define which residue kinds can be inserted or deleted from the subvolume. - Value 1: Integer - Sub-volume id. - Value 2: String - Residue kind inserted/deleted from subvolume - Value .: String - Residue kind inserted/deleted from subvolume - Value .: String - Residue kind inserted/deleted from subvolume - Value N: String - Residue kind inserted/deleted from subvolume

SubVolumeRigidSwap

Define whether molecules are held rigid or the geometry is sampled per the coupled-decoupled CBMC scheme. - Value 1: Integer - Sub-volume id. - Value 2: Boolean - If true the molecule is held rigid. If false, geometry is sampled when inserting in the subvolume.

SubVolumeChemPot

Define the chemical potential of a residue kind in the subvolume. Only used in TargetedSwap, not IntraTargetedSwap. - Value 1: Integer - Sub-volume id. - Value 2: String - Residue kind - Value 3: Double - Chemical potential

SubVolumeFugacity

Define the fugacity of a residue kind in the subvolume. Only used in TargetedSwap, not IntraTargetedSwap. - Value 1: Integer - Sub-volume id. - Value 2: String - Residue kind - Value 3: Double - Chemical potential

######################################################################
# TARGETED SWAP (Static subVolume)
######################################################################
SubVolumeBox                  1       0
SubVolumeCenter               1       25.0 25.0 25.0
SubVolumeDim                  1       35 35 5
SubVolumeResidueKind          1       TIP3
SubVolumeRigidSwap            1       false
SubVolumeChemPot              1       TIP3    -800
######################################################################
# TARGETED SWAP (Dynamic subVolume)
######################################################################
SubVolumeBox                  1       0
SubVolumeCenterList           1       1-402
SubVolumeDim                  1       35 35 5
SubVolumeResidueKind          1       TIP3
SubVolumeRigidSwap            1       false
SubVolumeChemPot              1       TIP3    -800
useConstantArea

For Isobaric-Isothermal ensemble and Gibbs ensemble runs only: Considers to change the volume of the simulation box by fixing the cross-sectional area (x-y plane).

  • Value 1: Boolean - If true volume will change only in z axis, If false volume will change with constant axis ratio.

Note

By default, useConstantArea will be set to false if no value was set. It means, the volume of the box will change in a way to maintain the constant axis ratio.

FixVolBox0

For adsorption simulation in NPT Gibbs ensemble runs only: Changing the volume of fluid phase (Box 1) to maintain the constant imposed pressure and temperature, while keeping the volume of adsorbed phase (Box 0) fix.

  • Value 1: Boolean - If true volume of adsorbed phase will remain constant, If false volume of adsorbed phase will change.

CellBasisVector

Defines the shape and size of the simulation periodic cell. CellBasisVector1, CellBasisVector2, CellBasisVector3 represent the cell basis vector \(a,b,c\), respectively. This tag may occur multiple times. It occurs once for NVT and NPT, but twice for Gibbs ensemble or GC ensemble.

  • Value 1: Integer - Sets box number (first box is box ‘0’).

  • Value 2: Double - x value of cell basis vector \(Å\).

  • Value 3: Double - y value of cell basis vector \(Å\).

  • Value 4: Double - z value of cell basis vector \(Å\).

Note

If the number of defined boxes were not compatible to simulation type, the program will be terminated.

Example for NVT and NPT ensemble. In this example, each vector is perpendicular to the other two (\(\alpha = 90, \beta = 90, \gamma = 90\)), as indicated by a single x, y, or z value being specified by each and making a rectangular 3-D box:

############################################
# BOX DIMENSION #, X, Y, Z
############################################
CellBasisVector1  0   40.00   00.00   00.00
CellBasisVector2  0   00.00   40.00   00.00
CellBasisVector3  0   00.00   00.00   80.00

Example for Gibbs ensemble and GC ensemble ensemble. In this example, In the first box, only vector \(a\) and \(c\) are perpendicular to each other (\(\alpha = 90, \beta = 90, \gamma = 120\)), and making a non-orthogonal simulation cell with the cell length \(a = 36.91 Å, b = 36.91 Å, c = 76.98 Å\). In the second box, each vector is perpendicular to the other two (\(\alpha = 90, \beta = 90, \gamma = 90\)), as indicated by a single x, y, or z value being specified by each and making a cubic box:

############################################
# BOX DIMENSION #, X, Y, Z
############################################
CellBasisVector1  0   36.91   00.00   00.00
CellBasisVector2  0   -18.45  31.96   00.00
CellBasisVector3  0   00.00   00.00   76.98

CellBasisVector1  1   60.00   00.00   00.00
CellBasisVector2  1   00.00   60.00   00.00
CellBasisVector3  1   00.00   00.00   60.00

Warning

If Restart was activated, box dimension does not need to be specified. If it is specified, program will read it but it will be ignored and replaced by the printed cell dimensions and angles in the restart PDB output file from GOMC (OutputName_BOX_0_restart.pdb and OutputName_BOX_1_restart.pdb).

CBMC_First

Number of CD-CBMC trials to choose the first atom position (Lennard-Jones trials for first seed growth).

  • Value 1: Integer - Number of initial insertion sites to try.

CBMC_Nth

Number of CD-CBMC trials to choose the later atom positions (Lennard-Jones trials for first seed growth).

  • Value 1: Integer - Number of LJ trials for growing later atom positions.

CBMC_Ang

Number of CD-CBMC bending angle trials to perform for geometry (per the coupled-decoupled CBMC scheme).

  • Value 1: Integer - Number of trials per angle.

CBMC_Dih

Number of CD-CBMC dihedral angle trials to perform for geometry (per the coupled-decoupled CBMC scheme).

  • Value 1: Integer - Number of trials per dihedral.

#################################
# CBMC TRIALS
#################################
CBMC_First  10
CBMC_Nth    4
CBMC_Ang    100
CBMC_Dih    30

Next section specifies the parameters that will be used for free energy calculation in NVT and NPT ensembles.

FreeEnergyCalc

For NVT and NPT ensemble only: Considers to calculate the free energy data (the energy different between current lambda state and all other neighboring lambda states, and calculate the derivative of energy with respective to current lambda) or not. If it is set to true, the frequency of free energy calculation need to be set. The free energy data will be printed into Free_Energy_BOX_0_ OutputName.dat.

  • Value 1: Boolean - True enabling free energy calculation during the simulation, false disabling the calculation.

  • Value 2: Ulong - The frequency of calculating the free energy.

MoleculeType

Sets the solute molecule kind (residue name) and molecule number (residue ID), which absolute solvation free will be calculated for.

  • Value 1: String - The solute name (residue name).

  • Value 2: Integer - The solute molecule number (residue ID).

InitialState

Sets the index of the LambdaCoulomb and LambdaVDW vectors, to determine the simulation lambda value for VDW and Coulomb interactions.

  • Value 1: Integer - The index of LambdaCoulomb and LambdaVDW vectors.

Note

  • Multiple initial states need to be run to perform a free energy analysis on the system.

LambdaVDW

Sets the intermediate lambda states to which solute-solvent VDW interaction to be scaled.

  • Value 1: Double - Lambda values for VDW interaction in ascending order.

Warning

All lambda values must be stated in the ascending order, otherwise the program will terminate.

LambdaCoulomb

Sets the intermediate lambda states to which solute-solvent Coulombic interaction to be scaled.

  • Value 1: Double - Lambda values for Coulombic interaction in ascending order.

Warning

All lambda values must be stated in the ascending order, otherwise the program will terminate.

Note

  • By default, the lambda values for Coulombic interaction will be set to zero if ElectroStatic or Ewald is deactivated.

  • By default, the lambda values for Coulombic interaction will be set to Lambda values for VDW interaction if ElectroStatic or Ewald is activated.

-The LambdaVDW and LambdaCoulomb lists must be equal in length.

ScaleCoulomb

Determines to scale the Coulombic interaction non-linearly (soft-core scheme) or not.

  • Value 1: Boolean - True if coulombic interaction needs to be scaled non-linearly, False if coulombic interaction needs to be scaled linearly.

Note

By default, the ScaleCoulomb will be set to false.

ScalePower

Sets the \(p\) value in soft-core scaling scheme, where the distance between solute and solvent is scaled non-linearly.

  • Value 1: Integer - The \(p\) value in the soft-core scaling scheme.

Note

By default, the ScalePower will be set to 2.

ScaleAlpha

Sets the \(\alpha\) value in soft-core scaling scheme, where the distance between solute and solvent is scaled non-linearly.

  • Value 1: Double - \(\alpha\) value in the soft-core scaling scheme.

Note

By default, the ScaleAlpha will be set to 0.5.

MinSigma

Sets the minimum \(\sigma\) value in soft-core scaling scheme, where the distance between solute and solvent is scaled non-linearly.

  • Value 1: Double - Minimum \(\sigma\) value in the soft-core scaling scheme.

Note

By default, the MinSigma will be set to 3.0.

Note

Scaling the distance between solute and solvent using soft-core scheme:

\[r_{sc} = \bigg[\alpha {\big(1 - \lambda \big)}^{p}{\sigma}^6 + {r}^6 \bigg]^{\frac{1}{6}}\]

Here is the example of solvation free energy of CO2, in intermediate state 3.

#################################
# FREE ENERGY PARAMETERS
#################################
FreeEnergyCalc true   1000
MoleculeType   CO2   1
InitialState   3
ScalePower     2
ScaleAlpha     0.5
MinSigma       3.0
ScaleCoulomb   false
#states        0    1    2    3    4
LambdaVDW      0.00 0.50 1.00 1.00 1.00
LambdaCoulomb  0.00 0.00 0.00 0.50 1.00

Output Controls

This section contains all the values that control output in the control file. For example, certain variables control the naming of files outputed of the block-averaged thermodynamic variables of interest, the PDB files, etc.

OutputName

Unique name with no space for simulation used to name the block average, PDB, and PSF output files.

  • Value 1: String - Unique phrase to identify this system.

#################################
# OUTPUT FILE NAME
#################################
OutputName  ISB_T_270_K
CoordinatesFreq

Controls output of PDB file (coordinates). If PDB outputing was enabled, one file for NVT or NPT and two files for Gibbs ensemble or GC ensemble will be outputed into OutputName_BOX_n.pdb, where n defines the box number.

  • Value 1: Boolean - “true” enables outputing these files; “false” disables outputing.

  • Value 2: Ulong - Steps per dump PDB frame. It should be less than or equal to RunSteps. If this keyword could not be found in configuration file, its value will be assigned a default value to dump 10 frames.

Note

  • DCDFreq should be used unless the low precision and slower PDB trajectory is needed, perhaps beta and occupancy values are desired. The PDB trajectory is much larger and consumes more disk space.

  • The PDB file contains an entry for every ATOM, in all boxes read. This allows VMD (which requires a constant number of atoms) to properly parse frames, with a bit of help. Atoms that are not currently in a specific box are given the coordinate (0.00, 0.00, 0.00). The occupancy value corresponds to the box a molecule is currently in (e.g. 0.00 for box 0; 1.00 for box 1).

  • At the beginning of simulation, a merged PSF file will be outputed into OutputName_merged.psf, in which all boxes will be outputed. It also contains the topology for every molecule in both boxes, corresponding to the merged PDB format. Loading PDB files into merged PSF file in VMD allows the user to visualize and analyze the results.

DCDFreq

Controls output of DCD file (binary coordinates). If DCD outputing was enabled, one file for NVT or NPT and two files for Gibbs ensemble or GC ensemble will be outputed into OutputName_BOX_n.dcd, where n defines the box number.

  • Value 1: Boolean - “true” enables outputing these files; “false” disables outputing.

  • Value 2: Ulong - Steps per dump PDB frame. It should be less than or equal to RunSteps. If this keyword could not be found in configuration file, its value will be assigned a default value to dump 10 frames.

Note

  • The DCD file contains an entry for every ATOM, in all boxes read. This allows VMD (which requires a constant number of atoms) to properly parse frames, with a bit of help. Atoms that are not currently in a specific box are given the coordinate (0.00, 0.00, 0.00). The occupancy value corresponds to the box a molecule is currently in (e.g. 0.00 for box 0; 1.00 for box 1).

  • At the beginning of simulation, a merged PSF file will be outputed into OutputName_merged.psf, in which all boxes will be outputed. It also contains the topology for every molecule in both boxes, corresponding to the merged PDB format. Loading DCD files into merged PSF file in VMD allows the user to visualize and analyze the results.

  • “The DCD files are single precision binary FORTRAN files, so are transportable between computer architectures. The file readers in NAMD and VMD can detect and adapt to the endianness of the machine on which the DCD file was written, and the utility program flipdcd is also provided to reformat these files if needed. The exact format of these files is very ugly but supported by a wide range of analysis and display programs. The timestep is stored in the DCD file in NAMD internal units and must be multiplied by TIMEFACTOR=48.88821 to convert to fs. Positions in DCD files are stored in Å. Velocities in DCD files are stored in NAMD internal units and must be multiplied by PDBVELFACTOR=20.45482706 to convert to Å/ps. Forces in DCD files are stored in kcal/mol/Å.”

  • source : https://www.ks.uiuc.edu/Research/namd/2.9/ug/node11.html

RestartFreq

Controls the output of the last state of simulation at a specified step in

  • PDB files (coordinates)

  • PSF files (structure)

  • XSC files (binary box dimensions)

  • COOR files (binary coordinates)

  • If provided as input: VEL files (binary velocity)

  • Value 1: Boolean - “true” enables printing block average; “false” disables it.

  • Value 2: Ulong - Number of steps per checkpoint output file.

OutputName_BOX_n_restart.*, where n defines the box number. Header part of this file contains important information and will be needed to restart the simulation:

Restart PDB files, one file for NVT or NPT and two files for Gibbs ensemble or GC ensemble, will be outputed with the following information. - Simulation cell dimensions and angles.

  • Maximum amount of displacement (Å), rotation (\(\delta\)), and volume (\(Å^3\)) that used in Displacement, Rotation, and Volume move.

Note

  • The restart PDB/PSF/COOR/VEL files contains only ATOM that exist in each boxes at specified steps. These box restart files allows the user to load a box into NAMD and run molecular dynamics in Hybrid Monte-Carlo Molecular Dynamics (py-MCMD).

  • When restarting the GOMC simulation from two restart files, the order of the molecules in the trajectory may differ preventing trajectory concatenation, unless the CHK file is loaded.

  • Only restart files must be used to begin a GOMC simulation with Restart simulation active. The merged psf is NOT a restart file.

  • CoordinatesFreq must be a common multiple of RestartFreq or vice versa.

CheckpointFreq

Controls the output of the checkpoint file containing values neccessary to continue a simulation as if it never ended.

  • Value 1: Boolean - “true” enables printing block average; “false” disables it.

  • Value 2: Ulong - Number of steps per checkpoint output file.

Note

  • RestartFreq must also be enabled, and CheckpointFreq must equal RestartFreq.

ConsoleFreq

Controls the output to STDIO (“the console”) of messages such as acceptance statistics, and run timing info. In addition, instantaneously-selected thermodynamic properties will be output to this file.

  • Value 1: Boolean - “true” enables message printing; “false” disables outputing.

  • Value 2: Ulong - Number of steps per print. If this keyword could not be found in the configuration file, the value will be assigned by default to dump 1000 output for RunSteps greater than 1000 steps and 100 output for RunSteps less than 1000 steps.

BlockAverageFreq

Controls the block averages output of selected thermodynamic properties. Block averages are averages of thermodynamic values of interest for chunks of the simulation (for post-processing of averages or std. dev. in those values).

  • Value 1: Boolean - “true” enables printing block average; “false” disables it.

  • Value 2: Ulong - Number of steps per block-average output file. If this keyword cannot be found in the configuration file, its value will be assigned a default to dump 100 output.

HistogramFreq

Controls the histograms. Histograms are a binned listing of observation frequency for a specific thermodynamic variable. In this code, they also control the output of a file containing energy/molecule samples; it only will be used in GC ensemble simulations for histogram reweighting purposes.

  • Value 1: Boolean - “true” enables printing histogram; “false” disables it.

  • Value 2: Ulong - Number of steps per histogram output file. If this keyword cannot be found in the configuration file, a value will be assigned by default to dump 1000 output for RunSteps greater than 1000 steps and 100 output for RunSteps less than 1000 steps.

#################################
# STATISTICS Enable, Freq.
#################################
CoordinatesFreq   true 10000000
RestartFreq       true 1000000
CheckpointFreq    true 1000000
ConsoleFreq       true 100000
BlockAverageFreq  true 100000
HistogramFreq     true 10000

The next section controls the output of the energy/molecule sample file and the distribution file f or molecule counts, commonly referred to as the “histogram” output. This section is only required if Grand Canonical ensemble simulation was used.

DistName

Sets short phrase to naming molecule distribution file.

  • Value 1: String - Short phrase which will be combined with RunNumber and RunLetter to use in the name of the binned histogram for molecule distribution.

HistName

Sets short phrase to naming energy sample file.

  • Value 1: String - Short phrase, which will be combined with RunNumber and RunLetter, to use in the name of the energy/molecule count sample file.

RunNumber

Sets a number, which is a part of DistName and HistName file name.

  • Value 1: Uint – Run number to be used in the above file names.

RunLetter

Sets a letter, which is a part of DistName and HistName file name.

  • Value 1: Character – Run letter to be used in above file names.

SampleFreq

Controls histogram sampling frequency.

  • Value 1: Uint – the number of steps per histogram sample.

#################################
# OutHistSettings
#################################
DistName   dis
HistName   his
RunNumber  5
RunLetter  a
SampleFreq 200
OutEnergy, OutPressure, OutMolNum, OutDensity, OutVolume, OutSurfaceTension

Enables/Disables for specific kinds of file output for tracked thermodynamic quantities

  • Value 1: Boolean – “true” enables message output of block averages via this tracked parameter (and in some cases such as entry, components); “false” disables it.

  • Value 2: Boolean – “true” enables message output of a fluctuation into the console file via this tracked parameter (and in some cases, such as entry, components); “false” disables it.

The keywords are available for the following ensembles

Keyword

NVT

NPT & Gibbs

GC

OutEnergy

\(\checkmark\)

\(\checkmark\)

\(\checkmark\)

OutPressure

\(\checkmark\)

\(\checkmark\)

\(\checkmark\)

OutMolNum

\(\checkmark\)

\(\checkmark\)

OutDensity

\(\checkmark\)

\(\checkmark\)

OutVolume

\(\checkmark\)

\(\checkmark\)

OutSurfaceTension

\(\checkmark\)

Here is an example:

#################################
# ENABLE: BLK AVE., FLUC.
#################################
OutEnergy         true true
OutPressure       true true
OutMolNum         true true
OutDensity        true true
OutVolume         true true
OutSurfaceTension false false

Binary input file types

Binary representations of the system.
  • XSC

  • COOR

  • VEL

  • CHK

  • “NAMD uses a trivial double-precision binary file format for coordinates, velocities, and forces. Due to its high precision this is the default output and restart format. VMD refers to these files as the ``namdbin'' format. The file consists of the atom count as a 32-bit integer followed by all three position or velocity components for each atom as 64-bit double-precision floating point, i.e., NXYZXYZXYZXYZ... where N is a 4-byte int and X, Y, and Z are 8-byte doubles. If the number of atoms the file contains is known then the atom count can be used to determine endianness. The file readers in NAMD and VMD can detect and adapt to the endianness of the machine on which the binary file was written, and the utility program flipbinpdb is also provided to reformat these files if needed. Positions in NAMD binary files are stored in Å. Velocities in NAMD binary files are stored in NAMD internal units and must be multiplied by PDBVELFACTOR=20.45482706 to convert to Å/ps. Forces in NAMD binary files are stored in kcal/mol/Å.”

XSC (eXtended System Configuration file) File

GOMC allows the box dimensions to be defined in one of three ways:

  • In the control file

  • In the header of restart PDB file

  • In a binary XSC file

The XSC file contains the first step of the simulation, cell vectors, and cell origin. Currently, GOMC only uses the cell vectors.

COOR (binary coordinates) File

GOMC allows the box coordinates to be overwritten by a binary coordinates file. The COOR file should have the same number of atoms in it as the PDB file which it is overwriting. The actual coordinates can vary dramatically, which allows the user to sample the coordinates with other engines (MCMD), or transform it however one sees fit.

VEL (binary velocity) File

GOMC allows the velocities associated with each atom to be maintained and output for continuing MD simulations. In the event a molecule transfer occurs, all the atoms of the transferred molecule are given new velocities by Langevin dynamics. These VEL files must originate from NAMD, as GOMC will not produce them without first being provided them.

CHK (checkpoint) File

GOMC contains several variables which, if not accounted for, will produce different outputs even if the initial conditions are exactly the same. These variables are contained in the checkpoint file, and allow the user to pick up a GOMC simulation where it left off without altering the course of the simulation. Also, the checkpoint file is essential for MCMD as molecules are treated as distinguishable in molecular dynamics due to the fact that MD is a continuous trajectory through time. The checkpoint file contains the original atom order of the molecules, and coordinates and velocities are loaded into this order to ensure the trajectories are consistently arranged.

Checkpoint file contents:

  • Last simulation step that saved into checkpoint file (Start step can be overriden).

  • True number of simulation steps that have been run.

  • Maximum amount of displacement (Å), rotation (\(\delta\)), and volume (\(\AA^3\)) that used in Displacement, Rotation, MultiParticle, and Volume move.

  • Number of Monte Carlo move trial and acceptance.

  • Random number sequence.

  • Molecule lookup object.

  • Original pdb atoms object to reload new positions into.

  • Original molecule setup object generated from parsing first PSF files.

  • Accessory data for coordinating loading the restart coordinates into the original ordering.

  • If built with MPI and parallel tempering was enabled: Random number sequence for parallel tempering.