From 4d36b7976ba278b2669276e31b87fa8e39811353 Mon Sep 17 00:00:00 2001 From: Eric Berquist Date: Sun, 30 Sep 2018 16:22:56 -0400 Subject: [PATCH 1/4] Initial implementation of G1 --- Composite-Methods/G1.py | 150 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 150 insertions(+) create mode 100644 Composite-Methods/G1.py diff --git a/Composite-Methods/G1.py b/Composite-Methods/G1.py new file mode 100644 index 00000000..628536d1 --- /dev/null +++ b/Composite-Methods/G1.py @@ -0,0 +1,150 @@ +""" +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.get_variable("ZPVE") +# du = psi4.get_variable("INTERNAL ENERGY CORRECTION") +du = psi4.get_variable("THERMAL ENERGY CORRECTION") +dh = psi4.get_variable("ENTHALPY CORRECTION") +dg = psi4.get_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)) From 62d6d64a1aec8a4ec8c2031b5099ce2bf064d8df Mon Sep 17 00:00:00 2001 From: Eric Berquist Date: Thu, 15 Oct 2020 00:03:54 -0400 Subject: [PATCH 2/4] get_variable -> core.scalar_variable --- Composite-Methods/G1.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/Composite-Methods/G1.py b/Composite-Methods/G1.py index 628536d1..54ae7b1e 100644 --- a/Composite-Methods/G1.py +++ b/Composite-Methods/G1.py @@ -40,11 +40,11 @@ psi4.frequencies("HF") T = psi4.core.get_global_option("T") -zpe = psi4.get_variable("ZPVE") -# du = psi4.get_variable("INTERNAL ENERGY CORRECTION") -du = psi4.get_variable("THERMAL ENERGY CORRECTION") -dh = psi4.get_variable("ENTHALPY CORRECTION") -dg = psi4.get_variable("GIBBS FREE ENERGY CORRECTION") +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() From 3c26be77fc06ba78fb7a913ce18afd6a267d48cb Mon Sep 17 00:00:00 2001 From: Eric Berquist Date: Thu, 15 Oct 2020 00:06:33 -0400 Subject: [PATCH 3/4] start testing composite methods under CI --- tests/test_RI_composite.py | 5 +++++ 1 file changed, 5 insertions(+) create mode 100644 tests/test_RI_composite.py diff --git a/tests/test_RI_composite.py b/tests/test_RI_composite.py new file mode 100644 index 00000000..46897cef --- /dev/null +++ b/tests/test_RI_composite.py @@ -0,0 +1,5 @@ +from utils import exe_py + + +def test_G1(workspace): + exe_py(workspace, tdir, 'G1') From 6c919c84cc4c261e01a23a528ef178525cb2a75d Mon Sep 17 00:00:00 2001 From: Eric Berquist Date: Tue, 20 Oct 2020 00:07:35 -0400 Subject: [PATCH 4/4] Add G2 and G2(MP2) --- Composite-Methods/{G1.py => Gaussian_n.py} | 62 ++++++++++++++++++++++ tests/test_RI_composite.py | 4 +- 2 files changed, 64 insertions(+), 2 deletions(-) rename Composite-Methods/{G1.py => Gaussian_n.py} (64%) diff --git a/Composite-Methods/G1.py b/Composite-Methods/Gaussian_n.py similarity index 64% rename from Composite-Methods/G1.py rename to Composite-Methods/Gaussian_n.py index 54ae7b1e..a1657030 100644 --- a/Composite-Methods/G1.py +++ b/Composite-Methods/Gaussian_n.py @@ -148,3 +148,65 @@ 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 index 46897cef..f11a8a0d 100644 --- a/tests/test_RI_composite.py +++ b/tests/test_RI_composite.py @@ -1,5 +1,5 @@ from utils import exe_py -def test_G1(workspace): - exe_py(workspace, tdir, 'G1') +def test_Gaussian_n(workspace): + exe_py(workspace, tdir, 'Gaussian_n')