Changed how ORCA specifications work

Author: RaphaelRobidas

ccinput/packages/gaussian.py

@@ -12,6 +12,7 @@
     get_npxyz,
     check_fragments,
     add_fragments_xyz,
+    parse_specifications,
 )
 from ccinput.constants import CalcType, ATOMIC_NUMBER, LOWERCASE_ATOMIC_SYMBOLS
 from ccinput.exceptions import InvalidParameter, ImpossibleCalculation
@@ -104,39 +105,7 @@ def add_commands(self, commands):
     def handle_specifications(self):
         s = self.clean(self.calc.parameters.specifications.lower())
 
-        # Could be more sophisticated to catch other incorrect specifications
-        if s.count("(") != s.count(")"):
-            raise InvalidParameter("Invalid specifications: parenthesis not matching")
-
-        _specifications = ""
-        remove = False
-        for c in s:
-            if c == " " and remove:
-                continue
-            _specifications += c
-            if c == "(":
-                remove = True
-            elif c == ")":
-                remove = False
-
-        for spec in _specifications.split(" "):
-            if spec.strip() == "":
-                continue
-            if spec.find("(") != -1:
-                key, options = spec.split("(")
-                options = options.replace(")", "")
-                for option in options.split(","):
-                    if option.strip() != "":
-                        self.add_option(key, option)
-            elif spec.find("=") != -1:
-                try:
-                    key, option = spec.split("=")
-                except ValueError:
-                    raise InvalidParameter(f"Invalid specification: {spec}")
-                self.add_option(key, option)
-            else:
-                self.add_option(spec, "")
-
+        parse_specifications(s, self.add_option)
         if self.calc.parameters.d3:
             self.add_option("EmpiricalDispersion", "GD3")
         elif self.calc.parameters.d3bj:

ccinput/packages/orca.py

@@ -13,6 +13,7 @@
     get_abs_basis_set,
     clean_xyz,
     warn,
