Vibration Analysis The vibrational analysis section of CHARMM has been designed to be a general purpose normal mode generation and analysis facility. Also included is an extensive set of vector analysis and comparison features. Support programs such as the iterative diagonalization program, or the restartable large matrix diagonalization program are compatable with this facility. Also included are routine to generate trajectories which can be used for examining modes on the picture system. In order to process commands with the vibrational analysis routines, The energy terms must all be defined, and the structure must be determined (see *note Needs: (energy.doc)Needs.). At present, SHAKE and images (see *note Images: (images.doc).) are not supported. Systems with atoms fixed by the CONStraint FIX command can be treated by using the REDU FIX option. Within the vibrational analysis command mode, all miscellaneous (MISCOM), coordinate manipulation (CORMAN), and internal coordinate (IC) commands are allowed. Keywords used to define Hydrogen bonds and nonbonded interactions may be included in the command that invokes VIBRAN. * Menu: * Syntax:: Syntax of the VIBRan command and all commands * Normal modes:: Description of normal modes * I/O:: Description of the read and write commands. * Diagonalization:: Description of the diagonalization command. * Quasiharmonics:: Description of the quasiharmonics command. * Reduce:: Reduced basis normal mode analysis * Dimb:: Iterative diagonalization (DIMB). * Explore:: Command to explore the energy hypersurface * Fluctuations:: Description of the fluctuation command. * Princ-Axis-Flucts:: Description of the PAFL command. * Projections:: Description of the projection command. * Mode-definition:: What vectors may be specified * Rayleigh:: Compute Rayleigh quotients for selected modes * PED:: Compute Potential Energy Distribution for a mode * Edit:: Description of the edit facility. * Basis:: Generate entire basis sets. * Fill:: Description of the fill command. * Second derivatives:: Generation, storage, and uses of second derivatives * Optimization:: Optimize parameter set to fit spectra

Syntax for vibrational analysis command and subcommands [SYNTAX VIBRational analysis] Main command: VIBRan [hbond-spec] [nbond-spec] [nmode-spec] [IHBFRQ 0] [INBFRQ 0] nmode-spec ::= NMODe integer (this specification allocates space on the heap and its ) (default is NATOM*3. For large systems, if normal modes are) (to be used, it should be set to the largest number needed.) (Its default value is set to 1 if NATOM is greater than 50) hbond-spec ! see *note hbonds:(hbonds.doc). nbond-spec ! see *note nbonds:(nbonds.doc). Subcommands: miscellaneous-command-spec ! see *note miscom:(miscom.doc). IC ic-subcommand-spec ! see *note intcor:(c22doc.intcor.doc). COOR coor-subcommand-spec ! see *note corman:(corman.doc). READ { NORMal-modes [FILE] unit-spec [APPEnd] [mode-spec] } WRITe { NORMal-modes [mode-spec] [CARD] } unit-spec { SECOnd-derivatives [RAISe] [MASS] [CARD [FINIt [STEP val] [TOL val] [atom-selection] ]] } { TRAJectory mode-spec magnitude-spec [PHAS real] [SHAKe] [SEQUential-files] [NCYC integer] [STEP real] [SUPErpose] [RANDom] [ISEEd integer] } PRINt NORMal-modes [mode-spec] [magnitude-spec] [INTDer [FINIt]] [VECTors] [DOTProducts] [DIPOles] [STATistics] [atom-selection] PROJect mode-definition [mode-spec] [magnitude-spec] DIAGonalize [NFREquencies integer] [NADD integer] [RAISe] [FINIte [STEP real]] QUASIharmonics [NUNIts integer] FIRStunit integer [NSKIp integer] [NFREquencies integer] [NADD integer] TEMP real REDUce IUNBas int [IUNTrans int] [BIG] [NFREq integer] [RAISe] [FIX] [CMPAct] RBQUasi IUNBas int [IUNTrans int] [NFREq integer] TEMP real [NUNIts integer] FIRStunit integer [NSKIp integer] DIMB IUNMode int [IUNRead int] [PARDim int] [NFREq int] [NADD int] [RAISe] [CUTF1 real] [ITERation int] [SAVFreq int] [TOLErance real] [STRENgth real] [DWINdow] [BIG] [SCILib] EXPLore [mode-spec] [magnitude-spec] unit-spec [GRID int] [COMP] [SHAKe] [ [ADJUst] [BOLTzmann temp] [GAUSs factor] ] FLUCtuations { ATOM [atom-selection] } [mode-spec] [magnitude-spec] { IC } [QUANtum] [VERBose] { USER [atom-selection] } PAFLuctuations { ATOM } [atom-selection] [mode-spec] [magnitude-spec] { GROUP [MASS] } [QUANtum] [VERBose] [COORdinate] [SAVE] { USER } [CONTinue] RAYLeigh [mode-spec] [SAVE] PED [mode-spec] [magnitude-spec] [TOL real] THERmodynamic [mode-spec] [TEMPerature real] [STEP real] EDIT { INCL mode-definition [ORTHog] [TO mode] } [atom-selection] { REMOve [mode-spec] mode-definition [NONOrm] } { DELEte [mode-spec] } { ORTHogonalize [PURGe] [mode-spec] [TOL real] } { SHAKe mode-spec } { ZERO mode-spec } { MOVE MODE n [TO m] [SCALe real] [NONOrm] } { ADD MODE n [TO m] [SCALe real] [NONOrm] } { MULT MODE n SCALe real } { SET MODE n SCALe real [NONOrm] } { COPY MODE n TO m } BASIS { IC { FIRSt BOND } } [NOORthonorm] { { FIRSt ANGLe } } { { DIHEdral } } { { SECOnd ANGLe } } { { SECOnd BOND } } { } { TR atom-selection [BALAnce] } FILL { DIFF } [mode-spec] [magnitude-spec] [APPE] { COMP } CORREL [ MAXTimesteps int ] [ MAXSeries int ] [ MAXAtoms ] [ COVAriance] {mode-spec DATUnit integer } OPTImize [PARUnit integer] {NOMOdes } - [NRSTep integer] [NDSTep integer] [ABNR] - [TOLSsq real] [TOLPar real] [MAXIt integer] [MINIt integer] - [DIFZ real] [NPAR integer] [PRLEvel integer] [PRUNit integer] - [WEIG] [MASS] [COORdinates] [SYMMetric-differences] - [USERobservations integer] [CHARges real] - [atom-selection] MASS ! turn on mass weighting flag (default) NOMAS ! turn off mass weighting flag END unit-spec ::= UNIT unit-number mode-spec ::= MODE integer [ THRU integer ] mode-definition ::= { { TRAN } { X } } [NONOrm] [NOTR] { { ROTA } { Y } } { { Z } } { } { SPHEre { X } [IX int] } { { Y } [IY int] } { { Z } [IZ int] } { { R } [IR int] } { { TX } } { { TY } } { { TZ } } { } { COMP } { DIFF } { FORC } { USER integer } { BOND atom atom } { ANGL atom atom atom } { DIHE atom atom atom atom } atom::= {residue-number atom-name} { segid resid atom-name } { BYNUm atom-number } magnitude-spec ::= { TEMP real TFRE real } [NONOrm] { KCAL real TFRE real } { RMS real } { MRMS real } { FACT real } miscellaneous-command-spec *note misc:(miscom.doc). ic-subcommand-spec *note ic:(intcor.doc). coor-subcommand-spec *note coor:(corman.doc).

