***,HF dimer mp2/CP optimization

basis=avdz

maxit=20                !max number of iterations

text,OPTIMIZED VALUES OF GEOMETRY VARIABLES 

RFF=      5.31431160
R1=       1.75768738 
R2 =      1.75298524  
THETA1 =  7.03780227 
THETA2 = 111.25930975

geometry={x;noorient    !noorient must be specified since gradients are added
      f1
      f2  f1  rff
      h1  f1  r1   f2  theta1
      h2  f2  r2   f1  theta2  h1  180.}

do iter=1,maxit         !optimization loop

text, CALCULATION AT LARGE SEPARATION

rff_save=rff            !save current rff distance
rff=1000                !calculation at large separation

text, HF1 MONOMER

dummy,f2,h2;            !second hf is now dummy
hf;                     !scf for first monomer
mp2;                    !mp2 for first monomer
ehf1inf=energy          !save mp2 energy in variable
forces,numerical,startcmd=hf

text, HF2 MONOMER

dummy,f1,h1;            !first hf is now dummy
hf;                     !scf for second monomer
mp2;                    !mp2 for second monomer
ehf2inf=energy          !save mp2 energy in variable
forces,numerical,startcmd=hf
add,1                   !add from previous gradient

rff=rff_save            !reset HF - HF distance to current value

text, HF1 CP MONOMER

dummy,f2,h2;            !second hf is now dummy
hf;                     !scf for first monomer
mp2;                    !mp2 for first monomer
ehf1=energy             !save mp2 energy in variable
forces,numerical,startcmd=hf
add,-1                  !subtract from previous gradient

text, HF2 CP MONOMER

dummy,f1,h1;            !first hf is now dummy
hf;                     !scf for second monomer
mp2;                    !mp2 for second monomer
ehf2=energy             !save mp2 energy in variable
forces,numerical,startcmd=hf
add,-1                  !subtract from previous gradient

dummy                   !reset dummies

text, DIMER CALCULATION
hf;                     !scf for dimer
mp2;                    !mp2 for dimer
edimer=energy           !save mp2 energy in variable
forces,numerical,startcmd=hf
add,1                   !add to previous gradient

decpc1=ehf1inf-ehf1     !counterpoise correction for first monomer
decpc2=ehf2inf-ehf2     !counterpoise correction for second monomer
decpc=decpc1+decpc2     !total counterpoise correction
de=edimer-ehf1-ehf2     !Interaction energy relative to unrelaxed monomers
etot=edimer+decpc       !total cpc corrected energy

opt;                    !find next geometry
show,optconv            !show convergence parameter
if(optconv.lt.0.0003) goto,lab:  !test for convergence



enddo            !end of optimzation loop

lab:

text, compute optimized monomer energy 

rhf=r1
geometry={h1
          F1,H1,rhf}
hf;                     !scf for relaxed monomer
mp2;                    !mp2 for relaxed monomer
ehf=energy              !save mp2 energy in variable
optg                    !optimize monomer structure

text, compute interaction energy including relaxation

erelax = 2*ehf - ehf1inf - ehf2inf   !relaxation energy
derel=de-erelax         !interaction energy relative to relaxed monomoers

text, results 

show,erelax,de,derel,r1,r2,rhf,rff,theta1,theta2

decpc=decpc*tocm        !total counterpoise correction in cm-1
erelax=erelax*tocm      !relaxation energy in cm-1
de=de*tocm              !interaction energy relative to unrelaxed monomers in cm-1
derel=derel*tocm        !interaction energy relative to relaxed monomers in cm'1


---;


