Emory University | Department of Chemistry | Bowman Group | Emerson Center


(Version 3.4 2/99)

Stuart Carter and Joel Bowman  


This manual describes version 3.4 of a code, MULTIMODE, which does Vibrational Self Consistent Field (VSCF) and several types of "CI" calculations for the rovibrational energies and wavefunctions of polyatomic molecules, transition states, clusters, etc. The code is based on the Watson Hamiltonian for non-linear molecules. (A version for linear molecules will be available soon.) The treatment is exact for zero and non-zero total angular momentum, J, and uses general potentials. For non-zero J, an approximate, adiabatic treatment of rotation can also be used. This approximate method, which is very efficient compared to an exact calculation, has been tested and should be very reliable for small values of J, except in those cases where there are strong Coriolis couplings.

The treatment of the potential is novel, and permits applications to many-mode systems. The potential is represented as an expansion in the number of simultaneously coupled modes. Currently this number can range from one to four. The computational effort increases significantly as this number increases; however, by doing calculations with increasing values of this number, convergence of this representation of the potential can be tested.

 Table of Contents

Bowman Page

Bowman Group Home Page






 We thank the National Science Foundation for Financial Support and S. Carter also thanks the Office of Naval Research for support. 




Brief description of theory

User -Specific Routines

User Guide 

Input Parameters

Sample Inputs and Outputs

Applications to Date

How to Order MULTIMODE

Bowman Page

Bowman Group Home Page

Back to Top











 Brief description of the theory


1.1 Non-rotating molecules

For a full account of the theory see References 1, 2 and 4 and references therein.

Here we give a brief recap of the theory. We begin with the Watson Hamiltonian for non-rotating, non-linear molecules:


which we assume is familiar to the reader. We represent the potential V by the following expansion:

where the one-mode representation of the potential contains only terms, the two-mode representation contains those terms plus the terms, etc. Thus, explicit coupling of up to a maximum of 4 modes can be done. To be more precise is given by. The potential is obtained from a user-supplied potential (given as a function of Cartesian coordinates) via the exact transformation between normal coordinates and Cartesian coordinates.


1.1.1 VSCF Theory

The many-mode vibrational wavefunction is written as a "Hartree product"

Applying the variational principle to the modals yields the VSCF equations

where V is the n-mode representation of the potential (n=1,2,3,4), Tc is the Coriolis coupling operator, and Ti is the Cartesian kinetic energy operator in the Watson Hamiltonian above. These equations are solved by a standard iteration method beginning with the zero-order normal-mode, harmonic oscillator approximation. Convergence is based on convergence of the eigenvalues. The approximate energy of the state is simply the expectation value of HJ=0.

The VSCF modals are expanded in a primitive basis of normal-mode harmonic-oscillator functions. Matrix elements over the (effective) potential are done by Gauss-Hermite quadrature. Optimized "HEG" quadrature points may be used to great advantage. The potential and Coriolis grids at the quadrature points are evaluated first and written to disk. The VSCF eigenvalue equations are solved by a standard diagonalization method. The order of the matrix to be diagonalized is typically quite small, e.g., of the order of 10 - 20. The CI equations are written using standard linear algebra methods and the diagonalization of the H-matrix is done using GIVENS/HOUSEHOLDER routines since the number of eigenvalues requested is usually much less than the order of the CI H-matrix. An optional interative diagonalization routine may be used. This is useful for very large matrices.

The code has the option to do a complete normal mode analysis and transformation to a principal axis system, starting with a general, user-supplied potential, and stationary geometry. The user may supply his/her own normal-mode eigenvectors and frequencies as input. In this case it is assumed that the eigenvectors are referenced to a principal axis system, at the reference, equilibrium geometry. The code will evaluate all the relevant Coriolis coupling coefficients, etc. from this input. The code also has an option to do a normal mode and principal axis analysis only.


1.1.2 VSCF+CI Theory

The many-mode wavefunction is expanded in a basis of VSCF states. These states are not orthogonal and this results in a generalized eigenvalue problem. This problem is solved by standard methods. The expansion can be done in several ways. An expansion in increasing energy is recommended to obtain true variational upper bound approximations to the exact eigenvalues and eigenfunctions.


1.1.3 VCI Theory

The many-mode wavefunction is expanded in terms of the eigenstates (also known as the virtual states) of a given VSCF Hamiltonian. Usually the ground state VSCF Hamiltonian is used. This basis is orthonormal and results in a standard eigenvalue problem. The resulting eigenvalues and eigenfunctions are true variational upper bound approximations to the exact eigenvalues and eigenfunctions. The code many options to develop this basis, and a basis in which selected modes are unexcited is an example of one such basis.


1.2.1 Exact treatment of rotating molecules


The Hamiltonian for non-zero total angular momentum J is,

where are the components of the total and vibrational angular momentum operators, respectively, µ is the effective reciprocal inertia tensor. The eigenvalues and eigenfunctions of this Hamiltonian are obtained as described in detail in ref. 4 below. These calculations are done using the N-mode representation of the potential as described above.

The computer time grows roughly quadratically with J, and so an approximate treatment of rotation, for which the computational effort grows linearly with J, can also be used, as described next.


1.2.2 Adiabatic approximation for rotating molecules

The adiabatic rotation Hamiltonian is given by

where is the Watson Hamiltonian for non-rotating molecules, given above, and is an eigenvalue of the asymmetric top Hamiltonian , where A, B, and C are rotation constants which depend on the nuclear configuration. The sum of V and form an effective potential. The eigenvalues and eigenfunction of HAR are obtained using the VSCF, VSCF+CI, or VCI methods.



(1). "Vibrational Self-Consistent Field Method for Many-Mode Systems: A New approach and application to the vibrations of CO adsorbed on Cu(100)", S. Carter, S. J. Culik and J. M. Bowman, J. Chem. Phys. 107, 10458 (1997).

