GMIN

File:Gmin3.jpg

300px

GMIN page - very much under construction!

Introduction

GMIN is a program that attempts to find the global potential energy minimum for a collection of atoms or molecules using the `basin-hopping' algorithm described by Wales and Doye. A constant temperature Monte Carlo (MC) run is performed on the transformed potential energy surface (PES), and the configuration point may either be reset to the latest minimum in the chain or vary freely. The program knows many different empirical potentials, and it is straightforward to add new systems. From version 2.2 basin-sampling thermodynamics has been added, and from version 2.3 parallel tempering basin-sampling and basin-hopping have been implemented with MPI.

To start a calculation you need a file called data in the current directory, along with a file called coords containing the initial coordinates, which can be random. If the SEED keyword is present in data you also need a file called seed containing the seed coordinates. Most output is written to stdout, although the file lowest is always created at the end of the run, containing the energies and geometries of the lowest few configurations found in the given run. The geometries are saved in XMakemol xyz format.

Source code

You can download the latest version of the source code for all our programs under the GPL here. This tarball does not contain the CHARMM and AMBER source. If you have a license and want to use them with GMIN, please contact Professor Wales.

Keywords

The data file

Input is keyword driven with sensible defaults in most cases. Free format may be used within each line. Blank lines are ignored.

The following keywords are recognised, where n and x are integer and real data, respectively.

• 2D: enforce two-dimensional `flatland'.

• A9DIHE resid atom1 atom2 atom3 atom4 step: specifies LeaP should be used to take a step in the dihedral for the atoms specified. It is limited to dihedrals within a residue. This keyword is not yet fully functional.

• A9INTE: specifies that after each quench that does not lead to an inversion of chirality, isomerisation of a peptide bond or cold fusion - the interaction enthalpy between a specified residue and the rest of the system should be calculated using the external script 'AMBGMINintE.sh', and read back into GMIN.

This is intended for use with protein/ligand systems where you are searching for low energy docked structures. As the total energy does not fully correlate with the protein/ligand interaction enthalpy, it is often useful to retain not only the lowest SAVE total energy structures, but also the lowest SAVEINTE interaction enthalpy structures. To use this keyword, the `AMBGMINintE.sh' script (contained in the SVN repository in the SCRIPTS directory) must be present in the GMIN working directory. You should ensure that you have edited it to match the residue numbering of your system. You also need a full AMBER9+ installation with access to the `sander' executable. When using this keyword, an interaction enthalpy dump file is produced very DUMPINT steps, and at the end of the run, structural output files are produced for the SAVEINTE lowest interaction enthalpy geometries. After each quench, the structure with the current lowest interaction energy is dumped in pdb and rst format prefixed with `bestint.' to allow monitoring.

• ACCEPTRATIO accrat: accrat is the required acceptance ratio for the MC exploration of the transformed surface. For fixed temperature runs (the default) the maximum step size is adjusted to try and meet the requested value of accrat for a fixed maximum step size the temperature is adjusted instead. The default value of accrat is a half.

• ACKLAND id: specifies an Ackland embedded atom metal potential coded by Dr Mihai-Cosmin Marinica. id specifies the particular metal: 1 is ?, 2 is ?, 3 is ?, 4 is ?, 5 is iron, 6 is a different iron, 7 is tungsten. Positive values for id specify periodic boundary conditions, where box lengths must be specified by the PERIODIC keyword. Negative values for id specify a cluster calculation. A CUTOFF value can also be used for clusters.

• ALGLUE: specifies a glue potential for aluminium.

• AMBER: specifies the AMBER force field (WARNING: old implementation). See also DIELEC.

• AMBER9 inpcrd inpcrdformat: specifies a calculation with the interfaced version of the Amber 9 program package.

From this package the Amber force fields are being used, with small modifications (e.g. smooth cut-offs). Starting coordinates are read from inpcrd instead (default coords.inpcrd), in Amber inpcrd file format specified by the second optional argument inpcrdformat. If the second argument is missing, it is assumed that inpcrd contains only three columns with the xyz coordinates of all atoms, in the same order as in the topology file. To start a run with this interface, several auxiliary files are required in the same directory: input coordinate file coords.inpcrd, parameter topology file coords.prmtop, input file to Amber containing force field specifications min.in, and, if desired, a coordinate file different from coords.inpcrd containing starting coordinates. To turn on smooth cutoffs for the Generalised Born force fields, the keyword ifswitch=1 has to be used in the \&cntrl namelist block of min.in. When using the AMBER9 keyword, any calculated second derivatives will be numerical. If one wants analytical second derivatives, the NAB keyword should be used instead, with the same syntax. Additional keywords for the AMBER 9 runs are DUMPSTRUCTURES, AMBERMDSTEPS, LIGMOVE (0.0-1.0) (x.x) and MOVABLEATOMS.

