Users' Guide to the Monte Carlo Program MMC

for Computer Simulation of Molecular Solutions.

Mihaly Mezei

Department of Pharmacological Sciences

Icahn School of Medicine at Mount Sinai

E-mail address: Mihaly.Mezei@mssm.edu

April 04, 2022

TABLE OF CONTENTS

I. GENERAL DESCRIPTION

The present document describes the Force-Biased Metropolis Monte Carlo program for computer simulation studies of molecular liquids and solutions in the canonical, isobaric-isothermal and grand-canonical ensembles and the analysis of such calculations based on the Proximity Criterion. This section provides an overview of the capabilities of the program.

I.1. Systems modeled

The present program was originally written for simulations on infinitely dilute solutions (one solute in (numsol) solvent molecules). As a result of this origin, the program has developed two different treatments of molecules. Solvent molecules have to be rigid and be of the same type (this allows the use of special, rapid energy calculation code). The collection of all other molecules (i.e., the one(s) not declared to be solvents) are called by the program the solute. Thus the solute can consist of a single rigid molecule, a collection of freely translating and rotating rigid molecules or a collection of freely translating and rotating molecules with certain (user specified) torsion angles freely rotating.

The program establishes the connectivity of the solute molecule(s) based on interatomic distances and built-in threshold values. By default H-H bonds and bonds between an H and a heavy atom on different groups (residues) are ignored, but this can be changed (see the keys BDHH and BDRH). Also, the bond thresholds can be changed with the MODA key. This sometimes leads to missing or extra bonds, especally if the input structure (given by the SLTA key) is obtained from a molecular dynamics simulation. Such defects can be identified by using the FCGA and BNDL keys and corrected by the keys MAKB or BRKB or forced to be ignored (ok for molecules without active torsions) with the MOLD key.

Free energy simulations can be performed on selected molecules of the solute (see Sec. I.5.).

The torsion angles to be changed on the solute are specified by the two atoms defining the bond of the torsion (with the key TORD). The torsions can be grouped and it is the group of torsions that is selected for change. The main reason for this is that if a torsion affects atoms that are on an other rotatable bond then it may be efficient to move the affected bond as well. These moves, however, are not correlated (i.e., each is generated by a different random number). There is also an option of changing the torsion angles of a loop segment (with the key LOOP). The torsion potential parameters can be obtained from a CHARMM or AMBER parameter file.

It is possible to impose periodicity on the solute, either along the X axis (with the key SLTA POLY) i.e., with period equal to the edge of the cell in the X direction - an infinite polymer, or in all three directions (with the key SLTA XTAL) with its unit cell (or a multiple of a unit cell) coinciding with the periodic cell of the simulation (i.e., a crystal). The periodic extension of the solute outside the periodic cell is done automatically by the program.

The intermolecular interactions are assumed to be of the usual 12-6-1 type, i.e., Lennard-Jones terms for the repulsion and dispersion plus electrostatics based on partial charges on atoms (or pseudo atoms). The parameters can be obtained from libraries stores by the program or defined via input. The interactions can be truncated (see keys SVVC, SUVC, SUUC) either by a spherical cutoff (each atom interacts with all others within a certain distance) or by the minimum image convention (each atom interacts with one periodic image of all the other atoms). For the purpose of cutoff or image calculation, distances usually refer to molecule-molecule or residue-residue distances.

When spherical truncation is used, the long-range electrostatic energy can be estimated by the reaction field formula (provided all residues are neutral) - see key RFCR.

I.2. Properties calculated during a simulation

The major outputs of a Monte Carlo simulation (besides the configurations saved for further analysis by the TRAJ key) are as follows.

I.3. Potential functions

Several potential functions are built into the program. For the solute-solvent nonbonded interactions (specified with the SUPT key). Clementi et al.'s library, EPEN/QPEN, Kollman et al.'s libraries (AMBER-94 and AMBER-2002), CHARMM22 library, Berendsen-Van Gunsteren library (GROMOS87 and GROMOS96); and Jorgensen's OPLS library. Note that for the GROMOS force field only the combination rule can be used - the non-bonded parameters for pairs are not implemented. Except for the EPEN form, parameters are stored internally, but they can be modified from input and allowance is made for inputting new atom types besides the stored ones. As these parameter libraries are continually updated and extended, it is advisable to check if the value stored is the best one (using the PLBP key). For aqueous solutions, the solvent-solvent interaction (specified with the SUPT key) can be either in the form of MCY-CI type or TIPS/TIP3P/SPC type or TIPS2/TIP4PB-F type or TIP5P type potential. For the MCY-CI type the parameters of both the MCY-CI-II and the Yoon-Morokuma-Davidson parametrizations are built in. For TIPS2/TIP4P/B-F types, the A, B and q parameters may also be provided by the user, but the TIPS, SPC, TIP3P, Bernal-Fowler, TIPS2, TIP4P and TIP5P water model parameters are built in. For solvent other than water, 12-6-1 type potential can be used with parameters from the AMBER, GROMOS or OPLS libraries.

There is a framework for adding field-dependent potentials to the currently implemented atom-atom potentials, but this requires some additional coding.

I.4. Sampling techniques

Various sampling techniques can be selected with the MOVE key. Translations and rotations can be sampled with the original Metropolis Method or with the Force-Biased algorithm of Pangali, Rao and Berne (not implemented for multimolecular solute moves, though). Force biasing, of course, necessitates the recalculation of energies and forces at every accepted step, thereby significantly increasing the computer time required to perform a move. The improvement in convergence, however, more than outweighs the added cost. By default a lambda value of 1/2 is used but the scaled force-bias technique where lambda is scaled down near the solute has also been implemented. Novel torsion-angle sampling techniques include the extension-biased torsion moves and the local moves with the reverse proximity criterion.

Molecules to be moved can be selected randomly, in a cycle or in a shuffled cycle or using the preferential sampling scheme of J.C. Owicki, which leads to improved statistics for the local solution environment of the solute (see key MOVE PRF*). The cavity-biased insertion technique for grand-canonical ensemble simulations and the virial-biased volume-change for the isobaric ensemble have also been implemented.

I.5. Free energy methods

Free energy simulations can be performed using a coupling parameter lambda (CP in the followings) in the energy expression. There are two fundamentally different approaches to generate a path: (A) simultaneous creation/annihilation and (B) continuous deformation. In the first the Hamiltonians of the initial and final states are used in unchanged and the coupling parameter is used to combine the two; in the second case there is a single hybrid Hamiltonian that is based on a parametrization that is a mixture of the two systems. For multimolecular solute, the free energy simulation can be applied to just one or two solute molecules, depending on the type of coupling.

I.5.A. Creation/annihilation path

I.5.A.1. Creation/annihilation path to be used in thermodynamic integration (FREE TICA). Here

E(CP)=CPk*E1+(1-CP)k*Eo

where Eo and E1 represent the Hamiltonians initial and final state, respectively. k=1 represents linear coupling, k>1 has been referred to as or 'almost linear' coupling (e.g., k=4 for 1/r12 repulsion). It is also possible to use different k values for the different types of energy terms (1/r12, 1/r6, and 1/r) with the so-called polynomial TI. Eo and E1 can differ in conformation, partial charges, and chemical identity (the present version requires equal number of atoms in the two states - this can always be achieved with he use of dummy atoms. It is also possible to provide two more coupling parameter values CPo and CP1 such that CPo < CP < CP1. The program will compute the free energy difference between states represented by E(CPo) and E(CP1) with the perturbation method as the sum of two free energy differences: A1-Ao = -kT ln<exp((E(CP1)-E(CP))/kT)>CP + kT ln<exp((E(CPo)-E(CP))/kT)>CP. <E(CPo)>CPo and <E(CP1)>CP1 are also computed for consistency check. The two molecules have to contain the same number of residues (groups), but can have different number of atoms.

I.5.A.2. Perturbation method based on a sequence of intermediate states (FREE TILI). This option performs a perturbation method calculation between two states generated using the deformation-type path described below in I.5.B. The reference state will be (E(CP1)-E(CPo))/2. In this case, the inital and final states have to have the same number of atoms also.

I.5.B. Deformation path

I.5.B.1. The program also provides option for a fundamentally nonlinear coupling where the coupling parameter changes both the solute coordinates and the potential parameters to generate intermediate states (continuous deformation). Besides creating intermediate Eo and E1 states for the linear or 'almost linear' couplings described above in I.5.A.2, this coupling can be used with the probability ratio method (FREE PMF1) or with the overlap ratio method (FREE PMNL), as described below in I.5.B.2 and I.5.B.3, respectively. The initial and final states may differ in position, partial charge or identity of atoms (there are limitations for EPEN).

In each of the following cases, the path involves one or two molecules of the solute. The following paths are supported:

I.5.B.1.1. Linear combination of two conformations of a single molecule (FREE PMF1 GENL). In this case all atoms of the molecule may change positions independently. The configurations at the beginning and at the end of the path are given on input and the intermediate configurations are convex linear combinations (i.e., R(CP)=CP*R1+(1-CP)*Ro) of these two.

I.5.B.1.2. Moving two molecules of the solute in a correlated fashion (FREE PMF1 WRMM). In this case, the intramolecular distances within the two parts should remain constant. Note, that if the orientation of either part is different in the initial and final states then at the intermediate states the molecule will be distorted.

I.5.B.1.3. Varyig the distance between two solute molecules of the solute but keeping the one of them fixed (FREE PMF1 WRFM). In this case the configurations at the two ends of the path differ only in the conformation of the moving part of the solute. This option may be more economical when one part is large and the other is small.

I.5.B.1.4. Moving and rotating one molecule while keeping the other one fixed (FREE PMF1 WRFR). In this case the center of mass of the first part is moved between the center of masses at the two endpoints and the orientation of the moving part is also changed to transform it from the initial to the final orientation. This option differs from I.5.B.1.3 in that at each point of the path the moving molecule has the same intramolecular conformation.

I.5.B.1.5. Correlated rotation of a set of torsion angles (FREE PMF1 TORS). In this case a set of torsion angles are changed in a user-specified interval using a common CP parameter: CP=0 sets all the selected angles to the beginning of their interval and CP=1 sets all of them to the end of their interval.

I.5.B.2. The solvent contribution to the free energy change between two solute conformations can also be determined by computing the ratio of the Boltzmann probabilities as sampled along a line in the conformation space of the solute (described by any of the paths described in Sec. I.5.B.1, FREE PMF1 **** *USS). This approach is called the probability ratio method. Calculations of this type usually require umbrella sampling techniques. The program allows either the use of a harmonic potential to be given from input or the self-consistent determination of a tabulated weighing function using the Adaptive Umbrella Sampling Method.

I.5.B.2. Any of the paths descibed in Sec. I.5.B.1 can be used to estimate free energy changes simultaneosly using the overlap ratio method of Jacucci and Quirke] and the perturbation method (FREE PMNL). Here the calculations are run at a given fixed E(CP) and two more coupling parameter values, CPo and CP1 are specified, CPo < CP < CP1. For the overlap ratio method, the program will calculate the CP averaged distribution of <E(CP)-E(CPo)>CP and <E(CP1)-E(CP)>CP. Separate calculations at CPo and CP1 will give the needed CPo and CP1 averaged distributions to be used with the overlap ratio method. Simultaneously, the program also calculates the perturbation method differences between the states CPo and CP as well as CP and CP1. Notice that this kind of perturbation method calculation is different from the calculation described with the creation/annihilation path above: there the middle state E(CP) was a linear combination of the initial and final states, E(CPo) and E(CP1). For consistency check, (similarly to the perturbation method calculations described above) the CPo average of E(CPo) and the CP1 average of E(CP1) are also calculated.

I.5.C. Chemical potential.

Free energy can be also obtained from the chemical potential. The excess chemical potential can be obtained either from grand-canonical ensemble simulations (key GCEN), from calculating the excess free energy of changing a molecule to ideal gas (with one of the methods described above, e.g., with key FREE TICA) or by the Widom insertion method (key FREE WIDO). The program implements a version that is based on insertions into cavities.

I.6. Solute move selection

The sampling of the moves is specified as follows. The overall frequency of solute moves is specified by input provided with the STEP key. As discussed above, the whole solute can be a) rotated and displaced; b) one constituent molecule of the solute can be rotated and displaced; c) two constituent molecules can be swapped; d) a group of torsion angles can be changed; and e) the coupling parameter used in a potential of mean force calculation can be changed. If moves of type b, c, d, and e are required, then the program will request a relative frequency (weight) of making that move. Regular solute moves and coupling parameter changes are both given a weight of 1.0 . Thus giving a weight of 2.0 for type 3 moves, out of 4 solute move attempts, 1 will change the solute position, 1 will change the coupling parameter, and 2 will change one of the randomly selected torsion angle groups (on the average, of course). Clearly, the more types of moves on is using, the more frequently has the solute to be moved.

For the free energy options, the treatment of moves of type 2 or 3 depends of the type of path. For the creation/annihilation path using TI (FREE TICA) the torsions are independently sampled on both copies, since these can be, and generally are, different molecules. For the continuous path (FREE PMF1 or FREE PMNL) or the cration/annihilation PM (FREE PMLI), the torsions are treated identically in each copy since here it is assumed that the molecule is changing smoothly. This means, that for a creation/annihilation TI these moves have to be sampled twice as frequently.

I.7. Initial configuration

The starting configuration for the Monte Carlo (see key CNFG) is either read from a previously prepared file, generated from the solute's description (see key SLTA) and randomly placed solvents or from a configuration of another system through various transformations. Grand-canonical ensemble simulation (see key GCEN) is also an efficient way of obtaining a random initial configuration.

I.8. Information needed for a simulation

Several steps in the preparation of the input to and the analysis of the results form MMC can be helped by the program Simulaid. The list of its features are described in a J. Comp. Chem. paper.

I.9.1. T, temperature at which the simulations are to be run (specified with the TEMP key).

I.9.2. The coordinates of the nuclei, electrons and pseudoatoms (if any) in a local coordinate framework, and the charges associated with these centers for the solute and (optionally) for the solvent (specified with the SLTA key). The program Simulaid can help convert structures in different file formats (e.g., PDB, CHARMM .CRD) to MMC format and can optimize the orientation of a solute within a periodic box (to maximize the smallest inter-image distance and thereby minimize the artifactual effect of the solute images), within a boundig box, or optimally center a solute in a sphere (for droplet-type calculations).

I.9.3. Potential specifications. For each solute atom the potential type has to be specified from input (except for EPEN where for each atom and pseudoatom the A, B, and C parameters have to be specified individually). The program Simulaid can help to convert from PDB, CHARMM, AMBER or Macromodel input files to MMC format. Several water potentials are implemented with specific codes, but solvent can be defined similarly to the solute as well.

I.9.4. numsol, the number of solvent molecules used in the simulation, specified with the NSLV key.

I.9.5. Size and shape of the periodic cell (specified with the PBCN key).

I.9.6. When generating a new configuration the program can determine nmolec and the simulation cell size parameters from the partial molar volumes of the solute and solvent and the required number of solvent shells around the solute.

I.9.7. For partial solute displacements the groups of solute atoms (independent molecules) to be moved separately (specified with the PARD key).

I.9.8. For solute torsions (specified with the TORD key). the atoms defining the torsion angle and the torsion potential parameters.

I.9. Self testing

There is an extensive self testing routine in the program that tests the integrity/consistency of the data. It protects against both corruption of data and program errors. It is invoked at several ways:

I.10. General analysis

Except for the the bulk distributions and most thermodynamic variables, all the the properties discussed above in Sec. I.2., can be calculated from the history of a simulation as well. Note, that such analyses can be done on simulation trajectories generated by other programs, e.g., by CHARMM or AMBER. In addition to the analyses described above, several other types of analyses are implemented as summarized below.

II. FILE ORGANIZATION.

Beside then standard input (Fortran unit 5) - for the information needed for the run - and the standard output file (unit 6) - printed output, diagnostics, described in Section V. - the program uses a number of different files. In the followings, writing to the standard output will be refererred to as printing.

File names are limited to 80 characters and are structured the following way:

<jobname>[.<numrun>][_<version>].<file type>.

where

III. INPUT DESCRIPTION

The program input is based on keyword driven commands. The structure of a command is as follows:

Of these only the leading keyword appears all the time. Each keyword is a four-character, upper case string. Keywords and unformatted data have to be on the same line. Commands requiring more than one line require the continuation character "~" at the end of the line that is not the end of the command. Several commands can be put on the same line if they are separated by the separator character "|". Text after the comment character "!" is neglected (i.e., considered just comment).

The type of items required are defined by the leading keyword. Auxiliary keywords and data may have default values. The order of keywords and data is fixed, so if there is a needed auxiliary keyword or data, all other items have to be specified that precede the needed one. Formatted data (when required) always start in a line after the keywords and it can not be omitted, but for items with default value, can be secified as zero.

While there is no rigid order of the leading keywords, there are certain rules. Certain leading keywords have prerequisite keywords (i.e., the prerequisite keywords all have to be inputted earlier) and potential precursor keywords that, if used, must be read before the current one or after the next RUNS) or SCAN). If KEY1 is a potential precursor of KEY2 then Key2 is a potential successor of KEY1. If the same leading keyword is used more than once before a RUNS command, a warning is given and the last occurrence will be used (except for TITL).

