molecule water_water {
0 1
    O           -1.478998795354    -0.052839545458     0.005102811541
    H           -1.833481358066     0.836940113830    -0.016170713458
    H           -0.520967139608     0.065783864884    -0.000246507515
--
0 1
    O            1.407306450244     0.049995986940    -0.004848420469
    H            1.748173821752    -0.409838411044     0.765769338201
    H            1.744084218674    -0.447756230480    -0.753389488157 
}

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, dertype=1)

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