The MOLVIB Module of CHARMM By K.Kuczera & J.Wiorkiewicz-Kuczera, May 1991 MOLVIB is a general-purpose vibrational analysis program, suitable for small to medium sized molecules (say of less than 50 atoms). For larger systems the detail of description may be too great. The main options are: - the vibrational problem in internal coordinates (GF) - the vibrational problem in cartesian coordinates (GFX) - analysis of GAUSSIAN program output (GFX,GAUS) - analysis of dependencies in internal coordinate sets (G) - canonic force field calculations (KANO) - crystal normal mode analysis for k=0 (CRYS) - generating cartesian displacements along some interesing directions (STEP) The different options use mostly the same package of subroutines called in different order. New applications may thus be easily added when necessary. Of special interest is the symbolic PED analysis package, enabling a clear and condensed overview of the usually complex PED contributions. * Menu: * Syntax:: Syntax of the MOLVIB command * Function:: Purpose of each of the keywords * Input:: MOLVIB Input Description

[SYNTAX MOLVib command] MOLVib NDI1 int NDI2 int NDI3 int [NATOm int] [MAXSymbol int] [NGMAx int] [NBLMax int] [IZMAx int] [NOTOpology] [SECOnd] [PRINt]

The following section describes the keywords of the MOLVib command. NDI1,NDI2,NDI3 are the MOLVIB variables NQ, NIC, NIC0; their definition here effectively replaces the MOLVIB 'DIM ' card. Two cases: a. Molecular vibrations NIC0 - number of primitive internal coordinates (PIC's), this must correspond to the number of entries following the 'IC' card NIC - number of IC's left after transformation by first U matrix If only one U matrix is used, this should be the same as NQ; if no U matrices used, NQ=NIC=NIC0 NQ - number of vibrational degrees of freedom. Usually this is the famous number 3*Nat-6 (3*Nat-5), but also separate symmetry blocks of the vibrational Hamiltonian may be entered. b. Crystal vibrations NIC0 - no. of primitive molecular coords (MC), i.e. external coords + primitive IC's NIC - no. of vibrational degrees of freedom = 3*NAT, where NAT is the total no. of atoms in unit cell NQ - here =NIC NATOm - defines the number of atoms in the system. To be used only in conjunction with the NOTO flag. If NOTO is not provided, the number of atoms from the PSF will be used in MOLVIB and will override any values provided here. IZMAx - needed for 'CRYS' option only - specifies maximum number of molecules in unit cell. Default is 10. MAXSymbol, NGMax, NBLMax - dimensions for PED analysis arrays. They specify the maximum number of symbols, coordinate groups and symmetry blocks, respectively. Defaults are NQ (NDI1) for all three. It is recommended not to modify these defaults. NOTOpology - if flag is present, CHARMM data structures will not be used, all information required is to be read in inside the module. If flag is absent, cartesian coordinates, atomic masses and cartesian force constants from CHARMM may be passed to MOLVIB, as needed. SECOnd - calculate second derivatives (force constants in cartesian coordinates) and pass them to MOLVIB. This is done through a call to the CHARMM routine ENERGY, so all preconditions for energy (and second derivative) calculations must be met. PRINT - flag for test printout of the CHARMM second derivatives being passed to MOLVIB.

Input Description This data is processed by subroutine MOLINP. As the CHARMM command parser is not used, this input does not conform to CHARMM standards, e.g. - parameter substitution will not work - the STREAM command will not work, all commands will be read from the current input stream - OPEN, READ, WRITE, etc. commands will not work - most entries are not free format [SYNTAX MOLVIB input] The MOLVIB input consists of a series of blocks; each block consists of a command and an (optional) data structure; i.e. it has the form: command-spec [data-struc] command-spec ::== keyword [<int1>] [<int2>] [<int3>] [<int4>] format: A4,6X,4I5 data-struc ::== one of the MOLVIB input data structures; defined by the keyword. The list of currently supported keywords folows. One of the first group of keywords must be used first in order to define type of calculation. Keyword Interpretation G - perform redundancy analysis GF - solve standard Wilson GF problem GAUS - choose GAUSSIAN analysis option GFX - vibrational problem in cartesian coordinates KANO - determine canonical force field CRYS - crystal vibrations for k=0; STEP - generate cartesian displacements in a given For the remaining keywords, the order is arbitrary: Keyword Interpretation CART - read in cartesian coordinates MASA - interpret fourth column of cartesian coord input as A numbers MASZ - interpret the above column as Z numbers UMAT - read in U Matrix for similarity transformation FMAT - read in F matrix LX - read in cartesian eigenvectors IC - read in internal/external coordinate definitions PRNT - set print level TEST - set print level NULL - control card for 'G ' option with IGLEV=2 PED - read in PED data structure SCAL - read in scale factor for F matrix TRRM - remove translational and rotational contributions to LX MNAT - read in the numbers of atoms for each molecule in unit cell IFTR - specifies the dimension (and type) of F matrix SYMM - read in symmetry blocking data EXPF - read in reference frequencies for the system END - end input section, perform MOLVIB calculations and This section gives a more detailed explanantion of the keywords and the assocaited data structures. keyword Interpretation G - perform redundancy analysis <int1> == IGLEV IGLEV=1 - diagonalize G and write out eigenvalues and eigenvectors IGLEV=2 - additionally generate a set of null and independent coordinates orthogonal to the initially specified ones GF - solve standard Wilson GF problem GAUS - choose GAUSSIAN analysis option GFX - vibrational problem in cartesian coordinates KANO - determine canonical force field <int1> == ICANO ICANO=1 - preliminary run, just to output the FR matrix; one of the other keywords must follow GF, GFX or GAUS - so that FR is evaluated or just read in as part of those processes. ICANO=2 - evaluation of the canonic force field FR* N.B. No U matrix allowed here. Give: DIMensions, CARTesian coords, IC's and FMAT. CRYS - crystal vibrations for k=0; <int1> == IZMOL, <int2> == IFCRYS IZMOL - no. of molecules in unit cell IFCRYS=0 (default) - calculation analogous to GFX STEP - generate cartesian displacements in a given direction. <int1> == IFSTEP <int2> == ISTCOR <int3> == IFFMAT <int4> == IFLMAT Additionally, the card following the 'STEP' card contains the value of STPSIZ (real,free format) IFSTEP=1 - cartesian eigenvector no. ISTCOR (IFSTEP=2 - internal eigenvector no. ISTCOR, not implemented) (IFSTEP=3 - internal coordinate no. ISTCOR, not implemented) STPSIZ - step size, e.g. the transformation is X(I)=X(I)+STPSIZ*LX(I,ISTCOR) for cart. eigenvectors where LX the columns of LX are normalized. IFFMAT,IFLMAT - determine the starting point of the calculation: IFFMAT=0 and IFLMAT =1 - start from LX =2 - start from LS IFLMAT=0 and IFFMAT =1 - start from FX =2 - start from FS CART - cartesian coordinates for NAT atoms will follow <int1> == unused <int2> == IFC In MOLVIB <int1> usually used to define the number of atoms NAT. In the CHARMM version, NAT is specified on the MOLVIB command line (if NOTO flag is used) or is read from the PSF (if NOTO is absent). IFC - specifies format for cart. coords: IFC=0 free format, four real numbers per line X, Y, Z, and MASS (see below). IFC=1 CHARMM format; only atom entry lines, no titles or NATOM field, mass information in WMAIN field. N.B. For the 'GAUS' option use GAUSSIANxx CMS coordinates. FOR THE 'GFX ' option use GAUSSIANxx coordinates in the Z-matrix orientation. Mass specification : (1) enter mass in amu as fourth real number in entry line for each atom. (2) instead of mass place atomic number Z or mass number A as fourth real number and subsequently use a 'MASZ' or 'MASA' control cards. NB. For 'CRYS' NAT should be equal to no. of atoms in unit cell. MASA - interpret fourth column of cartesian coord input as atomic mass numbers (A) ; useful for isotopes, e.g. a mass of 2.0 will designate D, mass of 15.0 - 15N etc. MASZ - interpret the above column as atomic numbers (Z) UMAT - read in U Matrix for similarity transformation <int1> == IFU <int2> == INU <int3> == IUU <int4> == IZU IFU - defines format =0 Schachtschneider/Snyder format only supported INU = 1/0 - normalize/dont normalize rows of U IUU - defines FORTRAN unit for U read if left blank, unit input stream will be used if >0 then the data should be provided in the correct FORTRAN file IZU - multiplicity; usually IZU=(no. of molecules in unit cell). IZU.GT.1 turns on autogeneration of U for whole unit cell from the provided values for the first asymmetric unit. (see SUBROUTINE RDMAT in MOLVIO) Details: Two, one or none U matrices may be supplied on input. These are (generally) rectangular matrices which perform linear transformations on internal coordinate sets, of the type S=U*R ( or S(i) = {sum over j} U(i,j)*R(j) ), with S - final, and R initial coordinate sets. The function of the U matrix is e.g. to transform from primitive IC' s (of which there are NIC0>NQ) to a set of independent IC's NQ in number, or to scale the IC's by a factor (useful when trying to reproduce vibrations reported in the literature, as different research groups use different definitons of angle or dihedral IC's). If two U matrices are given, then the IC's (and the B and G matrices) are sequentially transformed using first U1, then U2. The F matrix is assumed to be expressed in the final IC's on input, and is not transformed (except for the 'KANO' option - see 'IFTR'). FMAT - read in F matrix, (the second derivatives of energy wrt coordinates) <int1> == IFF <int2> == ISF <int3> == IUF IFF - specifies format = 0 - Schatschneider/Snyder format = 1 - GAUSSIANxx format N.B. remember to use 'Z matrix' oriented cartesian coords. = 2 - CHARMM formatted SECO file format IFS = 1/0 - symmetrize/dont symmetrize (upper triangle assumed on input) IUF - FORTAN unit no., as for 'UMAT' (see RDMAT, RFORC, RDSECO in MOLVIO) LX - read in cartesian eigenvectors <int1> == IFL <int2> == unused <int3> == IUL IFL - specifies format = 0 - GAUSSIANxx format (see SUBROUTINE REIGEN) all 3*NAT eigenvectors read in N.B. remember to use 'standard' (or 'CMS') oriented cartesian coordinates = 1 - CHARMM binary format (see SUBROUTINE REIGCH) only the NQ=3*NAT-6 "vibrational" eigenvectors are expected by REIGCH; use "WRITE NORM 7 THRU ..." command to achieve this. NB. Binary files are machine specific. IUL - FORTAN unit no. from which to read , aas in 'UMAT' IC - read in internal/external coordinate definitions; <int1> == IZIC Five integers will be read from NIC0 lines in free format; each line contains: ITYP,I,J,K,L - specify type and four atom numbers as defined in cartesian coordinates Note: it is necessary to add zeros in unused fields. IZIC - multiplicity, usually = no. of molecules in unit cell. IZIC.GT.1 turns on autogeneration of internal/external coordinates for unit cell from the ones provided for the first asymmetric unit. ITYP=1,2,3,4 - internal coordinates ITYP = 1 - I-J bond stretch I --- J ITYP = 2 - I-J-K angle bend J / \ / \ I K ITYP = 3 - I-L bond angle with J-K-L plane (Wilson wag) K / / I --- L \ \ J ITYP = 4 - angle between IJK abd JKL planes I \ \ J --- K \ \ L For details : a) see SUBROUTINE BMAT b) see Wilson,Decius,Cross section 4.1, substituting their atom numbers with: ITYP=1 (12) -> (IJ) 2 (123) -> (IKJ) ! thats right ! 3 (1234) -> (IJKL) 4 (1234) -> (IJKL) A good reference for standard definitions of independent internal coordinates for a wide selection of chemical groups is: P.Pulay,G.Fogarasi,F.Pang & J.E.Boggs, JACS 101, 2550 (1979) For the 'CRYS' option, the external coordinates are defined here; their codes: ITYP=11 - x translation ITYP=12 - y translation ITYP=13 - z translation ITYP=14 - x rotation ITYP=15 - y rotation ITYP=16 - z rotation In this case the I field should hold the consecutive number of the molecule in the unit cell (consistent with MNAT data). PRNT - set print level <int1> == IPRNT IPRNT=0 - minimal printout IPRNT=5 - maximum printout [default is 2] TEST - equivalent to 'PRNT' with IPRNT=4 NULL - control card for 'G ' option with IGLEV=2 <int1> == NULL <int2> == NSTRT NULL = the number of orthonormal vectors for the null space to be read from the U2 matrix NSTRT = the number of starting vectors for the Gram-Schmidt procedure in the vibrational space Note: If any null coordinates are known, they should be orthonormalized and placed in the first NULL rows of U2. The program will then write out the complete set of orthonormal coordinates spanning the null space, starting with the ones provided. If NSTRT.GT.0 a completely independent calculation will be performed in the vibrational space. In that case, the NULL+1,...NULL+NSTRT rows of U2 should contain the known coordinates of the vibrational space orthogonal to each other and the redundancies (null space vectors). The program will construct an orthonormal basis of the vibrational space which is orthogonal to the redundancies, starting with the provided vectors. PED - symbolic PED analysis will be performed <int1> == NGRUP <int2> == NCUTP NGRUP - number of coordinate groups to be defined NCUTP - cutoff level; PED contributions below NCUTP % will not be printed, for clarity (default is 3%). The following cards must contain: 1. for each group I=1,NGRUP: LGRUP(I),IGRUP(I,J), J=1,LGRUP(I) - the number of coords in group and their consecutive numbers (these are the final numbers, i.e. after all U matrix operations) (20I3) 2. for each coordinate : IS,SS - its consecutive number (after all U matrix operations) and the assigned symbol. 4(I3,2X,A8,2X) - zero to four entries per line; blank fields skipped, negative IS value to end this input section. Only the first coord of each group needs a symbol defnition, the rest are set to this string; contributions from the whole group are added up and printed beside the group symbol. SCAL - scale the F matrix Fij' = FACT*Fij; the real value of FACT will be read from next line (F10.6). TRRM - remove translational and rotational contributions from cartesian coordinate vibrational eigenvectors. (currently used only for GAUS) MNAT - lines following this card will contain the numbers of atoms of the individual molecules comprising the unit cell (or molecular aggregate) in 20I3 format. Application - makes possible external coordinate use in vibrational analysis of mixed crystals or molecular aggregates (use CRYS option in both cases). The value of IZMOL should already be defined for this card IFTR - specifies the dimension (and type) of F matrix <int1> == IFTRAN = 0 - F is in primitive ICs R, NIC0xNIC0 = 1 - F is in S1=U1*R, NICxNIC = 2 - F is in S2=U2*U1*R coords, NQxNQ If card is not given, default IFTRAN=NUMAT is assumed (works only for 'KANO' option) SYMM - use symmetry (in symbolic PED analysis only) <int1> == NBLOCK It is assumed that by use of similarity transformations (the U matrices), the vibrational problem has been transformed to such coordinates that the Hamiltonian (G and F) is block-diagonal. This usually happens if the coordinates form a basis for the irreducible representations of the molecular point group. The following cards should contain the data: IBLOCK(I),I=1,NBLOCK - sizes of consecutive blocks (coordinate numbering is as for PED analysis, i.e. after all U matrix transformations) SBLOCK(I),I=1,NBLOCK - block symbols (e.g. representation names) EXPF - read in reference frequencies for the system Frequencies should be in ascending order (if 'SYMM' is present, the ordering should be separate within each block). The frequencies from MOLVIB will be printed out side-by-side with the reference set, differences and an rms deviation will computed. (If 'SYMM' is present, a separate analysis will be performed for each block). Format: free, 1 real value per line. END - end input section, perform MOLVIB calculations and return to CHARMM. Note: the Schactschneider/Snyder format This format is very useful for i/o of sparse matrices (or small and not so sparse ones). The basic format is: 4(2I3,F12.6) The two integer fields specify the row and column number, the real field - the value of the array element. Any elements not explicitly specified are set to zero. Each line of input may contain 0-4 entries, blank lines are ignored, a negative value for the column number terminates input. See subroutine RDMAT in MOLVIO.

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