Table Of Contents

Previous topic

Combined Quantum Mechanical and Molecular Mechanics Method Based on MNDO97 in CHARMM

Next topic

Multi Scale Command: MSCALE

This Page


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 interesting directions (STEP)
  • the vibrational analysis in presence of Drude particles

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.


[SYNTAX MOLVib command]

    MOLVib NDI1 int NDI2 int NDI3 int
           [NATOm int] [MAXSymbol int] [NGMAx int] [NBLMax int]
           [IZMAx int] [NOTOpology] [SECOnd] [PRINt] [FINIte real]

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
 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.

FINIte - finite step size (default 0.01) for numerical calculation of
  second derivatives, which is automatically invoked when Drude particles
  are present in the system. Note, Drude coordinates are reset to atomic
  centers upon leaving the MOLVIB analysis.

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


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 ::== keyword [<int1>] [<int2>] [<int3>] [<int4>]
                  format: A4,6X,4I5

data-struc ::== one of the MOLVIB input data structures; defined by the

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
FINI step size for numerical icalculation of second derivatives (applicable to classical Drude oscillator polarizable model only) Note, Drude coordinates are reset to atomic centers upon leaving the MOLVIB analysis.
END end input section, perform MOLVIB calculations and

This section gives a more detailed explanation of the keywords and the associated 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
                     <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
              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.
                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.

   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
            <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
                = 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,5,6 - internal coordinates
           ITYP = 1 - I-J  bond stretch

                   I --- J

           ITYP = 2 - I-J-K  angle bend

                     / \
                    /   \
                   I     K

           ITYP = 3 - I-L bond angle with J-K-L plane (Wilson wag)

                 I --- L

           ITYP = 4 - angle between IJK abd JKL planes

                  J --- K

           ITYP = 5 - I-J-K linear bend in IJL plane
           ITYP = 6 - I-J-K linear bend perpendicular to IJL plane

               I --- J --- K
                 .       .
                   .   .

             Note: For ITYP = 5 and 6, atom L can be omitted, in which
                 case the reference plane will be defined arbitrarily
                 based on the cartesian basis.

           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)   ! that's 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
             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.

 FINI  - a finite difference step size can be specified (default is 0.01
             Angstrom). Applicable to classical Drude oscillator systems only.
             MOLVIB supports such systems by automatic switching to numerical
             second derivatives when coordinates of Drude particles are
             self-consistently adjusted to positions of real atoms. Note,
             Drude coordinates are reset to atomic centers upon leaving the
             MOLVIB analysis.

 END   - end input section, perform MOLVIB calculations and
                     return to CHARMM.


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.