This is an old revision of the document!

# Vibration correlation programs

VCI,options [vci]

VCI calculations account for vibration correlation effects and use potential energy surfaces as generated from the SURF program and a basis of VSCF modals. For each vibrational state an individual VCI calculation will be performed. As VCI calculations may require substantial computer resources, these calculations can be rather expensive. Currently, two different VCI programs (configuration selective and conventional) are available (see below). Moreover, VCI calculations can be performed using the grid-based version of the program or within an analytical representation. The latter is significantly faster and is thus recommended. The different versions of the configuration selection VCI program and the underlying configuration selection scheme are described in detail in:
M. Neff, G. Rauhut, Toward large scale vibrational configuration interaction calculations, J. Chem. Phys. 131, 124129 (2009).
M. Neff, T. Hrenar, D. Oschetzki, G. Rauhut, Convergence of vibrational angular momentum terms within the Watson Hamiltonian, J. Chem. Phys. 134, 064105 (2011).
The anharmonic frequencies and intensities calculated by the VCI program can be used to plot an IR spectrum, using the PUT command (see subsection writing files for postprocessing(PUT)) with the style IRSPEC.

The following options are available:

VCI solutions can be obtained using a potential in grid representation, i.e. POT=GRID (default), or in an analytical representation, POT=POLY, POT=GAUSS or POT=BSPLINE. In the latter cases the POLY program needs to be called prior to the VSCF and VCI programs in order to transform the potential.

By default, i.e. SADDLE=0, the VCI program assumes, that the reference point of the potential belongs to a local minimum. Once the PES calculation has been started from a transition state, this information must be provided to the VCI program by using SADDLE=1. Currently, the VCI program can only handle symmetrical double-minimum potentials.

Both, the grid-based and the analytical versions of the VCI programs offer 2 different kinds of VCI implementations:
VERSION=3 (which is the default) is a configuration selective and most efficient VCI program. VERSION=4 is a conventional VCI program without configuration selection. It is thus computationally extremely demanding.

CITYPE defines the maximum number of simultaneous excitations, i.e. Singles, Doubles, Triples, … and thus determines the kind of calculations, i.e. VCISD, VCISDT, … The default is CITYPE=4 (VCISDTQ) for 3D potentials and CITYPE=5 for 4D potentials, which appears to be a fair compromise between accuracy and computational speed. The maximum excitation level is currently limited to CITYPE=9.

LEVEX determines the level of excitation within one mode, i.e. $0\rightarrow 1$, $0\rightarrow 2$, $0\rightarrow 3$, … The default is LEVEX=5, which was found to be sufficient for many applications.

CIMAX is the maximum excitation level corresponding to CITYPE and LEVEX. In principle, a triple configuration $(1^42^43^4)$ would contribute to the VCI space. However, CIMAX=7 restricts this to $(1^42^3)$, $(1^32^33^1)$, $(1^32^23^2), \ldots$. The default is CIMAX=12 for 3D potentials and CIMAX=15 for 4D potentials.

THRSEL controls the determination of the iterative configuration selection scheme. By default the wavefunction is considered to be converged once energy changes drop below THRSEL=0.02 cm$^{-1}$.

THRCF is the threshold for selecting individual configurations. The default is given by THRCF=$5\cdot 10^{-10}$.

This thresholds controls the exclusion of selected configurations within the perturbative configuration selection criterion. The default is THREX=5.d-4.

Within the perturbative selection criterion, unselected configurations are excluded from the procedure based on a criterion, which checks on the selection of configuration, which are excited in the same modes. The default is THRUP=$10^{-9}$.

By default SELSCHEME=1, configurationis will be selected by a perturbative criterion. Alternative one may use a criterion based on 2$\times$2 VCI matrices (SELSCHEME=2). Usually the differences are extremely small and the matrix based criterion is slightly more time-consuming.

Within the evaluation of the VCI integrals contraction schemes are used to reduce the computational effort. In the analytical VCI program values 0, 1 and 2 are supported, while in the grid based version only the options 0 and 1 exist. Memory demands and CPU speed-ups increase with increasing values. The default is set to 1. On machines with limited memory a value of 0 is recommended for this keyword.