• NO KEYWORD INFO FOR AMH :( - Mike?

• AMCHNMAX: The maximum number of angles that will be changed by up to STEP during an AMBER dihedral step. If this is not set or is set to zero, cartesian steps of maximum size STEP are taken instead.

• AMCHNMIN: The minimum number of angles that will be changed during an AMBER dihedral step.

• AMCHPMAX: The maximum probability for a single angle to be twisted in an AMBER dihedral step.

• AMCHPMIN: The minimum probability for a single angle to be twisted in an AMBER dihedral step.

• ANGSTROM: specifies coordinates in $\mathrm {\AA}$ngstrom for the FRAUSI potential.

• ARGON: introduces a diatomics-in-molecules calculation for a neutral, cationic or electronically excited argon cluster. See also GROUND, PLUS, TWOPLUS and STAR.

• ARM arma armb: use the acceptance-ratio method (Bouzida et al., Phys.Rev.A, 45, 8894, 1992) to adjust the step size to achieve the requested acceptance ratio.

A scaling factor is calculated and applied to STEP, rotmax, and/or transmax. The scaling factor is calculated according to $\log(arma*P_{\rm {t}}+armb)/\log(arma*P_{0}+armb)$ where $P_{\rm {t}}$ defines the target acceptance ratio and $P_{0}$ the actual acceptance ratio. Both values arma and armb default to 0.4.

• ARNO: specifies a diatomics-in-molecules potential for Ar$_{N}$-NO clusters.

• AVOID dist maxsave: specifies that the geometry should be reseeded if the latest structure gets within a distance dist of the maxsave members of a cyclic list.

• AXTELL zstar: specifies an additive Axilrod-Teller term for certain diatomics-in-molecules potentials as well as the Pacheco-Ramelho intermolecular potential for C$_{60}$.zstar is the coefficient multiplying this term.

• BASIN bgmax: specifies a basin-hopping run (as opposed to standard MC on the untransformed surface). bgmax is the convergence threshold on the RMS force in the basin-hopping quenches. If this criterion is too strict then the run time will be greatly increased. If it is too sloppy then the performance of the algorithm is impaired. Different values are needed for different potentials. SLOPPYCONV can be used instead.

• BFGS: specifies that the full BFGS minimiser should be used. Inefficient compared to LBFGS.

• BHPT pttmin pttmax exchprob: specifies minimum (pttmin) and maximum (pttmax) temperatures for a parallel tempering basin-hopping run and the probability of attempting replica exchange (exchprob). Should be used together with the MPI keyword. (Only available if the source is compiled with MPI enabled.)

• BINARY ntypea epsab epsbb sigmaab sigmabb: specifies a binary Lennard-Jones system. ntypea is the number of type A atoms - the rest are assumed to be type B and appear at the end of the list of coordinates. $\epsilon _{\rm {AA}}=\sigma _{\rm {AA}}=1$ define the units of energy and length, and epsab=$\epsilon _{\rm {AB}}$, epsbb=$\epsilon _{\rm {BB}}$, sigmaab=$\sigma _{\rm {AB}}$, sigmabb=$\sigma _{\rm {BB}}$. The box parameters and cutoff should be specified with the PERIODIC keyword.

• BINSTRUCTURES SaveNth: requests that the geometry of every SaveNth new structure found during basin-sampling is recorded in binstructures.j, where j is the index of the bin to which a given minimum belongs. If this keyword is present then GMIN switches from plain PTMC to BSPT. Without BINSTRUCTURES the BSPT keyword will perform a standard PTMC run with no quenching.

• BLJCLUSTER ntypea epsab epsbb sigmaab sigmabb cutoff: specifies a binary Lennard-Jones cluster. The parameters are the same as for BINARY, above.

• BLN $k_{r}$ $k_{\theta }$: specifies a BLN off-lattice protein model with bond-length and bond-angle force constants $k_{r}$ and $k_{\theta }$. An auxiliary file BLNsequence is required.

• BSMIN: specifies a Bulirsch-Stoer minimisation scheme. Very inefficient compared to LBFGS.

• BSPT histmin histmax ptemin ptemax pttmin pttmax exchprob nequil ptsteps nquench nenrper hbins qfrq: requests a basin-sampling run to accumulate the quench probability for local minima as a function of potential energy using a parallel-tempering algorithm.

This keyword also specifies the energy range for the histogram of quench energies, histmin to histmax, the energy range for the histogram of instantaneous configurations, ptemin to ptemax, the temperature range (pttmin and pttmax), the probability of attempting an exchange exchprob, the number of equilibration steps, nequil, the number of parallel tempering MC steps without quenching, ptsteps, the number of parallel tempering MC steps with quenching, nquench, the number of bins for the histogram of instantaneous potential energy, nenrper, the number of bins for the histogram of quench energies, hbins, and the quench frequency, qfrq. Should be used together with the MPI keyword. (This option is only available if the source is compiled with an MPI enabled.)

• BSPTDUMPFRQ n: n is the interval at which intermediate statistics and BSPTRESTART files are dumped. If n is less than one these files will only be dumped at the end of a complete run. See also BSPTRESTART. Note BSPTRESTART is NOT documented.

• BSPTDUMPFRQ: restart a previous BSPT or PTMC run. The instantaneous and quench potential energy histograms are read from the last Visits.his and Visits2.his files, and the current state from bsptrestart files (one per node, numbered from zero). A finished run can be continued with more steps by changing the nquench or ptsteps parameters on the BSPT or PTMC line of the data file. Setting the interval for BSPTDUMPFRQ to

minus one will read the last set of dump files.

• BSWL: obsolete Wang-Landau basin-sampling; do not use.

• CAPSID rho epsilon radius height: specifies a coarse-grained potential to represent virus capsid pentamers with parameters $\rho$, $\epsilon _{2}$, $r$ and $height$, respectively. If $height$ is omitted the default is 0.5.

• CENTRE: if present the system will be translated so that the centre-of-mass lies at the origin after every quench.

• CENTREXY: if present the system will be translated so that the centre-of-mass lies at the centre of the xy plane after every quench. This is useful when using an implicit membrane like IMM1 where you have directionality only in the z-direction, so centreing in x and y should have no delaterious effect.

• CG: specifies a conjugate-gradient minimisation scheme. Inefficient compared to LBFGS.

• CHANGEACCEPT naccept: naccept is an integer which sets the interval at which the acceptance ratio is checked and possible adjustments are made to the maximum step size or the temperature. The default is naccept$=50$.

• CHARMM: specifies that a CHARMM potential should be used. See also keywords CHARMMTYPE, CHPMAX, CHPMIN, CHNMAX, CHNMIN, NOPHIPSI, TOMEGA, INTMIN, CHFREQ, CHRIGIDROT, CHRIGIDTRANS, and RMS. If CHNMAX is not specified, a cartesian displacement step taking scheme will be used. For cartesian steps, rings are moved as rigid bodies to avoid false knotted minima. See RINGROTSCALE. Finally, Molecular Dynamics (MD) can be employed to generate new geometries. See CHMD below.

• CHMD chmdfreq: Requests Molecular Dynamics (MD) runs to be performed every chmdfreq steps to generate new geometries. A chmdfreq setting of 20 will execute an MD run every 20$^{\mathrm {th} }$ step, while dihedral or cartesian moves are applied otherwise as specified in the data file. A CHARMM parameter file named chmd.par containing all relevant keywords for the CHARMM DYNA module has to be present in the working directory. All CHARMM keywords must be uppercase and given in the first line. A typical example is:

VERL NSTEP 500 TIMESTEP 0.002 TWINDH 10.0 IEQFRQ 200 ICHECW 1 IASORS 0 IASVEL 1 FIRS 500 FINA 500

Please consult the CHARMM manual for further details on the DYNA module. Currently, the length of the input string given in chmd.par is limited to 500 characters.

• CHARMMENERGIES: prints the components of the total CHARMM energy after each step.

• CHARMMTYPE topfile paramfile: topfile and paramfile are the common CHARMM top and param files, e.g. 'toph19\_eef1\_perm.inp' and 'param19\_eef1\_perm.inp'.

• CHFREQ nfreq: used with CHARMM keyword to specify that every nfreq basin-hopping steps dihedrals are twisted. Default is nfreq=1.

• CHNMAX: used with CHARMM keyword to specify the maximum allowed number of angles to be twisted. Specifies a dihedral angle step taking scheme.

• CHNMIN: used with CHARMM keyword to specify the minimum allowed number of angles to be twisted.

• CHPMAX: used with CHARMM keyword to specify the maximum allowed probability for twisting an angle.

• CHPMIN: used with CHARMM keyword to specify the minimum allowed probability for twisting an angle.

• CHRIGIDROT prot rotmax nrot: used with CHARMM keyword to support rigid body rotation every nrot basin-hopping steps with maximum allowed probability prot and maximum allowed rotation angle rotmax (in degrees). The keyword CHRIGIDROT requires a file segments.tomove, which specifies the segments for rigid rotation. The segments are numbered and each line contains only one number.

• CHRIGIDTRANS ptrans transmax ntrans: used with CHARMM keyword to support rigid body translation every ntrans basin-hopping steps with maximum allowed probability ptrans and maximum allowed translation transmax (in $\mathrm {\AA}$). The keyword CHRIGIDTRANS requires a file segments.tomove, which specifies the segments for rigid translation. The segments are numbered and each line contains only one number. CHRIGIDROT and CHRIGIDTRANS use the same segments.tomove.

• CISTRANS: disables all checks for cis or deformed amide/peptide bonds.

• COLDFUSION thresh: if the energy falls below threshold thresh then cold fusion is assumed to have occurred and geometry optimisation stops. The default value is $-10^{6}$.

• COMPRESS comp: add a harmonic compression potential with force constant comp using the centre-of-mass distance for each atom.

• COMMENT: the rest of the line is ignored.

• COOPMOVE n cut: specifies cooperative moves in the step-taking routine. An atom is selected at random, and the n nearest neighbours (default 5) that lie within a cutoff distance of cut (default 1.0) are moved by the same amount.

• CPMD sys: specifies that the CPMD program should be called for energies and gradients. Not tested!

• CUTOFF cutoff: sets a cutoff beyond which the potential is truncated. This only has an effect for tight-binding silicon at present.

• DBRENT specifies minimisation using Brent's method with first derivatives in the conjugate-gradient procedure. Inefficient compared to LBFGS.

• DEBUG: sets various debug printing options including the dumping of initial geometries and energies (to dump.X.xyz) if DUMP is also set.

• DECAY x: magnitude of random move decays according to parameter x with distance from a randomly chosen atom.

• DF1: specifies a binary 2D potential. The first $N/2$ atoms have unit radius and the rest have radius 1.4, with a cutoff for each pair type at the average radius. The keyword 2D must also be specified, along with a PERIODIC line to specify two box-lengths. Initial work uses box lengths of 3.31437171 for a number density of 0.9.

• DFTB: specifies a DFT-based tight-binding potential; the multiplicity is specified by keyword MULTIPLICITY.

• DGUESS dguess: initial guess for diagonal elements of the inverse Hessian, used whenever the LBFGS optimiser is reset. The default is dguess=0.1.

• DIELEC dparam: specifies dielectric constant for AMBER - old implementation.

• DIPOLES: causes the first order induction energy to be included in a diatomics-in-molecules calculation for Ne$_{n}^{+}$ or Ar$_{n}^{+}$. By default this term is neglected, although it may be significant.

• DUMP: if present will cause the energy and quench geometry for every step to be dumped into dump.X.xyz where X is an integer. The geometries are saved in XMakemol xyz format. If CHARMM is also specified, dump.pdb and dump.crd are produced containing each quench geometry in PDB and CHARMM CRD format.

• DUMPINT int: changes the default interval for dumping a restart GMIN.dump file from 1000 basin-hopping steps to int.

• DUMPQU: when using AMBER9, dumps each quench geometry in rst format to quenchX.rst and pdb format to quenchX.pdb. Dumping does not occur if a chirality check fails.

• DZUGUTOV dzp1 dzp2 dzp3 dzp4 dzp5 dzp6 dzp7: Dzugutov potential in a general form. The parameters are $m$, $A$, $aa$, $B$, $d$, $bb$ and $m2$.

• EAMAL: specifies an embedded atom model for aluminium.

• EAMLJ A0 beta Z0: specifies the EAMLJ potential (Baskes, Phys.Rev.Lett., 27, 2592, 1999) with parameters A0, beta and Z0.

• EDIFF econv: quench minima are only considered to be different if their energies differ by at least $econv$. This option mainly affects the lowest energy saved geometries. If the current quench energy is within $econv$ of a saved energy, but lies lower, then the saved energy and geometry are replaced. The default is $0.02$ but different values are appropriate for different potentials.

• EQUILIBRATION equil DumpEveryNthQuench: equil is the number of MC steps preceding the accumulation of the density of states histogram in a Wang-Landau basin-sampling run. The default is 0. DumpEveryNthQuench specifies how often the statistics are recorded into the output files.

• FAKEWATER: specifies a distance-dependent dielectric in AMBER - old implementation.

• FIXEDEND: requires documentation :(.

• FAL: specifies the Farkas potential for aluminium.

• FIXBOTH: both the temperature and maximum step size are fixed regardless of the calculated acceptance ratio.

• FIXCOM: fix centre of mass rather than centre of coordinates.

• FIXSTEP: the maximum step size is fixed and the temperature is varied to try and achieve the requested acceptance ratio.

• FIXTEMP: explicitly fixes the temperature. Only used if FIXSTEP is set, in which case using FIXTEMP gives a result equivalent to FIXBOTH.

• FNI: specifies the Farkas potential for nickel.

• FRAUSI: specifies a particular tight-binding potential for silicon. See also keyword ANGSTROM.

• FREEZE n1 n2 $\ldots$: freeze the coordinates of atoms n1, n2,$\ldots$}. Note that all the FREEZE keywords do NOT function when using AMBERMDSTEPS.

• FREEZERES n1 n2 $\ldots$: freeze the coordinates of all atoms in residues n1, n2,$\ldots$.

• FREEZEALL n1 n2 $\ldots$: freeze the coordinates of all atoms.

• FS gatom: specifies a Finnis-Sinclair potential using parameters from Finnis and Sinclair, Phil.Mag.A, 50, 45 (1984) and corresponding erratum Phil.Mag.A, 53, 161 (1986). gatom=1 for V, gatom=2 for Nb, gatom=3 for Ta, gatom=4 for Cr, gatom=5 for Mo, gatom=6 for W, gatom=7 for Fe (original parameters), gatom=8 for Fe (modified parameters in erratum). Subtoutine FS was coded by James Elliott in April 2009.

• GROUND: when combined with keywords NEON or ARGON uses an accurate (Aziz) potential to model the ground state neutral cluster.

• GUIDE guidecut: specifies the RMS force below which the real potential is used rather than a guiding potential. The systems affected are CPMD and WELCH, which are guided by AMBER and TOSI, respectively, and also PACHECO, where the Axilrod-Teller contribution is only included when the RMS force falls below guidecut. Default guidecut=0.0001. New guided potentials are ZETT1 and ZETT2 (guided by Morse with $\rho =5$) and NATB (guided by GUPTAT). Parameters for the guiding potential must also be specified in data.

• GUPTA gatom: specifies a Gupta potential using parameters from Cleri and Rosato, Phys.Rev.B, 48, 22 (1993). gatom=1 for Ni, gatom=2} for Cu, gatom=3 for Rh, gatom=4 for Pd, gatom=5 for Ag, gatom=6 for Ir, gatom=7 for Pt, gatom=8 for Au, gatom=9 for Al, gatom=10 for Pb, gatom=11 for Ti type 1, gatom=12 for Ti type 2, gatom=13 for Zr type 1, gatom=14 for Zr type 2, gatom=15} for Co, gatom=16 for Cd type 1, gatom=17 for Cd type 2, gatom=18 for Zn, gatom=19 for Mg, gatom=20 for V, gatom=21 for Na, gatom=22 for Sr (Wang and Blaisten-Barojas, J.Chem.Phys., 115, 3640 (2001)), gatom=? for Au as used by Garzon et al. The Gupta subroutine was recoded more efficiently by James Elliott in April 2009.