(2). "The adiabatic rotation approximation for rovibrational energies of many-mode systems: Description and tests of the method", S. Carter and J. M. Bowman, J. Chem. Phys. 108, 4397 (1998).

(3). "Vibrational spectrum of (CO)2 on Cu(100): Quantum calculations with 18 coupled modes", F. Dzegilenko, J. M. Bowman and S. Carter, J. Chem. Phys. 109, 7509 (1998).

(4) "Extensions and tests of 'multimode': A code to obtain accurate vibration/rotation energies of many-mode molecules", S. Carter, J. M. Bowman, and N. Handy, Theoretical Chem. Accnts. 100, 191 (1998).

(5) Non-separable transition state theory for non-zero total angular momentum: Implications for J-shifting and application to the OH+H2 reaction, J. M. Bowman and H. Shnider, J. Chem. Phys. 110, xxxx (1999).

(6) "Calculations of low-lying vibrational states of cis and trans-HOCO", J. M. Bowman, K. Christoffel, and G. Weinberg, Theochem, in press.

(7) "Non-separable RRKM theory applied to unimolecular dissociation of HCO2 H+CO2, J. Phys. Chem., K. Christoffel and J. M. Bowman, submitted.

(8) "Variational calculations of ro-vibrational energies of CH4 and isotopomers in full dimensionality using an ab initio force field", S. Carter, H. M. Shnider, and J. M. Bowman, J. Chem. Phys., submitted.


Table of Contents

Bowman Page

Bowman Group Home Page

Back to Top
















User-Specific Routines


The body of the code is `multimode.3.4.f'; user-specific codes are `memo.f' and `user.f'. A brief summary of the latter two codes is as follows:




The module `memo.f' sets the size of RAM to be used, and it also opens all input/output files. The RAM is set by the statement




If the current RAM specified by MAXSIZ is exceeded during a run of `multimode', the program stops with the message:




followed by the current and projected size of MAXSIZ. The user should recompile `memo.f' with the new value of MAXSIZ and try again.

There are several files used in `multimode'. All files are opened in `memo.f'; the user must modify these OPEN statements in order to write them onto a disc of his/her choosing. All files except units 1 and 2 are written UNFORMATTED.

The definitions of the files are as follows:

Unit 1: Contains the standard input data for `multimode' in FREE FORMAT


Unit 2: Contains the FORMATTED output from `multimode'


Unit 20: Lanczos file holding hamiltonian half-matrix


Unit 21: Basic rovibrational integrals for one-mode coupling


Unit 22: Basic rovibrational integrals for two-mode coupling


Unit 23: Basic rovibrational integrals for three-mode coupling


Unit 24: Basic rovibrational integrals for four-mode coupling


Unit 25: Lanczos work file


Unit 30: Lanczos work file


Unit 35: Optional Lanczos work file


Unit 40: Contains the final VSCF information in each UNFORMATTED record

State number, assignment, J, Ka, Kc, energy


Unit 45: Contains the final CI information in each UNFORMATTED record

State number, assignment, J, Ka, Kc, energy


For BOTH units 40 and 45:

  • state number is an (integer) quantity representing the numbering of the VSCF or CI state,
  • assignment is an (integer) sequence of NMODE numbers which represent the number of quanta in each Normal mode (see INPUT Section for definition of NMODE),
  • J is the (integer) total angular momentum quantum number,
  • Ka is the (integer) asymmetric top quantum number along the principal Z-axis,
  • Kc is the (integer) asymmetric top quantum number along the principal X-axis,
  • energy is the (double precision) energy of the VSCF or CI state in cm-1 relative to the relative to the reference energy EVLJO (see INPUT Section).
  • In the event that EVLJ0 is given a value of ZERO, the Band origin of the zero point level is written as its absolute energy in cm-1. In this case, all other energies are relative to this zero point energy.


Unit 50: Scratch file used to store VSCF wavefunctions if followed by an SCFCI calculation


Unit 55: Optional Lanczos work file


Unit 58: VCI dump scratch file


Unit 59: VCI dump scratch file


Unit 60: HEG work file


Unit 61: Rigid-rotor energies at each one-mode grid point


Unit 62: Rigid-rotor energies at each two-mode grid point


Unit 63: Rigid-rotor energies at each three-mode grid point


Unit 64: Rigid-rotor energies at each four-mode grid point


Unit 71: Potential energy at each one-mode grid point


Unit 72: Potential energy at each two-mode grid point


Unit 73: Potential energy at each three-mode grid point


Unit 74: Potential energy at each four-mode grid point


Unit 81: Vibrational coriolis coupling data at each one-mode grid point


Unit 82: Vibrational coriolis coupling data at each two-mode grid point


Unit 83: Vibrational coriolis coupling data at each three-mode grid point


Unit 84: Vibrational coriolis coupling data at each four-mode grid point


Unit 91: Rotational coriolis coupling data at each one-mode grid point


Unit 92: Rotational coriolis coupling data at each two-mode grid point


Unit 93: Rotational coriolis coupling data at each three-mode grid point


Unit 94: Rotational coriolis coupling data at each four-mode grid point


For the sets of files (units 21, 22, 23, 24), (units 61, 62, 63, 64), (units 71, 72, 73, 74), (units 81, 82, 83, 84), (units 91, 92, 93, 94), the sizes grow permutationally from one-mode to four-mode coupling.



The module user.f contains two routines which may be modified by the user in order to interface with his/her potential energy routines.






The calling sequences to these routines are:






where V is a DOUBLE PRECISION variable which must contain the potential energy in hartrees on exit. NATOM is an INTEGER variable denoting the number of atoms in the molecule, and XX, RR are DOUBLE PRECISION arrays of dimensions.