+    parse_specifications,
 )
 from ccinput.exceptions import (
     InvalidParameter,
@@ -69,6 +70,7 @@ def __init__(self, calc):
         self.specifications = {}
         self.solvation_radii = {}
         self.aux_basis_sets = {}
+        self.specifications_list = []
 
         if self.calc.type not in self.CALC_TYPES:
             raise ImpossibleCalculation(
@@ -98,7 +100,7 @@ def confirmed_specifications(self):
 
     def clean(self, s):
         WHITELIST = set(
-            "0123456789abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ/()=-,. "
+            "0123456789abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ/()=-,. {}_[]"
         )
         return "".join([c for c in s if c in WHITELIST])
 
@@ -108,12 +110,27 @@ def add_to_block(self, block, lines):
         else:
             self.blocks[block] = lines
 
+    def add_option(self, key, option):
+        if option == "":
+            if key[-2:] == "/c":
+                self.aux_basis_sets["C"] = get_basis_set(key[:-2], "orca")
+            elif key[-3:] == "/jk":
+                self.aux_basis_sets["JK"] = get_basis_set(key[:-3], "orca")
+            elif key[-2:] == "/j":
+                self.aux_basis_sets["J"] = get_basis_set(key[:-2], "orca")
+            elif key not in self.specifications_list:
+                self.specifications_list.append(key)
+        else:
+            self.add_to_block(key, [option.replace("=", " ")])
+
     def handle_specifications(self):
         _specifications = (
             self.clean(self.calc.parameters.specifications).lower().strip()
         )
 
-        specifications_list = []
+        parse_specifications(_specifications, self.add_option, condense=False)
+
+        """
         if _specifications != "":
             sspecs = _specifications.split()
             ind = 0
@@ -129,24 +146,15 @@ def handle_specifications(self):
                         raise InvalidParameter("Invalid specifications")
                     self.specifications["nimages"] = nimages
                     ind += 1
-                elif spec[-2:] == "/c":
-                    self.aux_basis_sets["C"] = get_basis_set(spec[:-2], "orca")
-                elif spec[-3:] == "/jk":
-                    self.aux_basis_sets["JK"] = get_basis_set(spec[:-3], "orca")
-                elif spec[-2:] == "/j":
-                    self.aux_basis_sets["J"] = get_basis_set(spec[:-2], "orca")
-                elif spec not in specifications_list:
-                    specifications_list.append(spec)
+            
 
                 ind += 1
+        """
 
         if self.calc.parameters.d3:
-            specifications_list.append("d3zero")
+            self.specifications_list.append("d3zero")
         elif self.calc.parameters.d3bj:
-            specifications_list.append("d3bj")
-
-        if len(specifications_list) > 0:
-            self.additional_commands = " ".join(specifications_list)
+            self.specifications_list.append("d3bj")
 
     def handle_command(self):
         if self.calc.type == CalcType.NMR:
@@ -244,13 +252,15 @@ def handle_command(self):
             self.command_line = "SP "
         elif self.calc.type == CalcType.MEP:  #### Second structure to handle
             self.command_line = "NEB "
-            if "nimages" in self.specifications:
-                nimages = self.specifications["nimages"]
+            self.add_to_block("neb", [f'product "{self.calc.aux_name}.xyz"'])
+            if "neb" in self.blocks:
+                for entry in self.blocks["neb"]:
+                    if entry.find("nimages") != -1:
+                        break
+                else:
+                    self.add_to_block("neb", ["nimages 8"])
             else:
-                nimages = 8
-            self.add_to_block(
-                "neb", [f'product "{self.calc.aux_name}.xyz"', f"nimages {nimages}"]
-            )
+                self.add_to_block("neb", ["nimages 8"])
 
         method = get_method(self.calc.parameters.method, "orca")
         if self.calc.parameters.theory_level not in [
@@ -446,7 +456,11 @@ def create_input_file(self):
 
         self.block_lines += f"""%MaxCore {self.mem_per_core}"""
 
+        if len(self.specifications_list) > 0:
+            self.additional_commands = " ".join(self.specifications_list)
+
         cmd = f"{self.command_line} {self.additional_commands}".replace("  ", " ")
+
         raw = self.TEMPLATE.format(
             cmd,
             self.calc.charge,

ccinput/tests/test_orca.py

@@ -414,13 +414,13 @@ def test_sp_DFT_multiple_specifications(self):
             "method": "M06-2X",
             "basis_set": "Def2-SVP",
             "charge": "-1",
-            "specifications": "TightSCF GRID6",
+            "specifications": "TightSCF DEFGRID3",
         }
 
         inp = self.generate_calculation(**params)
 
         REF = """
-        !SP M062X Def2-SVP tightscf grid6
+        !SP M062X Def2-SVP tightscf defgrid3
         *xyz -1 1
         Cl 0.0 0.0 0.0
         *
@@ -441,13 +441,13 @@ def test_sp_DFT_duplicate_specifications(self):
             "method": "M06-2X",
             "basis_set": "Def2-SVP",
             "charge": "-1",
-            "specifications": "tightscf TightSCF GRID6",
+            "specifications": "tightscf TightSCF DEFGRID3",
         }
 
         inp = self.generate_calculation(**params)
 
         REF = """
-        !SP M062X Def2-SVP tightscf grid6
+        !SP M062X Def2-SVP tightscf defgrid3
         *xyz -1 1
         Cl 0.0 0.0 0.0
         *
@@ -1372,6 +1372,78 @@ def test_ts_DFT_custom_bs_synonym(self):
 
         self.assertTrue(self.is_equivalent(REF, inp.input_file))
 
+    def test_ts_bond_following(self):
+        params = {
+            "nproc": 8,
+            "type": "TS Optimisation",
+            "file": "mini_ts.xyz",
+            "software": "ORCA",
+            "charge": "0",
+            "method": "B3LYP",
+            "basis_set": "6-31+G(d,p)",
+            "specifications": "geom(TS_Mode {B 0 108 } end)",
+        }
+
+        inp = self.generate_calculation(**params)
+
+        # One cannot calculate the Hessian for use in a TS optimization
+        # when using xtb as QM engine.
+        REF = """
+        !OPTTS B3LYP 6-31+G(d,p)
+        *xyz 0 1
+        N   1.08764072053386     -0.33994563112543     -0.00972525479568
+        H   1.99826836912112      0.05502842705407      0.00651240826058
+        H   0.59453997172323     -0.48560162159600      0.83949232123172
+        H   0.66998093862168     -0.58930117433261     -0.87511947469677
+        *
+        %geom
+        ts_mode {b 0 108 } end
+        Calc_Hess true
+        end
+        %pal
+        nprocs 8
+        end
+        %MaxCore 125
+        """
+
+        self.assertTrue(self.is_equivalent(REF, inp.input_file))
+
+    def test_ts_bond_following_alt(self):
+        params = {
+            "nproc": 8,
+            "type": "TS Optimisation",
+            "file": "mini_ts.xyz",
+            "software": "ORCA",
+            "charge": "0",
+            "method": "B3LYP",
+            "basis_set": "6-31+G(d,p)",
+            "specifications": "geom(TS_Mode={B 0 108 } end)",
+        }
+
+        inp = self.generate_calculation(**params)
+
+        # One cannot calculate the Hessian for use in a TS optimization
+        # when using xtb as QM engine.
+        REF = """
+        !OPTTS B3LYP 6-31+G(d,p)
+        *xyz 0 1
+        N   1.08764072053386     -0.33994563112543     -0.00972525479568
+        H   1.99826836912112      0.05502842705407      0.00651240826058
+        H   0.59453997172323     -0.48560162159600      0.83949232123172
+        H   0.66998093862168     -0.58930117433261     -0.87511947469677
+        *
+        %geom
+        ts_mode {b 0 108 } end
+        Calc_Hess true
+        end
+        %pal
+        nprocs 8
+        end
+        %MaxCore 125
+        """
+
+        self.assertTrue(self.is_equivalent(REF, inp.input_file))
+
     def test_opt_DFT_custom_bs_ecp(self):
         params = {
             "nproc": 8,
@@ -1677,7 +1749,7 @@ def test_NEB2(self):
             "file": "elimination_substrate.xyz",
             "auxiliary_file": "elimination_product.xyz",
             "software": "ORCA",
-            "specifications": "--nimages 12",
+            "specifications": "neb(nimages=12)",
             "charge": -1,
             "method": "gfn2-xtb",
         }
@@ -1698,8 +1770,8 @@ def test_NEB2(self):
         H         -1.82448        0.94856        3.28105
         *
         %neb
-        product "calc2.xyz"
         nimages 12
+        product "calc2.xyz"
         end
         %pal
         nprocs 8
@@ -1715,7 +1787,7 @@ def test_NEB_aux_name(self):
             "file": "elimination_substrate.xyz",
             "auxiliary_file": "elimination_product.xyz",
             "software": "ORCA",
-            "specifications": "--nimages 12",
+            "specifications": "neb(nimages=12)",
             "charge": -1,
             "method": "gfn2-xtb",
             "aux_name": "product",
@@ -1737,8 +1809,8 @@ def test_NEB_aux_name(self):
         H         -1.82448        0.94856        3.28105
         *
         %neb
-        product "product.xyz"
         nimages 12
+        product "product.xyz"
         end
         %pal
         nprocs 8
@@ -1756,7 +1828,7 @@ def test_hirshfeld_pop(self):
             "method": "M06-2X",
             "basis_set": "Def2-SVP",
             "charge": "-1",
-            "specifications": "--phirshfeld",
+            "specifications": "output(Print[ P_Hirshfeld] 1)",
         }
 
         inp = self.generate_calculation(**params)
@@ -1766,7 +1838,7 @@ def test_hirshfeld_pop(self):
         Cl 0.0 0.0 0.0
         *
         %output
-        Print[ P_Hirshfeld] 1
+        print[ p_hirshfeld] 1
         end
         %pal
         nprocs 8

ccinput/utilities.py

@@ -480,3 +480,40 @@ def check_fragments(counterpoise, fragments, xyz):
         raise InvalidParameter(
             "Counterpoise keyword must be equal to the total number of fragments"
         )
+
+
+def parse_specifications(specs, add_option_fn, condense=True):
+    # Could be more sophisticated to catch other incorrect specifications
+    if specs.count("(") != specs.count(")"):
+        raise InvalidParameter("Invalid specifications: parenthesis not matching")
+
+    _specifications = ""
+    remove = False
+    for c in specs:
+        if c == " " and remove:
+            if not condense:
+                _specifications += "&"  # Replace spaces temporarily
+            continue
+        _specifications += c
+        if c == "(":
+            remove = True
+        elif c == ")":
+            remove = False
+
+    for spec in _specifications.split(" "):
+        if spec.strip() == "":
+            continue
+        if spec.find("(") != -1:
+            key, options = spec.split("(")
+            options = options.replace(")", "")
+            for option in options.split(","):
+                if option.strip() != "":
+                    add_option_fn(key, option.replace("&", " "))
+        elif spec.find("=") != -1:
+            try:
+                key, option = spec.split("=")
+            except ValueError:
+                raise InvalidParameter(f"Invalid specification: {spec}")
+            add_option_fn(key, option.replace("&", " "))
+        else:
+            add_option_fn(spec, "")