• INTMIN: used with CHARMM keyword to specify minimisation in internal coordinates. This generally appears to be slower than using Cartesian coordinates.

• JC: Specifies Murrell's two-plus three-body potential. A file {\tt JMparams} must exist in the current directory containing the parameters $c_{0}$, $c_{1}$, $\ldots$, $c_{10}$, $r_{e}$, $D$, $a_{2}$ and $a_{3}$. An optional cutoff parameter can also be provided at the end of the JMparams file. Subroutines used: jmec, jm2c, jm3c.

• JUMPMOVE np1 np2 int: specify J-walking type attempts between parallel runs np2 and np1 at intervals of int steps.

• LB2: specifies the potential $V={\frac {\epsilon }{2}}\sum _{i<j}\left[\left({\frac {r_{ij}}{\sigma }}\right)^{2}+\left({\frac {\sigma }{r_{ij}}}\right)^{2}\right]$, where $\epsilon$ and $\sigma$ are set to unity.

• LIGMOVE ligrotscale ligcartstep: used with AMBER9 and MOVABLEATOMS. Specifies ligand only rotation and cartesian perturbations. The ligand is defined by atom index in the file movableatoms. Setting ligrotscale less than 1.0 limits the ammount of rotation possible - this may be required to prevent cold fusion with non-spherical ligands. ligcartstep defines the maximum size (in angstroms) of the random cartesian perturbations applied to the ligand. Both rotation and cartesian ligand moves are applied AFTER any MD if AMBERMDMOVES is on to prevent the MD exploding.