Normal Modes There are two ways that normal modes are stored internally in CHARMM. The most common usage is as one double precision mass weighted array. A series of such arrays usually span an orthonormal basis (as would be the case upon diagonalization). The second method is to represent a normal mode as three non mass weighted coordinate displacement arrays, stored in double precision. The program automatically converts between them as necessary. Whenever interconversion is to be done, a "magnitude specification" may be given. This specification requires a step type and step length. The valid step types are; FACT (simple factor) TEMP (put mode at this temperature) KCAL (put in the specified Kcals) MRMS (step along until this mass weighted RMS is obtained) RMS (step along until this RMS is reached). When specifying TEMP or KCAL, a terminal frequency (cm-1) may be specified (default TFRE is 5.0) which prevents excessive stepping along very low frequency or translation-rotation modes. The procedure used in going from double precion to coordinate displacement arrays is; 1) Normalize double precision vector (unless NONOrm keyword is used) 2) Mass weight by dividing by root(mass) 3) Multiply by appropriate scale factor from step type and length To convert from single precision into double precision, the procedure is; 1) Save inner product (as initial step length) 2) Mass weight by multiplying by root(mass) 3) Normalize vector (unless NONO keyword is used) 4) Compute scale factor using initial step length and step type Whenever a magnitude-specification is called for, some interconversion will take place. The conversion from double precision to coordinate displacements takes place in the subcommands; PRINt NORMal-modes WRITe TRAJectory PROJect FILL EXPLore FLUCtuations The interconversion from coordinated displacements to double precision takes place in the subcommands; PROJect EDIT NORMal-modes The Normal Mode data structure of double precision arrays is local to the Vibrational analysis section of CHARMM, and the storage space for these arrays is released when exiting to CHARMM via the END command.

I/O For Normal modes The VIBRAN section supports its own I/O commands. Commands to read, write and print coordinates (see *note corman:(corman.doc).) and internal coordinates (see *note intcor:(intcor.doc).) are identical with those in the main program. This section can read and write the Normal Mode data structure and write out the second derivative matrix in several ways (for external use) and normal mode trajectories (for the movie programs). Also, useful information about normal modes may be printed using the PRINT NORM command. READ { NORMal-modes [CARD] unit-spec [APPEnd] [mode-spec] } WRITe { NORMal-modes [CARD] [mode-spec] } unit-spec { SECOnd-derivatives [RAISe] [MASS] - [CARD [FINIt [STEP val] [TOL val] [atom-selection] ]] } { TRAJectory mode-spec magnitude-spec [PHAS real] [SHAKe] - [SEQUential-files] [NCYC integer] [STEP real] [SUPErpose] [RANDom] [ISEEd integer] } PRINt NORMal-modes [mode-spec] [magnitude-spec] [INTDer [FINIt]] - [VECTors] [DOTProducts] [DIPOles] [STATistics] - [atom-selection] By default, normal mode vectors are read and written in binary (FILE) format, but ascii (CARD) can be specified. When writing, a unit must be specified, and a contiguous subset of modes may be specified using the mode-spec. When reading modes, all selected modes (from mode-spec) of the modes in the normal mode file will be read (assuming there is enough space). Note, when modes are selected on input, the selection is relative to the ordering of modes in the input file and does NOT correspond to the destination of these modes. If the available space is exhausted, a warning is issued, and further reading stops. Existing modes will be deleted when the READ NORM command is executed unless the append option is used, in whichcase, the new modes are added sequentially at the end. No modification of modes is done upon reading (i.e. normalization, or orthogonalization). When printing Normal Modes, a variety of options may be specified. A contiguous subset of modes may be specified (the default is all modes), and an appropriate magnitude may be specified (see *note modes:Normal Modes.). For each specified mode, the frequency (cm**-1), eigenvalue (Kcal/gram/A**2), force projection (Kcal/mole/A), % translation-rotation, and magnitude information will be printed. In addition, internal derivatives (or optionally by finite differences with the FINIte keyword) of the internal coordinate data structure (INTD keyword), displacements in coordinate space (VECTors keyword), dipole derivatives (DIPOles keyword), dotproducts with other modes can be printed (DOTProducts keyword), and some vector statistics (STATistics keywords). The Second Derivative matrix may be written out in binary or card format for input to other programs. In particular, there is a program that extracts normal modes from a large secular equation (DIAGIT). There is also a restartable diagonalization program for medium sized systems. The RAISe keyword will shift the translation-rotation modes to high frequency (currently 5000 cm-1). This option only work with the card output format. The MASS keyword will calculate and write out M**(-1/2) * H * M**(-1/2) which is the matrix used in the eigensolver with the DIAG command. Finally, this command may be used to test the second derivative energy surface determination by a comparison of finite derivatives and analytic derivatives. This is done with the FINIte and CARD keywords. Two values are also looked for, the finite difference step size (STEP, def=0.005), and the difference tolerence for printing (TOL, def=0.0001). Trajectory files may be written out for a set of modes with a given magnitude factor. The modes for all specified modes may be written out in one file, or in separate files for different modes where sequential unit numbers are used starting with the specified unit number. The SEQU keyword will cause sequential files to be written. By default, one cycle of 12 frames will be written. The keywords NCYC and PHASe can be specified to alter the number of cycles and the phase angle between frames within a cycle. The PHASe keyword should be an integer divisor of 360.0 ( a PHASe value of 7.2 will result in 50 frames per cycle). The total number of frames is given by NCYC * 360/PHASe. If it is desired to specify a particular time step between frames, the STEP keyword will specify this value in picoseconds. The total number of frames will remain unchanged will remain unchanged. If a movie of 10 modes and displaying one picosecond of each mode is desired with 500 frames per picosecond (about 25 seconds of film time), then the appropriate input could be: WRIT TRAJ MODE 7 THRU 16 STEP 0.002 PHAS 3.6 NCYC 5 TEMP 2000.0 It is also possible to write out a trajectory that consists of a superposition of a number of normal modes. The STEP parameter needs to be specified, together with PHAS and NCYC, to set the length of the trajectory in ps. By default, all normal modes in the trajectory have an initial phase of 0. Random initial phases can be specified by using the keyword RANDom and specifying a seed for the random number generator (default ISEED = 314159). A movie of 50 ps (10 frames per ps) for modes 7 to 100 with random initial phases could be generated by: WRIT TRAJ MODE 7 THRU 100 SUPE STEP 0.1 PHAS 3.6 NCYC 5 RAND TEMP 300.0

