From 2dd1769c34cd2a96c62a8abedc42eacf83726020 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rapha=C3=ABl=20Robidas?= Date: Sun, 6 Nov 2022 15:07:37 -0500 Subject: [PATCH] Added support for special basis sets in Gaussian (like SDD) --- ccinput/constants.py | 2 +- ccinput/packages/gaussian.py | 46 ++++++---- ccinput/tests/test_gaussian.py | 148 +++++++++++++++++++++++++++++++++ 3 files changed, 178 insertions(+), 18 deletions(-) diff --git a/ccinput/constants.py b/ccinput/constants.py index 0eb840f..8a75204 100644 --- a/ccinput/constants.py +++ b/ccinput/constants.py @@ -850,7 +850,7 @@ class CalcType(Enum): "cep4g": [], "cep31g": [], "cep121g": [], - "sdd": [], + "sdd": [], # D95, also named DZ or DZP (Dunning-Hay) in the BSE up to Ar, then "Stuttgart/Dresden ECPs" (called Stuttgart RLC in the BSE) for the rest of the periodic table "sddall": [], "dgdzvp": [], "dgdzvp2": [], diff --git a/ccinput/packages/gaussian.py b/ccinput/packages/gaussian.py index 3d8ae08..19a8279 100644 --- a/ccinput/packages/gaussian.py +++ b/ccinput/packages/gaussian.py @@ -233,23 +233,35 @@ def parse_custom_basis_set(self, base_bs): f"Invalid atom in custom basis set string: '{el}'" ) - bs = basis_set_exchange.get_basis( - bs_keyword, fmt="gaussian94", elements=[el_num], header=False - ) - if bs.find("-ECP") != -1: - ecp = True - sbs = bs.split("\n") - ecp_ind = -1 - for ind, line in enumerate(sbs): - if sbs[ind].find("-ECP") != -1: - ecp_ind = ind - break - bs_gen = "\n".join(sbs[: ecp_ind - 2]) + "\n" - bs_ecp = "\n".join(sbs[ecp_ind - 2 :]) - to_append_gen.append(bs_gen) - to_append_ecp.append(bs_ecp) + try: + bs = basis_set_exchange.get_basis( + bs_keyword, fmt="gaussian94", elements=[el_num], header=False + ) + except: + # Some basis sets are built-in, but use different names as the BSE does (e.g., SDD) + # In this case, just feed the user keyword in and hope it works. + # The basis set string has been recognized by ccinput, so it should exist in the program. + # ECP is added if Z > 18 (Ar) + to_append_gen.append(f"{el} 0\n{bs_keyword}\n****\n") + if el_num > 18: + ecp = True + to_append_ecp.append(f"{el} 0\n{bs_keyword}\n") + else: - to_append_gen.append(bs) + if bs.find("-ECP") != -1: + ecp = True + sbs = bs.split("\n") + ecp_ind = -1 + for ind, line in enumerate(sbs): + if sbs[ind].find("-ECP") != -1: + ecp_ind = ind + break + bs_gen = "\n".join(sbs[: ecp_ind - 2]) + "\n" + bs_ecp = "\n".join(sbs[ecp_ind - 2 :]) + to_append_gen.append(bs_gen) + to_append_ecp.append(bs_ecp) + else: + to_append_gen.append(bs) if len(custom_atoms) > 0: if ecp: @@ -264,7 +276,7 @@ def parse_custom_basis_set(self, base_bs): custom_bs += base_bs + "\n" custom_bs += "****\n" - custom_bs += "".join(to_append_gen) + custom_bs += "".join(to_append_gen) + "\n" custom_bs += "".join(to_append_ecp).replace("\n\n", "\n") return gen_keyword, custom_bs diff --git a/ccinput/tests/test_gaussian.py b/ccinput/tests/test_gaussian.py index 90ed74c..135cb94 100644 --- a/ccinput/tests/test_gaussian.py +++ b/ccinput/tests/test_gaussian.py @@ -2267,6 +2267,154 @@ def test_multiple_ecp(self): self.assertTrue(self.is_equivalent(REF, inp.input_file)) + def test_builtin_custom_bs(self): + params = { + "nproc": 8, + "mem": "10000MB", + "type": "Geometrical Optimisation", + "file": "ethanol.xyz", + "software": "Gaussian", + "charge": "0", + "method": "B3LYP", + "basis_set": "6-31+G(d,p)", + "custom_basis_sets": "O=SDD;", + } + + inp = self.generate_calculation(**params) + + REF = """ + %chk=ethanol.chk + %nproc=8 + %mem=10000MB + #p opt B3LYP/Gen + + File created by ccinput + + 0 1 + C -1.31970 -0.64380 0.00000 + H -0.96310 -1.65260 0.00000 + H -0.96310 -0.13940 -0.87370 + H -2.38970 -0.64380 0.00000 + C -0.80640 0.08220 1.25740 + H -1.16150 1.09160 1.25640 + H -1.16470 -0.42110 2.13110 + O 0.62360 0.07990 1.25870 + H 0.94410 0.53240 2.04240 + + C H 0 + 6-31+G(d,p) + **** + O 0 + sdd + **** + + """ + + self.assertTrue(self.is_equivalent(REF, inp.input_file)) + + def test_builtin_genecp_bs(self): + params = { + "nproc": 8, + "mem": "10000MB", + "type": "Geometrical Optimisation", + "file": "Ph2I_cation.xyz", + "software": "Gaussian", + "charge": "+1", + "method": "B3LYP", + "basis_set": "6-31+G(d,p)", + "custom_basis_sets": "I=SDD;", + } + + inp = self.generate_calculation(**params) + + REF = """ + %chk=Ph2I_cation.chk + %nproc=8 + %mem=10000MB + #p opt B3LYP/GenECP + + File created by ccinput + + 1 1 + C -3.06870 -2.28540 0.00000 + C -1.67350 -2.28540 0.00000 + C -0.97600 -1.07770 0.00000 + C -1.67360 0.13090 -0.00120 + C -3.06850 0.13080 -0.00170 + C -3.76610 -1.07740 -0.00070 + H -3.61840 -3.23770 0.00040 + H -1.12400 -3.23790 0.00130 + H 0.12370 -1.07760 0.00060 + H -1.12340 1.08300 -0.00130 + H -4.86570 -1.07720 -0.00090 + I -4.11890 1.94920 -0.00350 + C -4.64360 2.85690 -1.82310 + C -3.77180 3.76300 -2.42740 + C -5.86360 2.55380 -2.42750 + C -4.12020 4.36650 -3.63560 + H -2.81040 4.00240 -1.95030 + C -6.21180 3.15650 -3.63650 + H -6.55070 1.83950 -1.95140 + C -5.34050 4.06290 -4.24060 + H -3.43340 5.08120 -4.11170 + H -7.17360 2.91710 -4.11310 + H -5.61500 4.53870 -5.19320 + + C H 0 + 6-31+G(d,p) + **** + I 0 + sdd + **** + + I 0 + sdd + """ + + self.assertTrue(self.is_equivalent(REF, inp.input_file)) + + def test_multiple_builtin_genecp(self): + params = { + "nproc": 8, + "mem": "10000MB", + "type": "Geometrical Optimisation", + "file": "AuI.xyz", + "software": "Gaussian", + "charge": "0", + "method": "B3LYP", + "basis_set": "6-31+G(d,p)", + "custom_basis_sets": "I=SDD;Au=SDD;", + } + + inp = self.generate_calculation(**params) + + REF = """ + %chk=AuI.chk + %nproc=8 + %mem=10000MB + #p opt B3LYP/GenECP + + File created by ccinput + + 0 1 + Au -9.27600 -1.06330 0.00000 + I -6.60600 -1.06330 0.00000 + + I 0 + sdd + **** + Au 0 + sdd + **** + + I 0 + sdd + Au 0 + sdd + """ + + self.assertTrue(self.is_equivalent(REF, inp.input_file)) + def test_global_specification(self): params = { "nproc": 8,