• LJCOUL nc $q'$ f $T_{\rm {swap}}$: specifies a cluster of Lennard-Jones particles in which the first nc particles carry identical reduced charges $q'$ in addition to the Lennard-Jones interaction. The parameter f specifies what fraction of the Monte Carlo steps should be swaps between the positions of a charged and a neutral particle, rather than a conventional step. $T_{\rm {swap}}$ is the temperature to be used in the acceptance criterion for swap moves, overriding that specified using the TEMPERATURE keyword. Generally, a lower temperature is more effective at finding the lowest-energy permutation of charges. The default value of $T_{\rm {swap}}$ is zero. The reduced charge $q'$ is related to the actual charge $q$ by $q'=q/(4\pi \epsilon _{0}\epsilon \sigma )^{1/2}$, where $\epsilon$ and $\sigma$ are the Lennard-Jones well depth and length parameter respectively. This way, the reduced energy of two charges is $E'=q'^{2}/r'$, where $r'=r/\sigma$ is the reduced distance between the charges.

• LOCALSAMPLE abthresh acthresh: Keyword currently under construction! For three groups of atoms defined in movableatoms (A,B,C), a step is quenched when the centre of coordinates of A->B is less than abthresh AND A->C is less than acthresh. If this condition is broken AFTER the quench, it is automatically rejected.

• MAXBFGS max: max is the largest permitted LBFGS step.

