/usr/share/psi4/samples/freq-isotope/input.dat 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 | #! 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)
entropy = 1000 * psi_hartree2kcalmol * (get_variable('enthalpy correction') - get_variable('gibbs free energy correction')) / get_global_option('t')
clean()
set t 400.0
#freq('hf', molecule=h2o)
thermo(wfn, wfn.frequencies())
entropy = 1000 * psi_hartree2kcalmol * (get_variable('enthalpy correction') - get_variable('gibbs free energy correction')) / get_global_option('t')
clean()
set t 298.15
freq('hf', molecule=d2o)
entropy = 1000 * psi_hartree2kcalmol * (get_variable('enthalpy correction') - get_variable('gibbs free energy correction')) / get_global_option('t')
clean()
e, wfn = freq('hf', molecule=hdo, return_wfn=True)
entropy = 1000 * psi_hartree2kcalmol * (get_variable('enthalpy correction') - get_variable('gibbs free energy correction')) / get_global_option('t')
# 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())
entropy = 1000 * psi_hartree2kcalmol * (get_variable('enthalpy correction') - get_variable('gibbs free energy correction')) / get_global_option('t')
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)
entropy = 1000 * psi_hartree2kcalmol * (get_variable('enthalpy correction') - get_variable('gibbs free energy correction')) / get_global_option('t')
|