XX contains the cartesian coordinates of the NATOM atoms, where XX(NATOM,1), XX(NATOM,2), XX(NATOM,3) are the x-, y-, z-coordinates respectively. RR is reserved space for the internuclear distances between the atoms. If bond distances are required, the user may call the utility routine BONDS from within GETPOT.




The routines USERIN and GETPOT must be present ALWAYS.


USERIN can be used to read (for example) any data that the user may require to define his potential. If the user wishes to append the standard input data file to `multimode' with his own data, and/or output any information to the standard output file produced by `multimode', he may use INP and/or IOUT in the COMMON block




respectively. If no such data are required, then USERIN may be treated as a dummy routine which immediately exits, but it is ALWAYS called by `multimode'.


GETPOT is used as an interface with the user's own routine(s) which evaluate the potential at the cartesian coordinates held in XX(NATOM,3) which are passed to GETPOT (or at the internuclear bond distances held in RR(NATOM,NATOM) which are returned from a user-initiated call to BONDS, if required).

If the arrays XX or RR are required by the user in his own routine(s), they must be passed, along with NATOM, in the calling sequence(s) to his routine(s), AND MUST BE DIMENSIONED: XX(NATOM,3),RR(NATOM,NATOM) IN HIS ROUTINE(S).


The COMMON block




contains some useful fundamental constants:

WAVENM: Conversion from hartrees to wavenumbers

ATTOJ: Conversion from atojoules to hartrees

BOHR: Conversion from bohr to Angstroms

ELMASS: Conversion from Unified atomic mass units to electron mass units

RAD: Conversion from radians to degrees


which have been initialised to the latest IUPAC standards. These are for use by the user, if required. These constants require, for example, that

hartrees*WAVENM = wavenumbers, etc.

The code `multimode' contains several machine-specific calls to timing routines. Currently computers supported are CRAY, IBM, and SGI. The user must comment out or delete timing routines to machines not being used. Also, there are calls to 'FLUSH' that are currently commented out. This routine flushes the output buffer.


