[molpro-user] Potential Energy Surface Diabatization

Shapera, Ethan eshaper at sandia.gov
Tue Oct 3 23:53:16 CEST 2017


We are attempting to reproduce Figures 2 and 3 of http://dx.doi.org/10.1063/1.479214 showing the CASSCF and MRCI potential energy surfaces (PESs) of H2S.  Our input files are closely based on the h2s_diab1 example.  The multi and adiabatic MRCI PES match the literature well.  However, the diabatic MRCI calculations do not match the paper and produce inconsistent results.  With the script provided below, we obtain a conical intersection at 2.7 a0 with the MRCI calculation of the diabatized PES.  Figure 2 of the included paper shows a lack of intersection between the diabatic surfaces.  Further, by reversing the order in which we take the non-fixed H-S bond length values we obtain a conic intersection at ~3.5 a0.



Similar issues arise with recreating Figure 3 of PES versus bond angle.  Starting with small bond angles produces a conical intersection at ~110 degrees; starting with large bond angles has a conical intersection at ~95 degrees, relatively close to the literature.


How should we fix this issue and consistently produce reliable diabatized MRCI surfaces?



***,h2s Diabatization
memory,3,m

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=aug-cc-pVQZ                  !This basis is too small for real application

r1=2.7                      !Reference geometry
theta=92

r=[ 3.6,3.55,3.5,3.45,3.4,3.35,3.3,3.25,3.2,3.15,3.1,3.05,3.0,2.95,2.9,2.85,2.8,2.75,2.7,2.65,2.6,2.55,2.5,2.45,2.4,2.35,2.3,2.25,2.2,2.15,2.1,2.05,2.0]

reforb=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=3.6

{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,reforb               !Save reference orbitals on reforb
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,reforb              !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
r2=r(i)                     !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;             !Dump record for orbitals
diab,reforb                 !Generate diabatic orbitals relative to reference geometry
noextra}                    !Dont use extra symmetries

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

e1(i)=energy(1)             !Save adiabatic energies
e2(i)=energy(2)

{ci;trans,savci,savci        !Compute transition densities at R2
dm,7000.2}                  !Save transition densities on this record
{ci;trans,savci,refci;       !Compute transition densities between R2 and R1
dm,7100.2}                  !Save transition densities on this record

{ddr
density,7000.2,7100.2       !Densities for <R2||R2> and <R2||R1>
orbital,3140.2,2140.2       !Orbitals for <R2||R2> and <R2||R1>
energy,e1(i),e2(i)          !Adiabatic energies
mixing,1.2,2.2}             !Compute mixing angle and diabatic energies

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)
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)


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

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

enddo                       !end loop over i

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.molpro.net/pipermail/molpro-user/attachments/20171003/9e56e110/attachment.html>


More information about the Molpro-user mailing list