• MAXERISE maxez: specifies the largest rise in energy permitted during an LBFGS minimisation, default $10^{-10}$. Useful for potentials with imprecise derivatives.

• MAXIT maxit maxit2: maxit and maxit2 are integers specifying the maximum number of iterations allowed in the conjugate gradient quenches. maxit applies to the 'sloppy' quenches of the basin-hopping run and maxit2 to the final quenches that are used to produce the output in file lowest.

• MGGLUE: specifies a glue potential for magnesium.

• MORSE rho: specifies a Morse potential with range parameter rho.

• MPI: specifies an MPI parallel job. (only available if the source is compiled with MPI enabled).

• MSORIG: specifies a particular tight-binding potential for silicon.

• MSTRANS: specifies an alternative tight-binding potential for silicon.

• MULLERBROWN: specifies the 2D Muller-Brown potential.

• MULTIPLICITY xmul: specifies the multiplicity of the electronic state in DFTB calculations.

• NATB: specifies the sodium tight-binding potential of Calvo and Spiegelmann. This potential can be guided by also specifying GUPTA 21 in the data file.

• NEON: introduces a diatomics-in-molecules calculation for a neutral, cationic or electronically excited neon cluster. See also GROUND, PLUS, TWOPLUS and STAR.

• NEWJUMP prob: for (serial) 'parallel' runs specifies a jump probability between runs (parallel tempering) of prob. See the BHPT keyword for a better alternative.

• NEWRESTART nrelax nhs MD newrestemp: reseed runs if the energy does not decrease within nrelax steps. nhs is the number of hard sphere moves used to produce the new starting configuration. If nhs=0 (the default) then the geometry is changed by reseeding. If MD is present, then a short AMBER or CHARMM MD run is performed at temperature newrestemp to generate the new configuration.

• NMAX nmax: nmax is an integer that specifies the maximum number of dihedral angles to be twisted in AMBER - old implementation.

• NMIN nmax: nmin is an integer that specifies the minimum number of dihedral angles to be twisted in AMBER - old implementation.

• NOCHIRALCHECKS: disables checks for inversion of CA atoms and chiral side-chains for ILE and THR.

• NOCISTRANS minomega: set on by default with a threshold minomega of 150 degrees. If an amide bond is deformed to a angle below the specified threshold, the structure is discarded. i.e. with $|\omega |<$minomega. minomega defaults to 150 degrees. For proline every $\omega$ is allowed. To enable cis-trans isomerisation, the CISTRANS keyword should be used.

• NOCISTRANSDNA minomega: should be specified when working with DNA in AMBER to ensure the correct bonds are checked. As above, the deformation threshold minomega can be set. It is defaulted to 150 degrees.

• NOCISTRANSRNA minomega: should be specified when working with RNA in AMBER to ensure the correct bonds are checked. As above, the deformation threshold minomega can be set. It is defaulted to 150 degrees.

• NOFREEZE: don't freeze the core atoms when doing the initial geometry optimisations in a run where SEED is specified.

• NOPHIPSI: used with the CHARMM keyword to specify twisting of sidechain dihedrals only.

• NORESET: by default the configuration point is set to that of the quench minimum in the Markov chain during a basin-hopping simulation. This keyword turns off the resetting so that the geometry varies continuously.

• NOTE: the rest of the line is ignored.

• ODIHE: order parameter - requires documentation.

• OEINT: interaction energy between 2 peptides will be used as an order parameter - requires further documentation.

• ORGYR: radius of gyration will be calculated as an order parameter - requires further documentation.

• OSASA: order parameter - requires documentation.

• P46: specifies a 46-bead three-colour model polypeptide. See also the BLN keyword, which implements this potential in a more general way and uses unit bond lengths.

• PACHECO: specifies the intermolecular Pacheco-Ramelho potential for C$_{60}$. The Axilrod-Teller contribution, specified with the AXTELL keyword, is included when the RMS force falls below the value entered with guidecut.

• PAH: specifies a polycyclic aromatic hydrocarbon potential.

• PARALLEL npar: npar is the number of parallel runs within GMIN.

• PBGLUE: specifies a glue potential for lead.

• PERIODIC boxlx boxly boxlz: specifies periodic boundary conditions for potentials which understand such a directive (such as tight-binding silicon). The three double precision variables are the box lengths. If only one box length is given the others are set to the same value to give a cube.

• PERMDIST: minimise distances between the coordinates in files coords and the fixed coordinates in file finish with respect to permutational isomerisation. Requires the auxiliary file perm.allow to specify permutable atoms, otherwise all atoms are assumed to be permutable. The absence of a perm.allow file is considered a mistake for CHARMM runs. The first line of the perm.allow file must contain an integer that specifies the number of groups of interchangeable atoms. The groups then follow, each one introduced by a line with two integers that specify the number of permutable atoms and the number of other pairs of atoms that must swap if the first pair are permuted. The following line contains the numbers of the permutable atoms and then the numbers of the swap pairs, in order. For CHARMM runs the perm.allow file can be generated automatically from the input.crd file using the python script perm.py written by Dr Mey Khalili. Note that it is then essential to use CHARMM topology files that include symmetrised definitions of all the relevant amino acids if

energies are actually calculated.

\item {\it PLUS\/}: when combined with keywords {\it NEON\/} or {\it ARGON\/} uses a diatomics-in-molecules potential for the singly charged cation.

\item {\it PMAX max\/}: {\it max\/} is the maximum probability of dihedral angle twisting for {\it AMBER\/}.

\item {\it PMIN min\/}: {\it min\/} is the minimum probability of dihedral angle twisting for {\it AMBER\/}.

\item {\it POWER ipow\/}: {\it ipow\/} is the initial power for shifts in the old line minimisation routine for conjugate gradient. LBFGS minimisation should be used instead.

\item {\it PROJI\/}: turns on projection operator to enforce $I$ point group symmetry in {\bf mylbfgs.f}. The geometry is projected after every proposed step.

\item {\it PROJIH\/}: turns on projection operator to enforce $I_h$ point group symmetry in {\bf mylbfgs.f}. The geometry is projected after every proposed step.

