Non adiabatic coupling matrix elements
Non-adiabatic coupling matrix elements can be computed by finite differences for MCSCF
or CI
wavefunctions using the DDR
program. For state-averaged MCSCF
wavefunctions, they can also be computed analytically (cf. section difference gradients for SA-MCSCF).
Note that present numerical procedure has been much simplified relative to Molpro96
. No GEOM
and DISPL
input cards are needed any more, and the three necessary calculations can be done in any order.
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|R)$, $\gamma(R|R+DR)$, and $\gamma(R|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. 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$_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.
- examples/lif_nacme.inp
***,lif non-adiabatic coupling basis,f=avdz,li=vdz !define basis r=[10.0,10.5,11.0,11.5,12.0] !define bond distances dr=0.01 !define increment geometry={li;f,li,rlif} !define geometry rlif=3 !first calculation at R=3 {hf;occ,4,1,1} !SCF {multi;closed,3; !CASSCF, 3 inactive orbitals wf,12,1;state,2; !Two 1A1 states orbital,2140.2} !dump orbitals to record 2140.2 do i=1,#r !loop over geometries rlif=r(i) !set bond distance {multi;closed,3; !CASSCF, 3 inactive orbitals wf,12,1;state,2; !Two 1A1 states orbital,2140.2} !Overwrite previous orbitals by present ones {ci;state,2;noexc; !CI for 2 states, no excitations save,6000.2; !save wavefunction to record 6000.2 dm,8000.2} !save (transition) densities to record 8000.2 rlif=r(i)+dr !increment bond distance by dr {multi;closed,3; !same CASSCF as above wf,12,1;state,2; !Two 1A1 states start,2140.2; !start with orbitals from reference geometry orbital,2141.2; !save orbitals to record 2141.2 diab,2140.2} !generate diabatic orbitals by maximizing the !overlap with the orbitals at the reference geometry {ci;state,2;noexc;save,6001.2} !CI for 2 states, wavefunction saved to record 6001.2 {ci;trans,6000.2,6001.2; !Compute overlap and transition density <R|R+DR> dm,8100.2} !Save transition density to record 8100.2 rlif=r(i)-dr !repeat at r-dr {multi;closed,3; !same CASSCF as above wf,12,1;state,2; !Two 1A1 states start,2140.2; !start with orbitals from reference geometry orbital,2142.2; !save orbitals to record 2142.2 diab,2140.2} !generate diabatic orbitals by maximizing the !overlap with the orbitals at the reference geometry {ci;state,2;noexc;save,6002.2} !CI for 2 states, wavefunction saved to record 6002.2 {ci;trans,6000.2,6002.2; !Compute overlap and transition density <R|R-DR> dm,8200.2} !Save transition density to record 8200.2 {ddr,dr,2140.2,2141.2,8100.2} !compute NACME using 2-point formula (forward difference) nacme1p(i)=nacme !store result in variable nacme1p {ddr,-dr,2140.2,2142.2,8200.2} !compute NACME using 2-point formula (backward difference) nacme1m(i)=nacme !store result in variable nacme1m {ddr,2*dr !compute NACME using 3-point formula orbital,2140.2,2141.2,2142.2; !orbital records for R, R+DR, R-DR density,8000.2,8100.2,8200.2} !transition density records for R, R+DR, R-DR nacme2(i)=nacme !store result in variable nacme2 end do !end of loop over differend bond distances nacmeav=(nacme1p+nacme1m)*0.5 !average the two results forward and backward differences table,r,nacme1p,nacme1m,nacmeav,nacme2 !print a table with results title,Non-adiabatic couplings for LiF !title for table
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.
Please also see diabatic_orbitals.