molecule water_water {
0 1
O   -0.066999140   0.000000000   1.494354740
H    0.815734270   0.000000000   1.865866390
H    0.068855100   0.000000000   0.539142770
--
0 1
O    0.062547750   0.000000000  -1.422632080
H   -0.406965400  -0.760178410  -1.771744500
H   -0.406965400   0.760178410  -1.771744500 
}

set basis  aug-cc-pVTZ
set opt_coordinates both
set geom_maxiter 2000
# set disp_size 0.001
set g_convergence interfrag_tight
memory 100000mb
# set molden_write true

# optimize('mp2', molecule=water_water)

set freeze_core true

energy('mp2', bsse_type='cp', molecule=water_water)
print("CP Corrected MP2 Binding Energy")
print("CP_Corrected_2_Body_Interaction_Energy|{}|kcal/mol".format(psi_hartree2kcalmol*get_variable('cp-corrected 2-body interaction energy')))


energy('sapt2+(3)dmp2')
print("SAPT2+(3)dmp2 Binding Energy")
print("SAPT_Disp_Energy|{}|kcal/mol".format(psi_hartree2kcalmol*get_variable('sapt disp energy')))
print("SAPT_Elst_Energy|{}|kcal/mol".format(psi_hartree2kcalmol*get_variable('sapt elst energy')))
print("SAPT_Exch_Energy|{}|kcal/mol".format(psi_hartree2kcalmol*get_variable('sapt exch energy')))
print("SAPT_Ind_Energy|{}|kcal/mol".format(psi_hartree2kcalmol*get_variable('sapt ind energy')))
print("SAPT_Total_Energy|{}|kcal/mol".format(psi_hartree2kcalmol*get_variable('sapt total energy')))

energy('sapt2+(3)-ct')
print("SAPT2+(3) Charge Transfer Estimate|{}|kcal/mol".format(psi_hartree2kcalmol*get_variable('sapt ct energy')))

set t 298.15
energy, wavefunction = frequencies('MP2', molecule=water_water, return_wfn=True)

electronic_energy = psi_hartree2kcalmol*(get_variable('enthalpy')-get_variable('enthalpy correction'))
entropy = psi_hartree2kcalmol * (get_variable('enthalpy correction') - get_variable('gibbs free energy correction')) / get_global_option('t')

print("Thermodynamic Analysis at 298.15K")
print("Electronic_Energy|{}|kcal/mol".format(electronic_energy))
print("Zero_Point_Vibrational_Energy|{}|kcal/mol".format(psi_hartree2kcalmol*get_variable('zpve')))
print("Enthalpy_Correction|{}|kcal/mol".format(psi_hartree2kcalmol*get_variable('enthalpy correction'))) 
print("Enthalpy|{}|kcal/mol".format(psi_hartree2kcalmol*get_variable('enthalpy')))
print("Thermal_Energy_Correction|{}|kcal/mol".format(psi_hartree2kcalmol*get_variable('thermal energy correction')))
print("Thermal_Energy|{}|kcal/mol".format(psi_hartree2kcalmol*get_variable('thermal energy'))) 
print("Gibbs_Free_Energy_Correction|{}|kcal/mol".format(psi_hartree2kcalmol*get_variable('gibbs free energy correction')))
print("Gibbs_Free_Energy|{}|kcal/mol".format(psi_hartree2kcalmol*get_variable('gibbs free energy')))
print("Entropy|{}|kcal/mol.K".format(entropy))