**This is an old revision of the document!**

# Vibrational SCF programs

## The VSCF program (VSCF)

`VSCF`

,*options* [vscf]

The `VSCF`

program is exclusively based on the Watson Hamiltonian
\begin{equation}\label{eq:1}
\hat{H} = \frac{1}{2} \sum_{\alpha\beta} ( \hat{J}_\alpha - \hat{\pi}_\alpha) \mu_{\alpha\beta}
(\hat{J}_\beta - \hat{\pi}_\beta)
-\frac{1}{8}\sum_\alpha \mu_{\alpha\alpha} -\frac{1}{2}\sum_i \frac{\partial^2}{\partial q_i^2} + V(q_1,\dots,q_{3N-6})
\end{equation}
in which the potential energy surfaces, $V(q_1,\dots,q_{3N-6})$, are provided by the `SURF`

module. The Watson correction term and the 0D term of the vibrational angular momentum terms are by default (`VAM=2`

) included. Within the grid-based version of the program the one-dimensional Schrödinger equation is solved by the DVR procedure of Hamilton and Light. Note that, the number of basis functions (e.g. distributed Gaussians) is determined by the grid points of the potential and cannot be increased without changing the PES grid representation. In contrast to that, the number of basis functions can be modified without restrictions in the version based on an analytical representation of the potential (polynomials, B-splines, Gaussians). As VSCF calculations are extremely fast, these calculations cannot be restarted. For details see:

G. Rauhut, T. Hrenar, *A Combined Variational and Perturbational Study on the Vibrational Spectrum of P$_2$F$_4$*, Chem. Phys. **346**, 160 (2008).

The anharmonic frequencies and intensities calculated by the `VSCF`

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:

VSCF solutions can be obtained using a potential in grid representation, i.e.`POT`

=*variable*`POT=GRID`

, or in an analytical representation,`POT=POLY`

,`POT=BSPLINE`

,`POT=GAUSS`

. In the latter case the`POLY`

program needs to be called prior to the`VSCF`

program in order to transform the potential.By default, i.e.`SADDLE`

=*n*`SADDLE=0`

, the`VSCF`

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`VSCF`

program by using`SADDLE=1`

. Currently, the`VSCF`

program can only handle symmetrical double-minimum potentials.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 by default (`VAM`

=*n*`VAM=2`

) included.

`VAM=1`

adds the Watson correction term (see eq. \eqref{eq:1}) as a pseudo-potential like contribution to the fine grid of the potential.

`VAM=2`

allows for the calculation of the integrals of the `VAM`

operator using the approximation that the $\mu$ tensor is given as the inverse of the moment of inertia tensor at equilibrium geometry.

When using `VAM=4`

the expansion of the effective moment of inertia tensor will be truncated after the 1D terms (rather than the 0D term in case of `VAM=2`

). Note that values higher than 2 are only active for non-linear molecules. `VAM=5`

truncates the series after the 2D term. In almost all cases `VAM=2`

is fully sufficient. Vibrational angular momentum terms are accounted for in a perturbational manner and do not affect the wavefunction.

Plots all $\mu$-tensor surfaces up to`MUPLOT`

=*n**n*D and a corresponding Gnuplot script in a separate subdirectory (`plots`

). This option works only in combination with`POT=POLY`

. The`VAM`

option has to be set accordingly.By default the`COMBI`

=*n*`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`

. See also the`VIBSTATE`

program (section the VIBSTATE program (VIBSTATE)) for even more possibilities of defining vibrational states.Once vibrational states have been defined with the`USERMODE`

=*n*`VIBSTATE`

program (section the VIBSTATE program (VIBSTATE)), the VSCF 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`

=*n*`UBOND`

, 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}$.For solving the one-dimensional Schrödinger equation within a grid representation two different algorithms can be used. The default, i.e.`SOLVER`

=*n*`SOLVER=1`

, calls the discrete variable representation (DVR) as proposed by Hamilton and Light. Alternatively, the collocation algorithm of Young and Peet can be used (`SOLVER=2`

).The initial guess for the VSCF programs is by default generated from the uncoupled one-dimensional potentials, i.e.`GUESS`

=*n*`GUESS=1`

. Alternatively, one may start within the calculation of excited vibrational states from the solution of the vibrational ground state,`GUESS=2`

.`BASIS`

=*variable*`BASIS=DGB`

(default) defines a mode-specific basis of distributed Gaussians and distributes the Gaussians in a way, that the overlap integral between two functions is always the same (controlled by`THRBASOVLP`

). This guarantees that an increasing number of basis functions will always lead to an improvement.`BASIS=HO`

defines a harmonic oscillator basis. Using this basis together with`MAXITER=0`

and`VAM=0`

provides a simple harmonic oscillator basis to be used in all subsequent programs, e.g. in the VCI program.`BASIS=SIN`

uses a basis of sine functions. This is not a fully implemented feature, but primarily available for experimental purposes.Overlap between two Gaussian basis functions, once`THRBASOVLP`

=*value*`BASIS=DGB`

has been chosen. The default is 0.75.Determines the type of orthogonalization within the VSCF program.`ORTHO`

=*n*`ORTHO=1`

invokes a symmetrical orthogonalization,`ORTHO=2`

a canonical one and`ORTHO=3`

uses a canonical one together with an elimination of linear dependencies (see also keyword`THRLINDEP`

. The default is`ORTHO=1`

.Threshold for eliminating linear dependencies in the VSCF procedure (see keyword`THRLINDEP`

=*value*`BASIS=DGB`

). The default is`THRLINDEP=1e-8`

.This key sets the maximum number of iterations to be performed in the VSCF program.`MAXITER`

=*n*By default VSCF calculations will be performed for non-rotating molecules, i.e. J=0. Rovibrational transitions can be computed for arbitrary numbers of J$=n$ within the adiabatic rotation approximation.`JMAX`

=*n*`THERMO`

=*n*`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`VSCF`

