Added support for special basis sets in Gaussian (like SDD)

cyllab · Nov 6, 2022 · 2dd1769 · 2dd1769
1 parent 0e3124b
commit 2dd1769
Show file tree

Hide file tree

Showing 3 changed files with 178 additions and 18 deletions.
diff --git 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
@@ -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
@@ -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,