/usr/share/psi4/samples/freq-isotope/test.in is in psi4-data 1:1.1-5.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 | #! Vibrational and thermo analysis of several water isotopologs.
#! Demonstrates Hessian reuse for different temperatures and pressures
#! but not for different isotopologs.
# all reference values from NWChem (see other files in this dir) unless otherwise notated
molecule h2o {
units au
O 0.00000000 0.00000000 0.00000000
H 0.00000000 1.93042809 -1.10715266
H 0.00000000 -1.93042809 -1.10715266
}
molecule d2o {
units au
O 0.00000000 0.00000000 0.00000000
H@2.014101779 0.00000000 1.93042809 -1.10715266
H@2.014101779 0.00000000 -1.93042809 -1.10715266
}
molecule hdo {
units au
O 0.00000000 0.00000000 0.00000000
H@2.014101779 0.00000000 1.93042809 -1.10715266
H 0.00000000 -1.93042809 -1.10715266
}
molecule dto {
units au
O 0.00000000 0.00000000 0.00000000
H@2.014101779 0.00000000 1.93042809 -1.10715266
H@3.01604927 0.00000000 -1.93042809 -1.10715266
}
set basis sto-3g
set e_convergence 9
set g_convergence gau_verytight
set scf_type pk
optimize('hf', molecule=h2o)
# can't currently update masses with set_mass() in the mols and redetect the right symmetry
# so have to create as as above and update the geometry
ogeo = h2o.geometry()
d2o.update_geometry()
d2o.set_geometry(ogeo)
hdo.update_geometry()
hdo.set_geometry(ogeo)
dto.update_geometry()
dto.set_geometry(ogeo)
e, wfn = freq('hf', molecule=h2o, return_wfn=True)
compare_values(-74.96590119, get_variable('current energy'), 6, 'H2O E0') #TEST
compare_values(0.024367, get_variable('ZPVE'), 4, 'H2O ZPVE') #TEST
compare_values(0.027199, get_variable('thermal energy correction'), 4, 'H2O dE') #TEST
compare_values(0.028143, get_variable('enthalpy correction'), 4, 'H2O dH') #TEST
entropy = 1000 * psi_hartree2kcalmol * (get_variable('enthalpy correction') - get_variable('gibbs free energy correction')) / get_global_option('t')
compare_values(45.283, entropy, 2, 'H2O S') # molpro #TEST
clean()
set t 400.0
#freq('hf', molecule=h2o)
thermo(wfn, wfn.frequencies())
compare_values(-74.96590119, get_variable('current energy'), 6, 'H2O E0 @400K') #TEST
compare_values(0.024367, get_variable('ZPVE'), 4, 'H2O ZPVE @400K') #TEST
compare_values(0.028170, get_variable('thermal energy correction'), 4, 'H2O dE @400K') #TEST
compare_values(0.029436, get_variable('enthalpy correction'), 4, 'H2O dH @400K') #TEST
entropy = 1000 * psi_hartree2kcalmol * (get_variable('enthalpy correction') - get_variable('gibbs free energy correction')) / get_global_option('t')
compare_values(47.603, entropy, 1, 'H2O S @400K') #TEST
clean()
set t 298.15
freq('hf', molecule=d2o)
compare_values(-74.96590119, get_variable('current energy'), 6, 'D2O E0') #TEST
compare_values(0.017731, get_variable('ZPVE'), 4, 'D2O ZPVE') #TEST
compare_values(0.020566, get_variable('thermal energy correction'), 4, 'D2O dE') #TEST
compare_values(0.021510, get_variable('enthalpy correction'), 4, 'D2O dH') #TEST
entropy = 1000 * psi_hartree2kcalmol * (get_variable('enthalpy correction') - get_variable('gibbs free energy correction')) / get_global_option('t')
compare_values(47.525, entropy, 2, 'D2O S') # molpro #TEST
clean()
e, wfn = freq('hf', molecule=hdo, return_wfn=True)
compare_values(-74.96590119, get_variable('current energy'), 6, 'HDO E0') #TEST
compare_values(0.021103, get_variable('ZPVE'), 4, 'HDO ZPVE') #TEST
compare_values(0.023935, get_variable('thermal energy correction'), 4, 'HDO dE') #TEST
compare_values(0.024878, get_variable('enthalpy correction'), 4, 'HDO dH') #TEST
entropy = 1000 * psi_hartree2kcalmol * (get_variable('enthalpy correction') - get_variable('gibbs free energy correction')) / get_global_option('t')
compare_values(47.828, entropy, 2, 'HDO S') # molpro sym=cs #TEST
# For a symmetry-lowering isotopic substitution like HDO, psi4 recomputes
# the rotational symmetry number with the lower point group. This only
# affects rotational entropy. That's why the above (non-default) molpro
# value from the above was computed with sym=cs. To replicate molpro &
# qchem's default, have to set rotational_symmetry_number to H2O value.
set rotational_symmetry_number 2
thermo(wfn, wfn.frequencies())
compare_values(-74.96590119, get_variable('current energy'), 6, 'HDO E0') #TEST
compare_values(0.021103, get_variable('ZPVE'), 4, 'HDO ZPVE') #TEST
compare_values(0.023935, get_variable('thermal energy correction'), 4, 'HDO dE') #TEST
compare_values(0.024878, get_variable('enthalpy correction'), 4, 'HDO dH') #TEST
entropy = 1000 * psi_hartree2kcalmol * (get_variable('enthalpy correction') - get_variable('gibbs free energy correction')) / get_global_option('t')
compare_values(46.450, entropy, 2, 'HDO S') # molpro default #TEST
clean()
# Could just reset the symmetry number to 1 since that's the correct
# value for DTO. But the below will calculate the default.
psi4.revoke_global_option_changed('rotational_symmetry_number')
freq('hf', molecule=dto)
compare_values(-74.96590119, get_variable('current energy'), 6, 'DTO E0') #TEST
compare_values(0.016317, get_variable('ZPVE'), 4, 'DTO ZPVE') #TEST
compare_values(0.019154, get_variable('thermal energy correction'), 4, 'DTO dE') #TEST
compare_values(0.020098, get_variable('enthalpy correction'), 4, 'DTO dH') #TEST
entropy = 1000 * psi_hartree2kcalmol * (get_variable('enthalpy correction') - get_variable('gibbs free energy correction')) / get_global_option('t')
compare_values(49.603, entropy, 2, 'DTO S') # molpro sym=cs #TEST
|