calculation. Default:`THERMO=0`

.`DIPOLE`

=*n*`DIPOLE=1`

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`

=*n*`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`NDIM`

=*n*`VSCF`

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.Term after which the $n$-body expansions of the dipole surfaces are truncated. The default is set to 3. Note that`NDIMDIP`

=*n*`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`

=*n*`NDIMPOL`

has to be lower or equal to`NDIM`

and must be smaller than 4.The number of basis functions (distributed Gaussians) to be used for solving the VSCF equations can be controlled by`NBAS`

=*n*`NBAS`

=*n*. The default is`NBAS=20`

. This option is only active once an analytical representation of the potential has been chosen, see the option`POT`

and the`POLY`

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

=*n*`NVARC=2`

. This may be altered by the`NVARC`

keyword.This option provides an extended output.`PRINT`

=*n*`PRINT=1`

prints the vibrationally averaged rotational constants for all computed states and the associated vibration-rotation constants $\alpha$. Moreover, the temperature dependence of bond lengths will also be printed, when the potential is represented by a linear combination of basis functions.`PRINT=2`

prints the effective 1D polynomials in case that the potential is represented in terms of polynomials, see the option`POT=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`

.`INFO`

=*n*`INFO=1`

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

The following input example for a grid based calculation of anharmonic frequencies and intensities on the `VSCF`

level (1) optimizes the geometry of water, (2) computes the harmonic frequencies,(3) generates a potential energy surface around the equilibrium structure and (4) computes the nuclear wave function and the infrared intensities at the `VSCF`

level. Vibrational angular momentum terms (`VAM`

) are included. Note, that it is recommended to perform a `VCI`

calculation after a `VSCF`

calculation. The details of the `VCI`