During input and initializations the data read is cheked for consistency and for current program limits. This may result in messages prefixed by >>>>> OVERRIDE, ----- WARNING, ===== STRONG WARNING or ***** ERROR. At the conclusion of the run, the number of each type of messages is printed. ERROR messages prevent a simulation or an analysis run but don't abort immediately the input. When the input has more than one ERROR messages, some of the additional ERROR messages may only be the consequence of the first one. WARNING message indicates that under some circumstances the input may be incorrect and STRONG WARNING message indicates that the input is likely but not necessarily incorrect or that the error is not fatal. OVERRIDE messages are generally harmless, but they can indicate that the input was not what the user intended.

This section describes the input for each leading keyword and specifies their prerequisites and potential precursors. The description of the commands is given by leading keywords. For example CNFG Key1 [Key2 Key3 Data] Fdata means that the leading keyword CNFG has to be followed by at least one other auxiliary keyword and optionally two more ("Key 2 Key 3"), followed by free-format numerical data ("Data") (all on the line containing the leading keyword), followed by formatted data, started in a new line. Note also that the type of (or the need for) auxiliary keys, free-format and formatted data may depend on the actual choice of auxiliary keys. The possible keywords to be used for Key1, Key2 and Key3 and their effects are explained subsequently, as well as the nature and (for Fdata) the format of the data. Keywords are always given with upper case BOLD characters while variable (data) names are given in lower case bold characters. Input item that is optional is given in [ brackets ] and the default values are given in { braces }. The default value for an auxiliary keyword is the first one on the list.

A record type number (RECTYPE) is assigned to every formatted data record. It is used in error messages on the program's output to identifying the type of data read when the error was encountered.

III.1. File handling keywords

  1. FILE: Job name
  2. WCNF: Configuration saving
  3. WCKP: Creation of a checkpoint file
  4. SCKP: Periodic archiving of checkpoint files
  5. RMCK: Checkpoint file removal
  6. TRAJ: History file specification
  7. MINE: Minimum energy extraction into history file
  8. IDAC: Proximity insertion/deletion acceptance rate
  9. IDLG: Insertion/deletion log file specification
  10. PXPL: Proximity analysis distribution plot file writing
  11. PXWR: Proximity information file specification
  12. WTRA: Trajectory conversions
  13. STIR: Contact elimination
  14. LPAT: Path of local disk directory
  15. BUFF: History file buffer save

III.01.01. Job name

FILE Data

Data: jobname[, numrun]

III.01.02. Configuration saving

WCNF

Prerequisite: CNFG, SLTA
Related key: PDBT.

Key1 (coordinate file format):

III.01.03. Creation of a checkpoint file

WCKP [Data]

Prerequisite: FILE

Data: numrunw

III.01.04. Periodic archiving of checkpoint files

SCKP [Data]

Data: nsavckpf

Each checkpoint file will have the runnumber of the current run and consecutive version numbers, starting at 2. This is option allows the later calculation of proximity analysis averages over different segments of the run (see key CMPC), provides a backup procedure in case the program crashes while writing the checkpoint file, and is useful for debugging.

III.01.05. Checkpoint file removal

RMCK [Data]

Data: numrunp

This option is useful for preventing the disk to become too full when the calculation is broken up into several RUNS or SCAN steps.

III.01.06. History file specification

TRAJ Key1, [Key2,Key3] [Data]

Can not be followed by: MINE
Potential successors: RCKP
Prerequisite of: DENF, WTRA
Variable to set to 'T' in the preprocessor: ME
Related keys: MINE, WTRA, BUFF.

For MC runs this key sets up the saving or reading of the run history on file <jobname>.<runnum>_<runvers>.hst (For runnum=1 the run number part is omitted and for runvers=1 the version part is omitted). See also the key BUFF.

For analysis (SCAN) or filtering (FILT **** TRAJ) this command specifies the format of the history file. An input history can consist of several files with the same runnumber and with consecutive versionnumbers (i.e., the second, third, etc. piece should have names <jobname>.<runnum>_<runvers+1>.hst , <jobname>.<runnum>_<runvers+2>.hst , etc.). Note, that this naming convention also applies to non-MMC formats that have customary extensions (e.g., .pdb, .DCD) requiring either renaming these files or setting up aliases for them.

Key1 (trajectory file format - if any):

When running a MC simulation, full configurations (Key1=ALL*) will be written on the history file at every nmcrec-th MC step (read with the RUNS key).

III.01.07. Minimum energy extraction

MINE Key1

Related key: TRAJ

Key1 (trajectory or overall):

III.01.08. Proximity insertion/deletion acceptance rate

IDAC [Data]

Data: rpxidlim, nrange, irangestart, irangestopnrange times)

III.01.09. Insertion/deletion log file specification

IDLG Key1 [Data]

Prerequisite: GCEN.
Prerequisite of: IDAG
Related key: IDAG.

Key1 (ins/del log file creation):

III.01.10. Proximity analysis distribution plot file writing

PXPL If this keyword is present, the program will write an ASCII file containing all the calculated proximity distributions to the .pxp file. The distributions written on the .pxp file can be plotted with the companion (Fortran-77) program pxps.f. pxps.f produces Postscipt files of plots, 8 plots per page. It can put two distributions on the same plot, but both have to use the same dependent variable. The right-side vertical axis will belong to the second plot.

A typical run of pxps.f is as follows (user input is displayed in bold face): 

% pxps
 
Proximity plot generator Version:/07/22/98
Data file name root=wt1
 File wt1.pxp opened                             <-  the file wt1.pxp
                                                     was generated by
                                                     the PXPL key
Use @ for: @f=Greek fnt, @n=Normal fnt, @+=superscript
 @-=subscript, @m=newline, @#n#=font size n, @@=@symbol
Title of graph page:
GnRH Wild Type, Extended Conformation
Do you want grid lines drawn (y/n)? n   <- If yes, only one plot
                                                  can be plotted
Plot types:
 1  Solute-solvent radial distribution functions       <- from key PXGR
 2  Solute-solvent pair energy distribution functions  <- from key PXBE
 3  Mean solvent orientation as a function of R        <- from key PXDP
 4  Distribution of first shell solvent orientations   <- from key PXDP
 5  Solvent-solvent radial distribution functions      <- from key PXWW
Plot type number=1
Data colums: 1 R        2 gpx(R)   3 Kpx(R)   4 gt(R)    5 Kt(R)
Data column number=2
Do you want a second function also plotted (y/n)? y
Plot types:
 1  Solute-solvent radial distribution functions
 2  Solute-solvent pair energy distribution functions
 3  Mean solvent orientation as a function of R
 4  Distribution of first shell solvent orientations
 5  Solvent-solvent radial distribution functions
Plot type number=1
Data colums: 1 R        2 gpx(R)   3 Kpx(R)   4 gt(R)    5 Kt(R)
Data column number=3
Number of graphs to print (0: will ask for range)=0
First and last distribution function=1,18
Data set not found (skipped): gkr is=   1    <- empty proximity
regions
Data set not found (skipped): gkr is=  10
 File wt1_2_9_gpx_Kpx.ps opened
 File wt1_11_18_gpx_Kpx.ps opened
More plots (y/n)? n   <- allows to plot other distributions
%

The format of the .pxp file is as follows:

III.01.11. Proximity information file specification

PXWR Key1 [Key2 Data]

Prerequisite: PXCR,

Key1 (file format - if any):

The proximity information file has extension .pxi.

III.01.12. Converting/shifting a trajectory

WTRA Key1 Key2 [Data]

Prerequisite: TRAJ

Key1 (transformation type):

This key allows the conversion of the input trajectory to a new format and gives the option of centering it in two different ways.

III.01.13. Contact elimination

STIR [Data]

Prerequisite: CNFG

Data: eijmin, disp, dispa

This option is useful for 'cleanng up' a starting configuration generated with the RANC option.

III.01.14. Local disk directory

LPAT [Data]

Data: scratchpath

III.02. System descriptor keywords

  1. NSLV: Number of solvent molecules
  2. TEMP: Simulation temperature
  3. TITL: Description
  4. PBCN: Periodic boundary conditions
  5. GCEN: Grand-canonical (T,V,mu) ensemble
  6. STPX: Stop GCE run at nx
  7. BTUN: Tune GCE B parameter
  8. LIMG: Limit the cavity grid
  9. IBEN: Isothermal isobaric (T,P,N) ensemble
  10. LCMP: Virial biasing factor
  11. CNFG: Initial configuration
  12. SCRM: Scramble torsion angles
  13. SCAL: Simulation cell rescaling
  14. OVST: Overlap of the second solute copy with the first
  15. SPRD: Spread out solute molecules
  16. PBGR: Solute reset with grids
  17. LIGA: Multi-copy ligand
  18. STUN: Stepsize tuning
  19. TAC0: Stepsize accumulator reset
  20. SANN: Simulated annealing
  21. SACP: Simulated annealing of the chemical potential

III.02.01. System size

NSLV Data

Can not be followed by: HRDW
Potential successors: CNFG, RCKP

Data: numsol, nsolfix

nsolfix > 0 requires that the solute does not undergo translation/rotation (cedslt and rtxslt should both be 0.0 - see the STEP key).

If the input has no NSLV key then the number of solvent molecules will be established from the .crd or .hst files.

III.02.02. Simulation temperature

TEMP Data

Prerequisite of: RUNS, SANN

Data: T

III.02.03. Title

TITL Data

Potential successor: RCKP

Data: description of the calculation.

The program is prepared to read at most two title lines.

III.02.04. Periodic boundary conditions and cell dimension.

PBCN Key1 [Key2 Key3] Data [Fdata]

Prerequisites: SUVC, SVVC
Potential successor: RCKP
Prerequisite of: SLTA, GCEN, MOND, LIMG, VORO, FREE, PRPB, GENS

Key1 (cell shape):

Note: Input boundary conditions work only if the image cell can be found by comparing the distances of a point from the various cell centers, the smallest distance will belong tot he wanted cell number. It works when the periodic system can be described by an orthogonal crystal coordinate system. For other cases, more complicated calculations are required - to start with, a transformation to non-orthogonal coordinate system - but it is not implemented yet.

For pure liquids, FCC is the best since this maximizes the distance of a molecule from its image. On the other hand, RECT is conceptually the simplest. Furthermore, for RECT the dimensions of the cell can be different along the three axes, allowing for various types of crystals. Solutes that are far from spherical can also be wrapped around by waters more efficiently. (In this case, however, the solute rotation should be prevented - see key STEP.) HEXG have been developed for long solutes or solutes representing an infinite polymer chain, to maximize the distance of the chain from its image. SPHR simply encloses the system into a sphere of inputted radius. Note that the FCC calculations are only partially vectorized.

For the rectangular, face-centered cubic, hexagonal prism and sphere boundary conditions there are special algorithms in the program, one for each. When boundary conditions are inputted, the cell number of the periodic image cell where a molecule is located is determined simply by comparing the distance of the molecule from the various cell centers read in.

The unit cell parameters are determined from the bulk density or molar volume of the solution, Nmol, and the choice of boundary conditions. If the density of the system is d g/ml , then

d (g/ml) = weight in g of Nmol / volume of Nmol in ml =

(Nmol/L)*Mol.wt / (Vol/1024)

where L is Avogadro's number. The cell parameter(s) can be obtained from the resulting Vol value.

III.02.05. Grand-canonical ensemble selection

GCEN[ Key1 Key2 Key3 Key4 Key5 Data ] Fdata

Prequisite: PBCN
Can not be followed by: LIMG
Prerequisite of: PRTG, PRFI, BTUN, STVG, STPX, SACP
Related keys: IDLG, IDAG
Variable to set to 'T' in the preprocessor: CB

Grand-canonical ensemble simulation using the cavity-biased method with grid-insertion or the original, random insertion method.

Key1 (GCE insertion strategy):

Other aspects of a GCE simulation are controlled by the keys PRFI, IDLG, IDAG, STPX, LIMG, and BTUN.

III.02.06. Stop GCE run at nx

STPX Data

Prerequisite: GCEN,

Data nslvxp

When this key is present the run will stop if the number of solvent molecules reached nslvxp.

III.02.07. Tune the GCE B parameter

BTUN Key1 [Data]

Prerequisite: GCEN, MOND or BLKW

Key1 (tuning algortithm)

The program will periodically change the GCE B parameter to achieve the targeted density in the 'bulk' region - defined as the outside region by the MOND key.

Note that the keys NONE and AVRG are only meaningful after having restored a checkpoint file (key RCKP).

III.02.08. Limit the cavity grid

LIMG [Key1,Data]

Prerequisite: PBCN
Potential successor: GCEN

Key1 (used only for GCE simulations to specify treatment of boundary crossings):

III.02.09. Isothermal isobaric ensemble simulation (TPN)

IBEN Key1, Key2, Data

Can not be followed by: VCHA
Potential successors: CNFG, VCHA
Prerequisite of: LCMP
Related key: LCMP.

Variable to set to 'T' in the preprocessor: IB

Key1 (volume sampling technique):

Note, that putting the key IBEN after CNFG will result in ignoring the cell size on the .crd file and the vaulues inputted with the PBCN key will be used instead.

III.02.10. Virial-biasig factor

LCMP Data

Prerequisite: IBEN

Data:

III.02.11. Initial configuration

CNFG Key1 [ Key2 Key3 Data ] [Fdata]

Prerequisites: SLTA, FILE
Can not be followed by: FREE, TORD, CMPC, NSLV, IBEN, VCHA, PXYZ
Potential successors: RCKP, SPST
Prerequisite of: RUNS, SCAL, PRTG, REGE, SPRD, VORO, STIR, SCRM, TORT, LPST, GENT
Related keys: GENT, SPRD, SCRM, STIR.

Key1 (data source):

III.02.12. Scramble torsion angles

SCRM

Prerequisite: CNFG

When this key is present, all torsion angles are given a random value and the solute is regenerated in this random conformation.

III.02.12. Simulation cell rescaling

SCAL Data

Prerequisite: CNFG

Data: e1, e2, e3

III.02.13. Overlap of the second solute copy with the first

OVST Data

Data: i1, i2, i3

III.02.14. Spread out the solute molecules for display

SPRD [Data]

Prerequisite: CNFG

Data:

This key generates a PDB file containing the solute molecules spread out in the x-y plane so that they don't overlap.

III.02.15. Force solute reset with grids

PBGR

When this key is present, the grid calculations will reset each solute atom to the periodic cell before calculating the cover list. This is a safety option, since it is automatically set when the inital structure contains solute atoms outside the simulation cell.

III.02.16. Multi-copy ligand

LIGA Data

Data:

III.02.17. Stepsize tuning

STUN Key1, Key2, Key3, Data

Prerequisite: MVRT
Can not be followed by: MVRT
Prerequisite of: TAC0

Related key: TAC0.

Key1 (move type):

Note that the solute sampling should be set to alternating translation and rotation with the MVRT key. For cloned solute molecules (see the key CLON) the tuning will generate identical stepsizes for each instance of a torsion angle.

III.02.19. Stepsize accumulator reset

TAC0 Key1, Data

Prerequisite: STUN

Key1 (move type): same as Key1 of STUN

Data: nmc_zeroacc

III.02.19. Simulated annealing

SANN Key1, Data

Prerequisite: TEMP

Key1 (annealing schedule):

III.02.20. Simulated annealing of the chemical potential

SACP Key1, Data

Prerequisite: GCEN

Key1 (iteration length)

Data: bincr, nmcsacpstep, fract_lim

When this option is used and the trajectory writing (key TRAJ) is not used, at each B parameter change the solvent coordinates will be saved on the history file as PDB MODELs.

III.3. Free-energy calculation keywords

  1. FREE: Free energy options
  2. RAUS: AUS parameter change
  3. WMAT: PMF matching
  4. WPLT: AUS iteration plot
  5. TIQU: TI quadrature
  6. OVRA: Overlap ratio
  7. FETK: Temperature scaling for TI evaluation
III.03.01. Free energy options

FREE Key1, [Key2-4 Data Fdata]

Prerequisite: PBCN
Potential successors: CNFG, SLTA, TORD, RCKP, PARD
Related keys: TIQU, WMAT, WPLT.

Key1 (free-energy method choice):

III.03.02. Adaptive Umbrella Sampling parameter change

RAUS Key1, Fdata

Key1 (initialization option):

Note that cplmin and cplmax have to stay at their original value.

III.03.03. PMF matching

WMAT Key1, Key2, Key3, Key4, Data, Fdata

This function joins the adaptive umbrella sampling estimates of the potential of mean force obtained from the checkpoint files of calculations on overlapping coupling parameter ranges (windows).

Key1 (checkpoint list specification):

III.03.04. AUS iteration plot

WPLT Key1, Data, [Data]

Key1 (scale for data points):

This key will produce a family of curves showing the PMF calculated from each iteration as lines, with circles representing the extent of sampling at each coupling parameter grid. The matched PMF is shown with large filled circles.

III.03.05. TI quadrature

TIQU Key1,[Key2] Data, [Fdata]

Key1 (data source specification):