\item {\it PTMC histmin histmax ptemin ptemax pttmin pttmax exchprob nequil ptsteps nenrper hbins\/}: requests a standard parallel tempering MC run. This keyword also specifies the energy range for the histogram of quench energies, {\it histmin\/} to {\it histmax\/}, the energy range for the histogram of instantaneous configurations, {\it ptemin} to {\it ptemax}, the temperature range ({\it pttmin} and {\it pttmax}), the probability of attempting an exchange {\it exchprob}, the number of equilibration steps, {\it nequil}, the number of parallel tempering MC steps without quenching, {\it ptsteps}, the number of bins for the histogram of instantaneous potential energy, {\it nenrper}, and the number of bins for the histogram of quench energies, {\it hbins}. Should be used together with the {\it MPI\/} keyword. % and {\it BINSTRUCTURES\/} keywords. (This option is only available if the source is compiled with an MPI enabled.)

\item {\it PULL a1 a2 f\/}: apply a static force to the potential, equivalent to adding the term $V_{\rm pull}=-f(z_{a1}-z_{a2})$. Here $z_{a1}$ and $z_{a2}$ are the $z$ coordinates for atoms $a1$ and $a2$, and $f$ specifies the force. This potential is designed to simulate a pulling experiment with static force where a molecule is pulled along the $z$ axis from atoms $a1$ and $a2$.

\item {\it QCUTOFF qcut\/}: {\it qcut\/} is a distance cut-off for Coulomb interactions in {\it AMBER\/}.

\item {\it QMAX cgmax\/}: {\it cgmax\/} is the tolerance for the RMS force in the final set of quenches that are used to produce the output for file {\tt lowest}. The default is {\it cgmax\/}$=10^{-3}$, but the appropriate value depends upon the system in question. {\it TIGHTCONV} can be used instead.

\item {\it QUAD\/}: requires documentation.

\item {\it QUCENTRE\/}: sets the centre of coordinates to the origin (0,0,0) before each MC step is taken (so after each quench), but not during the minimisation itself unlike {\it CENTRE}.

\item {\it RADIUS radius\/}: sets the radius of the container that prevents particles evaporating during quenches. If unset the program calculates an appropriate value based upon the volume per particle for close-packed material and the known pair equilibrium distance for the given potential. The formula employed is $$ RADIUS=r_e\left[1 + \left(3 n \over 4\pi\sqrt{2}\right)^{1/3}\right], $$ where $n$ is the number of atoms and $r_e$ is the pair equilibrium separation.\cite{kittel76} The `1' in this formula is to allow some extra space for more open structures.

\item {\it RANDOMSEED\/}: specifies that the random number generator should be seeded with system time after each quench, allowing simple parallel use. Currently functional only for the CHARMM and AMBER potentials.

\item {\it RANSEED i\/}: integer seed for the random number generator.

\item {\it RCUTOFF rcut\/}: {\it rcut\/} is a distance cut-off in {\it AMBER\/}.

\item {\it RESIZE resize\/}: all the coordinates are multiplied by {\it resize\/} after they have been read in, before any other operations are performed. This command is useful for scaling results obtained with one potential for a system with a different pair equilibrium distance.

\item {\it RESTART\/ nrelax nhs\/}: reseed runs if a step in not accepted within twice {\it nrelax} steps. {\it nhs} is the number of hard sphere moves used to produce the new starting configuration. If {\it nhs=0} (the default) then the geometry is changed by reseeding.

\item {\it RESTORE\/ dumpfile intEdumpfile\/}: restore a previous {\tt GMIN} run from a {\it dumpfile}. The number of basin-hopping steps performed will be the difference between the number requested for the run that produced the dumpfile, minus the number that were completed at the point the dumpfile was created. This option is not available before version 2.3. If you are using the {\it A9INTE\/} keyword, you can specify the interaction enthalpy dump file to restore from as a second arguement.

\item {\it RGCL2\/}: specifies a DIM rare gas-Cl$_2$ potential.

\item {\it RINGROTSCALE factor\/}: when applying cartesian moves with CHARMM, amino acid rings are moved as rigid units. Setting {\it factor} (default 0.0) between 0.0 and 1.0 will apply a random rotation to these rings during step taking. The suggested value is 0.1 to prevent the regular formation of high energy structures.

\item {\it RKMIN\/}: specifies a Runga-Kutta minimisation scheme. Very inefficient.

\item{\it RMS rmslimit rmstol rmssave lca best}: used with {\it CHARMM} keyword to specify that the RMSD compared to a reference geometry is calculated. The reference geometry must be given in xyz-format in an additional file {\tt compare}. {\it rmssave} is an integer that specifies the number of lowest energy geometries and RMSD $\le$ {\it rmslimit} to save. Geometries are only saved if their RMSD's are more than {\it rmstol} different. The flag {\it lca} controls whether the all-atom RMSD ({\it lca}=0) or the $C_{\alpha}$-RMSD ({\it lca}=1) is calculated. The flag {\it best} determines which structure is compared to the reference after each quench. {\it best}=0 implies the current quench minimum and {\it best}=1 implies the current best (lowest energy) minimum. If {\it TRACKDATA} is also specified, the RMSD calculated after each quench is produced in the file `rmsd' in gnuplot readable format.

\item {\it SAVE nsave\/}: {\it nsave\/} is an integer that specifies the number of lowest energy geometries to save and summarise in the file {\tt lowest}. Arrays are now dynamically allocated, so any positive integer can be specified.

\item {\it SAVEINTE nsaveinte\/}: {\it nsaveinte\/} is an integer that specifies the number of lowest interaction enthalpy geometries to save and summarise in the file {\tt intelowest}. See {\it A9INTE\/}.

\item {\it SETCENTRE x y z\/}: Sets the centre of mass/coordinates (before the initial quench) to ({\it x,y,z\/}). For example, {\it SETCENTRE 0.0 0.0 0.0\/} would translate the centre of mass to the origin.