input are described in the next chapter the VCI program (VCI).

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 {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 put,irspec,irspec.gnu !writes a gnuplot file to plot an IR !spectrum of the VSCF calculation

### Record handling

`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:

Surface information shall be read from the specified record. This must correspond to the record, to which the potential has been dumped in the`STARTSURF`

=*record*`SURF`

program.Polynomial and other information shall be read from the specified record. This must be the same record, to which the polynomials have been dumped in the`START`

=*record*`POLY`

program.This specifies the record, where to dump the VSCF information. Usually this is the same record as specified in the`SAVE`

=*record*`START`

option. Note that the VSCF information is currently stored in the same record as the polynomial information.Rather than using the options`AUTO`

=*n*`START`

and`SAVE`

one may simply assign a label*n*to a a certain PES and all the records will be set automatically.

## The VMCSCF program (VMCSCF)

`VMCSCF`

,*options* [vmcscf]

The `VMCSCF`

program is still under development and has not yet been fully optimized with respect to speed and efficiency. It is strongly recommended to run a VSCF calculation prior to the VMCSCF calculation and thus to use VSCF modals as an initial guess. By default configuration-selective VMCSCF calculations will be performed. The VMCSCF calculation based on analytical representations of the potential is significantly faster than the grid-based version and should be used whenever possible. By default the active space will be determined automatically, but it can also be controlled by the input stream. For details see:

S. Heislbetz, G. Rauhut, *Vibrational multiconfiguration self-consistent field theory: Implementation and test calculations*, J. Chem. Phys. **132**, 124102 (2010).

S. Heislbetz, F. Pfeiffer, G. Rauhut, *Configuration selection within vibrational multiconfiguration self-consistent field theory: Application to bridged lithium compounds*, J. Chem. Phys. **134**, 204108 (2011).

P. Meier, D. Oschetzki, F. Pfeiffer, G. Rauhut, *Towards an automated and efficient calculation of resonating vibrational states based on state-averaged multiconfigurational approaches*, J. Chem. Phys. **143**, 244111 (2015). [2ex]

The following *options* are available:

Is the number of active modals for each mode. The smallest meaningful value for `NACT`

is 3, which is the current default. Configurations will only be generated within this space. Note that there is no equivalent in the VMCSCF program to the closed and core orbitals in electronic structure theory.

Is the number of modals on top of the number of active modals. The default is `NVIRT=3`

. The virtual modals are needed for the modal rotations. Modals above the virtual modals are entirely neglected.

By default state-specific VMCSCF calculations will be performed. This may be altered by `AVERAGE=1`

, which calls the state-averaged VMCSCF program. The states, which will automatically be chosen for averaging, are displayed in th output.

Specifies, if only active-virtual modal rotations shall be considered (`ROT=0`

, default) or if active-active modal rotations shall be considered as well (`ROT=1`

). Formally one would need active-active modal rotations for configuration-selective VMCSCF calculations, but the effects are usually extremely small.

The procedure how to determine the rotational angles between the modals within the Jacobi rotations can be specified by this keyword. `CALCANGLE=1`

uses the standard quadratic procedure, which is the current default. `CALCANGLE=2`

uses a cubic equation instead. `CALCANGLE=3`

switches to a full numerical determination of the rotational angle.

Once vibrational states have been defined with the `VIBSTATE`

program (section the VIBSTATE program (VIBSTATE)), the VMCSCF 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.

By default configuration-selective VMCSCF calculations will be performed (`VERSION=3`

). Note that the selection is always performed prior to the VMCSCF iterations but not within them. The selection of configurations can be switched off by (`VERSION=4`

).

Defines the maximum number of simultaneous excitations within the configurations generated from modals of the active space, i.e. Singles, Doubles, Triples, … The maximum excitation level is limited to `CITYPE=9`

, the default is set to 4 for 3D potentials and to 5 for 4D potentials.

`REF=0`

specifies ground-state based VMCSCF calculations, i.e. the initial guess of the VMCSCF calculation is given by the VSCF wave function of the vibrational ground state. In contrast to that, `REF=1`

allows for a state-specific initial guess (default).

By default all VMCSCF calculations are based on ground-state based VSCF modals, `GSMODALS=1`

. `GSMODALS=0`

uses state-specific modals for all VMCSCF calculations.

Controls the maximum number of macroiterations of the VMCSCF program. The default is 25.

Controls the convergence threshold for the macroiterations of the VMCSCF program. The default is 1.d-2.

Controls the maximum number of microiterations of the VMCSCF program. The default is 50.

Controls the convergence threshold for the microiterations of the VMCSCF program. The default is 1.d-5.

Provides additional information within the VMCSCF iterations, once a value larger than 0 (default) is used.

Specifies the record and file, on which the VMCSCF wave function shall be stored. The default is 5950.2.

Once set to 1, this option allows to store the VMCSCF modals in the record of the VSCF modals (which will be overwritten). This allows for VCI calculations with VMCSCF modals.

Besides these VMCSCF specific keywords, a number of option can be used, which are identical with those provided for the VCI program. These keywords are `NDIM, VAM, POT, NBAS, DIPOLE, NDIMDIP, COMBI, INFO, MPG, DIAG, CONT, THRSEL, THRCF, ANALYZE`

.

### Explicit definition of active spaces

`NACT`

,*options*

Within the VMCSCF program the active space can be specified in a general manner (option `NACT`

), which means that the same number of active modals for each mode will be used. Alternatively, one may use the `NACT`

directive. This allows to specify active spaces for the individual modes. Typically such definition refer to just one vibrational state and thus the program should be used within the `USERMODE`

mode.

By default, the active space will be determined by a VMP2 based criterion (`AUTO`

=*on / off*`AUTO=on`

). This can be switched off by`AUTO=off`

. The options exists only for analytical representations of the potential.The number of active modals for mode`MODE(n)`

=*m**n*is set to*m*. This always includes the vibrational ground state. In order to resolve Fermi resonances, one would need to choose*m*at least to be 3, which corresponds to the modals 0, 1 and 2.

## The VIBSTATE program (VIBSTATE)

`VIBSTATE`

,*options* [vibstate]

The `VIBSTATE`

program allows to specify the occupation number vectors for individual vibrational states to be calculated in the following vibrational SCF and vibration correlation programs. Within the input stream, the `VIBSTATE`

program needs to be called prior to the first call of the VSCF program. Note that, the `VIBSTATE`

program needs only to be called, if a limited number of states shall be computed. On the contrary, the `COMBI`

keyword within the `VSCF`

and `VCI`

programs generate large lists of vibrational states, which may result in significant computational effort. The `VIBSTATE`

program necessarily requests the `VSTATE`

directive to be called as described below. If only those vibrational states shall be computed, which are defined by the `VIBSTATE`

program, all subsequent programs need to call the `USERMODE`

option.

### Definition of vibrational states

`VSTATE`

,*options*

This directive specifies the occupation number vector of the vibrational state to be calculated.

The $n$th mode of the molecule is supposed to have the quantum number`MODE(n)`

=*m**m*. If several modes are excited, the option can be repeated accordingly, e.g.`VSTATE,MODE(1)=3,MODE(3)=2,MODE(6)=2`

. All other modes are not excited and thus have the quantum number 0.