This function performs the numerical (Gaussian) quadrature on thermodynamic integration runs. It prepares a plot of the polynomial fitting the quadrature ponts and, optionally, performs other numerical integrations.

III.03.06. Overlap ratio evaluation

OVRA Data

This function obtains the free energy from two overlap ratio method runs.

Data: numrun1, numrun2

III.03.07. Temperature scaling for TI

Prerequisite: TEMP

FETK Data

This option scales the temperature by dividing it with cplparkexp to represent a run where the total system is made to turn into ideal gas

Data: kexp, nquad, iquad

III.4. Molecule descriptor keywords
  1. SLTA: Solute description
  2. SLVA: Solvent description
  3. CLON: Solute cloning
  4. MAKB: Make bonds
  5. BRKB: Break bonds
  6. BDHH: Forcing H-H bonds
  7. BDRH: Forcing bonds of an H to a different group
  8. MOLD: Define solute molecule limits
  9. RDBD: Bond list read
  10. TORD: Torsion angle definitions
  11. CSEG: CHARMM segment id handling
III.04.01. Solute description

SLTA Key1 Key2 Key3, Data, [Data, Fdata]

Prerequisites: FILE, PBCN, SVPT, SUPT
Can not be followed by: SLVA, PMOD, SLVA, PFRD, FREE, PROT, p FLXR, CNST, CHRG, CLON, MAKB, BRKB, RFSL, MIXR, SETC, MODA, MOLD, STEP, RDBD
Potential successors: RCKP, PRMF
Prerequisite of: CNFG, RUNS, TORD, WCNF, FCGA, PARD, FCGD, PMFM, AROM, ENGL, STVG, CMPC, BNDL, HBMO, SVDP, GENV, PXCR, SCAN, HBBR
Related keys: CLON, MAKB, BRKB, MOLD, RDBD.

Key1 (possible solute periodicity):

The program is currently prepared to handle atoms up to atomic number 86. The atomic number for lone-pair is 89 (in case of EPEN, it is disregarded) and atomic number 90 is an electron. Atomic numbers 91-93 have been reserved for the united atoms CH1-CH3, respectively. The various atomtypes for the different potential libraries implemented and the actual functional form of the potentials are specified in a separate document.

Important: For EPEN, the different types of centers have to be grouped together. Lone-pairs should follow nuclei and electrons should follow nuclei and lone-pairs.

When the coupling parameter is used to determine conformational transition (FREE PMF1), three sets of solute coordinates are used: the two conformations Ro and R1 and their combination CP*R1+(1-CP)*Ro. The solute description thus has to contain the coordinates, charges of the two conformations as well as room for the mixture. In this third set, only the atomtypes have to be given. When CNFG PMFN, however, only the Ro conformation is required here since the R1 conformation will be given at a different part of the input. It is important that the different descriptions of the conformation give the solute atoms in identical order. Solute atom type change during the transition can not involve EPEN.

III.04.02. Solvent description

SLVA [Data Fdata]

Prerequisite: SVPT
Potential successors: SLTA, RCKP
Prerequisite of: SVDP

Data: nslv, islvrep, label, islv4, filename

Pseudo-atoms are to be placed after the regular atoms.

Note: This key is ignored if the solvent is not GENL type.

Note: the dipole of the solvent molecule should be parallel to the X axis for the purpose of the near-neighbor solvent dipole correlation function computation (DSTC ALL).

III.04.03. Solute molecule cloning

CLON Data

Potential successor: SLTA

Data: nclone, iclonef(1), iclonel(1), ncopy(1), ... iclonef(nclone), iclonel(nclone), ncopy(nclone)

III.04.04. Making additional bonds

MAKB Data

Potential successor: SLTA

Data: nmake, imake1(1), imake2(1), ..., imake1(nmake), imake2(nmake)

Note that this key is not implemented for the case when cloning (CLON) is in effect.

III.04.05. Breaking bonds

BRKB Data

Potential successor: SLTA

Data: nbreak, ibreak1(1), ibreak2(1), ... ibreak1(nbreak), ibreak2(nbreak)

Note that this key is not implemented for the case when cloning (CLON) is in effect.

III.04.06. Forcing H-H bonds

BDHH

By default, the topology search excludes bonds between hydrogens. If this key is present, these bonds will be kept.

III.04.07. Forcing bonds of an H to a different group

BDRH

By default, the topology search excludes bonds of hydrogens to heavy atoms of a different group/residue. If this key is present, these bonds will be kept.

III.04.08. Defining solute molecule limits

MOLD Data

Potential successor: SLTA

Data: molslt, ilastm(1), ilastm(2), ...

This key overrides the molecular topology for defining the molecules within the solute (or just a part of it). This is useful when the molecule is distored and too many incorrect bonds are found or too many bonds are missing. When applied to cloned (CLON) solute atoms, the definitions have to be repeated for all clones. The solute atoms not covered by this operation have to be placed contiguously at the end of the solute definition read with the (SLTA) key.

III.04.09. Bond list input instead of calculating

RDBD Key1 [Data]

When this key is present, the solute bond list is read from a file.

Key1 (bond list type):

Data file name

III.04.10. Torsion angle definitions

TORD Key1 Key2 Data Fdata

Prerequisite: SLTA
Can not be followed by: PRMF, SETC, LOOP, FREE, PARD, SKWT
Potential successors: CNFG, SPST
Prerequisite of: REGE, TORT, LPST, GENT, TAND
Related keys: LOOP, PART, TORT, LPST.

Key1 (torsion definiton source):

III.04.11 CHARMM segment id handling

CSEG

When this key is present, the last 4 characters of spgroup (read with SLTA **** MMC) will be used for the segment id when CHARMM CRD format is written with the WCNF CHRM or DENF CHRM keys.

III.5. Potential descriptor keywords

  1. SUPT: Solute potential type
  2. SVPT: Solvent potential type
  3. MIXR: LJ parameter mixing rule
  4. PRMF: Torsion potential specification
  5. PMOD: Potential library modifications
  6. MODA: Element modification
  7. PFRD: Potential library input
  8. SUVC: Solute-solvent cutoff
  9. SUUC: Solute-solute cutoff
  10. SVVC: Solvent-solvent cutoff
  11. INCT: Solvent inner cutoff
  12. PROT: Minimum repulsion for general solvent
  13. RFCR: Reaction-field correction
  14. SETC: Set various constants
  15. CHRG: Solute charge manipulations
  16. SPST: Field-dependent potentials
  17. SPPS: Solute molecule-dependent potentials
  18. AROM: List of aromatic carbons
  19. SVIN: Solvation parameter input
  20. ENHB: Hydrogen-bond potential
  21. HBMO: Non-default hydrogen-bond
  22. FLXR: Residues to sample
  23. CNST: Constraint potential
  24. VVNE: No solvent-solvent elecrostatic in Solvation parameter input
  25. IGTT: Solute-solute potential is set to zero
III.05.01. Solute-solvent potential

SUPT : Key1 [Key2]

Potential successor: RCKP
Prerequisite of: SLTA, PROT, RFSL

Related keys: PMOD, MODA, PFRD, MIXR.

Key1 (solute PF library source):

III.05.02. Solvent-solvent potential

SVPT Key1 Key2 [Data]

Potential successor: RCKP
Prerequisite of: SLTA, SLVA

Key1 (solvent PF type):

Key2=GENL also requires the key SLVA READ

III.05.03. Lennard-Jones parameter mixing rules

MIXR Key1

Potential successors: SLTA, RCKP

Key1:

ELJ(r)=4*EPS[(SIG/r)12-(SIG/r)6] = c12/r12-c6/r6. Relevant only for AMBER, CHARMM, GROMOS and OPLS potentials.

III.05.04. Torsion potential parameter input

PRMF Key1 Data

Can not be followed by: SLTA
Potential successor: TORD
Related key: PFRD

Key1: Same as Key1 of SUPT - the name of the potential library.
Data: Name of the file containing the parameters.

The input syntax for the parameter files is defined by the respective libraries. At present, only AMBER, CHARMM, OPLS (in CHARMM format) and GROMOS libraries can be read. Furthermore, for GROMOS input the file name has to be of the form <name>bon.itp and the file <name>.rpt also has to be present.

III.05.05. Potential (re)definition

PMOD Key1 Data Fdata

Potential successor: SLTA

Key1:

III.05.06. Atom (element) descriptor modifications

MODA Data Fdata

Potential successor: SLTA

Data: nmod

Using this key the built in constants for each chemical element can be modified or new ones can be defined. Fields left blank or given as zero will retain their original value

III.05.07. Potential library input

PFRD Key1 Data

Potential successor: SLTA
Related key: PRMF

Key1 (library type):

III.05.08. Solute-solvent cutoff

SUVC Key1 [Data]

Potential successor: RCKP
Prerequisite of: PBCN

Key1 (cutoff type):

For the particular case of a single solute molecule in nmolec-1 solvent it is advisable to use the minimum-image cutoff for solute-solvent interaction to avoid the single discontinuity at the boundary of the cutoff sphere. Multimolecular solutes or for solutes with torsional degree of freedom require group-based cutoffs.

In the minimum-image cutoff for each pair of molecules only the interaction between the nearest images are taken into account, irrespective of the intermolecular separation. This is equivalent to the use of a cutoff cube, prism, etc., identical in size, shape and orientation with the periodic cell.

The cutoff radius may be based on the centers of solute groups (see SPGC or MIGC above) - useful for large solutes. Also, interionic PMF calculations can use the iso-energy cutoff introduced by Mezei.

III.05.09. Solute-solute cutoff

SUUC Key1

Potential successor: RCKP
Prerequisite of: PARD, PART, ENGL

Key1 (cutoff type):

Note, that the **GC options will update the energies of all atoms that are in a residue when a torsion move moves the group center atom. To reduce the extra calculations, a) you should not use too large groups when this option is active and b) whenever possible, group torsions so that they don't leave too many group members unchanged.

III.05.10. Solvent-solvent cutoff

SVVC Key1 [Data]:

Potential successors: RCKP, PART
Prerequisite of: PBCN

Key1 (cutoff type):

III.05.11. Inner cutoff on the solvent

INCT Key1 [Data]

Potential successor: RCKP

Key1:

III.05.12. Minimum repulsion.

PROT Data

Prerequisite: SUPT
Potential successor: SLTA

Data: c12prot

III.05.13. Reaction-field correction

RFCR [Data]

Variable to set to 'T' in the preprocessor: RF

Data: epsrf

The use of the reaction field correction requires the use of identical cutoffs for the solvent-solvent (SVVC), solute-solvent (SVVC) and solute-solute (SUUC) interactions. Also, the cutoffs are to be group based and each group has to be electrically neutral.

III.05.14. Set various constants

SETC Key1 Data

Potential successors: SLTA, TORD

III.05.15. Solute charge manipulations

CHRG Key1

Potential successor: SLTA

Key1:

III.05.16. Field-dependent solute potential

SPST Key1 Data

Prerequisite: FLXR

Can not be followed by: AROM, SVIN, ENHB, PARD, TORD, CNFG, HBMO
Related keys: AROM, ENHB
Variable to set to 'T' in the preprocessor: GM

Key1:

For Key1=EMAP, after establishing the grid constants the program reads the energy maps from the files prot.pi, lp.map, pa.map, pc.map, pd.map, ph.map, pn.map, po.map, pp.map, ps.map. for the electrostatic energy, van der Waals energy for C, A, N, O, P, S and H atoms, respectively (A is aromatic C) and, whenever needed, the file lp.map containing the contribution of the (Autodock 4) solvation term and the file ov.map containing the O-H VdW term. These are binary files containig the 3-dimensional array e(ningridmap,ningridmap,ningridmap) in a single record. The electrostatic map (read first) also contains a second record containig title, nx, ny, nz, ep, es, micro, stern, radd, ho, hi2, akp, aks where map(hi2*[(ningridmap-1)/2+1], hi2*[(ningridmap-1)/2+1], hi2*[(ningridmap-1)/2+1]) is the energy at <0.0,0.0,0.0>; map(1,1,1) is the energy at <-hi2*(ningridmap-1)/2, -hi2*(ningridmap-1)/2, -hi2*(ningridmap-1)/2>

The solvation term Esolv(i,j) contribution from atoms i and j is of the form
Esolv(i,j) = (ViSj + VjSi)*exp{r(i,j)/gaussdist2}
Si=Ai + qk_par*abs(q(i))
where gaussdist2, and qk_par are constants Vi, and Ai are atom type-dependent constants and r(i,j) is the distance between atoms i and j and q(i) is the partial charge on atom i. All constants have built-in values that can be changed with the SVIN key (the built-in values serve as defaults).

Furthermore, if the map file hb.map is present then the program will treat all fixed oxygens as hydrogen bond donors and include a 12-10 hydrogen-bond potential. Hydrogen-bond acceptors can not be included in such a map since the H-X bond's orientation is conformation-dependent. It can be calculated explicitly during the simulation with the ENHB key

Note that on restart the energy maps have to be read again the same way as for startup.