Diagonalization of the second derivative matrix The DIAGonalize command will generate the second derivative matrix, mass weight this matrix, and then diagonalize to generate normal modes. There are some needs and restrictions for this command. They can be summarized; 1) SHAKE may not be used 2) ST2 waters may not be used 3) Periodic boundaries and images may not be used without FINIte (see crystl.doc) 4) Group electrostatics and extended electrostatics may not be used without FINIte 5) All coordinates and energy lists (nonbond, hbond) must be set up 6) The number of atoms does not exceed the limit (300) (unless BOMLEV is reduced) Once normal modes have been generated, they may be saved (WRITE NORM...), or analyzed directly. DIAGonalize [NFREquencies integer] [NADD integer] [RAISe] [FINIte [STEP real]] The RAISe keyword, will cause the normal modes corresponding to translation-rotation to have a very high frequency (currently 5000 cm-1). This option is intended for calculations where rotational coupling terms are to be removed. The NADD option will cause the the specified number of lowest modes to be skipped in the evaluation. For example, if "NADD 100" is specified, then the first mode of the result will correspond to mode 101 in the actual matrix. This option has been added, so that modes of moderatly sized systems (200-400 atoms) can be found in groups when the memory requrement for a full calculation are prohibitive. The FINIte keywords causes the Hessian to be generated from the finite differences of the forces. This option requires 6*N+2 energy determinations, so it should be reserved for smaller systems. The step length for finite difference may be specified with the STEP keyword (default 0.005 Angstroms). A STEP value that is too large will cause errors due to anharmonicity. A STEP value that is too small will result in inaccurate hessian elements (force differences are divided by the step length).

Quasiharmonic dynamics from molecular dynamics QUASIharmonics [NUNIts integer] FIRStunit integer [NSKIp integer] [NFREquencies integer] [NADD integer] TEMP real For quasiharmonic dynamics, a dynamics trajectory file(s) is read. During this read, atom position fluctuation tensors are generated. The reference coordinates must be present in the main coordinate arrays. The COOR DYNA command can preceed the QUASi command if the average coordinates are to be used as a reference. Another good choice is to use an energy minimized coordinates set as the reference. The resulting fluctuation matrix is mass weighted and diagonalized. The resulting modes are the quasiharmonic modes of the system, and should roughly match those from a strightforward normal mode calculation. The significant differences of these methods are that anharmonic terms are included here. This method can yield misleading results when there are long lived transitions, or when there is very slow interchange of energy between degrees of freedom. The estimated frequency for a given mode will strongly depend on the average energy in this mode. For the case of slow energy transfer, this may be significantly different from kT. If internal motion is to be studied, it is strongly recommended that translation and rotation be removed from the dynamics trajectory. This is done with the MERGe command using all atoms and mass weighting. Without mass weighting, the translation-rotation motion is not projected out. In order for this method to estimate frequencies from fluctuations, the average temperature (TEMP real) must be specified. The dynamics trajectory is read in the usual manner. Specifications are; NUNIts integer - number of I/O units FIRStu integer - first I/O unit NSKIP integer - integration step modulo to use As with the DIAG command, this command will also accept the keywords; NFREquencies integer, and NADD integer. Reduced basis quasiharmonic calculations are also possible (RBQUas). See the section for reduced basis normal mode analysis for information on setting up a basis set. The method is exactly the same as above except the fluctuation matrix is generated in the reduced basis instead of the mass weighted cartesian basis. RBQUasi IUNBas int [IUNTrans int] [NFREq integer] TEMP real [NUNIts integer] FIRStunit integer [NSKIp integer]

