***,h2s Diabatization and NACME calculation

gprint,orbitals,civector

symmetry,x
orient,noorient             !noorient should always be used for diabatization
geometry={
         s;
         h1,s,r1;
         h2,s,r2,h1,theta}

basis=avdz                  !This basis is too small for real application

r1=2.5                      !Reference geometry
theta=[92]

r=[2.55,2.60]               !Displaced geometries
dr=[0,0.01,-0.01]           !Samll displacements for finite difference NACME calculation

reforb1=2140.2              !Orbital dumprecord at reference geometry
refci=6000.2                !MRCI record at reference geometry
savci=6100.2                !MRCI record at displaced geometries

text,compute wavefunction at reference geometry (C2v)
r2=r1

{hf;occ,9,2;wf,18,2,4;orbital,2100.2}

{multi;occ,9,2;closed,4,1;
wf,18,2;state,2;            !1B1 and 1A2 states
natorb,reforb1              !Save reference orbitals on reforb1
noextra}                    !Dont use extra symmetries

{ci;occ,9,2;closed,4,1;      !MRCI at reference geometry
wf,18,2,0;state,2;          !1B1 and 1A2 states
orbital,reforb1             !Use orbitals from previous CASSCF
save,refci}                 !Save MRCI wavefunction

Text,Displaced geometries

do i=1,#r                   !Loop over different r values
data,truncate,savci+1       !truncate dumpfile after reference
reforb=reforb1

do j=1,3                    !Loop over small displacements for NACME
r2=r(i)+dr(j)               !Set current r2

{multi;occ,9,2;closed,4,1;
wf,18,2,0;state,2;          !Wavefunction definition
start,reforb                !Starting orbitals
orbital,3140.2+j;           !Dumprecord for orbitals
diab,reforb                 !Generate diabatic orbitals relative to reference geometry
noextra}                    !Dont use extra symmetries

reforb=3141.2               !Use orbitals for j=1 as reference for j=2,3

{ci;occ,9,2;closed,4,1;
wf,18,2,0;state,2;
orbital,diabatic            !Use diabatic orbitals
save,savci+j}               !Save MRCI for displaced geometries

eadia=energy                !Save adiabatic energies for use in ddr
if(j.eq.1) then
e1(i)=energy(1)             !Save adiabatic energies for table printing
e2(i)=energy(2)
end if

{ci;trans,savci+j,savci+j;  !Compute transition densities at R2+DR(j)
dm,7000.2+j}                !Save transition densities on this record
{ci;trans,savci+j,refci;    !Compute transition densities between R2+DR(j) and R1
dm,7100.2+j}                !Save transition densities on this record
{ci;trans,savci+j,savci+1;  !Compute transition densities between R and R2+DR(j)
dm,7200.2+j}                !Save transition densities on this record

{ddr
density,7000.2+j,7100.2+j   !Densities for <R2+DR||R2+DR> and <R2+DR||R1>
orbital,3140.2+j,2140.2     !Orbitals for <R2+DR||R2+DR> and <R2+DR||R1>
energy,eadia(1),eadia(2)    !Adiabatic energies
mixing,1.2,2.2}             !Compute mixing angle and diabatic energies

if(j.eq.1) then             !Store diabatic energies for R2  (DR(1)=0)
  mixci(i)=mixangci(1)      !Mixing angle obtained from ci vectors only
  h11ci(i)=hdiaci(1)        !Diabatic energies obtained from ci vectors only
  h21ci(i)=hdiaci(2)        !HDIA contains the lower triangle of the diabatic hamiltonian
  h22ci(i)=hdiaci(3)
  mixtot(i)=mixang(1)       !Mixing angle from total overlap (including first-order correction)
  h11(i)=hdia(1)            !Diabatic energies obtained from total overlap
  h21(i)=hdia(2)
  h22(i)=hdia(3)
end if

mix(j)=mixang(1)            !Store mixing angles for R2+DR(j)

enddo                       !End loop over j

dchi(i)=(mix(3)-mix(2))/(dr(2)-dr(3))*pi/180  !Finite difference derivative of mixing angle

{ddr
density,7201.2,7202.2,7203.2   !Compute NACME using 3-point formula
orbital,3141.2,3142.2,3143.2
states,2.2,1.2}
nacmeci(i)=nacme

{table,r,mixci,mixtot,dchi,nacmeci
Title,Mixing angles and non-adiabatic coupling matrix elements for H2S
format,'(f10.2,4f14.4)'
sort,1
}

{table,r,e1,e2,h11ci,h22ci,h21ci
Title,Diabatic energies for H2S, obtained from CI-vectors
format,'(f10.2,5f16.8)'
sort,1}

{table,r,e1,e2,h11,h22,h21
title,Diabatic energies for H2S, obtained from CI-vectors and orbital correction
format,'(f10.2,5f16.8)'
sort,1}

enddo