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 ,
, and
, 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.
In that case, and also when there are just two states whose spatial
symmetry is not 1, it is
necessary to specify for which states the coupling is to be computed using
the `STATE` directive:

`STATE`,*state, state*

where *state* 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.

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.

molpro@molpro.net 2019-06-17