diff --git a/Composite-Methods/Gaussian_n.py b/Composite-Methods/Gaussian_n.py new file mode 100644 index 00000000..a1657030 --- /dev/null +++ b/Composite-Methods/Gaussian_n.py @@ -0,0 +1,212 @@ +""" +A reference implementation of Gaussian-n composite methods. +""" + +__authors__ = "Eric J. Berquist" +__credits__ = ["Eric J. Berquist"] + +__copyright__ = "(c) 2014-2018, The Psi4NumPy Developers" +__license__ = "BSD-3-Clause" +__date__ = "2018-09-30" + +import numpy as np +np.set_printoptions(precision=5, linewidth=200, suppress=True) +import psi4 + +psi4.set_memory("2 GB") +psi4.core.set_output_file("output.dat", False) + +print("\n Performing a G1 calculation.\n") + +# The starting geometry is not critical, since all compound methods include +# one or more geometry optiimizations. +mol = psi4.geometry(""" +O 0.000000000000 0.000000000000 -0.075791843589 +H 0.000000000000 -0.866811828967 0.601435779270 +H 0.000000000000 0.866811828967 0.601435779270 +symmetry c1 +""") + +psi4.set_options({ + "basis": "6-31G*", + "scf_type": "direct", + "df_scf_guess": False, + "g_convergence": "gau", +}) + +print("1a. Geometry optimization: HF/6-31G*") +psi4.optimize("HF") +print("1b. Harmonic frequencies: HF/6-31G*") +psi4.frequencies("HF") + +T = psi4.core.get_global_option("T") +zpe = psi4.core.scalar_variable("ZPVE") +# du = psi4.core.scalar_variable("INTERNAL ENERGY CORRECTION") +du = psi4.core.scalar_variable("THERMAL ENERGY CORRECTION") +dh = psi4.core.scalar_variable("ENTHALPY CORRECTION") +dg = psi4.core.scalar_variable("GIBBS FREE ENERGY CORRECTION") + +psi4.core.clean() + +print("2. Geometry optimization: MP2/6-31G* [no frozen core]") + +psi4.set_options({ + "basis": "6-31G*", + "scf_type": "direct", + "df_scf_guess": False, + "mp2_type": "conv", + "freeze_core": "false", + "g_convergence": "gau", +}) + +psi4.optimize("MP2") +psi4.core.clean() + +print("3. Single-point energy: MP4(SDTQ)/6-311G** [frozen core]") + +psi4.set_options({ + "basis": "6-311G**", + "scf_type": "direct", + "df_scf_guess": False, + "freeze_core": "true", +}) + +e_step3, wfn = psi4.energy("MP4", return_wfn=True) +psi4.core.clean() + +print("4. Single-point energy: MP4(SDTQ)/6-311+G** [frozen core]") + +psi4.set_options({ + "basis": "6-311+G**", + "scf_type": "direct", + "df_scf_guess": False, + "freeze_core": "true", +}) + +e_step4 = psi4.energy("MP4") +psi4.core.clean() + +e_step4 -= e_step3 + +print("5. Single point energy: MP4(SDTQ)/6-311G**(2df) [frozen core]") + +psi4.set_options({ + "basis": "6-311G(2df,p)", + "scf_type": "direct", + "df_scf_guess": False, + "freeze_core": "true", +}) + +e_step5 = psi4.energy("MP4") +psi4.core.clean() + +e_step5 -= e_step3 + +print("6. Single-point energy: QCISD(T)/6-311G** [frozen core]") + +psi4.set_options({ + "basis": "6-311G**", + "scf_type": "direct", + "df_scf_guess": False, + "freeze_core": "true", +}) + +e_step6 = psi4.energy("QCISD(T)") +psi4.core.clean() + +e_step6 -= e_step3 + +e_combined = e_step3 + e_step4 + e_step5 + e_step6 + +print("7. High-level correction (HLC)") + +nfc = wfn.nfrzc() +nalpha = wfn.nalpha() - nfc +nbeta = wfn.nbeta() - nfc +e_hlc = (-0.19 * nalpha) - (5.95 * nbeta) +e_hlc /= 1000 + +e_e = e_combined + e_hlc + +print("8. ZPE addition") + +scale_fac = 0.8929 +e_zpe = scale_fac * zpe +e_thermal = du - (1 - scale_fac) * zpe +g1_0k = e_e + e_zpe +g1_energy = g1_0k - zpe + du +g1_enthalpy = g1_0k - zpe + dh +g1_free_energy = g1_0k - zpe + dg + +print("E(ZPE)= {:10.6f}".format(e_zpe)) +print("E(Thermal)= {:10.6f}".format(e_thermal)) +print("E(QCISD(T))= {:10.6f}".format(e_step6 + e_step3)) +print("E(Empiric)= {:10.6f}".format(e_hlc)) +print("DE(Plus)= {:10.6f}".format(e_step4)) +print("DE(2DF)= {:10.6f}".format(e_step5)) +print("G1(0 K)= {:10.6f}".format(g1_0k)) +print("G1 Energy= {:10.6f}".format(g1_energy)) +print("G1 Enthalpy= {:10.6f}".format(g1_enthalpy)) +print("G1 Free Energy= {:10.6f}".format(g1_free_energy)) + +print("\nAdding components for G2 and G2(MP2) calculation.\n") + +print("Single-point energy: MP2/6-311+G(3df,2p) [frozen core]") + +psi4.set_options({ + "basis": "6-311+G(3df,2p)", + "scf_type": "direct", + "df_scf_guess": False, + "mp2_type": "conv", + "freeze_core": "true", +}) + +psi4.energy("MP2") +e_mp2_plus_3df2p = psi4.get_variable("MP2 TOTAL ENERGY") +psi4.core.clean() + +# e_mp2_base = -76.263654158359 +# e_mp2_plus = -76.274546774908 +# e_mp2_2df = -76.298942654928 +# This calculation can be avoided by rearranging the expressions. +# e_mp2_plus_2df +# e_mp2_plus_3df2p = -76.318108464496 + +# de_mp2_plus_2df +de_mp2_plus = e_mp2_plus - e_mp2_base +de_mp2_2df = e_mp2_2df - e_mp2_base + +# de_1 = de_plus_2df + de_plus + de_2df +# de_2 = e_mp2_plus_3df2p - e_mp2_plus_2df +# de = de_1 + de_2 +de = e_mp2_plus_3df2p - e_mp2_2df - e_mp2_plus + e_mp2_base + +# Calculate the number of valence pairs. +npairs = min(nalpha, nbeta) +e_hlc_2 = 1.14 * npairs / 1000 + +g2_0k = g1_0k + de + e_hlc_2 +g2_energy = g2_0k - zpe + du +g2_enthalpy = g2_0k - zpe + dh +g2_free_energy = g2_0k - zpe + dg + +# G2(MP2) bits + +de_mp2 = e_mp2_plus_3df2p - e_mp2_base +g2mp2_0k = (e_step6_qci + e_step3) + de_mp2 + e_hlc_1 + e_hlc_2 + e_zpe + +g2mp2_energy = g2mp2_0k - zpe + du +g2mp2_enthalpy = g2mp2_0k - zpe + dh +g2mp2_free_energy = g2mp2_0k - zpe + dg + +print("e(Delta-G2)= {:10.6f}".format(de)) +print("E(G2-Empiric)= {:10.6f}".format(e_hlc_2)) +print("G2(0 K)= {:10.6f}".format(g2_0k)) +print("G2 Energy= {:10.6f}".format(g2_energy)) +print("G2 Enthalpy= {:10.6f}".format(g2_enthalpy)) +print("G2 Free Energy= {:10.6f}".format(g2_free_energy)) +print("DE(MP2)= {:10.6f}".format(de_mp2)) +print("G2MP2(0 K)= {:10.6f}".format(g2mp2_0k)) +print("G2MP2 Energy= {:10.6f}".format(g2mp2_energy)) +print("G2MP2 Enthalpy= {:10.6f}".format(g2mp2_enthalpy)) +print("G2MP2 Free Energy= {:10.6f}".format(g2mp2_free_energy)) diff --git a/tests/test_RI_composite.py b/tests/test_RI_composite.py new file mode 100644 index 00000000..f11a8a0d --- /dev/null +++ b/tests/test_RI_composite.py @@ -0,0 +1,5 @@ +from utils import exe_py + + +def test_Gaussian_n(workspace): + exe_py(workspace, tdir, 'Gaussian_n')