VAM=0: switches off all vibrational angular momentum terms and the Watson correction term.
VAM=1: adds the Watson correction term (see eq. (1) in vibrational SCF programs) as a pseudo-potential like contribution to the fine grid of the potential. In the analytical representation of the potential this will already be done in the POLY program.
VAM=2: (default) the 0D terms of the vibrational angular momentum terms, i.e. $\frac{1}{2} \sum_{\alpha\beta} \hat{\pi}_\alpha\mu_{\alpha\beta} \hat{\pi}_\beta$, and the Watson correction term are included. The VAM-terms will be added to the diagonal elements of the VCI-matrix only. This approximations works rather well for many applications.
VAM=3: again, the $\mu$ tensor is given as the inverse of the moment of inertia tensor at equilibrium geometry, but is added to all elements of the VCI matrix.
VAM=4: extends the constant $\mu$-tensor (0D) by 1D terms and is added to all elements of the VCI matrix. A prescreening technique is used for the 1D terms, in which the convergence of the VAM operator will be checked for each VCI matrix element.
VAM=5: includes 0D, 1D and 2D terms of the $\mu$-tensor, which are added to all elements of the VCI matrix. Prescreening is used for the 1D and 2D terms.
VAM=6: includes 0D and 1D terms of the $\mu$-tensor without prescreening.
VAM=7: includes 0D, 1D and 2D terms of the $\mu$-tensor, which are added to all elements of the VCI matrix. Prescreening is used for the 2D terms only.
VAM=8: includes 0D, 1D and 2D terms of the $\mu$-tensor without any prescreening.
Note that the 1D and 2D corrections increase the computational cost considerably and are only available for non-linear molecules.

By default the VSCF program calculates the fundamental modes of the molecule only. However, choosing COMBI=$n$ allows for the calculation of the vibrational overtones and combination bands. The value of $n$ controls the excitation level, i.e. the number of states to be computed increases very rapidly for large values of $n$. Therefore, by default the upper limit is set to 5000 cm$^{-1}$, but this cutoff can be changed by the option UBOUND.

Once vibrational states have been defined with the VIBSTATE program (section the VIBSTATE program (VIBSTATE)), the VCI program can be forced to compute just these states by the option USERMODE=1. Note that the vibrational ground state will always be computed and needs not to be specified explicitly.

Once overtones and combination bands shall be computed, the upper energy limit is controlled by the keyword UBOUND, i.e. states, for which the harmonic estimate is larger than $n$, will not be computed. the default is set to $n$=5000 cm$^{-1}$.

Rovibrational energy levels can be computed within the adiabatic rotation approximation (ARA). All energy levels for J-values up to $n$ will be calculated.

THERMO=1 allows for the improved calculation of thermodynamical quantities (compare the THERMO keyword in combination with a harmonic frequency calculation). However, the approach used here is an approximation: While the harmonic approximation is still retained in the equation for the partition functions, the actual values of the frequencies entering into these functions are the anharmonic values derived from the VCI calculation. Default: THERMO=0.

DIPOLE=1 (default) allows for the calculation of infrared intensities. Calculation of infrared intensities requires the calculation of dipole surfaces within the SURF program. By default the intensities will be computed on the basis of Hartree-Fock dipole surfaces.

POLAR=1 allows to compute Raman intensities in addition to infrared intensities, but of course requires polarizability tensor surfaces from the SURF program. By default Raman intensities are switched off.

The expansion of the potential in the VCI calculation can differ from the expansion in the SURF calculation. However, only values less or equal to the one used in the surface calculation can be used. Default: NDIM=3.

Term after which the $n$-body expansions of the dipole surfaces are truncated. The default is set to 3. Note that NDIMDIP has to be lower or equal to NDIM.

Term after which the $n$-body expansions of the polarizability tensor surfaces are truncated. The default is set to 0. Note that NDIMPOL has to be lower or equal to NDIM and must be samller than 4.

By default the symmetry of the molecule will be recognized automatically within the VCI calculations. MPG=1 switches symmetry off.

This keyword specifies the reference for the definition of the configurations. By default, REFERENCE=0 the reference for all state-specific calculations is the vibrational ground-state configuration. This leads to a violation of the Brillouin condition, but often to also to faster convergence. REFERENCE=1 uses the VSCF configuration as reference for generating all excited configurations. This is the proper way of doing it, but usually requests higher excitation levels.

