
Stuart Carter and Joel Bowman
Introduction:
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 nonlinear molecules. (A version for linear molecules will be available soon.) The treatment is exact for zero and nonzero total angular momentum, J, and uses general potentials. For nonzero 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 manymode 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.
We thank the National Science Foundation for Financial Support and S. Carter also thanks the Office of Naval Research for support.
Brief description of the theory
1.1 Nonrotating 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 nonrotating, nonlinear molecules:
which we assume is familiar to the reader. We represent the potential V by the following expansion:
where the onemode representation of the potential contains only terms, the twomode 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 usersupplied potential (given as a function of Cartesian coordinates) via the exact transformation between normal coordinates and Cartesian coordinates.
1.1.1 VSCF Theory
The manymode vibrational wavefunction is written as a "Hartree product"
Applying the variational principle to the modals yields the VSCF equations
where V is the nmode representation of the potential (n=1,2,3,4), T_{c} is the Coriolis coupling operator, and T_{i} is the Cartesian kinetic energy operator in the Watson Hamiltonian above. These equations are solved by a standard iteration method beginning with the zeroorder normalmode, harmonic oscillator approximation. Convergence is based on convergence of the eigenvalues. The approximate energy of the state is simply the expectation value of H^{J=0}.
The VSCF modals are expanded in a primitive basis of normalmode harmonicoscillator functions. Matrix elements over the (effective) potential are done by GaussHermite 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 Hmatrix is done using GIVENS/HOUSEHOLDER routines since the number of eigenvalues requested is usually much less than the order of the CI Hmatrix. 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, usersupplied potential, and stationary geometry. The user may supply his/her own normalmode 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 manymode 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 manymode 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.
The Hamiltonian for nonzero 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 Nmode 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 nonrotating 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 H^{AR} are obtained using the VSCF, VSCF+CI, or VCI methods.
(1). "Vibrational SelfConsistent Field Method for ManyMode 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 manymode 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 manymode molecules", S. Carter, J. M. Bowman, and N. Handy, Theoretical Chem. Accnts. 100, 191 (1998).
(5) Nonseparable transition state theory for nonzero total angular momentum: Implications for Jshifting and application to the OH+H2 reaction, J. M. Bowman and H. Shnider, J. Chem. Phys. 110, xxxx (1999).
(6) "Calculations of lowlying vibrational states of cis and transHOCO", J. M. Bowman, K. Christoffel, and G. Weinberg, Theochem, in press.
(7) "Nonseparable RRKM theory applied to unimolecular dissociation of HCO2 H+CO2, J. Phys. Chem., K. Christoffel and J. M. Bowman, submitted.
(8) "Variational calculations of rovibrational 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.
UserSpecific Routines
The body of the code is `multimode.3.4.f'; userspecific codes are `memo.f' and `user.f'. A brief summary of the latter two codes is as follows:
memo.f
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
PARAMETER (MAXSIZ=XXXX)
If the current RAM specified by MAXSIZ is exceeded during a run of `multimode', the program stops with the message:
`MEMORY EXHAUSTED'
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 halfmatrix
Unit 21: Basic rovibrational integrals for onemode coupling
Unit 22: Basic rovibrational integrals for twomode coupling
Unit 23: Basic rovibrational integrals for threemode coupling
Unit 24: Basic rovibrational integrals for fourmode 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:
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: Rigidrotor energies at each onemode grid point
Unit 62: Rigidrotor energies at each twomode grid point
Unit 63: Rigidrotor energies at each threemode grid point
Unit 64: Rigidrotor energies at each fourmode grid point
Unit 71: Potential energy at each onemode grid point
Unit 72: Potential energy at each twomode grid point
Unit 73: Potential energy at each threemode grid point
Unit 74: Potential energy at each fourmode grid point
Unit 81: Vibrational coriolis coupling data at each onemode grid point
Unit 82: Vibrational coriolis coupling data at each twomode grid point
Unit 83: Vibrational coriolis coupling data at each threemode grid point
Unit 84: Vibrational coriolis coupling data at each fourmode grid point
Unit 91: Rotational coriolis coupling data at each onemode grid point
Unit 92: Rotational coriolis coupling data at each twomode grid point
Unit 93: Rotational coriolis coupling data at each threemode grid point
Unit 94: Rotational coriolis coupling data at each fourmode 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 onemode to fourmode coupling.
user.f
The module user.f contains two routines which may be modified by the user in order to interface with his/her potential energy routines.
1. USERIN
2. GETPOT
The calling sequences to these routines are:
CALL USERIN
CALL GETPOT(V,NATOM,XX,RR)
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.
DIMENSION XX(NATOM,3),RR(NATOM,NATOM)
XX contains the cartesian coordinates of the NATOM atoms, where XX(NATOM,1), XX(NATOM,2), XX(NATOM,3) are the x, y, zcoordinates 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.
CALL BONDS(NATOM,RR,XX)
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
COMMON/FILASS/IOUT,INP
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 userinitiated 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
COMMON/FUNDAM/WAVENM,ATTOJ,BOHR,ELMASS,RAD
contains some useful fundamental constants:
WAVENM: Conversion from hartrees to wavenumbersATTOJ: 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 machinespecific 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).
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 setup 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 (nonorthogonal) VSCF wavefunctions as a CI basis.
4..VCI: This category demonstrates the use of (orthonormal) virtual zpe
wavefunctions in nearvariational CI calculations.
5..LAN: This category demonstrates the use of an iterative
diagonalisation technique for use with very large VCI hamiltonian
matrices.
*************************************************************
Introduction
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 'UserSpecific 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 wellbehaved, 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 'UserSpecific 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 4atom (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 Setup 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 nonlinear 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 cm1) 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 Gausshermite 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 onedimensional 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 (nonequilibrium) geometry as input, and inspection of the output reveals that, after fitting to the maximumallowed 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 totallysymmetric 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 SelfConsistent 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.D2 cm1. 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 usersupplied 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 NSTAT1 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. Nonorthogonal 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 Kdiagonal (JMAX > 0) energies required, subject to the print cutoff parameter in CUT. All CI energies are printed in cm1 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 cutoff of 1.0D+4 cm1, 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 nonorthogonal 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 symmetryblocked, 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 symmetryblocked. 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 onemode, twomode, threemode terms in the basis. The second (MAXBAS) is the maximum quantum which appears (for all NMODE virtual functions) for onemode, twomode, threemode 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 Kdiagonal '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 Kdiagonal '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 Kdiagonal 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 nearvariational 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 incore methods may well be generated. To allow such matrices to exist in the MULTIMODE method, a Lanczos/Davidsontype 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 halfmatrix of order LANMAX is used to set up the complete hamiltonian halfmatrix before writing it to disc. During the course of the diagonalisation, ISCFCI eigenvalues (for vibrational or Kdiagonal 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 incore 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 onedimensional arrays and requires extra disc space to store intermediate arrays; the second (NCYCLE < 0) involves twodimensional arrays to hold these internal arrays, and therefore requires more memory but less disc space.
Test program LAN.1 diagonalises alternate Kdiagonal matrices and a final rovibrational matrix for J=5 by the highdisc/lowmemory iterative algorithm (NCYCLE = 10), setting up the matrices in a halfmatrix 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 halfmatrix of insufficient size (LANMAX = 100) to hold the complete hamiltonian halfmatrix.
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 'oneoff'. 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, Gausshermite 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 nearexact 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 onedimensional 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 onedimensional contraction, the integration scheme switches from Gauss to optimised HEG quadrature. The number of points required by this scheme is taken to be NVF + (MBFNBF); 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 Qexpectation matrix of the primitives, and so must involve NBF = NVF + (MBFNBF) 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 + (MBFNBF) can not be met, then MULTIMODE will redefine 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.
Definitions of input parameters on fort.1
C****************************************************************************
C**TITLE
C**TITLE: Any suitable title...maximum 80 characters
C****************************************************************************
C**NATOM,NSTAT,NPOT,CONV,JCOUPL,ISCFCI,IWHICH,IDISC,NROTTR,KMAX,MCHECK,INORM
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 usersupplied 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 cm1
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** Rigidrotor 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*NATOM6+NROTTR)
C**JMAX: Total angular momentum quantum number
C** If POSITIVE, perform exact WhiteheadHandy 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****************************************************************************
C**ICI,NMAX,CUT,EVLJ0,NVALR,KSTEP,IPRINT,MATSIZ
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 (cm1) for CI printout
C**EVLJ0: Reference energy (cm1) 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 cm1
C**NVALR: Number of rovibrational levels required
C**KSTEP: Increment between successive K (0 to 2*J+1)
C** For KSTEP = 0, only one Kdiagonal 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 Kdiagonal 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****************************************************************************
C**OMIT FOLLOWING INPUT IF SCF ONLY (ISCFCI.LE.0)
C**OMIT FOLLOWING INPUT IF SCFCI (ICI > 0)
C**OMIT FOLLOWING INPUT IF UNRESTRICTED OR RESTRICTED VCI (NMAX > 0)
C**
C**MAXSUM(I), I = 1,NMAX
C**MAXSUM: Maximum sum quanta in Jmode basis ... restricted only
C** Input single record of NMAX integers (maximum 4)
C**
C** MAXSUM(1) is maximum sum quanta in onemode basis sets
C** MAXSUM(2) is maximum sum quanta in twomode basis sets
C** MAXSUM(3) is maximum sum quanta in threemode basis sets
C** MAXSUM(4) is maximum sum quanta in fourmode basis sets
C****************************************************************************
C**OMIT FOLLOWING INPUT IF SCF ONLY (ISCFCI.LE.0)
C**OMIT FOLLOWING INPUT IF SCFCI (ICI > 0)
C**OMIT FOLLOWING INPUT IF UNRESTRICTED OR RESTRICTED VCI (NMAX > 0)
C**
C**MAXBAS(MODE,NMAX), MODE=1,NMODE
C**MAXBAS: Defines selective quanta in VCI basis
C** Input NMAX records (maximum 4) of NMODE integers
C**
C** MAXBAS(MODE,1) is maximum quanta for mode MODE in onemode basis sets
C** MAXBAS(MODE,2) is maximum quanta for mode MODE in twomode basis sets
C** MAXBAS(MODE,3) is maximum quanta for mode MODE in threemode basis sets
C** MAXBAS(MODE,4) is maximum quanta for mode MODE in fourmode basis sets
C****************************************************************************
C**NCYCLE,TOLLAN,LANMAX,IGIV
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 halfmatrix used to build Lanczos matrix
C**LGIV: 0: Use QL algorithm in Lanczos
C** 1: Use GIVENS in Lanczos
C****************************************************************************
C**NVSYM,NWSYM
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 fouratom or larger Cs molecules NVSYM=NWSYM=2 (A',A")
C** For triatomic C2v molecules NVSYM=NWSYM=2 (A1,B2)
C** For fouratom 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****************************************************************************
C**NSYM(NWSYM),ISYM(NWSYM,NSYM(NWSYM)) (order of NWSYM: A1,B2,B1,A2)
C** Input NWSYM records: NSYM ISYM(I), I=1,NSYM
C**
C** For triatomic C2v molecules, NWSYM=1 corresponds to A1,
C** NWSYM=2 corresponds to B2
C** For fouratom 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**
C**NSYM: Number of modes of current (NWSYM)symmetry
C**ISYM: NSYM specific modes of current (NWSYM) symmetry
C**
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****************************************************************************
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**
C** For J = 0 get vibrational energies of specified symmetry.
C** For J > 0 all vibrational symmetries are included in Kblocks, but
C** final rovibrational calculations are done for symmetries specified.
C**
C**Note: For a triatomic molecule, the number of possible rovibrational
C** species is twice the number of possible vibrational species
C********************************************************************
C** (Currently not yet implemented for J>0)
C********************************************************************
C****************************************************************************
C**IDUMP
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********************************************************************
C** (Currently internal use only........set IDUMP=0)
C********************************************************************
C****************************************************************************
C**NDUMP(I,NVSYM), I=1,IDUMP
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********************************************************************
C** (Currently internal use only........set IDUMP=0)
C********************************************************************
C****************************************************************************
C**NBF(NMODE),MBF(NMODE),NVF(NMODE)
C** Input NMODE records of three integers
C**NBF(K): Number of quanta in harmonicoscillator 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****************************************************************************
C**IPOT(NPOT),JPOT(NPOT),CPOT(NPOT)
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**
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 cm1
C****************************************************************************
C**NDEF
C** SCF state definition
C**
C**
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**
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**
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****************************************************************************
C**XM(I), I=1,NATOM
C**XM: Nuclear masses (U)
C****************************************************************************
C**X0(ATOM,I), I=1,3
C**X0: Input NATOM records of equilibrium x,y,z coordinates
C****************************************************************************
C**OMIT FOLLOWING INPUT IF INORM = 1
C**
C**XL(NATOM,NMODE,I), I=1,3
C** Massweighted Normal Coordinate displacement
C**XL: Input NMODE groups of NATOM records of x,y,z massweighted Normal
C** Coordinate displacements (INORM=0)
C****************************************************************************
C**USERDEFINED POTENTIAL
At this point, control passes to the user to input any data he requires in the
routine USERIN (see Manual)
Molecules
HCO(1), HO2 (2),H2CN and H2CS(4)
HOCO(6,7), CH4(8)
Transition states
H2OH(5)HCO2(7)
Adorbates
COCu (9 modes)(1)(CO)2Cu (18 modes)(3)
Van der Waals systems
Water dimer (12 modes) (unpublished)
HOW TO ORDER MULTIMODE
Please contact Stuart Carter directly. His email address is S.Carter@reading.ac.uk