From 2806dcd61eb0417ca495ce37d3880ee09e6dfea0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rapha=C3=ABl=20Robidas?= Date: Thu, 29 Aug 2024 11:15:24 -0400 Subject: [PATCH] Changed how ORCA specifications work --- ccinput/packages/gaussian.py | 35 +------------- ccinput/packages/orca.py | 56 ++++++++++++++-------- ccinput/tests/test_orca.py | 92 ++++++++++++++++++++++++++++++++---- ccinput/utilities.py | 37 +++++++++++++++ 4 files changed, 156 insertions(+), 64 deletions(-) diff --git a/ccinput/packages/gaussian.py b/ccinput/packages/gaussian.py index 29000c4..efd4103 100644 --- a/ccinput/packages/gaussian.py +++ b/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: diff --git a/ccinput/packages/orca.py b/ccinput/packages/orca.py index d8a31a6..fd8aea1 100644 --- a/ccinput/packages/orca.py +++ b/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, diff --git a/ccinput/tests/test_orca.py b/ccinput/tests/test_orca.py index 37f5f3d..4215a52 100644 --- a/ccinput/tests/test_orca.py +++ b/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 diff --git a/ccinput/utilities.py b/ccinput/utilities.py index 2965b71..b5552a0 100644 --- a/ccinput/utilities.py +++ b/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, "")