By default all VCI calculations will employ ground-state based VSCF modals, GSMODALS=0. GSMODALS=0 uses the state-specific VSCF modals the subsequent VCI calculations. In any case, individual VCI calculations will be performed for each vibrational state (in contrast to just one VCI calculations from which all solutions will be retrieved) and thus the final VCI wave functions may not be strictly orthogonal to each other once the VCI space is incomplete.

In the analytical configuration selective VCI program different diagonalization schemes can be used. DIAG=TRS uses our residual-based algorithm for the calculation of eigenpairs (RACE), which is the default. DIAG=CON specifies a conventional non-iterative diagonalization as used in the grid-based versions, which will be used once the ANALYZE keyword has been provided. DIAG=JAC uses a Jacobi-Davidson scheme. DIAG=HJD denotes a disk-based Jacobi-Davidson algorithm.

In case of resonances or strongly mixed states in general (i.e. low VCI coefficients) a multi-state analysis can be performed, which prints major contributions of the VCI-vectors for all states in a certain window around the state of interest. Typically a window between 10 and 20% (i.e. ANALYZE=0.1 or ANALYZE=0.2) provides all the information needed. As this analysis requires a conventional diagonalization (see DIAG), the CPU time may increase significantly.

By default the expansion of the $\mu$-tensor for calculating the vibrationally averaged rotational constants is truncated after the 2nd order terms, i.e. NVARC=2. This may be altered by the NVARC keyword.

This option provides an extended output. PRINT=1 prints the vibrationally averaged rotational constants for all computed states and the associated vibration-rotation constants $\alpha$. PRINT=2 prints the effective 1D polynomials in case that the potential is represented in terms of polynomials, see the option TYPE=POLY and the POLY program. In addition the generalized VSCF property integrals, i.e. $\left < VSCF \left | q_i^r \right | VSCF \right >$ are printed. These integrals allow for the calculation of arbitrary vibrationally averaged properties once the property surfaces are available. Default: PRINT=0.

If variable is set to FCON, important VCI information will be passed to the Franck-Condon calculation. Within Franck-Condon calculations this option has to be used.

INFO=1 provides a list of the values of all relevant program parameters (options).

LEVEX,options
Within the VCI program the correlation space can be specified in a general manner (keyword LEVEX), which means that the same number of modals for each mode will be used. Alternatively, one may use the LEVEX directive. This allows to specify the correlation spaces for the individual modes.

• MODE(n)=m The number of correlating modals for mode n is set to m.

The following input example for a grid based calculation of anharmonic frequencies and intensities (1) optimizes the geometry of water, (2) computes the harmonic frequencies, (3) generates a potential energy surface around the equilibrium structure, (4) computes the vibrational wave function and the infrared intensities at the VSCF level, and finally (5) a VCI calculation will be performed. Vibrational angular momentum terms (VAM) are included even for the non-diagonal elements of the VCI matrix.

memory,20,m
basis=vdz
orient,mass
geometry={
3
Water
O          0.0675762564        0.0000000000       -1.3259214590
H         -0.4362118830       -0.7612267436       -1.7014971211
H         -0.4362118830        0.7612267436       -1.7014971211
}

mass,iso

hf
mp2
optg                                     !(1) optimizes the geometry
frequencies,symm=auto                   !(2) compute harmonic frequencies

label1
int
{hf
start,atden}
{mp2
cphf,1}

{surf,start1D=label1,sym=auto            !(3) generate potential energy surface
intensity,dipole=2}
vscf                                     !(4) do a VSCF calculation
vci,vam=3                                !(5) do a VCI calculation
put,irspec,irspec.gnu                    !writes a gnuplot file to plot an IR
!spectrum of the last VCI calculation

The following input example for aan analytical calculation of anharmonic frequencies and intensities (1) optimizes the geometry of water, (2) computes the harmonic frequencies,(3) generates a potential energy surface around the equilibrium structure, (4) converts the potential energy surface into an analytical representation (5) computes the nuclear wave function and the infrared intensities at the VSCF level, and finally (6) performs a VCI calculation. Vibrational angular momentum terms (VAM) are included even for the non-diagonal elements of the VCI matrix.