Reduced basis normal mode analysis With the REDUce command, it is possible to do reduced basis normal mode analysis. This can be affected to constrain certain degrees of freedom (leave them out of the basis), or to reduce the size of a large problem. I can also be used to remove translation and rotation from a calculation. The same restrictions for a full normal mode calculation apply here as well. By defining the keyword CMPAct the initial second derivative matrix in Cartesian coordinates will be stored in a compressed form. This means that only non-zero elements are stored. This option is very useful if you deal with very large molecules, where the initial second derivative matrix may be too large to fit into memory. For instance, in a system with 2000 atoms, the regular Cartesian Hessian takes up about 144 Mbytes, whereas the compressed matrix takes up an absolute maximum of 43 Mbytes (usually less, depending on nonbonded cutoff distances). Currently, this option allows atoms to interact with up to 300 other atoms, on average. This means that if large cutoff distances (>15 A) are used, the command may exit with an error message. This command requires an external basis set. The form of the basis set is a CHARMM normal mode file of the correct dimension. The basis will include all of the vectors in that file. The basis is expected to span an orthonormal space. These basis sets may be made from the BASIS command followed by the EDIT ORTHog command if appropriate. The reduced basis may also be a subset of normal modes. For large systems, this is the refinement step. By repeating the calculation in a reduced basis, most roundoff effects from the tridiagonalization are removed. This command removes all existing vectors and replaces them with the result of the diagonalization in the original basis (i.e. the new vectors are linear combinations of the basis vectors). The number of vectors to compute may be specified with the NFREquencies keyword. The default is to compute all frequencies (for the entire reduced basis). REDUce IUNBas int [IUNTrans int] [BIG] [NFREquencies integer] [RAISe] [FIX] [CMPAct] The IUNBas keyword points to the basis vector file. If the IUNTrans value is specified, the transformation matrix (NDIMxNFREQ) relating the basis vectors and the final normal modes will be saved to that unit. The number of frequencies to be computed may be specified. The default is to compute all of them. The raise option raises the translation/rotation degrees of freedom to a high frequency. This keyword will have no affect if the net translation/rotation degrees of freedom are not in the basis. The BIG keyword signifies that the entire set of basis vectors is not to be stored in memory. This is essential for large calculations (i.e. when HEAP allocation errors occur without this option). In this case, VIBRAN ... NMODE 1 ... should be specified to save space, since no space is needed to store the basis vectors or the results. The BIG option will cause the backtransformation step to be suppressed, thus the IUNTrans unit should be specified (or all you will get is the frequencies). The will be a special purpose support program to backtransform this eigenvactors (normal modes) to the original basis. If the BIG option is not specified, there must be enough space allocated to store the entire basis in memory (VIBRAN NMODe int). The keyword FIX should be used when part of the system is fixed by the CONStraint FIX command. In the modes that are calculated, the fixed atoms will not move, but their presence is taken into account for the motions of the unfixed atoms. This option does not require a basis vector file, i.e., the IUNBas parameter should not be specified. For large systems, the use of a compressed second derivative matrix by specifying the keyword CMPAct is recommended. The BIG option is not allowed.

Iterative Mixed-Basis Diagonalization With the DIMB command, one can do an iterative diagonalization, which will consume a lot less memory space than a regular (DIAG) diagonalization, but will still result in exactly the same normal modes. The method (see L.Mouawad and D.Perahia (1993), Biopolymers, 33, 599) does repetitive reduced-basis diagonalizations. These reduced bases are constructed partially from the not yet converged eigenvectors and from regular Cartesian coordinates. A unit IUNMode for writing out the eigenvectors has to be specified. For restarts from existing eigenvector files, the unit IUNRead should be specified. The NFREq, NADD, and RAISe keywords have the same effect as in the regular DIAG command. The maximum size (in atoms) for a block to be diagonalized is defined with the PARDim (default=200) keyword, and can be tailored to the size of the available memory. This affects the number of blocks the system will be divided into. Make sure that the requested number of modes (lesser of NMODEs and NFREq) is smaller than PARDim*3. If no IUNRead unit is specified, initial basis vectors are calculated by diagonalizing main diagonal blocks of the full matrix. This initial calculation involves the use of residual (Lanczos) vectors, and will result in only half the number of requested modes. The full number of requested modes can be obtained by using the following sequence of commands (in this example NMODes is set to 200 and PARDim to 100, but these values can be changed): Create the initial basis and write out on unit 20. This results in (NMODes+6)/2 = 103 basis vectors. The usage of BIG reduces memory requirements by writing intermediate vectors to disk temporarily. Specification of ITERations 0 exits the routine after the calculation of the basis vectors: OPEN WRITe FILE UNIT 20 NAME initial.bas VIBRan NMODes 200 DIMB ITERations 0 PARDim 100 IUNMode 20 BIG END Do the iterative diagonalization, using the precalculated initial basis of 103 vectors. 103 converged modes will be written to unit 10: OPEN READ FILE UNIT 20 NAME initial.bas OPEN WRITe FILE UNIT 10 NAME modes.mod VIBRan NMODes 103 DIMB ITERations 100 TOLErance 0.05 PARDim 100 - IUNMode 10 IUNRead 20 DWIN END The usage of the keyword BIG is only allowed when ITERations 0 is specified. It will use unit IUNMode to temporarily store intermediate vectors. Make sure that there is enough available disk space. The size of the intermediate file will be about (NMODes*NAT3*8) bytes, where NAT3 = 3 * number of atoms. A rough estimate of the required memory space is: { 4*(PARDim3*PARDim3) + (2*NMODes+10)*NAT3 } * 8 bytes (without BIG) { 4*(PARDim3*PARDim3) + (1*NMODes+10)*NAT3 } * 8 bytes (with BIG) Usually convergence will be faster with larger PARDim and with smaller number of modes. Too small a number of modes will slow down converge due to a limited basis set of eigenvectors. With CUTF1 (default=50.0) a cutoff frequency (cm-1) can be defined. Only normal modes with frequencies below CUTF1 will be calculated. The maximum number of iterations can be specified with the ITERations (default=10) parameter, and the desired accuracy of the eigenvectors with TOLErance (default=0.05). TOLErance values range between 0 and 1, where 0 is the highest accuracy. The standard DIMB algorithm uses a single "window" of Cartesian coordinates that is added to the basis set. When the DWINdow keyword is specified, two different "windows" are added. The DWINdow method is more efficient, because it uses windows with strongly interacting atoms more often than those with weakly interacting atoms. When the DWINdow option is chosen, the parameter STREngth (default=0.0) defines which atom sets are considered to be interacting at all. If the sum of absolute values of the second derivatives between two particular atom sets is lower than STREngth, those second derivatives will not be considered. SAVFreq (default=number of blocks) defines the frequency of saving the eigenvectors to disk for the DWINdow option. When DIMB calculations are carried out on a Cray vector computer, it is advantageous to use the general EISPACK routines do diagonalize the submatrices. This can be defined by adding the keyword SCILib.