\item {\it SC nn mm sig sceps scc\/}: specifies a Sutton-Chen potential\cite{suttonc90} with parameters $n=${\it nn\/}, $m=${\it mm\/}, $a$={\it sig\/}, $\epsilon=${\it sceps\/} and $c=${\it scc\/}.

\item {\it SEED nsstop\/}: if the {\it SEED\/} keyword appears then the program looks for a file {\tt seed} containing coordinates, which are used to `seed' the new run. The number of coordinates given in this file should be no more than one less than the number given in {\tt coords}. The specified coordinates are frozen from the first step until step {\it nsstop\/}.

\item {\it SHIFTCUT\/}: specifies a shifted-truncated potential for bulk binary Lennard-Jones.

\item {\it SIDESTEP smax\/}: specifies the maximum step in Cartesian coordinates for side-chains in {\it AMBER\/}.

\item {\it SLOPPYCONV bgmax\/}: specifies a basin-hopping run (as opposed to standard MC on the untransformed surface). {\it bgmax\/} is the convergence criterion for the RMS force in the basin-hopping quenches. If this criterion is too strict then the run time will be greatly increased. If it is too sloppy then the performance of the algorithm is impaired. Different values are needed for different potentials. {\it BASIN} can be used instead.

\item {\it SORT}: for pairwise potentials the atoms can be sorted from most to least strongly bound. The {\it SORT} keyword enables this sorting for the coordinates printed in file {\tt lowest}. This can be useful for seeding subsequent runs by removing the most weakly bound atoms. This sort is not set by default and is meaningless if the pair energies are not computed.

\item {\it STAR}: specifies an excited state calculation for Ar$^*_n$ or Ne$^*_n$ for a diatomics-in-molecules potential when used with {\it NEON\/} or {\it ARGON\/}.

\item {\it STEP step astep ostep block\/}: specifies the maximum step sizes. {\it step\/} is for the maximum change of any Cartesian coordinate and {\it astep\/} specifies a tolerance on the binding energy of individual atoms (if available, i.e.~for Morse and LJ) below which an angular step is taken for that atom. See the following section for more details. {\it ostep\/} is the maximum displacement of an axis-angle coordinate for a rigid body system and {\it block\/} (an integer) is the block size for which separate translational and orientational displacements will be made for rigid bodies. Omitting {\it block\/} or using a value of zero results in translational and orientational steps being taken simultaneously for rigid bodies. The default values for {\it step\/}, {\it astep\/} and {\it ostep\/} are all 0.3 and the default value of {\it nblock\/} is zero.

\item {\it STEEREDMIN smink sminkinc smindiststart smindistfinish sminatoma sminatomb\/}: specified steered minimisation should be performed (must be used with {\it AMBER9}). For a protein/ligand system, this adds a translation to the MC move. The vector between the centre of coordinates of groups A and B (as defined in the file movableatoms) is calculated and set to {\it smindiststart}. During the following minimisation, a restoring force is applied to the ligand. The harmonic force constant is initially zero, and rises by {\it sminkinc} every LBFGS step up to a maximum of {\it smink}. The force is applied until the A->B distance is less than {\it smindistfinish}.

\item {\it STEPS mcsteps tfac\/}: determines the length of the basin-hopping run through the integer {\it mcsteps\/} and the annealing protocol through the real variable {\it tfac\/}. The temperature is multiplied by {\it tfac\/} after every step in each run.

\item {\it STICKY nrbsites, sigma\/}: specifies a `sticky patch' potential with {\it nrbsites} sites in the rigid body reference and a value of {\it sigma} for the $\sigma$ parameter.

\item {\it STOCK mu lambda}: specifies a Stockmeyer potential with parameters $\mu$ and $\lambda$, respectively.

\item {\it STRAND}: specifies a system of $\beta$ strands coded using the rigid body formalism.

\item {\it SW\/}: specifies the Stillinger-Weber Si potential.

\item {\it SYMMETRISE int tol1 tol2 tol3 tol4 tol5 qmax mdiff d}: specifies that the symmetrisation routine should be called every {\it int} steps. The five {\it tol} parameters are tolerances for various parts of the routine: {\it tol1} is used in {\bf ptgrp.f} in defining orbits; {\it tol2} is the distance tolerance used in {\bf ptgrp.f} to define point group symmetry operations; {\it tol3} is the maximum relative difference in principal moments of inertia used to diagnose point groups with degenerate irreducible representations in {\bf ptgrp.f}; {\it tol4} is the distance cutoff used to determine if a symmetry element has been lost in {\bf symmetry.f}. Since we are dealing with approximate symmetries, this parameter may be larger than {\it tol2}. It is compared to the largest atomic displacement divided by the corresponding radius for the closest permutation. {\it tol4} is also used to test whether atoms lie on a given symmetry element, and in testing whether orbits generated from `floaters' are actually contained in the core. {\it tol5} is generally to check for atom clashes in {\bf symmetry.f}, including analysis of missing sites in orbits, as well as overlap between orbits generated from `floaters' and previous core or new orbit sites. {\it qmax} is the maximum number of quenches allowed for each call to {\bf symmetry.f}. {\it mdiff} is used to test whether a generated symmetry operation is new. If any of the nine components of the corresponding $3\times3$ matrix differs by more than {\it mdiff} from an existing matrix then the operations are considered to be different. {\it d} is the exponential factor used in constructing a centre of mass that is biased towards core atoms. The contribution of each atom is weighted by $\exp(-dx(i))$, where $x$ is the centre of mass distance of atom $i$ on the previous cycle.

\item {\it TABOO nlist\/}: specifies a taboo list of the {\it nlist\/} lowest minima should be maintained.

\item {\it TARGET target1 target2 $\cdots$\/}: specifies any number of target energies. The current run stops in an orderly fashion if the current quench energy is within {\it econv\/} of any target (see {\it EDIFF\/}).