For Key1=CMAP or Key1=CMPW, the program needs the description of the atoms contributing to the grid potentials. These are assumed to be in a file <jobname>_2.slt that is in the same format as the records read by the key (SLTA, but only the fields tslt(i), (cslt(j,i),j=1,3), qslt(i), labslt2(i), indxrdf(i) are read. There should be a record for each atom that is not explicitly represented.

III.05.17. Molecule-dependent solute potential

SPPS Key1 Data

Key1:

This key allowed the use of potential terms depending on the field or on the molecule's center, orientation, etc. To use this option, the user also have to provide the corresponding code to the body of the subroutine esltmolec.

III.05.18. List of aromatic carbons

AROM Key1 Data

Prerequisite: SLTA
Potential successor: SPST
Prerequisite of: SVIN
Related key: SPST.

Key1: Solute potential type see key SUPT.

This information is needed for deciding what type of energy maps to use with these atoms - see, e.g., SPST EMAP.

III.05.19. Desolvation parameter input

SVIN Data

Prerequisite: AROM
Potential successor: SPST
Related keys: SPST, AROM.

Key1:

III.05.20. Hydrogen bond potential

ENHB Key1 Key2 Data

Prerequisite: SLTA
Potential successor: SPST
Prerequisite of: HBMO
Related keys: SPST, HBMO, SSVIN.

Key1 (solute-solute HB potential type):

When this key is used a hydrogen-bond term will be used between polar hydrogens and oxygens as specified by Key1 and Key2, subject to the filters specified by rmaxhb, rminhb, qmaxhbdon, qminhbacc. This term may be used between non-bonded solute atoms and - if used - the grid-based energy term (SPST EMAP).

When the grid energy is used, and there are hydrogen bonds between the explicitly represented atoms and the atoms contributing to the grid then the atoms (that are otherwise fixed) involved in hydrogen bonds have to be specified. The hydrogen-bond donors and acceptors will be extracted from the non-felxible part of the solute as specified by the FLXR key.

III.05.21. Non-default hydrogen bond

HBMO Key1 Data

Prerequisites: ENHB, SLTA
Potential successor: SPST
Related key: SPST, ENHB.

Key1 (type of action):

III.05.22. Residues to sample

Potential successor: SLTA

FLXR Key1 Data

Related key: TORD

Prerequisite of: SPST

Key1 (input syntax):

For calculations using the energy map this key defines the residues to be sampled. All other residues will contribute to the energy maps.

III.05.23. Constraint potential

Potential successor: SLTA

CNST Key1 Key2 Data

Key1 (new or modify)

III.05.24. No solvent-solvent electrostatic interaction

Potential successor: SLTA

VVNE

When this key is present solvent-solvent interactions will mot include the electrostatic terms.

III.05.25. Ignore solute-solute energy

IGTT

When this key is present all solute-solute interactions are set to zero.

III.6. Sampling-related keywords

  1. MOVE: Move selection strategy
  2. SAMP: Sampling technique
  3. NFBU: Solute force bias
  4. STEP: Stepsizes
  5. STPS: Stepsize scaling
  6. PARD: Solute molecule displacement
  7. SWAP: Solute molecule swap
  8. MV2S: Correlated 2-solute move
  9. PART: Solute torsion
  10. MVRT: Alternating translation and rotation
  11. PRFI: Insertion/deletion pref. sampl.
  12. SEED: Random number seed input
  13. NONB: Non-Boltzmann solute samplings
  14. STSC: Temporary stepsize scaling
  15. LOOP: Loop torsion moves
  16. IGJA: Ignore Jacobian
  17. SKWT: Skewed torsion sampling
  18. IGSV: Non-interactive solvent
  19. RNDG: Random number generator type
III.6.01. Molecule selection strategy specification

MOVE Key1 [Data] [Fdata]

Potential successor: RCKP
Prerequisite of: PRFI

Key1 (molecule selection strategy):

III.6.02. Sampling technique selection

SAMP Key1 [Data Fdata]

Variable to set to 'T' in the preprocessor for Key1=FBPR or FBSC: FR

Key1 (possible displacement biasing strategy):

III.6.03. Solute sampling technique selection

NFBU

If present, the Metropolis method will be used for for moves translating/rotating the whole solute (se key STEP), even when the key SAMP specifies force-biased sampling. In this case leaving the C@TS lines comments saves some computer time.

III.6.04. Solute and solvent stepsizes

STEP Data

Potential successor: SLTA
Prerequisite of: RUNS

Data: cedslt, rtxslt, cedslv, rtxslv, nsltfreq

The above variables are used to perturb a chosen molecule, solute or any of the solvent during the Metropolis walk. As a simple rule of thumb, these variables are chosen by trial and error to obtain 30-40% acceptance rate. For water the recommended values are: 0.30 Å for the COM of the solvent and 30 degrees for the rotation angle for standard Metropolis and 0.55 Å and 40 degrees for Force Bias. The stepsize of the solute should be adjusted to provide similar acceptance rate as the solvent. In general, the larger the solute the smaller stepsizes are to be used.

III.6.05. Stepsize scaling

STPS Data

Data: notmovcyc

III.6.06. Partial displacement of the solute

PARD Key1, Key2, Data

Prerequisites: SLTA, SUUC
Can not be followed by: FREE
Potential successors: TORD, SPST
Prerequisite of: SWAP, MV2S
Related key: SWAP.

Key1 (rotation stepsize strategy):

III.6.07. Solute molecule swap

SWAP Data

Prerequisite: PARD

Data: wswap

III.6.08. Correlated 2-solute molecule move

MV2S Data

Prerequisite: PARD

Data: wmv2s, r2scut, rtxcslt, iaxis2s

Key1 (stepsize selection strategy): See Key1 of key PART

III.6.09. Partial torsion of the solute

PART Key1, Key2, Data

Prerequisites: SLTA, SUUC
Can not be followed by: SVVC
Prerequisite of: TAUC
Related keys: LOOP, TORD.

Key1 (stepsize selection strategy):

III.6.10. Alternating translation and rotation

MVRT Key1

Potential successor: STUN
Prerequisite of: STUN

Key1 (solute and/or solvent moves)

III.6.11. Insertion/deletion preferential sampling

PRFI Key1 Key2 [Data] [Fdata]

Prerequisites: GCEN, MOVE

The syntax is the same as for the preferential sampling for the solvent moves. Thus, except for not allowing RAND, CYCL, SHCY, CYCI for Key1, the input is the same as for the key MOVE.

III.6.12. Random number seed initialization

SEED Data

Data: ix1, ix2, ix3, ix4, ixscr

The random number generator will be reinitialized with these seeds.

III.6.13. Non-Boltzmann solute samplings

NONB Key1 Data

Key1 (sampling type):

III.6.14. Temporary stepsize scaling

STSC Key1 Data

Key1 (move type to scale):

The scaled stepsizes will only be used for this run - the next calculation (initiated by a new RUNS key) will use the originally inputted values.

III.06.15. Loop torsion moves

LOOP Key1, Key2, [ Data]

Potential successor: TORD
Prerequisite of: LPST

Variable to set to 'T' in the preprocessor: LO

Related key: IGJA

When this key is present the program will analyze the torsion angle list and determine which torsions can be the first torsion of a loop move. Loop moves change six consecutive torsions in such a way that the bond of the first and last torsion remain fixed, thereby limiting the change in the molecular geometry due to the change in the torsionangle(s).

The formulae giving the loop-closing torsion angles were obtained based on a distance geometry approach. The current version is limited to loops where the torsions are on consecutive bonds (-------), or on consecutive bonds broken by a rigid bond (-=-----, --=----, ---=---, ----=--, -----=-) or on a peptide backbone with rigid peptide bonds (-[=--]3, -[-=-]3).
Key1 (solution selection algorithm):

III.06.16. Ignore Jacobian

IGJA

Related key: LOOP,

When this key is present, the acceptance probability of a local move (see key LOOP) will ignore the Jacobian. This will allow skipping the calculation of the loop-closure problem for the reverse move and thereby speeding up the calculation. The speedup, however, comes at the expense of the loss of Boltzmann sampling and thus this option is to be used for search or minimization runs.

III.06.17. Skewed Torsion sampling

SKWT Key1 Data

Related key: TORD,

Key1 (type of skewing)

When this key is present then the formatted data after the key TORD will specify a target angle for some of the torsions.

III.06.18. Non-intractive solvent

IGSV

When this key is present then solvent-solvent interactions will be ignored.

III.06.19. Random number generator information

RNDG Key1 Data

Key1 (generator type)

Data for Key1=LCG : lcg_fac, lcg_add, lcg_exp, seed

Data for Key1=MRTW: seed

Data for Key1=INPT or Key1=OUTP: rndgfile, nrand

III.7. Calculation flow keywords

  1. RUNS: Run MC simulation
  2. SCAN: Scan trajectory for proximity analysis
  3. PXAN: Proximity analysis frequencies
  4. RCKP: Checkpoint file read
  5. STOP: Stop the run

III.07.01. Execution of a MC simulation

RUNS Data

Prerequisites: SLTA, CNFG, TEMP, FILE, STEP

Data: nmcmax, nmcrep, nplt, nmcrec, ncntin, nmcadp

After reading a checkpoint file (see key RCKP) the rest of the data is ignored. III.07.02. Proximity analysis from history

SCAN [Key1] Data

Prerequisite: SLTA

Key1 (source of data):

For Key1=TRAJ or TRNC:
Data: nmcpxmax, navgpx, nranpx, nmcpxdsc, npxres, npxcntin, lumppr, nsymdim, ifrsdim For Key1=CONF:
Data: nranpx, lumppr, nsymdim, ifrsdim - see above.

Structures in history files in CHARMM (TRAJ CHRM) or AMBER (TRAJ AMBR) trajectory syntax are labeled 1,2,... while structures in all other syntaxes are labeled by the MC stepnumber they were saved at.

III.07.03. Proximity analysis frequency parameters (concurrent to MC run)

PXAN Data

Can not be followed by: PXGR, PXDP, PXBE, PXWW, TAND, ATFR, FLDG, RTIM, DIFC, VORO, CHKP
Related keys: PXLM, PCPA, PXYZ.

Data: navgpx, nranpx, nmcpxdsc, npxres, npxcntin, lumppr, nsymdim, ifrsdim

These variables are the same as defined for the keyword SCAN.

Using this keyword the proximity analysis will be performed concurrently to the MC simulation (eliminating the need of saving the history file).

III.07.04. Continue from a checkpoint file

RCKP [Key1 Key2 Data]

Prerequisite: FILE
Can not be followed by: TITL, PBCN, SUPT, SVPT, TRAJ, INCT, MIXR, SUVC, SVVC, SUUC, NSLV, MOVE, FREE, SLTA, SLVA, CNFG
Potential successor: VOLE
Related keys: CHKP, WCKP, RMCK, SCKP.

The program will read the content of the current checkpoint file, .ckp and restore from there the state of the variables, including the file unit numbers and the run number numrun (!), allowing the continuation of an earlier run as if it never stopped. If you want to reinput some of the data used by the run the run may abort as certain inputs are forbidden in different context. In this case either use the RUNS key with zero steps to run or the FILE key. Use with caution, though!

Key1 (checkpoint file checking procedure):

If the checkpoint file was found too short, the program will assume that the last configuration was still properly read (since that is the first item in the file) and write it on a new file <jobname>.99.crd in binary form (to be read with CNFG READ BNRY NOFX 99).

III.07.05. Termination of the calculation

STOP [Key1, Key2, Data]

Key1 (possible exit self test request):

III.8. Output related keywords
  1. OUTP: Print various system charateristics
  2. PRNT: Input echo level
  3. DSTP: Bulk solvent distribution print
  4. PRAC: Acceptance rate print detail
  5. BNDL: Listing of bonds found
  6. PRCO: Compilation option listing
  7. PLBP: Potential parameter echo
  8. KMNP: PDB file property threshold
  9. EPLT: Switch between Cv and Utot in convergence plot
  10. PRPL: Print plotted data
  11. NMVP: 'Not-moved' threshold to print solvents
  12. PXPA: Empty proximity region print
  13. PXPR: Proximity analysis printing options
  14. PDBT: PDB atomname convention
  15. PLCV: Postscript convergence plots
III.08.01. Print various system characteristics

OUTP Key1

Key1 (type of output):

III.08.02. Input echo level selection

PRNT Key1

Key1 (level of echoing input):

III.08.03. Bulk solvent distribution function print level specification

DSTP Key1

Related key: DSTC.

Key1:

III.08.04. Acceptance rate listing detail

PRAC Key1[,Key2]

Key1:

III.08.05. Bond list print

BNDL [Key1 Data]

Prerequisite: SLTA

Key1:

Note: atoms with one neighbor will not be listed separately, only as neighbors.

III.08.06. Compilation option print

PRCO [Key1 Data]

Key1:

The program prints all compilation options and array sizes as well as version dates used either by the currently running ptogram or the program that created the checkpoint file with jobname.numrunr.ckp. If requested by Key1, a file pre_jobname.f will be generated that can be inserted into the preprocessor pre.f to compile the program the same way. Furthermore, if the file pre.f is present in the current directory then the file pre_jobname.f will contain the full source code of the modified preprocessor.

III.08.07. Potential parameter listing

PLBP Key1 [Key2]

Key1:

III.08.8. Coordination number threshold for property output

KMNP [Data]

Data: rkpdbmin

III.08.9. Switch between Cv and Utot in the convergence plot

EPLT Key1

Key1:

III.08.10. Listing of function values for plots

PRPL Key1

Key1:

III.08.11. 'Not-moved' threshold

NMVP Data

Data: notmoved

III.08.12. Empty proximity region print

PXPA Key1

Key1:

III.08.13. Proximity analysis print levels

PXPR [Key1 Key2 Key3 Key4]

Key1 (distribution function print level):

The proximity analysis results consist of two parts: distribution functions (and their characteristics) and tables of calculated avareges. The distribution printing is governed by Key1 and Key 2. The program can print the the averages arranged in three different tables: III.08.14. PDB atomname convention

PDBT Key1

Key1 (name conversion):

When printing a structure in PDB format (e.g., keys WCNF, GENS, FILT) this option allows the choice of PDB conventions.

III.08.15. Postscript convergence plots

PLCV Key1, Data

Key1 (Property):

This key can be used to generate a Postscript plot of the block averages and cumulative averages of the energy, volume, number of molecules or the free energy calculated with thermodynamic integration. The name of the plot file will be of the form <jobname>.<ext>.<runnumber<.ps where <ext> is eng, vol, nml, bpr, fen for Key1= ENGR, VOLU, NUMB, BPAR, FETI, resp. If the .ps file already is present, a new file will be opened vith an added (or incremented) version number.

III.9. General analysis keywords

  1. DSTC: Bulk solvent distribution calculation
  2. FCGA: Functional group analysis
  3. ATFR: Calculation of mean force on solute atom
  4. FLDG: Field gradient calculation on solvents
  5. SENS: Sensitivity analysis
  6. ENGL: Energy contribution calculation
  7. FILT: Solvent filtering
  8. HBBR: Calculate hydrogen-bond bridges
  9. CMPC: Configuration comparison
  10. DENF: Combined solvent conformation
  11. PRTG: Print cavity grids
  12. STVG: Estimate solute and cavity volumes with multiple grids
  13. CPOC: Calculate cavity/pocket occupancie
  14. GENS: Calculate generic solvent sites
  15. IWSL: Prepare site assignment file, entropy input
  16. GSAN: Generic site analysis choice
  17. GSAO: Generic site analysis options
  18. GABF: Generic site graph analysis by frame
  19. GENT: Configuration with input torsion angles
  20. TAND: Calculate torsion angle distribution
  21. TORT: Energy analysis along a torsion
  22. LPST: Energy analysis of backbone loop solutions
  23. IDAG: Aggregate insertion and deletion sites
  24. SUPI: Trajectory superimposition
  25. SVDP: Solvent dipole distribution
  26. MOND: Bulk density monitoring
  27. BLKW: Bulk density monitoring
III.09.01. Bulk solvent distribution function calculation level specification

DSTC Key1 [Key2]

Related key: DSTP.

Key1 (level of distribution function calculations):

III.09.02. Functional-group analysis of the solute.

FCGA [Key1]

Prerequisite: SLTA
Potential successor: FCGD

Key1 (specifies groups to be printed):

  • ALL - All functional groups will be listed
  • OERR - Only the ???? and VERR functional groups will be listed. Atoms in these groups are likely to have missing or superfluous bonds.

    The solute is broken up into functional groups and a list of bonded atoms is given (based on atom-atom distances). It is also useful for the detection of errors in the connectivity generated by the program. The errors found can be corrected with the keys MAKB or BRKB

    III.09.03. Calculation of average force on the solute atoms

    ATFR

    Potential successor: PXAN

    Also requires the PXAN or the SCAN key.

    III.09.04. Calculation of average field gradient at the solvents

    FLDG Key1, Key2, Data

    Prerequisite: PXCR
    Potential successor: PXAN

    Variable to set to 'T' in the preprocessor: FG

    Key1 (specifies the type of file to write the gradients):

    Only works for general solvent. Also requires the PXAN or the SCAN key.

    III.09.05. Sensitivity coefficent calculation

    SENS Data

    Data: logfreq

    III.09.06. Calculate energy contributions

    ENGL Key1 Key2 [Data]

    Prerequisites: SLTA, SUUC

    Key1 (source):

    For a configuration analyzed, the program prints

    ia, ig,E(es), E(d), E(r), E14(es), E14(d), E14(r) (2i10,3e12.5,3e12.5)
    ia, Esw(es), Esw(d), Esw(r)(i10,3e12.5)
    ia, ER(es), ER14(d), ERsw(r), sigma, Nnb, N14, Nw (i10,3e12.5,3i6)

    where E(*) refers to different terms of energy between solute atom ia and the rest with (es), (d) and (r) indicating electrostatic (1/r), dispersion (1/r6) and repulsion (1/r12) terms, respectively; 14 refers to 1-4 interactions; w refers to solute-water interactions; R refers to the sensitivity coefficient of the energy with respect to the Lennard-Jones sigma of solute atom ia; Nnb, N14, Nw are the number of non-bonded, 1-4 and solute-water terms contributing, respectively.

    The calculated contributions are also calculated for each group (residue) as defined in the .slt file.

    For trajectory analysis,a the calculated properties are averaged over the trajectory; the standard deviations are also calculated.

    Also requires the PXAN or the SCAN key.

    III.09.07. Solute and/or solvent filtering

    FILT Key1 Key2 Key3 Key4 Key5 Data :

    Prerequisite: CNFG or TRAJ

    Key1 (basis of filtering):

    For Key1=SOLV: For Key1=GEOM: For Key1=NMBR: For Key1=ENRG: For Key1=FRAM: ADDITIONAL INPUT:
    Data: incvers ADDITIONAL INPUT for Key2=TR**:
    Data: nfrmmax, nfrmfreq, nfrmskip III.09.08. Calculate hydrogen-bond bridges

    HBBR Key1 Data

    Prerequisite: SLTA

    Key1 (anchor atoms specification)

    III.09.09. Configuration comparison

    CMPC Key1, Key2, Key3, Data

    Prerequisite: SLTA
    Potential successor: CNFG

    Key1 (input file format):

    This operation reads in two configurations, calculates the average, minimum and maximum RMS between the solute molecules (possibly subject to the operation specified by Key2) and spreads out the compared molecules of the solute in the ix1-ix2 plane for comparison. The centers of the molecules will be shifted to this plane in a rectangular arrangement.

    III.09.10. Combining solvents from a trajectory

    DENF Key1 [Data]

    Prerequisite: TRAJ

    Key1 (output format):

    Data: nmcmax, nmcfreq, nmcskip This operation scans the trajectory specified and puts the first atom of all the solvent molecules into a single configuration onto a file <jobname>_den.dat (with the run and version number taken from the trajectory file's name) in InsightII free format) as He atoms If the output file already exists, a file with incremented version number will be created. If parts of the solute were allowed to move, the affected solute atoms are also collected. The resulting configuration is a simple approximation to the three-dimensional density.

    III.09.11. Cavity grid listing

    PRTG Key1 Key2 Key3 [Data]

    Prerequisites: CNFG, GCEN
    Potential successor: MOND or BLKW
    Variable to set to 'T' in the preprocessor: PG

    Grid points representing cavities around the solute will be enumerated and separated into mutually disconnected clusters. The first cluster is the set of gridpoints surrounding the solute and the remaining ones (if any) represent cavities inside the solute. The program writes out the solute and the internal cavities to the file <jobname>.<numrunw>.grd and, if requested, the solute and external gridpoints (in the same format) to <jobname>.<numrunw>_2.<ext> iwhere <ext> is pdb, CRD or grd for PDB, CHARMM, Insight free format, resp.. The gridpoints are written as helium atoms with zero charge and the group (residue) number is the group number of the last solute atom plus the cluster number. Note, that only a layer of gridpoints near the solute atoms are written providing a bounding shell for each cavity.

    Key1 (grid file format):

    III.09.12. Estimate solute and cavity volumes with multiple grids

    STVG [Data]

    Prerequisites: CNFG, GCEN, SLTA
    Variable to set to 'T' in the preprocessor: PG

    Data: diamslvnew, ngridrep

    The progam estimates the volume of the solute and the volume of each cavity based on the number of gridpoints covered/left uncovered, resp. The first 'cavity' is the volume outside the solute. The averaging assumes that the cavities occur in the same order in each sample.

    III.09.13. Calculate cavity occupancies

    CPOC Key1 Key2 Key3

    Prerequisite: CNFG,PRTG, GCEN

    This key calculates the (average) number of solvents in each cavity as determined by the key PRTG
    Key1 (select cavity/pocket)

    III.09.14. Calculate generic solvent sites

    GENS Key1 Key2 Key3 Key4 Key5 Data

    Prerequisite: PBCN
    Can not be followed by: IWSL Related keys: GABF, GSAN, GSAO, IWSL Variable to set to 'T' in the preprocessor: PG

    Key1 (type of calculation)

    For Key1=CALC or Key1=CCAL:

    Both GENS CALC and GENS COMP are based on using the so-called Hungarian method to obtain an optimal matching of solvent coordinates to (generic) site coordinates. With GENS CALC the calculation scans a simulation trajectory and dedudeces iteratively a set of generic sites. For each configuration, the solvent coordinates are matched to the site estimates, allowing for the same solvent to contribute to different sites during the simulation.

    When GENS CALC is used, the program will write several files:

    Depending on Key4, the sites or solvents may be preceded by the full solute. Representative and the composite configurations list the solvents in the order of the sites they match to, but may be offset due to missing matches. Also, the standard output has a combined list of the sites and their representatives with their properties.

    The run and version number for these files are taken from the trajectory file and extension defined by Key3. If the output file already exists, a file with incremented version number will be created.

    III.09.15. Prepare site assignment file, entropy calculation, clustering sites

    IWSL Key1, Data

    Potential successor: GENS

    Key1 (ask for entropy input)

    III.09.16. Generic site analyses

    GSAN Key1 [Key2 [Data]]

    Prerequisite: IWSL
    Can not be followed by: GENS
    Related key: GSAO

    Key1 (analysis type)

    This key activates the selected analysis options after a generic site calculation (with the key GENS)

    III.09.17. Generic site analysis options

    GSAO Key1 Key2 Key3 [Data]

    Prerequisite: IWSL
    Can not be followed by: GENS
    Related key: GSAN

    Key1 (matrix sort option)

    III.09.18. Generic site graph analysis by frame

    GABF Key1 [Data]

    Key1 (print option)

    Calculate the hydrogen-bond matrix for each frame, find the connected clusters and accumulate (and average) the matrix where all cluster members are connected. The matrix of the connectivity probabilities will be written to a file with extension .hbf.

    III.09.19. Configuration with input torsion angles

    GENT Data

    Prerequisites: CNFG, TORD.

    Data ntorch

    Fdata: ntorch records If the atoms ia1 and ia2 are listed in the reverse order on the torsion list, the torsion angle used will be -ang.

    III.09.20. Calculate torsion angle distribution

    TAND

    Prerequisite: TORD
    Potential successor: PXAN

    When this key is present, analysis of a configuration or trajectory (activated by the SCAN key) will include the calculation of torsion angle distributions for torsions defined by the TORD key. The number of angular grids will be determined by the value of #TD set by the preprocessor. Besides the distribution, the program prints the standard deviation (meaningful for not too wide distribution only), the circular variance (on a scale of 0 to 1) and the number of grids sampled.

    III.09.21. Energy analysis along a torsion

    TORT Key1 Data

    Prerequisites: CNFG, TORD.

    Key1 (reference energy)

    This key calculates the various non-bonded and torsion terms as the selected torsion angle is driven through the full circle. The values printed are relative to the initial torsion angle.

    If a conformation was generated previously with the GENT key then the torsion scan will start at those angles (otherwise at 0).

    III.09.22. Energy analysis of backbone loop solutions

    LPST Key1 Data

    Prerequisites: CNFG, LOOP, TORD

    Key1 (reference energy)

    This key calculates all possible backbone solutions for the six torsions following the selected torsion angle, and calculates the same energy contributions as the TORT key.

    HBBR Key1 Data

    III.09.23. Aggregate insertion and deletion sites

    IDAG Key1, Key2, Data

    Prerequisite: IDLG
    Related key: GCEN

    Key1 (select insert/delete/both)

    This option will prepare a PDB file containing the solute and the sites of insertions and/or deletions as read from the insertion deletion log file (created by the IDLG key during a grand-canonical ensemble run (see key GCEN)). Insertion sites will be represented as oxygens while deletion sites as nitrogens.

    For solute atom records, the occupancy and temperature factor columns will contain the distance of the closest deletion and/or insertion sites (as specified by Key1) and the number of such sites, resp. For the deletion and/or insertion site records, the occupancy and temperature factor columns will contain the index of the proximal solute atom and the distance from it, resp.

    III.09.24. Trajectory superimposition

    SUPI Key1, Key2, Data

    Prerequisite: TRAJ

    When filtering a trajectory, the filtered configuration will be superimposed to reference configuration.
    Key1 (reference structure)

    III.09.25. Solvent dipole distribution

    SVDP Data

    Prerequisites: SLVA, SLTA

    Data exa, eyb, ezc

    The program will calculate the distributions of the angles between the solvent's dipole vector and the three coordinate axes.

    III.09.26. Monitoring bulk solvent density

    MOND [Data]

    Prerequisite: PBCN
    Prerequisite of: BTUN
    Can not be followed by: PRTG

    Data: glimminx, glimmaxx, glimminy, glimmaxy, glimminz, glimmaxz

    III.09.27. Monitoring bulk solvent density

    BLKW [Data]

    Prerequisite: PBCN
    Prerequisite of: BTUN
    Can not be followed by: PRTG

    Data: wlayer

    III.10. Proximity analysis keywords

    1. PXYZ: Permute coordinates read
    2. PXCR: Proximity criterion choice
    3. FCGD: Functional group definition
    4. PXLM: Omitting atoms from proximity considerations
    5. RFSL: Set first shell radii
    6. RMOD: First shell radius modifications
    7. DBLG: Double proximity radial gridsize
    8. DIPC: Dipole fluctutations
    9. GENV: Proximity volume element generation
    10. GRAV: Proximity RDF averageing
    11. PXTD: Proximity table difference calculation
    12. PXDP: Solute-solvent proximity dipole distr. calculation
    13. PXGR: Solute-solvent proximity g(r) and coordination number calculation
    14. PXBE: Solute-solvent proximity binding energy calculation
    15. PXWW: Solvent-solvent proximity distribution calculation
    16. VOLE: Volume element regeneration
    17. VORO: Voronoi analysis
    18. RTIM: Solvent residence time
    19. DIFC: Solvent diffusion constant
    20. TAUC: Torsion angle autocorrelation

    III.10.01. Permute coordinates read from history file

    PXYZ Data

    Potential successor: CNFG

    Data: ix, iy, iz

    The coordinates will be permuted before analysis or before writing the new trajectory. This functionality allows the use of a trajectory file where the hexagonal prism's orientation is different from the convention MMC is using. Note, that it is only meaningful to use this key with hexagonal prism cell.

    Note also that the change is applied only to the configurations analyzed. Thus, depending on whether the last configuration read was analyzed, the configuration current after a SCAN may or may not have the coordinates permuted.

    III.10.02. Proximity criterion choice

    PXCR Key1 Key2 Key3 [Data]

    Prerequisite: SLTA
    Prerequisite of: GENV, PXGR, PXWW, PXBE, PXLM, TIMM, FLDG, RTIM, DIFC, PXDP, PXTD, PXWR

    Key1 (solvent density source):

    III.10.03. Functional group definition

    FCGD Key1 Data Fdata

    Prerequisite: SLTA
    Can not be followed by: FCGA

    Key1 (definition type):

    Note, that care must be exercised when this option is used after restoring a checkpoint file. If the INTG key was used before resulting in the non-default grouping of rdf's then using it again or using the STND key (to revert to the default) will abort the run since the original atom by atom information is lost.

    III.10.04. Omitting atoms from proximity considerations

    PXLM Data

    Prerequisite: PXCR

    Data: nfadel, nladel

    III.10.05. First shell radius definition

    RFSL Key1 [Data]

    Prerequisite: SUPT
    Potential successor: SLTA

    Can not be followed by: SLTA
    Key1 (source of radii):

    III.10.06. First shell radius modification

    RMOD Data

    Data: nmod, nmod times : pfll, pfl, rfsl

    III.10.07. Radial gridsize doubling

    DBLG

    This key will double the gridsize for the proximity analysis radial grid, modify all the accumulators involved and print the full result pages with the modified gridsize. If the PXPL key was used earlier, the proximity distributions will also be written on the .pxp file. Note, that the change is permanent for the remainder of the calculations.

    III.10.08. Dipole fluctuation calculation

    DIPC

    The fluctuation of the dipole moments will be calculated in the proximity regions to give estimates of the dielectric constant.

    III.10.09. Monte Carlo volume element estimation

    GENV Key1 Data [Data]

    Prerequisites: SLTA, PXCR, PXGR

    Key1 (possible visualization):

    The calculation of the proximity RDF's require the volume element of shells in the various proximity regions. This volume is calculated by a Monte Carlo method: random points of uniform distributions are generated and the volume of each shell is assumed to be proportional to the number falling into them.

    Generating a certain number of random points after each snapshot analyzed (see keys SCAN or PXAN), instead of generating them in one step at the beginning or end, allows for an approximate accounting of the variation of these volume element as a result of conformational changes in the solute.

    III.10.10. Calculating RDF extrema averages

    GRAV Key1 [Data] Prerequisite: PXGR

    Key1 (averaging strategy):

    If primary or total RDF's were calculated, the program averages the location and values of the first and second maxima, and the first minima for atoms that are sufficiently solvated as defined by the threshold values inputted above.

    III.10.11. Proximity table difference generation

    PXTD Data Prerequisite: PXCR

    Data: filename1, runnum1, filename2, runnum2

    The program reads both checkpoint files and calculated the differences of the table entries (shell volumes, coordination numbers, densities, energies, etc.). The two checkpoint files should refer to the same system.

    III.10.12. Solute-solvent dipole distribution function calculation

    PXDP

    Prerequisites: PXCR, PXGR
    Potential successor: PXAN

    This key will result in the calculation of two functions for the solute atoms, based on the angle THETA between the line drawn from the solute atom to the first solvent atom and the bisector of the angle S2-S1-S3 (i.e., assuming that the solvent is water and the first atom is the oxigen, the dipole direction):

    III.10.13. Solute-solvent proximity radial distribution function and coordination number calculation

    PXGR Key1 Key2 Key3 Data

    Prerequisite: SLTA
    Prerequisite of: GENV, PXGR, PXWW, PXBE, PXLM, TIMM, FLDG, RTIM, DIFC, PXDP, PXTD, PXWR

    Key1 ('total' solute-solvent g(r) calc.):

    The total g(r) is calculated for the solute atoms by considering all solvents. The proximity g(r) is calculated for solute atoms by considering only solvents proximal to it. Coordination numbers (and their distributions) will be calculated if either total g(r) or proximity g(r) calculation was requested.

    III.10.14. Solute-solvent proximity binding energy calculation

    PXBE Data

    Prerequisite: PXCR
    Potential successor: PXAN

    Data: epmink, epdivk

    This calculation will generate the distribution of the solute solvent energies for solvents in the first shell, as well as the sum of solute solvent energies, summed both for all solvent is the first shell and for all solvents in the whole proximity region.

    III.10.15. Solvent-solvent proximity distribution calculation

    PXWW Key1 Key2 Data

    Prerequisite: PXCR
    Potential successor: PXAN

    Key1 (coordination number and pair energy distr. calc.):

    The solvent-solvent coordination number is the mean number of solvents within a sphere of radius rfsww, averaged over all solvents considered. The solvent-solvent pair energy is the average of water-water energies, calculated for waters within a distance of rfsww and at least one of the pair being considered for the solvent-solvent analysis. The solvent-solvent binding energy is the sum of water-water energies, averaged fover all solvents considered.

    III.10.16. Volume element preservation/regeneration

    VOLE Key1

    Can not be followed by: RCKP Key1:

    III.10.17. Voronoi polyhedra analysis

    VORO Key1, Key2, Key3, Data

    Prerequisites: CNFG, PBCN
    Potential successor: PXAN

    Key1 (use of PBC):

    NOT IMPLEMENTED YET

    III.10.18. Solvent residence time calculation

    RTIM Key1, Data

    Prerequisite: PXCR
    Potential successor: PXAN
    Variable to set to 'T' in the preprocessor: DR

    Key1 (shell boundary)

    The program will calculate the autocorrelation function of the residence function (0: out; 1: in) and save it on a file <jobname>.rtm to be processed into actual times by a seperate program.

    III.10.19. Solvent diffusion constant calculation

    DIFC Data

    Prerequisite: PXCR
    Variable to set to 'T' in the preprocessor: DR

    Data: timestep, dcgrid, rdc1, rdc2

    The program will calculate the conformational average of the solvent displacements in each proximity region; within the first shell, the first two shells and the range defined by rdc1 and rdc2.

    III.10.20. Torsion angle autocorrelation calculation

    TAUC Data

    Prerequisite: PART

    Data: nmctorauc,tauc_timestep, tauc_min, ntaucprint

    III.11. Numerical precision related keywords

    1. HRDW: Computer architecture
    2. SLFT: Self test
    3. SLFT: Self test tolerances
    4. FIXD: Round-off error elimination
    5. REGE: Regenerate from torsion angles
    III.11.01. Computer architecture-dependent algorithm selection

    HRDW Key1

    Potential successor: NSLV

    Key1 (code-determining hardware selection):

    III.11.02. Periodic self test

    SLFT Key1 [Key2] [Data]
    Related key: PDBT.

    At each nrep-th (see key RUNS) step the program performs a consistency check as specified by Key1.

    Key1 (self test level):

    III.11.03. Input self-test tolerances

    STOL Data

    Prerequisite: CNFG

    Data: engtol, virtol, tortol, comtol, zmattol, cslttol, d12tol, d13tol, wsumtol
    Self test tolerances for various quantities:

    The default of engtol is relaxed by a factor of 10 for (TPN) ensemble run (see key IBEN), for runs with more than 1000 solute atoms or 1000 solvent molecules, for systems with active torsions (see key PART), or when the C@G7 lines are active; the default of tortol is relaxed to 20.0 when loop moves don't even touch the rest of the chain (see LOOP LIMM); the default of cslttol is relaxed to 0.001 for systems with more than 200 solute atoms.

    III.11.04. Energy contribution recalculation to eliminate round-off errors

    FIXD Key 1 Data

    Key1: (action on self-test failure):

    There are several quantities that the program calculates by accumulating the changes to them. Round-off errors may cause some self tests to fail. This option eliminates periodically these round-off errors.

    III.11.05. Regenerate solute from torsion angles to eliminate round-off errors

    REGE

    Prerequisite: TORD, CNFG

    When this keyword is invoked, the program takes the torsion angles read from the .crd file and applies them to the solute molecule given as input after the SLTA key. The rigid core of the molecule will be taken from the .crd file. This way the input molecule's bond lengths and bond angles will be made to agree with the input definition.

    III.12. Developpers' keywords

    1. TEST: Special code-testing
    2. OPTN: Direct option input
    3. DBUG: Input debug keys
      listing
    4. CHKP: Checkpoint file saving frequency
    See also the key OUTP PROP.

    III.12.01. Special code testing

    TEST Key1 [Data]

    Key1: (test type):

    III.12.02. Direct option input

    OPTN Data

    Data: nopt, nopt times: index, value

    This key allows to set the internal option arrays without the keyword mechanism normally used.

    III.12.03. Input debug keys

    DBUG Data

    Data: nkeys, nkeys times: index, value

    numbers and the file unit numbers will be also listed. If

    III.12.04. Checkpoint file saving frequency

    CHKP Data

    Potential successor: PXAN

    Data: nrecd

    III.13. Keyword summary

    IV. EXAMPLES

    In this section a number of commented input files are presented illustrating several different options.

    1. Glycine in water - minimal input; restart from checkpoint
    2. Ethanol in ethanol (neat liquid)
    3. Grand-canonical ensemble
    4. Constant pressure ensemble
    5. Dimethylphosphate-Sodium complex with freely moving sodium
    6. Torsion angle sampling
    7. Potential of mean force between dimethyl phosphate and sodium ions with adaptive umbrella sampling
    8. Creation/annihilation polynomial path thermodynamic integration
    9. Finite difference thermodynamic integration
    10. Widom insertion method
    11. Proximity analysis, trajectory filtering
    12. Conformation filtering
    13. Generic site calculation
    14. Primary hydration shell calculation
    15. Cavity and pocket determination

    IV.1. Glycine in water - minimal input; restart from checkpoint file

    This is a simple canonical ensemble MC run of a simple solute in water. 

    FILE glycine ! File names start with glycine.
! Run # = 1
NSLV 215      ! One glycine and 215 waters
TEMP 298.0     ! 298 Kelvin
SVVC SPCC 7.75 ! Solvent-solvent cutoff=7.75 Å
SUVC MICC 0.0  ! MI on the solute
PBCN FCC  14.81526   ! Face-centered cubic cell
STEP 0.15 10.0 0.55 40.0 40
! Solute and solvent stepsizes, slt move frequency
SUPT AM94   ! Solute-solvent potential is AMBER
SVPT TIP4 TIP4!Solvent-solvent potential: TIP4P
SLTA SMPL MMC READ  10       ! 10 solute atoms
   35     2.131    -1.202    -0.377    -.278623
   34     2.096     1.058     0.348    -.260774
   33    -0.549    -1.370     0.0      -.385556
    1     1.526     0.0       0.0       .295148
   18     0.0       0.0       0.0      -.040730
   20    -1.602    -1.374     0.044     .147624
   20    -0.249    -1.903    -0.836     .158476
   22    -0.351     0.514    -0.895     .061766
   22    -0.381     0.536     0.868     .086565
   23     3.068    -0.947    -0.337     .216103
CNFG RANC ASCI SIZE !Random initial configuration
! Run 100000 steps
RUNS 100000 10000 25000 0100000  100000  100000
! Run  an other 100000 steps, starting from the end of the first run.
! Run # = 2
RUNS 100000 10000 25000 0100000  100000  100000
STOP
    Click HERE to view the output of this run

    Files needed for run: none
    Files created by the run:

    Suppose we decide to extend the run for an other 100000 steps. It can be done by reading in the checkpoint file created after the run: 
    ! Extend the previus run with an other 100000 steps
FILE glycine 2! Set Run # to 2
RCKP ! Restore the status from the checkpoint file
! glycine.2.ckp
RUNS 100000 10000 25000 0100000  100000  100000
STOP
    Click HERE to view the output of this run

    Files needed by the run:

    Files modified by the run:

    IV.2. Ethanol in ethanol (neat liquid).

    This input runs a MC simulation in liquid ethanol. One ethanol molecule is considered the 'solute' 

    FILE etnl
PRNT ECHO       ! Echo formatted input
TITL Ethanol
HRDW SCLR       ! Scalar machine
SVVC SPCC 10.0  ! Solvent-solvent cutoff
SUVC SPGC 10.0  ! Solute-solvent cutoff
PBCN RECT 26.0
! Unit cell is a cube with 26 Å egdes
NSLV 181     ! One solute and 181 solvents
TEMP 298.0   ! Simulation temperature: 298 K
MOVE SHCY
! Molecules are selected by the shuffled cyclic scheme
STEP  0.50 30.0   0.50  30.0 40
! Solute and solvent stepsize params
SUPT AM94        ! AMBER solute
SVPT GENL AM94 ! General solvent, AMBER library
SAMP METC        ! Metroplis sampling
SLVA 9      ! Read 9 solvent atoms
   18    3.108016  0.653187 -8.526237 -0.198000
   22    2.814817 -0.348721 -8.761003  0.066000
   22    2.517394  1.015322 -7.710808  0.066000
   22    2.956134  1.278341 -9.381230  0.066000
   18    4.596542  0.674197 -8.131964  0.131000
   22    4.748424  0.049042 -7.276970  0.066000
   22    5.187164  0.312062 -8.947392  0.066000
   36    4.988388  2.013194 -7.818209 -0.566000
   23    5.982338  2.026292 -7.548412  0.303000
SLTA  SMPL MMC READ 9   !read 9 solute atoms
   18    3.108016  0.653187 -8.526237 -0.198000    1
   22    2.814817 -0.348721 -8.761003  0.066000    1
   22    2.517394  1.015322 -7.710808  0.066000    1
   22    2.956134  1.278341 -9.381230  0.066000    1
   18    4.596542  0.674197 -8.131964  0.131000    1
   22    4.748424  0.049042 -7.276970  0.066000    1
   22    5.187164  0.312062 -8.947392  0.066000    1
   36    4.988388  2.013194 -7.818209 -0.566000    1
   23    5.982338  2.026292 -7.548412  0.303000    1
DSTC NONE ! No distribution function calc.
CNFG READ ASCI NOFX
! Read initial coordinates in ASCII from etnl.crd
RUNS 1000000 10000 500000 100000 10000
!Run 1M steps
RUNS 1000000 10000 500000 100000 10000
!The second run will start at the end of the first.
STOP
    Click HERE to view the output of this run

    Files needed for run:
    etnl.crd: initial configuration.
    Files created by the run:

    IV.3. Grand-canonical ensemble

    This input runs MC simulation in the grand-canonical ensemble. The solute is the protein BPTI (from the PDB). 

    FILE bpti
PRNT DETL !Echo formatted data, print extra inf
HRDW VC32          ! 32-bit vector
SVVC SPCC 7.75     ! Solvent cutoff
SUVC MIGC 0.0      ! MI on the solute
PBCN RECT 24.87  40.39 66.2 !Rectangular PBC
MOVE RAND          ! Random selections
TEMP 298
NSLV 10   !10 solvent molecules to start with
STEP  0.00   00.0   0.55    40.0  60
! slt, slv stepsize, slt move freq.
SVPT TIP4 TIP4
SUPT AM94
SAMP FBSC 0.5  ! Scaled force-biased sampling
    2
 4.00   0.2   7.0     1.0
PMOD AM94 1
   75   27
    8            1.6     0.2
SLTA  SMPL MMC READ 892
 N3      -2.55400   7.05800  12.80100  -0.37640    1  ARG N
   .....
   .....   891 more lines describing solute atoms ....
   .....
DSTC NONE
GCEN CAVB RSIG ALTI
! Grand-canonical ens., cavity biased,
! LJ sigma-based cavity radii, alternating i/d
    1.0    2.5      1000.0       30   30   30 1000    1    1   0000100     50000
CNFG RANC ASCI SIZE
!Random configuration of 10 (from NSLV) waters
! using the given cell size and writing ASCII file
RUNS 10000  100 25000 100000 100000 100000
STOP
    Click HERE to view the output of this run

    Files needed for run: none.
    Files created by the run:

    IV.4. Constant pressure ensemble run

    This input runs MC simulation in the isothermal-isobaric ensemble. The solute is the protein a Na+ - DMP- ion pair. 

    FILE dmpna
TITL run 2 slt molecs in water, only moves
TITL NA+DMP-: KOLLMAN-STRAATSMA(NA+)
HRDW VC32          ! 32-bit vector
SVVC SPCC 7.75     ! Solvent cutoff
SUVC MIGC 0.0      ! MI on the solute
PBCN RECT 32.031  21.354 21.354
TEMP 298
NSLV 484
STEP  0.10   00.0   0.55    40.0 05
SVPT TIP3 TIP3
SUPT AM94
SAMP FBSC 0.5  ! Scaled force-biased sampling
    2
 4.00   0.2   7.0     1.0
PMOD AM94 1   ! Potential library modification
   49   46
   11  1.59955E+00     0.478833E-01
SLTA  SMPL MMC READ 8
   49     2.70000   0.0       0.00000   1.00000    1
   38     0.0       0.0       0.0       0.91200    4
   35     0.74000   0.0       1.28171  -0.65500    4
   35     0.74000   0.0      -1.28171  -0.65500    4
   37    -0.99164  -1.25560   0.00000  -0.41000    4
   37    -0.99164   1.25560   0.00000  -0.41000    4
    4    -1.92075   1.43300   1.07251   0.10900    4
    4    -1.92075  -1.43300  -1.07251   0.10900    4
IBEN UNIF ISOT 1.0 250.0 500 100 !(TPN) ensemble
CNFG READ ASCI
RUNS 2000000 50000 1000000 0100000 100000
STOP SLFT
    Click HERE to view the output of this run

    Files needed for run:

    Files created by the run:

    IV.5. Dimethylphosphate-Sodium complex with freely moving sodium

    This input show how to allow tranlational/rotational degrees of freedom for molecules considered 'solute' by the program. 

    FILE dmpna
TITL Na+DMP-: Kollman-Straatsma(Na+)
TITL SPC water. H-R PBC
HRDW VC32          ! 32-bit vector algorithm
SVVC SPCC 7.75     ! Solvent spherical cutoff
SUVC MIGC 0.0
! MI on the solute based on group center distances
SUUC MIMC          ! Molec-based MI for intrasolute
PBCN RECT 32.031  21.354 21.354
! Rectangular PBC, long box
TEMP 298   |  NSLV 484
!Simulation at 298 Kelvin with 485 solvent molecules
STEP 0.1 0.0 0.55 40.0 10
! 0 slt stepsize, 0.55 Å, 40 deg solvent steps
SVPT TIP3 TIP3 ! Solvent-solvent pot: TIP3P
SUPT AM94      ! Solute-solvent pot: AMBER
SAMP FBSC 0.5  ! Scaled force-biased sampling
    2
 4.00   0.2   7.0     1.0
MOVE PRSP NRST
! Preferential sampling based on nearest slt atom dist
    3
 3.00   3.5   4.0    03.0    7.0    1.0
PMOD AM94 1    ! Potential library modification
   49   46
   11  1.59955E+00     0.478833E-01
! Solute's "own" solvent will be used for slt-slv
SLTA  SMPL MMC READ 8
! 8 solute atom - regular solute, read from input stream
   38     0.0       0.0       0.0       0.91200    1
   35     0.74000   0.0       1.28171  -0.65500    1
   35     0.74000   0.0      -1.28171  -0.65500    1
   37    -0.99164  -1.25560   0.00000  -0.41000    1
   37    -0.99164   1.25560   0.00000  -0.41000    1
    4    -1.92075   1.43300   1.07251   0.10900    1
    4    -1.92075  -1.43300  -1.07251   0.10900    1
   49     2.70000   0.0      -0.00000   1.00000    2
PARD UNIF SHCY 1.0 0.3 0.0 1.0 0.3 0.0 ~
1.0 0.3 0.0 1.0 -1 ! Last slt group is independently moved
DSTC NONE   ! No distribution function calc
SLFT BASC   ! Basic self tests performed
CNFG READ ASCI ! Input from dmpna.crd in ASCII
RUNS 2000000 50000 1000000 0100000 100000
! Run 2M steps
STOP
    Click HERE to view the output of this run

    Files needed for run:
    dmpna.crd: initial configuration.
    Files created by the run:

    IV.6. Torsion angle sampling

    This input show how to sample torsional degrees of freedom for molecules considered 'solute' by the program. 

    FILE aa3
TITL Three amino acids in vacuum
TITL Solute molecules move, swap and torsions
HRDW VC32          ! 32-bit vector
SVVC SPCC 12.00    ! Solvent cutoff
SUVC MIGC 0.0      ! MI on the solute
SUUC MIMC          ! Molec-based MI for intrasolute
PBCN HEXG 99.00  99.00  99.00    !hexagonal PBC
TEMP 998
NSLV 0   ! Only solute
STEP  0.00   00.0   0.55    40.0 10
SVPT TIP3 TIP3 ! Solvent-solvent pot: CHARMM
SUPT CHRM
SLTA  SMPL MMC READ 34
 NH1      1.24900   2.31500  -0.79400  -0.47000    1  ILE N
 CT1      1.55500   3.73000  -0.91700   0.07000    1  ILE CA
 C        1.31700   4.34000   0.45300   0.51000    1  ILE C
 O        2.09800   5.17100   0.91200  -0.51000    1  ILE O
 CT1      0.73100   4.43400  -2.05000  -0.09000    1  ILE CB
 H        0.48500   2.01200  -1.35900   0.31000    1  ILE HN
 HB       2.61200   3.83400  -1.12700   0.09000    1  ILE HA
 HA      -0.33600   4.20200  -1.80500   0.09000    1  ILE HB
 NH1      1.49900   2.78700   3.28800  -0.47000    2  LEU N
 CT1      2.55800   2.23900   4.11700   0.07000    2  LEU CA
 C        3.78000   3.10300   3.86100   0.51000    2  LEU C
 O        4.51100   3.44400   4.79000  -0.51000    2  LEU O
 CT2      2.81700   0.74400   3.78200  -0.18000    2  LEU CB
 CT1      3.92200   0.01400   4.58300  -0.09000    2  LEU CG
 CT3      4.04200  -1.44600   4.12400  -0.27000    2  LEU CD1
 CT3      3.68900   0.07700   6.10000  -0.27000    2  LEU CD2
 H        1.16000   2.16100   2.58800   0.31000    2  LEU HN
 HB       2.28700   2.35700   5.15800   0.09000    2  LEU HA
 HA       1.86200   0.19700   3.95500   0.09000    2  LEU HB1
 HA       3.03300   0.64900   2.69400   0.09000    2  LEU HB2
 HA       4.90300   0.49900   4.36300   0.09000    2  LEU HG
 HA       4.86700  -1.95800   4.66500   0.09000    2  LEU HD11
 HA       3.09800  -1.99800   4.32300   0.09000    2  LEU HD12
 HA       4.25600  -1.49900   3.03500   0.09000    2  LEU HD13
 HA       4.49200  -0.47200   6.63700   0.09000    2  LEU HD21
 HA       3.69100   1.12200   6.47100   0.09000    2  LEU HD22
 HA       2.71400  -0.38500   6.36400   0.09000    2  LEU HD23
 NH1      3.65600   6.57500   5.52300  -0.47000    3  GLY N
 CT2      3.90500   6.38700   6.93600  -0.02000    3  GLY CA
 C        5.38300   6.32500   7.16600   0.51000    3  GLY C
 O        5.89600   6.89700   8.12500  -0.51000    3  GLY O
 H        3.16200   5.83500   5.07000   0.31000    3  GLY HN
 HB       3.47900   5.44100   7.23900   0.09000    3  GLY HA1
 HB       3.51900   7.24100   7.47600   0.09000    3  GLY HA2
PRMF CHRM par_all22_prot_na.inp
!Torsion pot param file
PARD UNIF RAND 1.0 0.5 30.0 1.0 0.5 30.0 ~
 1.0 0.5 30.0 1.0 34 ! Note the continuation character ~
PART UNIF RAND RAND 1.0
!Torsion strategy and frequency
TORD INPT 7  !Definiton of active torsions
    5    2     180.0      30.0    1
    2    3     180.0      30.0    1
   10    9     180.0      30.0    2
   10   11     180.0      30.0    2
   14   15     180.0      30.0    3
   13   14     180.0      30.0    4
   14   16     180.0      30.0    4
CNFG RANC ASCI
SWAP 0.2 !Specify swapping frequency
RUNS 100000 100  50000 10000 5000
STOP SLFT
    Click HERE to view the output of this run

    Files needed for run:

    Files created by the run:

    IV.7. Potential of mean force between dimethyl phosphate and sodium ions with adaptive umbrella sampling

    This input performs one type of free energy simulation: it calculates the potential of mean force between two solute molecules using adaptive umbrella sampling. 

    FILE dmpna_pmf     ! SPC water. H-R PBC
TITL Na+DMP-: Kollman-Straatsma(Na+),
TITL #1   4 Å range  MI/RC G/A  Run A
HRDW VC32          ! 32-bit vector hardware
SVVC SPCC 7.75     ! Solvent-solvent cutoff
SUVC MIGC 0.0      ! MI on the solute
PBCN RECT 32.031 21.354 21.354 !Rectangular PBC
MOVE SHCY          ! Shuffled-cyclic selections
TEMP 298   |  NSLV 485
STEP  0.10 0.0 0.55  40.0  10
! 0.1 Åslt stepsize, no slt rotation, 0.55 Å,
!40 deg solvent steps; solute moves at every 10th MC step
SVPT TIP3 TIP3 ! Solvent-solvent pot: TIP3P
SUPT AM94      ! Solute-solvent pot: AMBER
SAMP FBSC 0.5  ! Scaled force-biased sampling
    2
 4.00   0.2   7.0     1.0
FREE PMF1 WRMM GEOC AUSE 2
!2-part PMF both parts of the slt moves
                      0.00      0.25      0.050     100.0     0.1
    2    0    3    1   0.01     17.      0.0       0.0      0.00002    0.9
    4    5    5   05    2.43      0.75
PMOD AM94 1    ! Potential library modification
   49   46
   11  1.59955E+00     0.478833E-01
SLTA  SMPL MMC READ 24 16 8
!3*8 solute atom - regular solute, last 1 atom is in the 2nd part
   38     0.0       0.0       0.0       0.91200    1
   35     0.74000   0.0       1.28171  -0.65500    1
   35     0.74000   0.0      -1.28171  -0.65500    1
   37    -0.99164  -1.25560   0.00000  -0.41000    1
   37    -0.99164   1.25560   0.00000  -0.41000    1
    4    -1.92075   1.43300   1.07251   0.10900    1
    4    -1.92075  -1.43300  -1.07251   0.10900    1
   49     2.70000   0.0       0.00000   1.00000    2
   38     0.0       0.0       0.0       0.91200    1
   35     0.74000   0.0       1.28171  -0.65500    1
   35     0.74000   0.0      -1.28171  -0.65500    1
   37    -0.99164  -1.25560   0.00000  -0.41000    1
   37    -0.99164   1.25560   0.00000  -0.41000    1
    4    -1.92075   1.43300   1.07251   0.10900    1
    4    -1.92075  -1.43300  -1.07251   0.10900    1
   49     6.70000   0.0      -0.00000   1.00000    2
   38     0.0       0.0       0.0       0.91200    1
   35     0.74000   0.0       1.28171  -0.65500    1
   35     0.74000   0.0      -1.28171  -0.65500    1
   37    -0.99164  -1.25560   0.00000  -0.41000    1
   37    -0.99164   1.25560   0.00000  -0.41000    1
    4    -1.92075   1.43300   1.07251   0.10900    1
    4    -1.92075  -1.43300  -1.07251   0.10900    1
   49     2.70000   0.0       0.00000   1.00000    2
DSTC NONE
CNFG READ ASCI NOFX
! input from dmpna_pmf.crd in ASCII
RUNS 5000000 20000 2500000 0  100000 200000
STOP
    Click HERE to view the output of this run

    Files needed for run:
    dmpna_pmf.crd: initial configuration.
    Files created by the run:

    IV.8. Creation/annihilation polynomial path thermodynamic integration

    This is an example of a different free energy simulation type: calculation the solvation excess Helmholtz free energy of NO in water with polynomial thermodynamic integration over a creation/annihilation path. 

    FILE NO_ti 10
TITL Excess Helmholtz free energy of solvation
of NO in water
TITL Calculated by three-point Gaussian
quadrature
PRNT ECHO
HRDW VC32          ! 32-bit vector
SVVC MINI         ! Solvent cutoff
SUVC MIGC 0.0      ! MI on the solute
PBCN RECT 14.74    !Rectangular PBC
MOVE RAND          ! Random selections
!Set creation/annihilation TI using the first quadrature point
!Lambda exponents: 4, 3, 2 (for 1/r^12, 1/r^6, 1/r terms, resp.)
FREE TICA  NOMX
 4.0   3.0  2.0     0.112702
TEMP 298
NSLV 108
STEP  0.00   00.0   0.40    35.0 50
SVPT TIP3 TIP3 ! Solv-solv potential is TIP3P
SUPT CHRM
!Add NO parameters
PMOD CHRM 2
  150   70    3.3    N/NO
    7 1.824            0.17
  151   70    3.3    O/NO
    8 1.751            0.159
SAMP METC 0.5
!Make sure the dummy atoms in the first copy are treated as one molecule
MOLD 2 2 4
!Solute has 4 atoms, 4 free energy atoms, and the first copy has 2 atoms
SLTA SMPL MMC READ 4 4 2    !Dummy solute
 DUM      -.570      000.0    0.00000   0.00000    1
 DUM      0.58       000.0    0.00000   0.00000    1
  150     0.57000    000.0    0.00000   0.02800    2
  151     -.58000    000.0    0.00000  -0.02800    2
DSTC NONE
CNFG READ ASCI NOFX
FIXD 2500000
!Equilibration
RUNS 2000000 100000 100000 500000 100000
RMCK !Remove unneded checkpoint file
!Production run, runnumber=11
RUNS 10000000 100000 100000 500000 100000
FILE  NO_ti  20
!Set lambda to the 2nd quadrature points
FREE TICA  NOMX
 4.0   3.0  2.0     0.500000
CNFG READ ASCI NOFX  12
RUNS 2000000 100000 100000 500000 100000
RMCK !Remove unneded checkpoint file
!Production run, runnumber=21
RUNS 10000000 100000 100000 500000 100000
FILE  NO_ti  30
!Set lambda to the 3rd quadrature points
FREE TICA  NOMX
 4.0   3.0  2.0     0.887298
CNFG READ ASCI NOFX  22
RUNS 2000000 100000 100000 500000 100000 00000
RMCK !Remove unneded checkpoint file
!Production run, runnumber=31
RUNS 10000000 100000 100000 500000 100000
!Evaluate the quadrature from the data on the three checkpoint files
TIQU REGL ALL  3 11 10
    0  001    0
STOP SLFT
    Click HERE to view the output of this run

    Files needed for run:

    Files created by the run:

    IV.9. Finite difference thermodynamic integration

    This is a example combines an other two free energy simulation types. It calculates the Helmholtz free energy difference between ethanol and acetone in water, using the finite differenece thermdynamic integration. The integrand at each quadrature point is obtained by the popular perturbation method. 

    FILE pm 1
TITL Perturbation method test ethane ->
acetaldehyde
TITL 3-point Gaussian quadrature (probably
inadequate)
HRDW VC32
SVVC SPCC 10.0
SUVC MICC
PBCN FCC    19.57563
TEMP 298 | NSLV 500
STEP    0.10   10.0  0.55      30.0  30
SVPT TIP3 TIP3
SUPT CHRM
!Lambda=0.112702 run
FREE PMLI CCMX
                    0.112702  0.102702  0.122702         0.0      1.0
SLTA SMPL MMC FILE 16 16 8 16 1
CNFG READ ASCI NOFX 50
!Equilibrate
RUNS 1000000 100000 10000 1000000 100000
RMCK !Remove equilibration checkpoint file to
save disk space
!Production run
RUNS 5000000 100000 10000 1000000 100000
!New coupling parameter - lambda=0.5
FREE PMLI CCMX
                    0.50      0.49      0.51        0.102702 0.122702
!Restore solute to original initial and final state
SLTA SMPL MMC FILE 16 16 8 16 1
!Read starting config from disk to prepare the new
!initial and final states
CNFG READ ASCI NOFX
!Equilibrate
RUNS 1000000 100000 10000 1000000 100000
RMCK
!Production run
RUNS 5000000 100000 10000 1000000 100000
!New coupling parameter - lambda=0.887298
FREE PMLI CCMX
                    0.887298  0.877298  0.897298    0.49      0.51
!Restore solute to original initial and final state
SLTA SMPL MMC FILE 16 16 8 16 1
CNFG READ ASCI NOFX
!Equilibrate
RUNS 1000000 100000 10000 1000000 100000
RMCK
!Production run
RUNS 5000000 100000 10000 1000000 100000
!Evaluate the TI quadrature
TIQU 2 2
    3    0    1    0    0    0
STOP SLFT FULL
    Click HERE to view the output of this run

    Files needed for run:

    Files created by the run:

    IV.10. Widom insertion method

    This input shows how to run a Widom insertion calculation based on a pure water trajectory run earlier by CHARMM. 

    FILE methane
TITL Calculate the chemical potential of Methane in Water
TITL using Widom cavity insertion method
HRDW VC32          ! 32-bit vector
SVVC MINI 10.0
SUVC MIGC 10.0      ! MI on the solute
PBCN RECT 15.5  15.5 15.5  !Rectangular PBC
TEMP 300
NSLV 124 !no of solvent molecs in the CHARMM run
SVPT TIP4 TIP4 ! Solvent-solvent potential is TIP4P
SUPT CHRM      !Solute-solvent potential is CHARMM
FREE WIDO 2.6 64 64 64 2 5
SLTA SMPL MMC READ 5 5
 CT3     -0.13265   0.14073   0.03004  -0.36000    1  CSP3C1
 HA       0.84914   0.15299   0.45527   0.09000    1  CSP3H11
 HA      -0.35123  -0.84076  -0.33574   0.09000    1  CSP3H12
 HA      -0.84914   0.40992   0.77775   0.09000    1  CSP3H13
 HA      -0.18078   0.84077  -0.77776   0.09000    1  CSP3H14
TRAJ CHRM RGFX  !Open trajectory in CHARMM format
SCAN TRAJ 1000 5 00 0 500 2 1  2 3 1
STOP
    Click HERE to view the output of this run

    File needed for run:

    Files created by the run:

    IV.11. Proximity analysis, trajectory filtering

    This input show how to set up proximity analysis of a trajectory given as a series of PDB files and how to use this analysis to filter this trajectory file. 

    FILE 5ht1
TITL GCE run with fixed solute
HRDW VC32        ! 32-bit vector
SVVC SPCC 7.75  ! Solvent-solvent cutoff
SUVC MIGC 0.0 ! Group-cent based MI on the slt
PBCN RECT 70.0 75.0 58.0  !Rectangular PBC
TEMP 298
SVPT TIP3 TIP3
SUPT CHRM  ! Solute-solvent potential is CHARMM
MOLD 1 3192 ! Bypass the topology code
SLTA SMPL MMC FILE  3192
!Solute description is read from 5ht1.slt in MMC format
TRAJ ALLP RGFX ALST !Open trajectory in PDB format
FILT SOLV TRAJ ALLP SHL1 RVDW 1.4 0 0 10
!filter out solvents outside the first shell
TRAJ ALLP RGFX ALST 0 1 11
!Open filtered trajectory
PXCR WDEN BISE
!select proximity criterion type
PXWR ASCI
!Create ASCII proximity information file
PXBE -100.0  2.0
!Calculate solute-solvent energy analysis
SCAN TRAJ 250000 1 5000 0 2000 1000 2  3 1
!Run analysis
FILE 5ht1 1    !Reset run number to one
TRAJ ALLP RGFX ALST 0 1 11 !Open filtered traj
FILT ENRG TRAJ ALLP -100.0 -1.0 10
!Filter solvent by the solute-solvent energy
TRAJ ALLP RGFX ALST 0 1 11
!Open filtered trajectory again
FILT ENRG TRAJ ALLP -1.0 100.0 11
!Filter solvent by solute-solvent energy using a different range
TRAJ ALLP RGFX ALST 0 1 21
!Open trajectory after 1st energy filter
PXWR NONE
DENF INSG 0 1 0
!Create single configuration aggregating all solvents
TRAJ ALLP RGFX ALST 0 1 22
!Open trajectory after 2nd energy filter
DENF INSG 0 1 0
!Create single configuration aggregating all solvents
STOP
    Click HERE to view the output of this run

    Files needed for run:

    Files created by the run:

    IV.12. Conformation filtering

    This example shows the elimination of waters that are not internal to the solute according to their circular variance. The filtering is done on a CHARMM CRD file. 

    FILE scednas
PRNT ECHO
SVVC SPCC 12.00
SUVC SPGC 20.0
PBCN RECT 99.00  99.00  99.00 ! Large box
NSLV 16426
SVPT TIP3 TIP3 ! Solvent-solvent potential is TIP3
SUPT ATNO ! No potential information
BRKB  32 ~
    92    97    92    99   231   236   231   238   292   297   292   299 ~

   323   328   323   330   404   409   404   411   593  3691   856   861 ~
   856   863  1149  1154  1149  1156  1784  1789  1784  1791  1825  1830 ~
  1825  1832  1882  2138  2030  2035  2030  2037  2090  2095  2090  2097 ~
  2299  2304  2299  2306  2318  2323  2318  2325  2467  2472  2467  2474 ~
  2792  2797  2792  2799
!Break spurious bonds that would confuse the routine establishing
!the molecular connectivities
!Alternative: MOLD 1 3697
SLTA SMPL CRD FILE 3697
CNFG READ CHRM NOFX
FILT SOLV CONF CHRM CRCV RSIG 6.0 0.5 0 11
CNFG READ CHRM NOFX
FILT SOLV CONF CHRM CRCV RSIG 6.0 0.4 0 12
STOP
    Click HERE to view the output of this run

    Files needed for run:

    Files created by the run:

    IV.13. Generic site calculation 

    FILE GENSfile
SVVC SPCC 10.0  ! Solvent cutoff
SUVC MIGC 0.0   ! MI on the solute
PBCN RECT 56.0011  67.5172  81.0524 !Rectangular PBC
SVPT TIP3 TIP3
SUPT CHRM
SLTA SMPL MMC FILE 5402  0 0 0 1
TRAJ ALLH FLEX ALST
GENS CALC FIXI GETH SPDB FSRT ~
 125 5000000 1 10000 6.5 0.01 0.10 6.0 0.2 100 10
!Number of sites requested (initially): 500
!Number of simulation steps to read: 5000000
!Use every 1-th config (counted in MC stepnumbers)
!Discard the first 10000 MC steps
!For the matching, consider only sites within 6.5 Å
!Convergence threshold for the Hungarian method matching: 0.01
!Convergence threshold for the RMSD between iterations: 0.1
!Add a new site if a match is farther than 6.0 Å
!Previous site estimate will be mixed in with a factor of 0.2
!Maximum number of iterations over the trajectory: 199
!Minimum number of solvents in a configuration to be of use: 10
STOP NOSF
    Click HERE to view the output of this run

    Files needed for run:

    Files created by the run: from the simulation history, in PDB format

    IV.14. Primary hydration shell calculation 

    FILE alaval
TITL compacted alaval
TITL PHS, 2 shells
HRDW VC32
NSLV 190       !190 waters
TEMP 298.0
SUUC SPGC 21.00
SVVC SPCC 21.00
SUVC SPGC 21.00
PBCN PHS RSIG UPML 0.15 3.0 8.0 1.0 21.0 0.997 100 1000.0 3.
!PHS boundary conditions
STEP    0.00   0.0  0.40      30.0  20
SUPT CHRM
SVPT TIP3 TIP3
MIXR ARGE
SLTA SMPL MMC FILE 56
SAMP METC MOVE SHCY   !Plain Metropolis moves
PART UNIF SHCY CYCI 0.1          !Solute torsions active
PRMF CHRM  par_all22_prot_na.inp
TORD ALL SING 30.0  ! All torsions moved with same stepsize
CNFG READ ASCI
RUNS 200000 5000  100000  90000  9000000  90000000  1
STOP
    Click HERE to view the output of this run

    Files needed for run:

    File created by the run:

    IV.15. Cavity and pocket determination 

    FILE f8
PRNT ECHO
TITL 1tf_8_13
NSLV 0
TEMP 298.0
SVVC SPCC 9.05     ! Solvent cutoff
SUVC SPCC 12.0
PBCN RECT  70.0 70.0 70.0
STEP    0.00   00.0  0.55      40.0  30
SUPT CHRM
SAMP METC 0.5  ! Scaled force-biased sampling
MOVE RAND
SVPT TIP3 TIP3 ! Solvent-solvent potential is TIP4P
SLTA SMPL PDB FILE  3928
GCEN CAVB RSIG ! Grand-canonical ensemble (establishes grid)
    2.0    2.5      1000.0      250  250  250 1000    1    1  010000    20000
CNFG READ PDB
PRTG PDB ALLG AGLS 0 0 0.8 10.0 0.2 1.0 3.3 0.9 0.5 10.0
STOP NOSF
    Click HERE to view the output of this run

    Files needed for run:

    Files created by the run:

    V. BRIEF DESCRIPTION OF THE OUTPUT.

    Since the output of the program is largely self-explanatory, only a relatively brief description is given.

    The output starts with information about the program: self identification, dimensions used - important for compiling other programs that read the interrupt file of the program - the dates the program and the common blocks were modified and the date of the run.

    By default, the command lines are echoed verbatim. The command PRNT ECHO allows for direct echoing of the formatted data read as well. At the start of each run the the options selected are explained, numbers read are printed with their interpretation and unit assumed. These explanations are preceded by the keyword that was the source of the information, allowing the easy interpretation of the messages.

    The solute description consists of three parts:

    Each atom record consists of the following items (several of these may be blank):

    A very important feature of the program is the automatic self test on continuing a previously interrupted calculation (RCKP) as well as at every 2.5 millionth Monte Carlo step. It can also be called more frequently with the help of SLFT. If the self test failed a message "Discrepancy found" appears and the program stops (unless RCKP IGND or RCKP FIXD was used). Discrepancies may arise when the checkpoint file (.ckp) is corrupted or in case of a program error. If the latter is the case, contact the author of the program. A relatively harmless error can be an accumulation of round-off errors to unacceptable levels. Such errors can be prevented using the FIXD FIXD keys.

    The initial total energy and its components, and virial sum header precedes the periodic output (at each nmcrep-th step) during the simulation. This periodic output (one to three lines each) contains the stepnumber (N) the energy of the currently accepted configuration (E), the cumulative total energy average (<E>), the minimum and maximum energy value sampled (Emn and Emx), with the configuration number/1000 or /1000000 where it occurred after them in parenthesis, the solute binding energy (Us) the acceptance rate (a) of the current move type (for solute molecule move or torsion the acceptance rate printed is for the solute molecule or torsion angle sampled in this step) and the index of the molecule selected for move (I) and the move type coded as follows.

    For torsion changes the torsion number is also printed. For torsion on cloned solute, the torsion number corresponding to the original torsion list is also given. For molecule swap, the second molecule's index is also printed. The line is ended by the symbol A or R for accepted and rejected moves, resp. The second line gives the solute intramolecular energy terms when the solute is not rigid. An additional line gives grand-canonical ensemble information: current and average number of solvent molecules, smallest and largest number of solvent molecules sampled (Nmn, Nmx), the acceptance rate of the insertion/deletion steps (id acc), the number of uncovered gridpoints (Ngf) and the probability of finding a cavity (pcav).

    Periodically (at each nmcplt-th step), and at the end of each run a complete report of the averages and distributions is computed and printed. This report starts with the thermodynamic averages: energy, energy square, standard deviation and range for energy, total virial sum, solute virial sum and the pressure.This is followed by the error estimates using the method of batch means for the total energy or pressure. The error estimates are given as two standard deviations - corresponding to ca 95% confidence intervals. Since he method of batch means assumes that the actual batch means are statistically independent the program performs a test for independence - the so-called run test. At the and of this section there are two tables containing the critical values for the run test: for independence, the number of runs found should be between the two numbers corresponding to the number of "ups" (n1) and "downs" (n2) in the two tables. The program prints CORRELATED or UNCORRELATED when the run test was failed or passed, respectively and print ??? or >>> when the number of "ups" or "downs" were either too low or too high for the table. The test is repeated with double the blocksize until the actual run lengths prevent further doubling.

    During the run, various informations are gathered relevant to the convergence of the calculation: Number of times the solute is attempted to move, accepted to move; average translational and rotational displacements; the minimum, maximum and average correlations between the initial and final orientations; the total displacement of the system (square root of the molecular displacement square sums), information on the average molecular displacements, the correlation between two successive moves of a molecule (averaged over all moves and all molecules) and convergence indicators involving average molecular displacement square, acceptance ratio and correlation between successive moves; if preferential sampling was used (MOVE PRS*) the relative frequencies the molecules were attempted to move (perturbed), acceptance rates for individual molecules (useful to detect "stuck" molecules).

    For the various free energy options the distributions and free energy values calculated are printed here. Also, calculations involving the perturbation method print the solute binding energy estimate on the two limiting state ensembles to be used for convergence test and the range of energy values that went into the exponential averages - for reliability these ranges should not exceed cca 10 kcal/mol.

    For adaptive umbrella sampling, the range covered by each iteration is plotted with letters U and D indicating that the probability distribution is increasing or decreasing. In the plot each line corresponds to an iteration and the position of the U's and D's refer to the coupling-parameter range sampled.

    The program calculates first order quantum corrections to the intermolecular free energy with the method of Powles and Rickayzen. This calculation is based on the computed configurational average of the force and torque components and on the mass and moment of inertia of the molecule involved. The solute contribution to the quantum correction is also computed and printed, along with the computed force and torque component squares, molecular mass and moments of inertia.

    The various distribution functions are printed only optionally and they are self-explanatory. Each distribution is followed by a summary locating the extrema and other pertinent characteristics.

    Two more pieces of information relevant to the convergence of the run are given near the end of the output: the closest solvent-solvent distance and the number of steps elapsed during which all molecules were moved at least once.

    Proximity analysis results are printed by the atoms on which the proximity analysis was requested. The output may consist of several parts (depending on key PXPR and, of course, on the keys invoking the calculation of the various propeties).

    The output concludes with averages of the total energy and solute binding energy over the run, as well as the coupling parameter range sampled in the run, if the variable coupling parameter option was used.

    Throughout the output, lines with periods are used to delineate certain segments of the output. Warning messages (indicating potential problem with the input) are preceded by 5 dashes, strong warnings (likely, but not certain problems) are preceded by 5 equal signs, and error messages are preceded by 5 asterisks. Messages announcing some specific action (not taken all the time) are preceded by 5 plus signs. Messages announcing that an input setting is overriden by the program are preceded with five 'less-than' signs. Before stopping the program prints the run time and the number of overrides, warnings, strong warnings and errors encountered. It is always a good practice to find these messages in the output file to make sure the program was running as intended.

    Minimum critical values for the run test: 

      n2:  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
 n1  ----------------------------------------------------------
  2 |                                2  2  2  2  2  2  2  2  2 | 2
  3 |              2  2  2  2  2  2  2  2  2  3  3  3  3  3  3 | 3
  4 |           2  2  2  3  3  3  3  3  3  3  3  4  4  4  4  4 | 4
  5 |        2  2  3  3  3  3  3  4  4  4  4  4  4  4  5  5  5 | 5
  6 |     2  2  3  3  3  3  4  4  4  4  5  5  5  5  5  5  6  6 | 6
  7 |     2  2  3  3  3  4  4  5  5  5  5  5  6  6  6  6  6  6 | 7
  8 |     2  3  3  3  4  4  5  5  5  6  6  6  6  6  7  7  7  7 | 8
  9 |     2  3  3  4  4  5  5  5  6  6  6  7  7  7  7  8  8  8 | 9
 10 |     2  3  3  4  5  5  5  6  6  7  7  7  7  8  8  8  8  9 | 10
 11 |     2  3  4  4  5  5  6  6  7  7  7  8  8  8  9  9  9  9 | 11
 12 |  2  2  3  4  4  5  6  6  7  7  7  8  8  8  9  9  9 10 10 | 12
 13 |  2  2  3  4  5  5  6  6  7  7  8  8  9  9  9 10 10 10 10 | 13
 14 |  2  2  3  4  5  5  6  7  7  8  8  9  9  9 10 10 10 11 11 | 14
 15 |  2  3  3  4  5  6  6  7  7  8  8  9  9 10 10 11 11 11 12 | 15
 16 |  2  3  4  4  5  6  6  7  8  8  9  9 10 10 11 11 11 12 12 | 16
 17 |  2  3  4  4  5  6  7  7  8  9  9 10 10 11 11 11 12 12 13 | 17
 18 |  2  3  4  5  5  6  7  8  8  9  9 10 10 11 11 12 12 13 13 | 18
 19 |  2  3  4  5  6  6  7  8  8  9 10 10 11 11 12 12 13 13 13 | 19
 20 |  2  3  4  5  6  6  7  8  9  9 10 10 11 12 12 13 13 13 14 | 20
     ----------------------------------------------------------
       2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20

    Maximum critical values for the run test: 

      n2:  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
 n1  ----------------------------------------------------------
  2 |                                                          | 2
  3 |                                                          | 3
  4 |           9  9                                           | 4
  5 |        9 10 10 11 11                                     | 5
  6 |        9 10 11 12 12 13 13 13 13                         | 6
  7 |          11 12 12 12 14 14 14 14 15 15 15                | 7
  8 |          11 12 13 13 14 15 15 16 16 16 16 17 17 17 17 17 | 8
  9 |             13 14 14 15 16 16 16 17 17 18 18 18 18 18 18 | 9
 10 |             13 14 15 16 16 17 17 18 18 18 19 19 19 20 20 | 10
 11 |             13 14 15 16 17 17 18 19 19 19 20 20 20 21 21 | 11
 12 |             13 14 16 16 17 18 19 19 20 20 21 21 21 22 22 | 12
 13 |                15 16 17 18 19 19 20 20 21 21 22 22 23 23 | 13
 14 |                15 16 17 18 19 20 20 21 22 23 23 23 23 24 | 14
 15 |                15 16 18 18 19 20 21 22 22 23 23 24 24 25 | 15
 16 |                   17 18 19 29 21 21 22 23 23 24 25 25 25 | 16
 17 |                   17 18 19 20 21 22 23 23 24 25 25 26 26 | 17
 18 |                   17 18 19 20 21 22 23 24 25 25 26 26 27 | 18
 19 |                   17 18 29 21 22 23 23 24 25 26 26 27 27 | 19
 20 |                   17 18 20 21 22 23 25 25 26 26 27 27 28 | 20
     ----------------------------------------------------------
       2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20

    VI. COMPILATION OF THE PROGRAM.

    VI.1. Preprocessor

    The program source file mmc.for contains a number of special symbols allowing the customization of the program. To produce a valid Fortran source code mmc.f a preprocessor (also a Fortran program) pre.f is provided. The preprocessor performs two functions:

    If the generated mmc.f is inadequate for a particular input, mmc will abort the run and give specific instructions what to change.

    The preprocessor has to be compiled and executed (after which you can compile the program):

    The prepocessor as provided has default choices and values that will create a valid source file. There are three ways to change the defaults: Instructions for changing the source code of pre.f are at the beginning of that file. The part of the prepocessor that may be changed by the user is between the comment lines 'START Customization' and 'END Customization'. In the first segment, calls to the subroutine setopt select functionalities to enable/disable. In the second (longer) segment, calls to the subroutine setdim set maximum values for the various aspects of the system size. Some array sizes are generated automatically depending on the code segments activated.

    The following list gives the special symbols related to array sizes appearing in mmc.for with the names used in the preprocessor and the quantity they limit.

    Since the increased memory usage taxes the computer resources (both memory and disk-space), it is usually worth recompiling the program separately for significantly different system sizes (see the key PRCO SVMN for preparing a minimum-sized executable). Remember, however, that the checkpoint file (.ckp) is only readable by a program that was compiled with the same dimensioning!

    Code segment that can be selectively activated contains lines with the comment C@**. The selection and removal of such lines is controlled by a number of character variables in pre.f that can be set during the customization of the preprocessor. Except for the variable SYSTEM, the value is either 'T' or 'F'. The following variables can be set:

    Activating some of the codes requires deactivating a complementary code. Such code pairs are:
    C@NL - C@NA
    C@VC - C@NC
    C@DM - C@ND
    C@1R - C@NR
    C@RF - C@1R
    The preprocessor takes care of it automatically.

    Besides the activating and inactivating lines prefixed with C@**, the preprocessor provides the option to change all double precision (real*8) variables to quadruple precision variables (real*16). This option can be useful for calculations using the LOOP key. Note, however, that not all compilers support quadruple precision variables.

    The defaults for variables FG, PG, HU, and DR are set to 'F' because the corresponding code requires large arrays and RF is set to 'F' because it adds code to the energy calculations.

    The key PRCO produces an annotated list of the values used for the currently executing version as well as of the code segment selection choices made. A quick listing on the terminal can be obtained by typing mmc.bin, followed by the command PRCO at the MMC> prompt (followed by STOP at the next prompt for a graceful exit).

    If a compilation is needed to match a given checkpoint file the command PRCO SAVE will read the information of the current run's checkpoint file and write Fortran data statements to be inserted into pre.f

    VI.2. Compilation

    Once the mmc.f file is prepared it should be compiled with the appropriate Fortran compilation command, to produce an executable, say mmc.bin, e.g.,

    f77 -o mmc.bin mmc.f.

    Compiler options:

    Possible problems:

    For some systems, the source code is too big to compile it in one step. In this case, the program splitmmc.f (part of the distribution) should be compiled and run. This will result in five files: mmc1.f, ..., mmc5.f. Each one should be compiled separately with the -c option, i.e.,
    > f77 -c mmc1.f
    followed by
    > f77 -o mmc.bin mmc*.o

    Some compilers fail due to a so-called 'relocation error' when optimizing at levels higher than one is asked. When using the Intel Fortran compiler (ifort), adding the compiler directives
    -mcmodel=medium -share_intel
    solved the problem. With some of the other compiler (but not the GNU compiler) the compilation key -fpic was found to solve the problem.

    VI.3. Parallelization

    There are a number of limitations in the parallelized version of MMC. In particular, it is extremely inefficient (due to latency issues) on distributed memory systems so it should be run either on shared-memory systems, or limit the number of CPU's to the nmumber of procesors on a single node and make sure the the queuing system places each job on the same node.

    The following has been parallelized: solvent-solvent energy calculations; solute-solvent energy calculations; Widom-insertion free-energy calculations, and grand-canonical ensemble (GCE) simulations. However, force-calculations have to be disabled. Also, use any solute move (keys PARD or PART) would reduce the parallel efficiency.

    VIII. TESTING THE PROGRAM

    A test suite has been compiled to allow a systematic check on the integrity of the program. The tests are located in the directory testsuite. There is a c-shell script test.csh. It runs all the tests and examines the output. For each test, it prints different messages if there were input errors, inconcistency errors or no errors. For runs with special solute moves, it also prints the acceptance rates. The program is considered passing the test if there are no errors and the acceptance rates printed are all positive. Furthermore, if there is a subdirectory oldout containing reference output files, the new output will be compared to the old one. Differences should only show in CPU times and dates. Clearly, this is only meaningful if the reference output was generated by the same compiler and compilation options. The table below gives a summary of the features tested by each test system. 

     Nmolslt Nwat  PARD PART SWAP TICA PRM PRM GCEN IBEN SUUCm CLON PMLI PROX WIDO
                                       Ft  Dr
1   1     215
2   1     215
3   3      0     x   x     x
4   7      0     x         x
5   2      15                   x
6   2     484    x   x
7   1      0         x
8   4     484                            x
9   6     484                            x
10  1      9                                 x
11  6      0     x   x     x                                 x
12  6      0     x                                           x
13  2      0         x
14  2      0         x          x
15  2     484    x   x
16  3     484    x         x
17  1      4                                 x
18  3     484    x
19  3     693                        x
20  1                x
21  2      50    x                           x
22  2      50    x                                x
23  2      50                                     x
24  2      50    x   x                            x
25  3      0     x   x     x                           x
26  3     484    x   x                     x
27  3     500                                                    x
28  2     484                                                         x
29  4     215    x   x         x
30  1      0     x   x
31  1     190        x
32  1     107                                                              x
  Nmlslt Nwat  PARD PART SWAP TICA PRM PRM GCEN IBEN SUUCm CLON PMLI PROX WIDO
                                       Ft  Dr

    IX. REFERENCES

    Grand-canonical ensemble formalism: D.J. Adams, Mol. Phys, 29, 307 (1975).

    QCDF formalism: A. Ben-Naim, Water and aqueous solutions, Plenum, New York (1974).

    Overlap ratio method: Jacucci and N. Quirke, Mol. Phys., 40, 1005 (1980)

    EPEN/QPEN potential: F. Marchese, P.K. Mehrotra and D.L. Beveridge, J. Phys. Chem., 85, 1 (1981).

    Proximity criterion formalism: P.K. Mehrotra and D.L. Beveridge, J. Am. Chem. Soc., 102,4287 (1980).

    Original description of the Metropolis method: N. Metropolis et al., J. Chem. Phys., 21, 1087 (1953),

    Shuffled cyclic selection of molecules to be moved: M. Mezei, J. Comp. Phys., 39, 128 (1981)

    Click HERE to download a copy

    Cavity-biased insertion - original formalism: M. Mezei, Mol. Phys., 47, 1307 (1982)

    Virial-biased volume change in the (T,P,N) ensemble: M. Mezei, Mol. Phys., 48, 1075 (1983)

    Description of self tests: par_all27_prot_na_crb.acly.crbn.prmM. Mezei, CCP5 Newsletter, Daresbury Laboratory, No. 23, (1986)

    Click HERE to download a copy

    Free energy methodology review: M. Mezei and D.L. Beveridge, Ann. Acad. NY Sci., 482, 1, (1986).

    Click HERE to download a copy

    Cavity-biased insertion - grid-based formalism: M. Mezei, Mol. Phys., 61, 565-582 (1987); Erratum, 67,1207-1208 (1989).

    Click HERE to download a copy

    Adaptive Umbrella sampling: M. Mezei, J. Comp. Phys. 68, 237 (1987)

    Click HERE to download a copy

    Finite-difference thermodynamic integration: M. Mezei, J. Chem. Phys. 86, 7084 (1987)

    Click HERE to download a copy

    Radical-plane based proximity criterion: M. Mezei, Molecular Simulation, 1, 327 (1988).

    Click HERE to download a copy

    Using a bitmap for keeping neighbour lists: M. Mezei, Molecular Simulation, 1, 169 (1988).

    Click HERE to download a copy

    Scaled force bias method: M. Mezei, Molecular Simulation, 5, 405 (1990)

    Click HERE to download a copy

    Polynomial Thermodynamic Integration for the solvation free energy of liquid water: M. Mezei, J. Comput. Chem., 13, (1992) 651;

    Iso-energy cutoff for ion-ion PMF: M. Mezei, Int. J. Quantum Chem., 52, 147 (1994).

    Preferential sampling near the solute: J.C. Owicki, Computer Modelling of Matter, P.G. Lykos, ed., ACS, Washington, D.C., (1978), p 159.

    Force-bias sampling: C. Pangali, Rao and B.J. Berne, Chem. Phys. Lett., 55, 413 (1978), Mol. Phys., 37, 1773 (1979)

    Quantum corrections: Powles and Rickayzen, Mol. Phys., 38, 1875 (1979).

    Polynomial TI for the solvation free energy of the alanine dipeptede: H. Resat and M. Mezei, J. Chem. Phys., 99, (1993) 6052.

    Click HERE to download a copy

    Smart MC (a force-biasing variant): P.J. Rossky, J.D. Doll and H.L. Friedman, J. Chem. Phys., 69, 4628 (1978).

    Generic solvation sites: M. Mezei and D.L. Beveridge, J. Comput. Chem., 5, 523 (1984). (1992) 651;

    Click HERE to download a copy

    Primary Hydration Shell: D. Beglov and B. Roux, Biopolymers., 35, 171 (1995).
    A. Kentsis, M. Mezei, and R. Osman, Biophys. J, 84, 805-815 (2003).

    Sampling with Tsallis statistics: I. Andricioaei and J.E. Straub, 53, 3055 (1996).

    Cavity-biased Widom insertions: P. Jedlovszky, and M. Mezei, J. Am. Chem. Soc., 122, 5125-5131 (2000).

    Extension-biased torsion moves: P. Jedlovszky, and M. Mezei, J. Chem. Phys., 111, 10770-10773 (1999).

    Click HERE to download a copy

    Local torsion moves: M. Mezei, J. Chem. Phys., 118, 3874-3879 (2003).

    Automatic tuning of the B-parameter: J.A. Speidel, J.R. Banfelder, and M. Mezei, J. Chem. Theory and Comp., , (2006).

    Automatic tuning of the stepsize parameters: J.A. Speidel, J.R. Banfelder, and M. Mezei, Algorithms, 2, 215-226 (2009).

    Simulaid: a simulation facilitator and analysis program J. Comput. Chem., 31, 2658-2668 (2010).

    Return to contents page