Exploring Energy Surfaces The explore command is used to search along a particular mode and to compute the energy at regular intervals. For now, this search is limited to one dimension, but with the EXPL COMP (using the comparison coordinates) option, a two dimensional search is possible. This is done by filling the comparison coordinate set with a structure perturbed along one mode, and then exploring along another. EXPLore [mode-spec] [magnitude-spec] unit-spec - [GRID int] [COMP] [SHAKe] - [ [ADJUst] [BOLTzmann temp] [GAUSs factor] ] For this command, a magnitude specification and grid selection are read. The grid selection determines the number of energy evaluations. An odd number will include the center point, and an even number will not. The default value is 3, one energy evaluation at each extreem (determined by the magnitude specification), and one at the center. When the energies are computed, they are stored in a temporary array which may be output (for plotting) by specifying a unit number. They may also be used to adjust the eigenvalue of the particular mode base on a least squares quardratic fit. This is done with the ADJUst keyword. The weighting for this fitting is currently unity for every point. The adjust feature is intended for relatively small displacements. Large "rigid" dislacements will necessarily come up against bad close contacts. Quadratic fittings for this type of interactions, will be misleading. For this reason, two weighting options are available, an energy Boltzman weighting (BOLTzmann keyword-value), and a gaussian displacment weighting (GAUSs keyword-value) which weights nearby points with a larger factor. When these keyword are used together, the overall weighting is a product of individual weightings. To avoid problems with excessive (quartic) bond stretching which occurs with large displacements, SHAKE may be invoked with the SHAKe keyword. In order to use this option, the SHAKe command in CHARMM must fist be invoked specifying which bonds (and angles) are to be maintained. A step value may be specified. By default, all shake comparisons are done with respect to the center coordinates. This sometimes leads to "DEVIATION IN SHAKE TOO LARGE" errors. The step option allows shake to be invoked in small steps. The step value gives the relative step for each intemediate SHAKE calculation. A value of 1.0 is the default, and a value of 0.2 will cause the extreem points to be computed in 5 steps.

Computing fluctuations from normal modes The fluctuation command computes atom, internal coordinate, or user specified fluctuations from normal modes. For atom fluctuations, the magnitude of the components for the overall fluctuation is stored in the comparison coordinate arrays. In addition, information giving the contribution of each mode is printed. For the IC (internal coordinate) option, the overall fluctuations are stored in the IC table. the "IC WRITE..." command may be used to save this data. FLUCtuations { ATOM [atom-selection] } [mode-spec] [magnitude-spec] { IC } [QUANtum] [VERBose] { USER [atom-selection] } For the user specified option, the user may supply the subroutine USEFLU which is call once for initialization (IMODE=0), once for each mode (IMODE=n), and once for termination (IMODE=-1). In preparing this routine, see the existing routine for interfacing requirements. If no USEFLU routine is provided, the equal time cross-correlation matrix for all selected atoms will be computed. For this command, a selection of modes may be made (default all), a selection of atoms may be made (default all), and a magnitude specification may be made. For this application, a temperature factor is usually used. The QUANtum keyword may be used if a quantum scaling factor is desired. For this option, higher frequency modes would have a diminished amplitude for a given temperature. The VERBose option causes the contribution for each selected atom with each selected mode to be printed. This is useful when a power spectrum is to be computed or prepared for plotting.

Calculation of anisotropic (and isotropic) fluctuations from Normal Modes The PAFL command sums over the current set of normal modes (default) or some subset of this (given by the [mode-spec]) to create a fluctuation matrix. The three-by-three fluctuation matrix for each atom may be diagonalized (default) to produce the principal axis flucutions given in terms of three mutually perpendicular unit principle axes. In general, there will be a different set of principle axes for each atomic center. Alternatively, the COORdinate option may be used, in which case the fluctuation matricies are not diagonalized. Instead, the fluctuations are given in term of the coordinate axes in which the normal modes were calculated. Since the size of the fluctuations are temperature-dependent, a [magnitude-spec] in terms of temperature (such as TEMP 300.0) should be given. PAFLuctuations { ATOM } [atom-selection] [mode-spec] [magnitude-spec] { GROUP [MASS] } [QUANtum] [VERBose] [COORdinate] [SAVE] { USER } [CONTinue] The summation over modes which creates the fluctuation matrix may be controlled in a number of ways. The easiest of these is to specify a TFREquency. Normal modes of frequency lower than TFREquency are simply ignored in the summation. By setting TFREquency above the rotation-translation modes, but below the lowest vibrational mode, one achieves what is the desired result (the fluctuations due to internal modes of vibration) for most applications. For applications requiring more specificty (such as calculating the fluctuations due to modes 8,9,10,13,14, and 15 only) one may use a more detailed [mode-spec] to select the desired modes over which the summation is to be performed). The second way of controling the summation is necessary when the normal modes for the desired system can not all be read into CHARMM at one time due to space allocation limitations. In this case one uses the SAVE and CONTinue options. SAVE tells the routine that you don't want to manipulate the fluctuation matrix once it is created. CONTinue indicates that you don't want to reset (zero) the fluctuation matrix before the present summation. Therefore, if the normal modes are stored on, let us say, four files, one reads in the first normal mode file and gives the PAFL SAVE command. One then reads in the second normal mode file (without the APPEnd option) and gives the PAFL SAVE CONTinue command. One then reads in the third normal mode file (without APPEnd) and gives the PAFL SAVE CONTinue command. Finally, one reads in the fourth normal mode file (without APPEnd) and gives the PAFL CONTinue command. Each of these four PAFL commands should have the same options (such as GROUP, ATOM, USER; same TEMP; same atom selection, MASS or not). Moreover, there should be no overlap between modes in each of the four sets. If the four (or any n) files have overlapping modes, then a selection must be made so that no mode is counted twice in the summation (see the test case). For the user specified option, the user may supply the subroutine USPAFL which is call once for initialization (IMODE=0), once for each mode (IMODE=n), and once for termination (IMODE=-1). In preparing this routine, see the existing routine (and also USEFLU) for interfacing requirements. The PAFL command is quite similar to the FLUCtuations command. VERBose prints out lots of stuff and QUANtum uses quantum scaling in calculating the contribution of each mode to the fluctuations. When ATOM is specified, the fluctuations are calculated individually for each atom currently selected. The fluctuations are given in angstroms, and so mass-weighting would make no sense. When GROUp is specified, the fluctuations for the movement of all specified atoms moving as a group are calculated. When MASS is also specified, the fluctuation of the center of gravity of that group of atoms is calculated. Otherwise, it is the fluctuation of the center of geometry which is calculated. At present USPAFL does nothing, and one would need to read the code to be able to interface this routine properly with the surrounding code.

