30.1 The DDR procedure

In order to compute the coupling matrix elements by finite differences, one has to compute and store the wavefunctions at two (first-order algorithm) or three (second-order algorithm) slightly displaced geometries. The order of these calculations is arbitrary.

The typical strategy is as follows:

1.) Compute the wavefunction at the reference geometry. The wavefunctions for both states have to be stored using the SAVE command of the CI program. If the matrix elements are computed for MCSCF wavefunctions, it is necessary to recompute the wavefunction with the CI program, using the NOEXC option. The transition density matrix is stored using the DM directive of the CI program.

2.) Compute the wavefunctions at the (positively) displaced geometry and store the CI wavefunction in a second record.

3.) If the second-order (three-point) method is used, step (2) is repeated at a (negatively) displaced geometry.

4.) Compute the transition density matrices between the states at the reference geometry and the displaced geometr(ies). This is done with the TRANS directive of the CI program.

5.) Finally, the DDR program is used to assemble the matrix element. Using the first-order two-point method, only a single input line is needed:

DDR, dr, orb1, orb2, trdm2

where dr is the geometry increment used as denominator in the finite difference method, orb1 is the record holding the orbitals of the reference geometry, orb2 is the record holding the orbitals of the displaced geometry, and trdm2 is the record holding the transition density matrix computed from the CI-vectors at R and R+DR.

If central differences (three points) are used, the input is as follows:

DDR,2*dr
ORBITAL,orb1,orb2,orb3
DENSITY,trdm1,trdm2,trdm3

where dr, orb1, orb2 are as above, and orb3 is the record holding the orbitals at the negatively displaced geometry.

trdm1, trdm2, trdm3 are the records holding the transition densities $\gamma(R\vert R)$, $\gamma(R\vert R+DR)$, and $\gamma(R\vert R-DR)$, respectively.

If more than two states are computed simultaneously, the transition density matrices for all pairs of states will be stored in the same record. It is then necessary to specify for which states the coupling is to be computed using the STATE directive:

STATE,state$_1$, state$_2$

where state$_i$ is of the form istate.isym (the symmetries of both states must be the same, and it is therefore sufficient to specify the symmetry of the first state).

As an example the input for first-order and second-order calculations is given below. The calculation is repeated for a range of geometries, and at the end of the calculation the results are printed using the TABLE command.

In the calculation shown, the "diabatic" CASSCF orbitals are generated in the two CASSCF calculations at the displaced geometries by maximizing the overlap with the orbitals at the reference geometry. This is optional, and (within the numerical accuacy) does not influence the final results. However, the relative contributions of the orbital, overlap and CI contributions to the NACME are modified. If diabatic orbitals are used, which change as little as possible as function of geometry, the sum of overlap and orbital contribution is minimized, and to a very good approximation the NACME could be obtained from the CI-vectors alone.

Input: lif_nacme.com

Output: lif_nacme.out

This calculation produces the following table:

 Non-adiabatic couplings for LiF

       R      NACME1P       NACME1M       NACMEAV       NACME2
    10.0   -0.22828936   -0.22328949   -0.22578942   -0.22578942
    10.5   -0.51777034   -0.50728914   -0.51252974   -0.51252974
    11.0    0.76672943    0.76125391    0.76399167    0.76399167
    11.5    0.42565202    0.42750263    0.42657733    0.42657733
    12.0    0.19199878    0.19246799    0.19223338    0.19223338

Note that the sign changes because of a phase change of one of the wavefunctions. In order to keep track of the sign, one has to inspect both the orbitals and the ci-vectors.