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 computed analytically (cf. section 16.9.2).
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.