Projection of normal modes onto vectors of interest The PROJect command will project the selected modes onto one of the definable vectors. PROJect mode-definition [mode-spec] [magnitude-spec] The information displayed is; MODE integer - mode number FREQUENCY - frequency of this mode (cm-1) NORMAL DOTPR - actual dotproduct of normalized vectors PROJECTION - ratio of guess vector projection to step length APPROX DEL E - estimate of energy increase along this mode TYPE - type of step (FACT, RMS, KCAL, TEMP) STEP - step length for this step type The NONOrm keyword may be used to prevent normalization of modes before projecting.

Several commands use a mode_definition, which allows the specification of a vector. There is a wide variety of possible vectors which may be specified, and there is also a user supplied routine, USERNM, which may be used to get any other desired motion (vector). These vectors may be used for the analysis of normal modes, or as additions to the basis of vectors (EDIT INCLude) for further analysis. mode-definition ::= { { TRAN } { X } } [NONOrm] [NOTR] { { ROTA } { Y } } { { Z } } { } { SPHEre { X } [IX int] } { { Y } [IY int] } { { Z } [IZ int] } { { R } [IR int] } { { TX } } { { TY } } { { TZ } } { } { COMP } { DIFF } { FORC } { USER integer } { BOND atom atom } { ANGL atom atom atom } { DIHE atom atom atom atom } atom::= residue-number atom-name The NONOrm keyword supresses the automatic normalization of the specified vector. This is desired for some applications, but is not normally needed. The NOTR keyword removes any new translation/rotation from the specified mode. This may be needed when only internal motions are to be analysed. The vectors which may be defined are; TRAN X - Translation along the X-axis TRAN Y - Translation along the Y-axis TRAN Z - Translation along the Z-axis ROTA X - Rotation along the X-axis ROTA Y - Rotation along the Y-axis ROTA Z - Rotation along the Z-axis COMP - Use the comparison coordinates (as a vector) DIFF - Use the difference between the MAIN and COMP coordinates FORC - Use the forces from the last energy evaluation USER integer - Use a user specified vector (defined in USERNM) BOND atom atom - Use vector which stretches the specified bond ANGL 3X(atom) - Use vector which bends the specified angle DIHE 4X(atom) - Use vector which twists the specified dihedral SPHEre ... - Appropriate homogeneous motion (spherical harmonics) For the TRANslate, ROTAte, COMP, DIFF, and FORCe options, an atom selection may be specified. Any nonspecified atoms will have a zero values in the vector. For the BOND, ANGLe, and DIHEdral options, the atoms specified don't have to be bonded together or have any special connectivity, but the first atom specified must not close back through a loop to the last atom specified. If this is a problem, then one of the bonds in the loop must be deleted for this calculation to work properly. For angle terms, the first two atoms specified (and anything they are connected to) move as one unit, and the third atom and anything its connected to will move as a second unit. For dihedrals, the first 2 atoms specified define one block, and the last 2 atoms defines the second block. The axis of rotation is about the middle 2 atoms. For example if one specifies PROJECT DIHE 1 N 1 CA 1 CB 1 CG1 for an isoleucine residue, the atom 1 CG2 will also rotate with CG1 because it is connected to the third atom (CB) which is part of the second unit.

Rayleigh Quotients The RAYLeigh command will compute the second derivative matrix, and project all selected modes onto the second derivative matrix after mass weighting. If the modes are eigenvectors, then the resulting quotients will correspond to the eigenvalues. The quotient values are given by; -1/2 -1/2 Q = < v | M H M | v > RAYLeigh [mode-spec] [SAVE] If the SAVE keyword is given, then the quotients and corresponding frequencies are saved and become part of the data for the selected modes. For example, a subsequent PRINT NORM command would then use these new values in place of the original frequencies.

Potential Energy Distribution This command is designed to work with one mode at a time. For a given magnitude specification, it computes the expectation value for the energy contribution change for each internal coordinate term (bond, angle, dihedral, and improper dihedral) and prints that term if the fluctuation is greated than the tolerence (default TOL 0.0001). PED [mode-spec] [magnitude-spec] [TOL real] For this method, it assumes that the structure is at a stationary point and that the potential is quadratic. In addition to printing out the individual energy terms, it also prints out the total for each class (bonds, angles...). These values can be used to determine the type of any given mode.