Note...`multimode' calls timer routines for CRAY, IBM, SGI (these are also available from us).


 Table of Contents

Back to Top


User Guide


This Section describes in some detail the facilities offered by MULTIMODE. It should be read in conjunction with the two Sections 'Input Parameters' and 'Sample Inputs and Outputs'. In the Section 'Input Parameters' all data records are specified in full. In the Section 'Sample Inputs and Outputs' various facilities of MULTIMODE are tested, and in particular, key input parameters are highlighted. It will be seen that the test programs fall into 'five' basic categories:


1..AUX: This category is mainly concerned with basic set-up facilities,

in particular, the features associated with the establishment of the

Normal Modes.


2..SCF: This category demonstrates how VSCF (Vibrational SCF)

calculations are performed.


3..SCFCI: This category demonstrates how CI calculations can be

performed using the (non-orthogonal) VSCF wavefunctions as a CI basis.


4..VCI: This category demonstrates the use of (orthonormal) virtual zpe

wavefunctions in near-variational CI calculations.


5..LAN: This category demonstrates the use of an iterative

diagonalisation technique for use with very large VCI hamiltonian











In order to run MULTIMODE, the user must supply a subroutine which returns the value of the potential energy of his molecule for any arbitrary values of the cartesian coordinates of the atoms associated with the distortions from cartesian coordinates of the atoms associated with the distortions from equilibrium (see the Section 'User-Specific Routines'). It is advised that a careful check of his routine is first carried out to establish the presence of any spurious 'holes' which will play havoc with numerical integration of the potential. If his potential appears to be well-behaved, then the routines USERIN and GETPOT can be compiled in his 'user.f' module.

The next step is to choose a convenient size for the work array in the module 'memo.f' (see the Section 'User-Specific Routines'). This is done by setting the size of MAXSIZ in the PARAMETER statement which appears in the first line of MAIN. A recommended value is around 3000000. The value given to MAXSIZ determines the size of the work array which is used to dynamically allocate memory. If it is set too low, MULTIMODE will halt with a suggested alternative value. But remember; if it is set far too high, space will still be required on disc to hold unnecessary array space, and a excessive amount of memory will be sought.

Start with MAXSIZ as low as possible.

Since GETPOT is called many times during a run of MULTIMODE, it is advised to use the highest possible level of optimisation supported by the FORTRAN compiler when compiling 'user.f'. Some SGI compilers produce incorrect results if 'memo.f' is compiled at -O2, however, and we strongly recommend that the option -O1 is used when compiling 'memo.f'. All remaining modules can be compiled at the highest level of optimisation, and when all modules are linked, the user is ready to start. We recommend that he carries out test runs using the data files given in the Section 'Sample Inputs and Outputs' and which are summarised below. This will enable him to gain experience in setting up data files, as those in 'Sample Inputs and Outputs' will have to be modified for use with his choice of molecule. The following categories have been designed to cover features that the user would naturally wish to address in order to gain familiarity with MULTIMODE. The test examples are all for the 4-atom (NATOM = 4) molecule H2CS, for which a potential in internal coordinates (IWHICH > 0) is supplied. In all cases, coupling of a maximum of three modes (ICOUPL = 3) in the potential and coriolis terms is required.


1. Basic Set-up Facilities (AUX)

These test examples are mainly concerned with features relating to the Normal Coordinates, and are terminated after initial tests have been performed by setting ISCFCI < 0. For these initial test examples, many of the input parameters are ignored. Only those relevant to the specific test under discussion are highlighted below.


MULTIMODE Version 3.4 employs the Watson Normal Coordinate hamiltonian for non-linear molecules, and therefore requires an input equilibrium geometry in cartesian coordinates, the harmonic force constants of the potential (omegas) and the corresponding displacements of the atoms (the eigenvectors of FG). It could well be that the user knows all of these quantities, and wishes to use them in a run of MULTIMODE. Test program AUX.1 is such an example. The user inputs the NUCLEAR MASSES and corresponding EQUILIBRIUM CARTESIAN COORDINATES, and indicates to MULTIMODE that he wishes to use his own Normal Coordinates by setting INORM = 0. Since H2CS has six Normal Coordinates, he must input NPOT = 6 values of omega (in cm-1) as IPOT,JPOT,CPOT. He must also input sets of NORMAL COORDINATE DISPLACEMENTS for each Normal Mode as x,y,z displacements for each atom in turn. Note: there MUST be consistency between the order of the NUCLEAR MASSES, EQUILIBRIUM CARTESIAN COORDINATES, NORMAL COORDINATE DISPLACEMENTS and the input omegas at IPOT,JPOT,CPOT. It is most common when running an FG analysis to choose a 'convenient' cartesian coordinate system to describe the equilibrium configuration (one axis along a bond with an atom at the origin is common). The Watson hamiltonian, however, requires that the input axes are the principal axes of rotation, and to instruct MULTIMODE to find these from the input geometry, set MCHECK > 0. In example AUX.1, the molecule has been moved out of the yz plane on input to demonstrate this. On inspection of the output from AUX.1, it will be seen that the principal y and z axes are correctly in the plane of the molecule.

In the vast majority of cases, the user will rely on MULTIMODE to carry out his Normal Coordinate analysis. In this case he sets INORM > 0. There are now NPOT = 0 values of omega to input at IPOT,JPOT,CPOT and the NORMAL COORDINATE DISPLACEMENTS are omitted, as seen in test program AUX.2.

Finally, the user may have a potential function for which the precise position of the equilibrium structure is not known exactly, but an approximate structure is known. To obtain the correct Normal Coordinates, the principal axes must correspond to the true equilibrium, and to get MULTIMODE to find it, set MCHECK < 0. To find the true position of the minimum, use is made of the integration points generated by the input parameters NBF,MBF,NVF. These may either be Gauss-hermite or optimised HEG points as described in the Section 'Input Parameters'. For each mode, in turn, the position of the minimum is moved to Q=0 by first fitting the one-dimensional potential to a polynomial, using the three central integration points and finding the position of its minimum. This process is repeated with five, seven, ... points until the integration points have been exhausted, in which case, the 'best' geometry is used to find the 'PRINCIPAL AXES COORDINATES'. Test program AUX.3 is an example of this, and is in three parts. AUX.3a uses a modified (non-equilibrium) geometry as input, and inspection of the output reveals that, after fitting to the maximum-allowed eleven points, the true minimum has almost (but not quite) been found (the definition of the equilibrium is that Q=0 for all modes). The 'PRINCIPAL AXES COORDINATES' from AUX.3a are now used as input to AUX.3b. Because of symmetry, the new EQUILIBRIUM CARTESIAN COORDINATES have been 'tidied' to give a totally-symmetric structure for H2CS. This process must be repeated until the true equilibrium structure has been found, as is the case on inspection of the output from AUX.3b, where Q=0 for all modes.

Once the true equilibrium has been found, the resulting 'PRINCIPAL AXES COORDINATES' are transferred to the input EQUILIBRIUM CARTESIAN COORDINATES, taking into consideration any strict symmetry requirements, and some other task can be performed by first setting MCHECK = 0. In AUX.3b, for example, a VSCF job is submitted (ISCFCI = 0).


2. Vibrational Self-Consistent Field (SCF)

ISCFCI = 0 is the key parameter for VSCF calculations, and each state is found to an energy tolerance given by CONV. The SCF examples here are all to CONV = 1.D-2 cm-1. When any new system is being investigated, the potential and coriolis terms must be established at the integration grids for all of the coupled modes demanded by ICOUPL. If ICOUPL > 0, REAL*8 grids are formed and written to disc if IDISC = 0. If disc space is short, REAL*4 grids can be used by setting ICOUPL < 0. For large systems, a large amount of time can be spent writing these grids to disc, but, providing that the final Normal Modes have been found (as in the AUX examples) and the optimum number of Gauss points and/or HEG points determined by the parameters NBF,MBF,NVF have been determined, these grids need to be written ONCE ONLY. Thereafter, setting IDISC > 0 will omit writing these grids, and use the existing ones already on disc.

When performing an SCF calculation, a discreet number of vibrational states are found. There are a variety of ways to indicate those states required. The simplest way is merely to determine the lowest -NSTAT vibrational levels, by setting NSTAT < 0. In this case, there are no SCF STATE DEFINITIONS, and use is made of the omegas to calculate the states required. Test programs SCF.1 and SCF.2 are examples for which ten such states are found. SCF.1 uses Normal Modes determined by MULTIMODE, and SCF.2 uses user-supplied Normal Modes. If a more specific set of VSCF states are required, set NSTAT > 0 to indicate how many, and set the SCF STATE DEFINITION parameter. Two types of setting are supported. The first, as in test program SCF.3, is initiated by setting the SCF STATE DEFINITION < 0. This is followed by NSTAT records which define the SPECIFIC STATES. Each record consists of NMODE = 3*NATOM - 6 integers which signify the specific number of quanta in ALL the modes for all the NSTAT states to be determined. The second, as in test program SCF.4, first sets the SCF STATE DEFINITION to a default value (usually 0), and then inputs NSTAT-1 records which indicate number of quanta in corresponding mode, respectively. In test program SCF.4, all modes are given a default setting of 0 quanta (this is taken to be the first VSCF state, in this case, the zero point level), and then the next six VSCF states have a single quantum in each of the six modes, respectively. The final three states have two quanta in each of the first three modes, respectively. This default method is really only desirable for setting up fundamentals and overtones.

Test programs SCF.1 to SCF.4 all determine the VSCF states for zero total angular momentum (JMAX = 0). Test program SCF.5 calculates J=0 and J=1 VSCF states by setting JMAX = 1; J>0 energies for SCF states are obtained by using the adiabatic rotation approximation only. Currently, all J > 0 levels are assigned according to asymmetric top definitions of Ka,Kc.


3. Non-orthogonal VSCF Functions in CI (SCFCI)

ISCFCI > 0 is the key parameter for all CI calculations, and signifies the number of vibrational (JMAX = 0) or rovibrational K-diagonal (JMAX > 0) energies required, subject to the print cut-off parameter in CUT. All CI energies are printed in cm-1 relative to the base energy given in EVLJ0, which is usually set to zero for J=0 calculations. In this case, the energies are printed out relative to the J=0 zero point. For J>0 calculations, EVLJ0 is usually input as the absolute zero point energy for J=0. In all the following examples, ten CI energies are required (ISCFCI = 10) subject to a cut-off of 1.0D+4 cm-1, and only the leading CI coefficient is given in a medium print output (IPRINT = 2). An SCFCI calculation is indicated by setting ICI > 0, and this is the number of VSCF states to mix in the non-orthogonal CI basis.

In test program SCFCI.1, a J=0 calculation is performed (JMAX = 0) for the lowest fifty SCF states (NSTAT = -50), of which the first forty are used to form the CI basis (ICI = 40). These basis functions are symmetry-blocked, since the vibrations of H2CS span four C2v species (NVSYM = 4), resulting from combinations of the three distinct mode vibrations (NWSYM = 3). There are three A1 modes (NSYM = 3; ISYM = 3 4 5), two B2 modes (NSYM = 2; ISYM = 1 6) and one B1 mode (NSYM = 1; ISYM = 2). ISYM indicates the specific mode and can be obtained from inspection of the Normal Mode displacement vectors. All four symmetry blocks are required (MVSYM = 4; MWSYM = 1 2 3 4). MWSYM indicates which specific species are required, in this example all of them. All CI energies are given relative to the J=0 zero point (EVLJ0 = 0).

Test programs SCFCI.2 and SCFCI.3 perform the same task, namely an initial VSCF run for J=0 and J=1 (JMAX = 1) followed by an SCFCI calculation using the lowest twenty VSCF functions as basis (ICI = 20). The J > 0 calculations all make use of the adiabatic rotation approximation. The only difference between he two runs is that SCFCI.2 uses REAL*8 disc grids (ICOUPL = 3) and SCFCI.3 >uses REAL*4 disc grids (ICOUPL = -3). All CI energies are given relative to the J=0 zero point energy obtained from the output of test program SCFCI.1 (EVLJ0 = 5375.4543), which is obtained from the first (A1) CI block.


4. Orthonormal Virtual zpe Functions in CI (VCI)

This is the most exact algorithm in the MULTIMODE package, and is initiated by setting ISCFCI > 0 (as above), and ICI < 0. If NMAX > 0 or NMAX = 0, then -ICI means the maximum quantum of excitation (in ALL modes) in the VCI basis. If NMAX > 0, the VCI basis is restricted to a sum of quanta of NMAX. Test program VCI.1 is such an example. Each individual mode basis function can have up to three quanta of excitation (ICI = -3), but the total hartree product CI expansion function is only allowed a total sum of three quanta (NMAX = 3). As in the SCFCI test programs, the resulting CI matrix is symmetry-blocked. It will be noticed that the parameter NSTAT = -50, from earlier tests. Recall that this means that the lowest fifty VSCF states would be determined in an SCF run of MULTIMODE, starting with the zero point level. The fact that a VCI run is in progress curtails the VSCF procedure after the initial state has been determined, and the zpe virtual functions are therefore used in the VCI. If, for some reason, the user thinks that some other VSCF would be more appropriate for his VCI basis, he should use the NSTAT > 0 facility (see SCF above), and input his preferred choice of VSCF state as the initial SPECIFIC STATE. In all of the VCI tests, the leading two terms are required in the output of the vibrational CI levels (IPRINT = -2), together with the values of the CI coefficients. There is a research facility for VCI calculations (ICI < 0) whereby it is possible to dump information about wavefunctions to a disc file for further scrutiny, for example the study of transition moments, etc.. This facility is not yet released, but is controlled by the parameter IDUMP. For the time being the user should set IDUMP = 0.

If NMAX = 0, the VCI basis is unrestricted, but this may lead to extremely large CI matrices, especially as the number of modes increases. It is advised to first run a job with MATSIZ > 0. This merely instructs MULTIMODE to work out the matrix size(s) and then STOP. If they appear to be manageable, he can continue by running the same data set with MATSIZ = 0. Test programs VCI.2a and VCI.2b do this. As a matter of interest, it will be noticed from the output of VCI.2a that six modes excited to a maximum of three quanta without restriction (NATOM = 4, ICI = -3, NMAX = 0), combine to give four CI matrices each of size 1024, even when utilising full C2v symmetry, so beware!

If NMAX < 0, then -NMAX indicates the maximum number of coupled modes in the CI basis, and this method of selection gives the user maximum control over the composition of his hamiltonian matrix. In test program VCI.3, the value of NMAX = -3, and this means that two further types of input are required for up to three coupled modes in the CI basis. The first (MAXSUM) is the total sum of quanta for one-mode, two-mode, three-mode terms in the basis. The second (MAXBAS) is the maximum quantum which appears (for all NMODE virtual functions) for one-mode, two-mode, three-mode terms in the basis, respectively. Currently, the maximum number of coupled modes that are allowed in the CI basis is four (NMAX = -4); test program VCI.4 is an example of such a CI basis.

Test program VCI.5 is an example of a VCI calculation for J>0, in this case JMAX = 5. Only a single K-diagonal 'vibrational' set of functions are determined for Ka=0 (KSTEP = 0), and in the final rovibrational matrix, ten functions (ISCFCI = 10) are combined for each of the 2*JMAX+1 combinations of Ka,Kc. In the final output, twenty rovibrational energies are required for each symmetry block (NVALR = 20). In test program VCI.6, the same calculation is repeated, but in this case the K-diagonal 'vibrational' functions are evaluated for alternate combinations of Ka,Kc (KSTEP = 2), commencing with Ka=0,Kc=J. It will be seen that the final rovibrational energies of these two calculations are almost identical, and for very large values of JMAX, it is advised not to calculate K-diagonal functions for the maximum 2*JMAX+1 times, but to select a value of KSTEP which reduces the effort substantially (for example, for JMAX = 20, a value of KSTEP = 5 might be considered).


5. Iterative Diagonalisation (LAN)

It has already been pointed out that for the near-variational VCI algorithm, the size of the CI matrix can explode very rapidly. To reduce its size, a careful selection of individual basis functions can be used via the MAXSUM, MAXBAS facility described above, and for large molecules, this method is recommended at 'all' times. Even with the most careful choice of basis function, however, a matrix sufficiently large to prohibit its diagonalisation by in-core methods may well be generated. To allow such matrices to exist in the MULTIMODE method, a Lanczos/Davidson-type algorithm exists, and is initiated by setting NCYCLE > 0. This represents the number of cycles required in the iterative diagonalisation scheme, and TOLLAN is the tolerance required in the resulting eigenvalues. A half-matrix of order LANMAX is used to set up the complete hamiltonian half-matrix before writing it to disc. During the course of the diagonalisation, ISCFCI eigenvalues (for vibrational or K-diagonal steps), or NVALR eigenvalues (for rovibrational steps) are determined each cycle and therefore an internal matrix of maximum order NCYCLE*ISCFCI or NCYCLE*NVALR will require diagonalisation by an in-core method. This can be either a QL (LGIV = 0) or a Givens (LGIV > 0) algorithm. There are two algorithms present for the iterative scheme; the first (NCYCLE > 0) involves only one-dimensional arrays and requires extra disc space to store intermediate arrays; the second (NCYCLE < 0) involves two-dimensional arrays to hold these internal arrays, and therefore requires more memory but less disc space.

Test program LAN.1 diagonalises alternate K-diagonal matrices and a final rovibrational matrix for J=5 by the high-disc/low-memory iterative algorithm (NCYCLE = 10), setting up the matrices in a half-matrix of order LANMAX = 3000. Test programs LAN.2 (J=0) and LAN.3 (J=5) both demonstrate how to regulate the amount of memory required by setting up the hamiltonian matrix in stages. Both test programs use a half-matrix of insufficient size (LANMAX = 100) to hold the complete hamiltonian half-matrix.


Numerical Integration

The entire concept embraced by MULTIMODE revolves around the integration of the potential and coriolis terms of the hamiltonian. Furthermore, due to the fact that the potential can be supplied in any arbitrary coordinates, it is not possible to establish any form of analytical integration that can be used for different molecular species; each molecule must be treated as a 'one-off'. Because of this, there is only one acceptable method of integration, and that is numerical quadrature, and because the MULTIMODE method is centred on the Watson hamiltonian, Gauss-hermite numerical quadrature is the automatic choice for such a quadrature scheme. All algorithms commence with Gauss quadrature integration, and the relevant parameters to MULTIMODE are input as MBF and NBF. Here, MBF is the number of Gauss quadrature points required, and NBF is the maximum number of quanta required in the primitive basis function associated with MBF. It is important to note that, in order to obey Gauss quadrature rules, setting MBF > NBF will guarantee exact integration of the kinetic energy integrals. However, in order to obtain exact integration over potential terms will depend, to a large extent, on the form of the potential, and whether or not it can be expanded as a simple polynomial in the Normal Coordinates Q. It is quite usual for the optimum value of MBF to be determined by first setting it equal to some number > NBF, proportional to the anticipated order of the potential in Q, and then comparing the results of (say) a VSCF calculation with those obtained using a second value of MBF, for numerical stability. In this way, an optimum value of MBF can be found that leads to near-exact integration.

Since MULTIMODE is designed for use with large molecular systems, the time and disc space required to integrate all of the various terms in the complete hamiltonian matrix can become excessive. Because of this, an initial one-dimensional contraction of the basis takes place, with the number of functions contracted from NBF to NVF. Depending on the algorithm required, NVF is either supplied on input (for VSCF calculations if ISCFCI=0, or for SCFCI calculations if ISCFCI > 0 and ICI > 0) or it is determined from the parameters ICI, or MAXBAS (if available) for VCI calculations. Much of this selection is transparent to the user, but a value of NVF 'must' be present in the input file, and it is a good idea to use a value which takes into account the following. After the one-dimensional contraction, the integration scheme switches from Gauss to optimised HEG quadrature. The number of points required by this scheme is taken to be NVF + (MBF-NBF); in other words, the 'same' difference between points and functions is maintained as in the Gauss integration scheme. The HEG points are obtained by diagonalising the Q-expectation matrix of the primitives, and so must involve NBF = NVF + (MBF-NBF) such primitives. For VSCF or SCFCI calculations, this algorithm must be taken into account when selecting NVF for input. For VCI calculations, if NMAX < 0 there will be an array MAXBAS supplied on input which contains the maximum quanta required in each mode for construction of the CI basis, and NVF will automatically be taken from this array. In other cases, -ICI will be the maximum quanta in 'all' modes, and this will be used for NVF. If, after one or other of the above selection criteria, the condition imposed above relating NBF to NVF + (MBF-NBF) can not be met, then MULTIMODE will re-define whichever of these parameters it feels necessary to be able to continue. The user must inspect the output carefully to ensure that MULTIMODE is giving him what he really wants.


Input Parameters

Table of Contents

Back to Top


Definitions of input parameters on fort.1



C**TITLE: Any suitable title...maximum 80 characters



C**NATOM: Number of atoms

C**NSTAT: Number of SCF states required

C** NSTAT +ve: Input NSTAT specific SCF states depending on NDEF (below)

C** NSTAT -ve: Program generates -NSTAT SCF states in increasing energy

C**NPOT: Number of user-supplied Normal coordinate force field terms

C** IPOT,JPOT,CPOT (below). If INORM (below) = 1 set NPOT=0

C** If IWHICH (below) = 1 set NPOT=NMODE (below)

C**CONV: Convergence threshold for SCF states in cm-1

C**ICOUPL: Number of modes coupled (1,2,3,4)

C** If POSITIVE, use REAL*8 grid

C** If NEGATIVE, use REAL*4 grid (half the disk storage as *8)

C**ISCFCI: Number of CI energies required (if positive) for SCFCI or VCI

C** (see ICI below)

C** ISCFCI +ve: SCF + CI calculation

C** ISCFCI=0 : SCF only

C** ISCFCI -ve: Terminates after initial coordinate checks (see MCHECK,

C** INORM below)

C**IWHICH: Type of input potential

C** IWHICH=0: Conventional Normal coordinate force field...Normal mode

C** analysis and force field are assumed to be supplied by user

C** IWHICH=1: Internal coordinate potential

C**IDISC: IDISC=0: Write potential and coriolis grid values to disc

C** IDISC=1: Assumes grid values already on disc

C** Rigid-rotor energies written to (61),(62),(63),(64) for

C** ICOUPL=1,2,3,4

C** Potential written to units (71),(72),(73),(74) for ICOUPL=1,2,3,4

C** Coriolis (vibration) written to units (81),(82),(83),(84) for

C** ICOUPL=1,2,3,4

C** Coriolis (rotation) written to units (91),(92),(93),(94) for

C** ICOUPL=1,2,3,4

C**NROTTR: Number of rotational and translational degrees of freedom to include

C** in 'vibrational' modes (number of modes = NMODE=3*NATOM-6+NROTTR)

C**JMAX: Total angular momentum quantum number

C** If POSITIVE, perform exact Whitehead-Handy if ICI < 0 (below)

C** If NEGATIVE, perform adiabatic rotation if ICI < 0 (below)

C** For VSCF only (ISCFCI = 0) or SCFCI (ICI > 0) perform adiabatic

C** rotation

C**MCHECK: MCHECK=0: Use input coordinate system as principal axis system

C** MCHECK=1: Transform equilibrium geometry to principal axis system,

C** and calculate rotational constants. Program continues with geometry

C** in principal axis system (terminate if ISCFCI < 0)

C** MCHECK=-1: Finds minimum of potential (terminate if ISCFCI < 0)

C**INORM: INORM=0: Do not perform Normal mode analysis

C** INORM=1: Perform Normal mode analysis with internal coordinate

C** potential (IWHICH=1) (terminate if ISCFCI < 0)

C** Important note: Normal modes are ordered in increasing energy.



C**ICI: CI basis definition

C** ICI +ve: Number of specific SCF states in SCFCI calculation

C** ICI -ve: -ICI quanta in each mode in VCI calculation (if NMAX.GE.0):

C**NMAX: NMAX > 0: Maximum sum of quanta in VCI basis

C** NMAX = 0: Unrestricted

C** NMAX < 0: Maximum number coupled modes in basis ... maximum -4

C**CUT: Cutoff energy (cm-1) for CI printout

C**EVLJ0: Reference energy (cm-1) for J>0 CI calculations. In general, this

C** will be the zero point energy for J=0, which is printed at the end

C** of the output for J=0. This value will depend on the value of the

C** user's potential at the minimum. The user may wish to scale his

C** potential such that the value at the minimum is zero, in which case

C** the value printed for J=0 is the true zero point energy in cm-1

C**NVALR: Number of rovibrational levels required

C**KSTEP: Increment between successive K (0 to 2*J+1)

C** For KSTEP = 0, only one K-diagonal step is carried out, for Ka = 0

C**IPRINT: Level of output

C** IPRINT=3: Print everything

C** IPRINT=2: Omit rotation matrix element output if J>0

C** IPRINT=1: Omit SCF output, and K-diagonal output if J>0

C** IPRINT=0: Omit integration grid and normal coordinate output

C** IPRINT.GE.0 Print ONE (largest) CI coefficient

C** IPRINT < 0 Print -IPRINT largest CI coefficients

C**MATSIZ: Generate VCI matrix size

C** MATSIZ=0: Normal run

C** MATSIZ=1: Calculate current size of VCI matrix and terminate







C**MAXSUM: Maximum sum quanta in J-mode basis ... restricted only

C** Input single record of -NMAX integers (maximum 4)


C** MAXSUM(1) is maximum sum quanta in one-mode basis sets

C** MAXSUM(2) is maximum sum quanta in two-mode basis sets

C** MAXSUM(3) is maximum sum quanta in three-mode basis sets

C** MAXSUM(4) is maximum sum quanta in four-mode basis sets







C**MAXBAS: Defines selective quanta in VCI basis

C** Input -NMAX records (maximum 4) of NMODE integers


C** MAXBAS(MODE,1) is maximum quanta for mode MODE in one-mode basis sets

C** MAXBAS(MODE,2) is maximum quanta for mode MODE in two-mode basis sets

C** MAXBAS(MODE,3) is maximum quanta for mode MODE in three-mode basis sets

C** MAXBAS(MODE,4) is maximum quanta for mode MODE in four-mode basis sets



C**NCYCLE: Number of Lanczos cycles (Lanczos ignored if NCYCLE=0)

C** +ve: use LANCZA (less memory usage - more disk)

C** -ve: use LANCZB (more memory usage - less disk)

C**TOLLAN: Tolerance for Lanczos eigenvalues

C**LANMAX: Maximum order of half-matrix used to build Lanczos matrix

C**LGIV: 0: Use QL algorithm in Lanczos

C** 1: Use GIVENS in Lanczos



C**NVSYM: Number of vibrational symmetry species (1 or 2 or 4)

C**NWSYM: Number of mode symmetry species (1 or 2 or 3 or 4)

C** For no symmetry NVSYM=NWSYM=1 (C1)...can be used for molecules of

C** any symmetry

C** For four-atom or larger Cs molecules NVSYM=NWSYM=2 (A',A")

C** For triatomic C2v molecules NVSYM=NWSYM=2 (A1,B2)

C** For four-atom C2v molecules NVSYM=4 (A1,B2,B1,A2),

C** NWSYM=3 (A1,B2,B1)

C** For larger molecules, where reduction of symmetry to C2v is possible

C** NVSYM=4 (A1,B2,B1,A2), NWSYM=4 (A1,B2,B1,A2)



C** Input NWSYM records: NSYM ISYM(I), I=1,NSYM


C** For triatomic C2v molecules, NWSYM=1 corresponds to A1,

C** NWSYM=2 corresponds to B2

C** For four-atom C2V molecules, NWSYM=1 corresponds to A1,

C** NWSYM=2 corresponds to B2, NWSYM=3 corresponds to B1

C** etc.

C** The user must know the symmetry of the Normal modes which are

C** ordered either in increasing energy (see INORM=1),

C** or according to the user's convention (see INORM=0)


C**NSYM: Number of modes of current (NWSYM)symmetry

C**ISYM: NSYM specific modes of current (NWSYM) symmetry


C** For example, in the case of H2O, there are 2 mode vibrations

C** of A1 symmetry (modes 1 and 2, say) corresponding to NWSYM=1.

C** There is a single mode vibration of B2 symmetry (mode 3, say)

C** corresponding to NWSYM=2.

C** Hence for NWSYM=1 (totally symmetric), the input required is:

C** NSYM=2, ISYM(I)=1 2

C** For I=2, the input required is:

C** NSYM=1, ISYM(I)=3


C**MVSYM,MWSYM(I), I=1,MVSYM (order of MWSYM: A1,B2,B1,A2)

C** Input a single record of MVSYM+1 integers:

C**MVSYM: Single integer corresponding to the number of actual vibrational

C** (JMAX=0) or rovibrational (JMAX>0) symmetry species required

C**MWSYM: MVSYM integers corresponding to the Specific vibrational (JMAX=0)

C** or rovibrational (JMAX>0) symmetry species required


C** For J = 0 get vibrational energies of specified symmetry.

C** For J > 0 all vibrational symmetries are included in K-blocks, but

C** final rovibrational calculations are done for symmetries specified.


C**Note: For a triatomic molecule, the number of possible rovibrational

C** species is twice the number of possible vibrational species


C** (Currently not yet implemented for J>0)




C** Number of VCI functions written to disk

C** Input only if VCI (ICI<0, ISCFCI>0)

C** IDUMP>0: IDUMP consecutive functions for each symmetry

C** IDUMP<0: NDUMP(IABS(IDUMP)) specific functions for each symmetry


C** (Currently internal use only........set IDUMP=0)




C** Input only if IDUMP<0

C** Input NVSYM records of -IDUMP integers

C**NDUMP: -IDUMP specific functions for symmetry species NVSYM to dump to disc


C** (Currently internal use only........set IDUMP=0)




C** Input NMODE records of three integers

C**NBF(K): Number of quanta in harmonic-oscillator basis for mode K

C**MBF(K): Number of Gauss Hermite integration points for mode K

C**NVF(K): Number of quanta of contracted numerical functions for mode K



C** Normal coordinate force field, defined by:

C** Jelski et al in J. Comput. Chem., 17, 1645 (1996).

C** A maximum of 6 coupled modes is allowed.

C**IPOT: (first 6 integers): defines the Normal mode(s) involved.

C**JPOT: (second 6 integers): defines the number of quanta in the

C** corresponding mode(s)

C**CPOT: defines the potential contribution F(ijkl..) in hartrees


C**These parameters are input under the following conditions:

C**IWHICH=0: Input NPOT terms IPOT,JPOT,CPOT where the first NMODE terms are the

C** NMODE harmonic force constants F(ii)=0.5*omega(i)**2

C**IWHICH=1: Input NMODE terms corresponding to the NMODE F(ii) if INORM=0

C** Input NOTHING if INORM=1

C** note F(ii) is omega(i) in cm-1



C** SCF state definition



C** NSTAT +ve and NDEF -ve:

C** Input NSTAT records ISTAT(NSTAT,NMODE) corresponding to NSTAT

C** specific SCF state definitions (see Jelski et al)


C** NSTAT +ve and NDEF=0 or NDEF +ve:

C** Input NSTAT records N1 N2 where N1 are the number of quanta for

C** mode N2

C** All remaining modes will take on NDEF quanta (see Jelski et al)


C** NSTAT -ve:

C** Input NOTHING. In this case NDEF is a dummy parameter.

C** It is recommended that NSTAT is set -ve for all VCI calculations

C** (see ICI -ve above)



C**XM: Nuclear masses (U)


C**X0(ATOM,I), I=1,3

C**X0: Input NATOM records of equilibrium x,y,z coordinates





C** Mass-weighted Normal Coordinate displacement

C**XL: Input NMODE groups of NATOM records of x,y,z mass-weighted Normal

C** Coordinate displacements (INORM=0)



At this point, control passes to the user to input any data he requires in the

routine USERIN (see Manual)



Table of Contents

Back to Top




 Applications to date


HCO(1), HO2 (2),

H2CN and H2CS(4)

HOCO(6,7), CH4(8)

Transition states




CO-Cu (9 modes)(1)

(CO)2-Cu (18 modes)(3)

Van der Waals systems

Water dimer (12 modes) (unpublished)


Table of Contents

Back to Top




Please contact Stuart Carter directly. His email address is stuartcrtr@gmail.com


Emory University | Department of Chemistry | Bowman Group | Emerson Center