1. Nmode Usage: nmode [-O] -i nmdin -o nmdout -c inpcrd -p prmtop -r restrt -ref refc -v vecs -l lmode -t tstate -e expfile -O: Overwrite output files if they exist. 1.1. Introduction This program performs molecular mechanics calculations on proteins and nucleic acids, using first and second derivative information to find local minima, transition states, and to perform vibrational analyses. It is designed to read the prm- top and inpcrd files from the Amber suite of programs. Both Additive and Non-additive Hamiltonians are available in this version. There are accompanying programs nmanal (normal mode analysis) and lmanal (Langevin mode analysis) that use the out- put of these programs to compute molecular fluctuations and time correlation functions. Nmode was originally written at the University of California, Davis, by D.T. Nguyen and D.A. Case, based in part on code in the Amber 2.0 package. Major revisions were made at the Research Institute of Scripps Clinic by J. Kottalam and D.A. Case. M. Pique has provided valuable advice and help in porting it to many different machines. J.W.Caldwell inplemented the non-additive capabilities. 1.2. References The second derivative routines are based on expressions used in the Consistent Force Field programs; [1] similar information is given by K.J. Miller, et al., [2] although these expressions were not actually used in writing this code. The code also contains routines to search for tran- sition state, starting (generally) from a minimum. This proce- dure uses a modification of the procedure of Cerjan and Miller ____________________ [1] S.R. Niketic and K. Rasmussen, The Consistent Force Field: A Documentation, Springer-Verlag, 1977. [2] R.J. Hinde and J. Anderson, J. Comput. Chem. 1989, 10, 63. [3] as described elsewhere. [4] Langevin modes are analogous to normal modes, but in the presence of a viscous coupling to a continuum solvent. The basic ideas are presented by Lamm and Szabo, [5] and were implemented in the Amber environment by us. [6] 1.3. General description This program performs five tasks, depending on the value of the input variable ntrun (see below): (1) Perform a normal mode analysis from starting coordi- nates. Requires an input structure that has already been minimized, from process (4), below, or by some other method. In addition to the computation of normal mode frequencies, thermodynamic parameters are calcu- lated. (2) Search for transition state, starting (generally) from a minimum. See the references above for a detailed description of the method. (3) Perform a conjugate gradient minimization from the starting coordinates. This routine uses an IMSL library routine for this purpose, which is not supplied with this program. Persons who do not have access to the IMSL library should probably use the AMBER "sander" pro- gram to carry out conjugate gradient minimizations. (Use the double precision version for best convergence.) (4) Does a Newton-Raphson minimization from starting coordi- nates. A constant (tlamba) is added to the diagonal elements of the Hessian matrix to make it positive defi- nite. Tlamba is chosen in a manner such that the step is always downhill in all directions. Whenever the change in energy is > emx or the rms of step length is > smx, the step length is scaled back repeatedly until the above two conditions are satisfied. Note that this rou- tine will not converge to a transition state. (5) Perform a langevin mode calculation, starting from a minimized structure. This option is similar to (1), but includes the viscous effects of a solvent in the calcu- lation. ____________________ [3] C. Cerjan and W.H. Miller, J. Chem. Phys. 1981, 75, 2800. [4] D.T. Nguyen and D.A. Case, J. Phys. Chem., 1985, 89, 4020. [5] G. Lamm and A. Szabo, J. Chem. Phys. 1986, 85, 7334. [6] J. Kottalam and D.A. Case, Biopolymers 1990, 29, 1409-1421. Input files for this program are the same as for the regu- lar AMBER minimization and molecular dynamics programs, with the exception of File 5, whose parameters are given below. The defaults have been carefully selected, so that for most pur- poses, few of them need to be changed. See the sample runs for more information. 1.4. Files nmdin : control input for the run nmdout : standard output file for print and error messages prmtop : parameter file as output by the AMBER program parm inpcrd : starting coordinates refc : input coordinates for constraints restrt : output coordinates at end of minimization prlist : file for reading or storing the non-bonded pair list vecs : file containing output normal mode frequencies and eigenvectors tstate : output coordinates at a transition state expfile: file to read exposed surface area for atoms lmode : file to write Langevin modes 1.5. Input description Input found on nmdin: You can use as many title cards as you want, followed by the namelist &data, which contains the following variables. _______________________________________________________________ General flags describing the calculation _______________________________________________________________ ntrun 1: do normal mode analysis (default) 2: search for transition states 3: do conjugate gradient minimization (requires IMSL library) 4: do Newton-Raphson minimization 5: do Langevin mode analysis ibelly 1: some atoms are to be held fixed (default=0) icons 1: do constrained minimization to initial coordi- nates specified in refc. (default=0) maxcyc max. number of cycles for minimization (default=100) drms rms gradient to stop minimization (default=1.e-5) nvect number of vectors for normal mode analysis (default=0) nsave for every nsave steps the coordinates are saved. (default= 20) nprint every nprint-th step the energy will be printed ilevel if .ne. 0, then adjust second derivative matrix to put rotation and translation vectors to a high fre- quency; this can be useful if you want to perform a normal mode analysis from a not-completely-minimized structure, so that rotations and translations don't mix with the low-lying modes (default=0). ivform 0 if the normal mode eigenvectors are to be written out in unformattted form; 1 to use the formatted option (default). ntx 0 if the input coordinates are to be read in unfor- mattted form; 1 to use the formatted option (default). ntxo 0 if the output (restart) coordinates are to be written out in unformattted form; 1 to use the for- matted option (default). _____________________________________________________ Control of certain force field parameters _____________________________________________________ cut radius for non-bonded cutoff (default=99.) scnb 1-4 nonbonded scale factor (default=2.0) scee 1-4 electrostatic scale factor (default=2.0) dielc dielectric constant (default=1.0) idiel 0 for r**2 dielectric dependence (default); 1 for constant dielectric. ipol = 0 no polarization (default) = 1 include polarization modules i3bod = 0 no three-body interactions(default) = 1 readin and use specified three body interactions (see the sander input area for details) iprr 1: read in a non-bonded pair list from prlist; (default = 0) iprw 1: write out non-bonded pair list to prlist; (default = 0) _____________________________________________________ control of Newton-Raphson and transition-state searches _____________________________________________________ smx maximum rms step length (default = 0.08) emx maximum energy change per step ( default =0.3) alpha scale factor for step length (default = 0.8) (See Nguyen and Case paper for description of smx, emx, and alpha.) bdwnhl constant to determine tlamba, the value to be sub- tracted from the diagonal elements of Hessian matrix for a downhill step. tlamba is chosen as min ( (low- est eigenvalue - bdwnhl) , 0.00d0). ( default bdwnhl = 0.25) ndiag for every ndiag steps, the matrix is diagonalized to calculate tlamba, when ntrun=4 dfpred a rough estimate of the expected reduction in energy for the initial step (only for ntrun = 3). (default = 0.01 kcal/mol) _____________________________________________________ parameters for running Langevin modes (set ntrun = 5) _____________________________________________________ eta viscosity in centipoise ioseen 0: Stokes Law used for hydrodynamic interaction 1: Oseen interaction included 2: Rotne-Prager correction included hrmax hydrodynamic radius for the atom with largest area exposed to solvent. If a file named 'expfile' is present, then the relative exposed areas are read from that file as a namelist namelist /exposure/ expr(natom) If 'expfile' does not exist, then all atoms are assigned a hydrodynamic radius of hrmax. _____________________________________________________ parameters for transition state search (when ntrun = 2) _____________________________________________________ istart 0: new calc. (default) 1: restart calc. iflag 0: search for transition state then minimum (default) 1: search for minimum from a transition state -1: search for a transition state, then stop ivect no. of eigenvectors wanted (default=2) ( ivect has to be >=isdir ) isdir eigenvector along which search for transition state is to be made. Note that translations and rotations are removed from the Hessian, so this number refers to the ordering of the "true" vibrational normal modes. (default=1) idir search direction: = 1 along isdir direction (default); = -1 opposite isdir direction isw no. of steps before switching to the lowest mode (default=40) hnot initial step length (default=0.1 Ang.) buphl switch to Newton-Raphson step when lowest eigenvalue is less than this for uphill walk. (default=-0.1) _____________________________________________________ Cards 3 group cards for the parts of the molecule that move, if ibelly.ne.0. See group documentation for format. Cards 4 group cards for the part of the molecule to be con- strained, along with the constraint weights, if icons.ne.0. See group documentation of format. 2. Nmanal Usage: nmanal [-O] -i nmdin -o nmdout -p prmtop -v vecs -r rvecs -O Overwrite output files if they exist. _______________________________________________________________ This is a general routine to do vibrational analysis by projecting cartesian normal mode eigenvectors (generated by the nmode program ) onto internal coordinates or onto "rigid groups." For each internal coordinate, the program will deter- mine the projection of each normal mode onto that coordinate, and the fraction of the total potential energy change along the normal mode that is contributed by that internal coordinate (the "potential energy distribution" for each mode.) You can also sum over all modes to obtain the rms thermal fluctuations for any particular internal coordinate. The program can also compute the rms thermal fluctuations of atoms in a cartesian coordinate frame, and will compute time correlation functions and fluctuations of internal coordinates. The original code was written by D. T. Nguyen and D. A. Case at U. C. Davis, 1985. The RMS analysis added by Barbara Rudolph. Capabilities for time correlation functions were added by J. Kottalam at Scripps Clinic. Responsibility for the final versions of the codes (and for any bugs) rests with Dave Case. Files used in the program: nmdin : control input for the run, described below. nmdout : standard output file for print and error messages prmtop : formatted parameter file vecs : formatted file containing output normal mode frequencies and eigenvectors. This file is generated by the program nmode. rvecs : formatted file containing the reference eigenvectors and associated frequencies _______________________________________________________________ Input found on nmdin: Card 1: Title of the run Cards 2: namelist /data/, which contains the following parame- ters: ntrun Values of ntrun from -1 to 3 are used to analyze modes in terms of internal coordinates (if ipro = 1), to compute thermal fluctuations in internal coordi- nates (if ifluc = 1), or to compute time correlation and cross-correlation functions (if ntrun = -1). Options 4 and 5 are present only for historical rea- sons (although the code might be a good starting point for some interesting calculations), and options 6 to 8 carry out some specialized tasks: (default=1) =-1 project eigenvectors onto those internal coordi- nates read in on subsequent cards (labelled "3c", below). = 0 project eigenvectors onto bonds only = 1 project eigenvectors onto all internal coordi- nates = 2 " " " angles and dihedral angles = 3 " " " dihedral angles only = 4 project eigenvectors onto "dynamics groups" = 5 project eigenvectors of the system('system vec- tors') onto reference eigenvectors. (This is a fairly specialized option, and you will probably have to look at the code to see what you really get. It has rarely been used.) = 6 just calculate rms fluctuations in cartesian coordinates. = 7 compute dipole-dipole correlation functions. In this case, prmtop is not read, and the &data namelist must be followed by cards that have two integers per card (free format), giving the atom numbers for each pair for which the correlation functions are desired. See subroutine "corf" for details of the calcula- tional procedure. = 8 project MD snapshots onto normal mode directions nvect = number of eigenvectors in file vecs to be read in (default=50) ivform = 0 if the input vectors are in unformatted form = 1 for input vectors in formatted form (default) ieff = 0 use true frequencies (default) = 1 use effective frequencies (not implemented!) pcut cutoff value for printing out projections: print will occur if the estimated contribution of an internal coordinate to the total potential energy distribution along this mode exceed pcut. (default = 0.02) ibelly = 0 (default) no belly = 1 belly calculation ibeg first eigenvector to be analyzed (default = 7) iend last eigenvector to be analyzed (default = 50) ifluc = 0 don't do the rms internal crds. fluctua- tion(default) = 1 do the rms internal crds. fluctuation for the internal coordinates selected by the "ntrun" variable ipro = 0 don't print out the projections onto internal crds. = 1 print out the projections onto int. crds (default) bose true. if quantum (Bose) statistics are to be used in populating the modes; .false. (default) if classical (Boltzmann) statistics are to be used. natom number of atoms; only needed if ntrun=7 or 8 ihsful = 0 if dipole-dipole correlations do not include contributions from distance fluctuations =1 (default) if both distance and angle fluctuations are to be included in dipole-dipole correlations. tmax maximum value for time correlation functions if ntrun = 7. Default is 0.0, i.e. no time correlations will be carried out. tintvl interval for time correlation functions (default = 1.0). _____________________________________________________ The following are only used if ntrun=8: first first snapshot from MD simulation be be projected onto normal mode directions, when ntrun = 8. Default = 1. last last snapshot to be projected. Default = 9999. iskip every iksip-th snapshot will be projected. Default = 1. _____________________________________________________ The following are only used if ntrun=5: nrgrp number of rigid groups (default = 0) nrvec number of reference eigenvectors to be read in from file "rvecs" (only if ntrun=5; default=0) nrat number of atoms in reference system (default=0) iat first atom number of the part of the system to be excised and compared to reference system (only if ntrun = 5; default=1) jat last atom number of the part of the system to be excised (only if ntrun = 5; default = natom) imov flag to rotate/translate eigenvectors of system and reference to principal axes. This essentially decou- ples translation and rotation of the excised part of the system as a rigid body from the rest of the vibrational motion (only if ntrun=5, default=0). _____________________________________________________ Cards 3a group cards for the parts of the molecule that move, only if ibelly.ne.0. See group documentation for format. This card is not needed when ibelly = 0. Cards 3b group cards for subdividing the molecule into rigid groups (if ntrun=4). See group documentation for format. Each rigid group will have its own set of cards 3b. Cards 3c (if ntrun .eq. -1) Quantities for which time corre- lation functions are to be calculated. These quanti- ties are of the form of internal coordinates. All input is free format, which means that you must enter all numbers -- blanks are ignored. TYPE, IAT, JAT, KAT, LAT INTNAME = identifier for this internal coordi- nate (character*8) IAT, JAT, KAT and LAT are atom numbers. Set LAT to zero for bonds and angle, KAT to zero for bonds. Repeat this card for as many internal coordinates as you are interested in, up the the value of MAXINT specified in the "sizes.h" header file. Input is terminated when the end-of-file is reached. 3. Lmanal Usage: lmanal [-O] -i lmdin -o lmdout -c inpcrd -l lmode -O Overwrite output files if they exist. _______________________________________________________________ This program will compute time correlation functions from Langevin modes. Note that since the time-independent aspects of the molecular normal mode description are independent of solvent viscosity, all of the equal-time correlations (such as rms fluctuations in cartesian or internal coordinates) will be the same as for the vacuum calculation. Hence the companion program nmanal should be used to compute those. Input description for the lmdin file: namelist default meaning &data ntrun 1 'type of run' flag 1: correlation function calculated is for the deviation of the length of the vector from the reference value in the minimum energy structure. i.e., / 2: for the orientation of the vector i.e., 3: the frequency distribution is plotted i.e., the imaginary parts of the eigenvalues are counted at every interval of 10 wavenumbers. Other input parameters are irrelevant. kup 1 lup 2 atom numbers defining a position vector nvect 12 number of langevin modes to be used tf 2.0 final time for correlation functions i.e., t ranges from 0.0 to tf picoseconds np 1000 number of points at which to calculate correlation function between t = 0.0 and t = tf ps. bose .false. .true. if quantum (Bose) statistics are to be used in populating the modes; =.false. (default) if classical (Boltzmann) statistics are to be used. &end The input file inpcrd is a standard Amber coordinate file; the file lmode is that created by the nmode program with ntrun = 5. Output files are CORF ( if ntrun = 1 or 2 ) and DENS, CUMU (if ntrun = 3.) All of the output files (except lmdout) are input files for the package, which will create plots of the correlation functions. You should(?) have little trouble converting them to some other plotting package in order to see the correlation functions. +-----------------------------------------------------------------+ | Sample input file for nmode with ntrun=1 | +-----------------------------------------------------------------+ | | | get vibrational modes for staph nuclease ternary complex | | &data | | ntrun = 1, do vibrational calculation | | cut=10.0, cutoff; use same value in minimization | | idiel=0, distance-dependent dielectric | | nvect=6753, write out all 3*N modes... | | ivform=0, ...in unformatted form... | | ilevel=0, ...with no removal of trans. & rotation | | drms = 0.0001, will complain if rms gradient is not | | less than this | | &end | +-----------------------------------------------------------------+ +------------------------------------------------------------------+ | Sample input file for nmanal with ntrun=7 | +------------------------------------------------------------------+ |# | | Get N-H S**2 values from quasi-harmonic modes | | &data | | ntrun=7, compute dipole-dipole correlation fns. | | nvect=4496, this many modes in input file | | ibeg=1, iend=4496, use all of the modes for the calculation | | ivform=0 unformatted modes files | | natom=1529, molecule has this many atoms | | ihsful=0, do no include distance fluctuations | | &end | |3 4 atom numbers for N and H of residue 1 | |11 12 | |20 21 | |28 29 | |38 39 | +------------------------------------------------------------------+