Editing the set of normal modes There are several commands that can modify the normal mode vector space. In addition to the obvious ones such as READ NORMal-modes and DIAGonalize, there is also an EDIT command which can be used to modify the normal mode data structure. EDIT { INCL mode-definition [ORTHog] [TO mode] } [atom-selection] { REMOve [mode-spec] mode-definition [NONOrm] } { DELEte [mode-spec] } { ORTHogonalize [PURGe] [mode-spec] [TOL real] } { SHAKe mode-spec } { ZERO mode-spec } { MOVE MODE n [TO m] [SCALe real] [NONOrm] } { ADD MODE n [TO m] [SCALe real] [NONOrm] } { MULT MODE n SCALe real } { SET MODE n SCALe real [NONOrm] } The EDIT DELETE command will delete specified modes from the date structure. The EDIT ORTHogonalize command is used to orthogononalize and optionally normalize a particular subset of normal modes. Additional modes can be added with the EDIT INCLude command. This command will add on the defined mode to the end of the normal mode vector, unless the destination is specified with TO n in which case it will be placed as mode number n. The translation and rotation options are important when setting up guess vectors for the iterative diagonalization program, or when individual coreolis coupling terms are needed. Single additional modes may be added with the EDIT INCLude DIFF or the EDIT INCLude FORCe or the EDIT INCLude USER or the EDIT INCL COMP commands. For the DIFF option, the difference between the comparison and main coordinate sets will be appended, for the FORCe option, the current values in the force arrays (from the last energy evaluation) will be appended. The COMP option will use the comparison coordinates as the appended mode. This is option is intended for use in the case that the comparison coordinate set is filled with displacements. There is also the ability to specify a user vector. This is done by the inclusion of the subroutine USERNM in your USERSB and using the USERLINK facility. For all of the EDIT INCLude options, the defaults are to root mass weight, normalize, and orthogonalize to the rest of the vector space. To skip any of these steps, the NONO, and NOOR keywords must be specified. See (*note define:Projections.) for a complete list of definable modes. The following commands (which all can take an optional atom-selection) provide a simple means of manipulating, and reshuffling modes, primarily intended for the setting up of guess-vectors for the OPTImize DOTProd command. EDIT ZERO sets the specified modes to zero. EDIT MOVE MODE n [TO m] [SCALe real] mode(m)=scale*mode(n) EDIT ADD MODE n [TO m] [SCALe real] mode(m)=mode(m)+scale*mode(n) EDIT SET MODE n SCALe real mode(n)=scale EDIT MULT MODE n SCALe real mode(n)=scale*mode(n) If the destination is not specified then the result will be appended as a new mode after the existing ones. The default value for the SCALe factor is 1.0.

Generate entire basis sets The BASIs command will append requested basis vectors to an existing (or null) set of orthonormal vectors. In this way, basis sets for reduced diagonalizations, constrained normal modes, or other calculations may be generated in a simple manner. BASIS { IC { FIRSt BOND } } [NOORthonorm] { { FIRSt ANGLe } } { { DIHEdral } } { { SECOnd ANGLe } } { { SECOnd BOND } } { } { TR atom-selection [BALAnce] } The BASIS command is similar to the EDIT INCLude ORTHog command, except that many new vectors may be added. There are two modes for this command. The IC mode will include a set of vectors based on the IC table. The selection of which section of the table to use is required. The choices are; FIRSt BOND FIRSt ANGLE DIHEdral SECOnd ANGLe SECOnd BOND All valid IC table entries (i.e. all atoms of specified set defined) will result in the addition of one vector to the normal mode basis, provided that this vector is not a linear combination of existing vectors. See the description for IC mode specification (*note define:Projections.) for a description of these vectors. The second mode is the TR (translation/rotation) mode which will usually add the 6 translation rotation degrees of freedom for the selected set of atoms. If the selected set of atoms is linear, then only 5 vectors will be added. If only one atom is selected, then just the 3 translation vectors are added (for this case all numvers but one will contain a zero). The BALAnce keyword is suggested, and its operation it to remove the net (whole system) rotation/translation components from the added vectors. If all atoms are selected, the BALAnce keyword will result in NO vectors being added. The purpose of this command is to facilitate setting a vector basis for a reduced normal mode calculation. Another purpose is in setting up translation/rotation modes as initial vectors for the iterative diagonalization procedure (no BALAnce option). This command does normalize the new basis vectors, but it will not modify any existing vectors. Each added vector is then orthogonalized from all exisiting vectors (unless the NOORthonormalize keyword is specified). If the vector has a zero norm following orthogonalization, it is rejected. This avoids any possibility of basis interdependancies (zero determinant). The norm after orthogonalization is saved in the eigenvalue array. Each added vector is then normalized to form an orthonormal basis. NOTE: this method will not work unless all exisiting vectors form an orthonormal subbasis (i.e. before the BASIS command is invoked, the PRINT NORM DOTProducts should give the unit matrix). The EDIT ORTHog command may be used to orthonormalize a set of vectors. The user is expected to process these modes with the REDUce command after saving them, or in an external program.

The Fill command The comparison coordinates can be modified with the FILL DIFF command. This command will copy the main coordinate set to the comparison set, and then step along the specified mode by the specified magnitude. When the append (APPE keyword) is used, the main coordinates are not first copied. FILL { DIFF } [mode-spec] [magnitude-spec] [APPE] { COMP } The FILL COMP command will fill the comparison coordinate displacement arrays with the specified vector. For this option, the comparison coordinate are zeroed before the displacement is added. The coordinates will then contain a displacement vector. The append option will prevent the zeroing of these arrays before stepping along the mode.

