20.5.8 Diabatic orbitals

In order to construct diabatic states, it is necessary to determine the mixing of the diabatic states in the adiabatic wavefunctions. In principle, this mixing can be obtained by integration of the non-adiabatic coupling matrix elements. Often, it is much easier to use an approximate method, in which the mixing is determined by inspection of the CI coefficients of the MCSCF or CI wavefunctions. This method is applicable only if the orbital mixing is negligible. For CASSCF wavefunctions this can be achieved by maximizing the overlap of the active orbitals with those of a reference geometry, at which the wavefunctions are assumed to be diabatic (e.g. for symmetry reasons). The orbital overlap is maximized using using the new DIAB command in the MCSCF program. Only the active orbitals are transformed.

This procedure works as follows: first, the orbitals are determined at the reference geometry. Then, the calculations are performed at displaced geometries, and the "diabatic" active orbitals, which have maximum overlap with the active orbitals at the reference geometry, are obtained by adding a DIAB directive to the input:

Old form (Molpro96, obsolete):

DIAB,orbref, orbsav, orb1,orb2,pri

New form:

DIAB,orbref[,TYPE=orbtype][,STATE=state] [,SPIN=spin][,MS2=ms2][,SAVE=orbsav]
[,ORB1=orb1, ORB2=orb2][,PRINT=pri][,METHOD=method]

Here orbref is the record holding the orbitals of the reference geometry, and orbsav is the record on which the new orbitals are stored. If orbsav is not given (recommended!) the new orbitals are stored in the default dump record (2140.2) or the one given on the ORBITAL directive (see section 16.5.3). In contrast to earlier versions of MOLPRO it is possible that orbref and orbsav are the same. The specifications TYPE, STATE, SPIN can be used to select specific sets of reference orbitals, as described in section 2.16. orb1, orb2 is a pair of orbitals for which the overlap is to be maximized. These orbitals are specified in the form number.sym, e.g. 3.1 means the third orbital in symmetry 1. If orb1, orb2 are not given, the overlap of all active orbitals is maximized. pri is a print parameter. If this is set to 1, the transformation angles for each orbital are printed for each Jacobi iteration. method determines the diabatization method. method=1 (default): use Jacobi rotations; method=2: use block diagonalization. Both methods yield very similar results. method=2 must only be used for CASSCF wavefunctions. method=-1 and method=-2: as the positive values, but AO overlap matrix of the current geometry is used. This minimizes the change of the MO coefficients, rather than maximizing the overlap to the neighbouring orbitals.

Using the defaults described above, the following input is sufficient in most cases:

DIAB,orbref

Using Molpro98 is is not necessary any more to give any GEOM and DISPL cards. The displacements and overlap matrices are computed automatically (the geometries are stored in the dump records, along with the orbitals).

The diabatic orbitals have the property that the sum of orbital and overlap contributions in the non-adiabatic coupling matrix elements become approximately zero, such that the adiabatic mixing occurs only through changes of the CI coefficients. This allows to determine the mixing angle directly from the CI coefficients, either in a simple way as described for instance in J. Chem. Phys. 89, 3139 (1988), or in a more advanced manner as described by Pacher, Cederbaum, and Köppel in J. Chem. Phys. 89, 7367 (1988). Recently, an automatic procedure, as described in J. Chem. Phys. 102, 0000, (1999) has been implemented into MOLPRO. This is available in Version 99.1 and later and is described in section 27.

Below we present an example for the first two excited states of H$_2$S, which have $B_1$ and $A_2$ symmetry in $C_{2v}$, and $A''$ symmetry in $C_S$. We first perform a reference calculation in $C_{2v}$ symmetry, and then determine the diabatic orbitals for displaced geometries in $C_S$ symmetry. Each subsequent calculation uses the previous orbitals as reference. One could also use the orbitals of the $C_{2v}$ calculation as reference for all other calculations. In this case one would have to take out the second-last input card, which sets reforb=2141.2.

Input: h2s_diab.com

Output: h2s_diab.out

See section 27 for the automatic generation of diabatic energies.