***,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