Second Derivatives The second derivatives are computed during energy determination and they are stored in temporary arrays that are allocated dynamically. Once obtained, they can be written out, or diagonalized internally. If 'QSECD=.TRUE.' the second derivatives of the energy are returned in the array DD1. This array contains the full upper half of the second derivative matrix, and contains (NAT3*(NAT3+1)/2) REAL*8 elements (NAT3 = 3 * number of atoms). Memory storage for the DD1 array may cause memory overflows, especially for large systems. If one wants to get the normal modes for large structures, it is often advantageous to use the DIMB or REDU CMPAct option. These options will fill a DD1CMP second derivative array with only the non-zero elements. The size of this array depends heavily on the non-bonded cutoff distance, since increasing that distance increases the number of non-zero elements. At this moment, the DD1CMP array will allow the average number of interactions per atom to be as high as 300. This number will allow non-bonded cutoffs up to 13A without problems, for all- hydrogen systems. A comparison between allocated memory space for DD1 versus the maximum memory space for DD1CMP is given below: Structure #atoms(explicit H) DD1 (MB) DD1CMP (MB) Deca-alanine 66 0.2 1.4 BPTI 568 11.6 12.3 Lysozyme 1264 57.6 27.4 Hemoglobin 5600 1,129.2 121.2

Parameter Optimization Optimization of a parameter set to obtain agreement between calculated (in VIBRAN) and experimental normal mode spectra is done using a non-linear least squares fitting procedure based on the Levenberg- Marquardt algorithm. (For references see IMSL routine ZXSSQ. Our code does not come from IMSL, but the algorithm is the same.) The optimization may also include a fit to a reference structure which should then be present in the comparison coordinate set.An atom selection can be used to specify which atoms should be included in the fit. Keyword COOR specifies that structural information is to be used. Keyword NOMOdes specifies that no normal modes are to be used in the fitting procedure. Ie OPTI NOMO COOR will use only a reference structure. It is assumed that the structure has been minimized BEFORE this module is invoked. This module REQUIRES a mode-specification (so that it knows which of the calculated frequencies to use) or NOMO and one/two input files, specified by keywords: PARU integer is a card format parameter file where all (not necesary parameters that are allowed to vary during the if CHARges) fitting are given a value 1.0. Non-free parameters may either be given a value 0.0 or left out. Parameters given integer values other than 0.0 or 1.0 will be grouped together and covaried during the fitting; all parameters with the same value in this file will be kept equal, so even if you only have a few free parameters, each may actually correspond to a whole group of related parameters. If charges are to be optimized the keyword CHARges must be given. The variable charges are flagged in the WMAIN array which must be set (by reading it in together with the coordinates or by using the SCALAR WMAIN .... command) BEFORE the optimizer is invoked. The same rules apply here as for the other parameters: WMAIN(i)=0 => fixed charge 1 variable n grouped variable The total charge is conserved during the fitting, and (for now) the sign of each individual charge is also kept. The <real> following CHARges is used as the weight QW in the conservation constraint: QW*(QTOT-QINIT) DATU integer is a formatted file with the observed frequencies (not necessary (cm**-1) one per line (a CHARMM title is allowed in the if NOMO) beginning of the file). The routine will vary the free parameters to get the calculated freqs. to match these as closely as possible. NOTE that the assignment will be the simplest possible, frequencies will be matched in the order they are given, unless keyword DOTP is used (see below). If keyword WEIG or MASS is given each mode may be given a weight as a second real number on each data line. If a weight is absent and WEIG or MASS is specified a default weight of 1.0 will be used. This puts the mode weights on level with the standard coordinate weights (which are SQRT(WW/3) for each X,Y&Z where WW=MASS or WMAIN or MASS*WMAIN, depending on MASS and WEIG). The optimization may be further controlled with the following optional keywords: KEYWORD DEFAULT Function -------------------------------------------------------------------------- TOLSsq 0.001 Stop iterating when the relative change of the sum-of-squares (SSQ) is less than TOLSsq. TOLPar 0.001 Stop iterating when the relative change in each of the parameters is less than TOLPar. MINIt 1 Minimum number of iterations MAXIt 99 Maximum number of iterations. DIFZ 0.001 Ratio used in computing finite differences: (F(par(j))-F(par(j)+h))/h, h=abs(difz*par(j)) PRLEvel 2 Amount of output: 0=No output, 1=Basic, 2=Adequate 3=Very much (NRAPH & NORMDS output in each step). NRSTep 5 Number of steps to take in the reminimization of the structure at new points in parameter space. NRTOlstp 0.0002 Convergence criterion for NRAPH/ABNR step length NDSTep 1 (NRAPH) Number of steps to minimize for derivatives 5 (ABNR) PRUNit 6 Unit number to which output goes WEIG FALSE Use WMAIN for coordinate weighting MASS FALSE Use AMASS as coordinate weights IF both WEIG and MASS: WMAIN*AMASS is used. IF neither: No weighting at all, not even for modes. SYMMetric FALSE Use symmetric instead of forward differences NOMOdes FALSE DO NOT USE mode data COORd FALSE Use coordinate data ABNR FALSE Use the ABNER minimizer (so that systems > 200 atoms can be fitted to observed coordinates) CHARges FALSE Vary charges on atoms for which WMAIN > 0 when routine invoked NPAR 50 Maximum number of free parameters for which space will be reserved USERobs 50 Maximum number of user observables for which spcae will be reserved. DOTProd FALSE Each calculated mode will be assigned to the observed mode with which its overlap |<Qobs|Qcal>| is largest. For this to work you allocate space for twice the number of modes you will use, and read in the observed eigenvectors in the upper half of the eigenvector matrix before invoking OPTI. If the i:th observed eigenvector is zero, then the calculated i:th mode will be assigned to it. Currently force constants and equilibrium values for: ANGLes, BONDs, DIHEdrals, HBONDs, IMPRopers and CHARges may be free. There is also a user-hook into the module through subroutine USERF(MODE,NP,PAR,NATOM,X,Y,Z,ISLCT,NU,UY,NW,UW). With this routine you may include all sorts of information other than mode frequencies and coordinates in the fitting procedure. Check the code of the dummy USERF in USERSB.FLX for details of the calling sequence.

Information and HTML Formatting Courtesy of:

NIH/DCRT/Laboratory for Structural Biology

FDA/CBER/OVRR Biophysics Laboratory

Modified, updated and generalized by C.L. Brooks, III

The Scripps Research Institute