memory,20,m
basis=vdz
orient,mass
geometry={
3
Water
O          0.0675762564        0.0000000000       -1.3259214590
H         -0.4362118830       -0.7612267436       -1.7014971211
H         -0.4362118830        0.7612267436       -1.7014971211
}

mass,iso

hf
mp2
optg                                     !(1) optimizes the geometry
frequencies,symm=auto                   !(2) compute harmonic frequencies

label1
int
{hf
start,atden}
{mp2
cphf,1}

{surf,start1D=label1,sym=auto            !(3) generate potential energy surface
intensity,dipole=2}
poly,dipole=1                            !(4) converts potential energy surface
!    to a polynomial representation
vscf,type=poly                           !(5) do a VSCF calculation
vci,type=poly,vam=3                      !(6) do a VCI calculation
put,irspec,irspec.gnu                    !writes a gnuplot file to plot an IR
!spectrum of the last VCI calculation

DISK,options

The DISK directive allows to specify explicitly, from where the potential information shall be taken and where it shall be stored to disk. This can also be accomplished in an automated manner. These features are only relevant for the simulation of vibronic spectra as one has to deal with several PESs in the same input. For simple VCI calculations, no information is needed here.

The following options are available:

• START=record This card specifies the record from where to read the VSCF information. As the VSCF information usually is stored in the same record as the polynomials, it is usually defined in the POLY program.
• SAVE=record This keyword specifies the record where to dump the VCI information.
• AUTO=n Rather than using the options START and SAVE one may simply assign a label n to a a certain PES and all the records will be set automatically.

VMP2,options [vmp2]

The VMP2 program allows to perform 2nd order vibrational Møller-Plesset calculations. The program has been implemented in a grid-based and an analytical version. Most of the keywords as described for the VCI program are also valid for the VMP2 program, i.e. TYPE, CITYPE, LEVEX, CIMAX, NGRID, NDIM, NBAS, VAM, COMBI, THERMO, DIPOLE, MPG and INFO.

examples/vmp2.inp
basis=vdz
orient,mass
geomtyp=xyz
geometry={
3
Water
O          0.0675762564        0.0000000000       -1.3259214590
H         -0.4362118830       -0.7612267436       -1.7014971211
H         -0.4362118830        0.7612267436       -1.7014971211
}

hf
mp2
optg
{frequencies,symm=auto
print,low=50}

label1
{hf
start,atden}
{mp2
cphf,1}

{surf,start1D=label1,sym=auto
intensity,dipole=2}
poly,pmp=1,dipole=1,show=1
vscf,type=poly,pmp=2,dipole=1
vmp2,type=poly,pmp=3,dipole=1

VMRCI,options [vmrci]

The vibrational multi-reference CI program requests a preceding VMCSCF calculation and accounts for correlation effects on top of those considered in the VMCSCF run. Note, in all VMRCI calculations, configurations being considered within the active space of the underlying VMCSCF calculation, will be relaxed. Configurations to be generated for the VMRCI correlation space refer to all reference configurations and thus the VMRCI correlation space increases significantly faster than the VCI space. Once parameters have been altered with respect to the defaults in the VMCSCF program, e.g. NGRID, the same settings must be used within the VMRCI program. Further details are described in:
F. Pfeiffer, G. Rauhut, Multi-reference vibration correlation methods, J. Chem. Phys. 140, 064110 (2014).

The following options are available:

• VSCFMODALS=n Instead of using VMCSCF modals (default), VSCF modals will be used (once a VSCF calculation has been performed prior to the VMCSCF calculation), i.e. VSCFMODALS=1.

Most of the keywords as described for the VCI program are also valid for the VMRCI program, i.e. TYPE, CITYPE, LEVEX, CIMAX, NDIM, NBAS, DIAG, VAM, COMBI, DIPOLE, CONT, ANALYZE, MPG and INFO.

LEVEX,options
Within the VMRCI program the correlation space can be specified in a general manner (keword LEVEX), which means that the same number of modals for each mode will be used. Alternatively, one may use the LEVEX directive. This allows to specify the correlation spaces for the individual modes.

• MODE(n)=m The number of correlating modals for mode n is set to m.