\item {\it TEMPERATURE temp\/}: defines the temperature, {\it temp\/}, at which the MC runs are conducted. Different values can be specified for serial `parallel' runs if {\it PARALLEL} is set. For true parallel basin-hopping use the {\it BHPT\/} keyword and omit {\it TEMPERATURE\/}.

\item {\it TETHER hdistconstraint hwindows ExtrapolationPercent lnHarmFreq}: requests a calculation of the vibrational density of states for a given minimum. {\it hdistconstraint} is the minimised average radius of the basin of attraction to which the minimum belongs, {\it hwindows} is the number of potential energy windows into which a WL simulation is split. {\it ExtrapolationPercent} is the percentage of the whole potential energy spectrum, for which the density of states is estimated from the harmonic approximation and not sampled. {\it lnHarmFreq} the log product of positive Hessian eigenvalues.

\item{\it THOMSON q\/}: specify the Thomson problem for unit charges on a sphere. If {\it q\/} is present it is taken to be the charge on one particle, which can therefore be different from all the other unit charges and is read as a real number.

% Doesn't appear to be coded? % \item {\it THRESHOLD \/}: specifies threshold acceptance of steps. The change in potential energy must be % less than the value of the {\it TEMPERATURE\/} variable for a step to be accepted.

\item {\it TIGHTCONV cgmax\/}: cgmax is the tolerance for the RMS force in the final set of quenches that are used to produce the output for file {\tt lowest}. The default is {\it cgmax\/}$=10^{-3}$, but the appropriate values depend upon the system in question. {\it QMAX} can be used instead.

\item {\it TIP n\/}: specifies a TIP{\it n\/}P intermolecular potential for rigid body water molecules. $\ \le n \le 5$.

\item {\it TOLBRENT tolb\/}: parameter for {\it DBRENT\/} minimisation. Inefficient compared to LBFGS.

\item{\it TOMEGA}: used with the {\it CHARMM} keyword to specify that peptide bonds will be twisted along with all other dihedrals.

% \item {\it TN\/}: specifies a truncated Newton minimisation scheme. % Inefficient compared to LBFGS.

\item {\it TOSI app amm apm rho\/}: specifies the Tosi-Fumi potential\cite{tosif64} with parameters $A_{++}$, $A_{--}$, $A_{+-}$ and $\rho$.

\item {\it TRACKDATA}: produces `energy.dat' and `markov.dat' containing the quench number and associated energy and markov energy in two columns and `best.dat', containing the current quench number and the current lowest total energy. If {\it RMS\/} is also specified, a file called `rmsd.dat' is produced containing the RMSD from a reference structure. See {\it RMS\/} for more information. This allows plotting with gnuplot to monitor convergence of multiple runs. If {\it A9INTE} is also specified, two additional output files are produced, `intE.dat' containing the quench number and associated interaction enthalpy, and `bestintE.dat' containing the quench number and current lowest interaction enthalpy. This keyword does not yet function for MPI runs.

\item {\it TSALLIS q\/}: specifies that steps are accepted/rejected using Tsallis statistics with the given value of {\it q\/}, rather than the usual Boltzmann condition.

\item {\it TWOPLUS\/}: when combined with keywords {\it NEON\/} or {\it ARGON\/} uses a diatomics-in-molecules potential for the doubly charged cation.

\item {\it UACHIRAL\/}: MUST be included when using ff03ua, the AMBER united atom forcefield unless you have disabled the checks for inverted chiral carbons

with {\it NOCHIRALCHECKS\/}. {\it UACHIRAL\/} ensures the correct impropers are used to define sidechain chirality when HB hydrogen is missing. 

• UNFREEZE n1 n2 $\ldots$: unfreeze the coordinates of atoms n1, n2,$\ldots$. Only functions in conjunction with FREEZEALL.

• UNFREEZERES n1 n2 $\ldots$: unfreeze the coordinates of all atoms in residues n1, n2,$\ldots$. Only functions in conjunction with FREEZEALL.

\item {\it UPDATES nup\/}: {\it UPDATES\/} is the number of previous steps saved in the LBFGS routine, default 4.

\item {\it VGW ljsigma ljepsilon taumaxsg taumaxfg}: Specifies use of VGW quantum quenching in place of classical minimization routines such as LBFGS. {\it ljsigma} and {\it ljepsilon} are the corresponding Lennard-Jones parameters that must be specified, and taumaxsg and taumaxfg are the maximum value of ``imaginary time $\tau$ (inverse tempertaure) for the propagation. The former pertains to the faster ``single-particle SP-VGW used for quenching during the MC runs, and the latter for the more accurate ``fully-coupled VGW used for the final quenching (analogous to the tight convergence of the LBFGS). A $\tau$ of at least 2.5 is recommended for the SP-VGW and 5.0 for the FC-VGW. A file {\it vgwdata} containing the masses (in a.m.u.) of all particles, in order of the location of their {\it xyz} coordinates in ``{\it coords} must be present (e.g. for a 38 atom Ne cluster, {\it vgwdata} will have 38 lines of ``{\it 20}). Different masses are permitted, though the current version allows for only one set of LJ parameters.

\item{\it VGWCPS on magnitude}: Specifies use of contraining potential for SP-VGW (sloppy convergence), as clusters expand during quantum quenching with decreasing mass. 1 or 0 for {\it on} corresponds to on/off, and magnitude should range from 1 to 1000, with 1 having minimal effect, 1000 being highly constrained. Default value is ``on, with magnitude 1.

\item{\it VGWCPF on magnitude}: Same as VGWCPS but for FC-VGW, used for the final, full quenching (tight convergence).

\item{\it VGWTOL magnitude}: Absolute tolerance parameter for differential equation solver used for VGW quenching. Default value is 0.0001. For highly quantum or ``stiff systems this may need to be increased, while it may be decreased for ``softer or less quantum systems to enhance speed.

\item {\it VISITPROP}: if specified the Wang-Landau convergence is governed by proportionality of visits to the current value of the modification factor, and not the histogram flatness criterion \cite{ZhouB03}.

\item {\it WELCH $A_{++}\ A_{--}\ A_{+-}\ \rho\ Q_+\ Q_-\ \alpha_+\ \alpha_-$\/}: specifies a Welch binary salt potential with the parameters indicated.

\item {\it ZETT1\/} and {\it ZETT2\/}: specify the Zetterling potentials. \end{itemize}