From ef3cc1a5b095618ca6b873abf3c27dcc10868624 Mon Sep 17 00:00:00 2001 From: Kristjan Eimre Date: Mon, 23 Sep 2024 16:29:57 +0300 Subject: [PATCH] modify a bit; add example --- aiida_gaussian/calculations/gaussian.py | 8 ++- aiida_gaussian/parsers/gaussian.py | 52 +++++++++++--- examples/example_04_nmr_nics.py | 90 +++++++++++++++++++++++++ examples/napthalene_nics.xyz | 22 ++++++ 4 files changed, 159 insertions(+), 13 deletions(-) create mode 100644 examples/example_04_nmr_nics.py create mode 100644 examples/napthalene_nics.xyz diff --git a/aiida_gaussian/calculations/gaussian.py b/aiida_gaussian/calculations/gaussian.py index 61114f9..04fea63 100644 --- a/aiida_gaussian/calculations/gaussian.py +++ b/aiida_gaussian/calculations/gaussian.py @@ -170,8 +170,12 @@ def prepare_for_submission(self, folder): input_string = GaussianCalculation._render_input_string_from_params( self.inputs.parameters.get_dict(), pmg_structure ) - # for NICS replace X0+ with Bq - if "nmr" in self.inputs.parameters.get_dict()["route_parameters"]: + + # Handle Ghost atoms (e.g. when doing NICS calculations) + # Atoms with symbol 'X' in input structure are considered Ghost atoms. + # Pymatgen converts these into `X0+` in the input script + # and Gaussian needs these to be called `Bq`. + if pmg_structure and "X0+" in pmg_structure.labels: input_string = input_string.replace("X0+", "Bq") with open(folder.get_abs_path(self.INPUT_FILE), "w") as out_file: diff --git a/aiida_gaussian/parsers/gaussian.py b/aiida_gaussian/parsers/gaussian.py index 36dabb9..842f885 100644 --- a/aiida_gaussian/parsers/gaussian.py +++ b/aiida_gaussian/parsers/gaussian.py @@ -176,8 +176,6 @@ def _parse_log(self, log_file_string, inputs): if property_dict is None: return self.exit_codes.ERROR_OUTPUT_PARSING - # parse nics - property_dict.update(self._parse_nics(log_file_string)) property_dict.update(self._parse_electron_numbers(log_file_string)) @@ -187,6 +185,8 @@ def _parse_log(self, log_file_string, inputs): # separate HOMO-LUMO gap as its own entry in property_dict self._extract_homo_lumo_gap(property_dict) + property_dict.update(self._parse_nmr(log_file_string)) + # set output nodes self.out("output_parameters", Dict(dict=property_dict)) @@ -201,7 +201,23 @@ def _parse_log(self, log_file_string, inputs): return None - def _parse_nics(self, log_file_string): + def _parse_nmr(self, log_file_string): + """Parse the NMR magnetic shielding tensors. + + Example log: + + SCF GIAO Magnetic shielding tensor (ppm): + 1 C Isotropic = 64.6645 Anisotropy = 153.1429 + XX= 2.6404 YX= 36.9712 ZX= -0.0000 + XY= 39.8249 YY= 24.5934 ZY= -0.0000 + XZ= -0.0000 YZ= -0.0000 ZZ= 166.7598 + Eigenvalues: -26.3192 53.5530 166.7598 + 2 C Isotropic = 64.6641 Anisotropy = 153.1416 + ... + """ + + if "Magnetic shielding tensor" not in log_file_string: + return sigma = [] @@ -211,19 +227,33 @@ def extract_values(line): lines = log_file_string.splitlines() i_line = 0 + sigma_block = False while i_line < len(lines): - if "Bq Isotropic" in lines[i_line]: - s = np.zeros((3, 3)) - s[0] = extract_values(lines[i_line + 1]) - s[1] = extract_values(lines[i_line + 2]) - s[2] = extract_values(lines[i_line + 3]) - sigma.append(s) - i_line += 4 + if "Magnetic shielding tensor" in lines[i_line]: + sigma_block = True + + if sigma_block: + + if "Leave Link" in lines[i_line]: + sigma_block = False + + if "Isotropic" in lines[i_line] and "Anisotropy" in lines[i_line]: + s = np.zeros((3, 3)) + s[0] = extract_values(lines[i_line + 1]) + s[1] = extract_values(lines[i_line + 2]) + s[2] = extract_values(lines[i_line + 3]) + sigma.append(s) + i_line += 4 + else: + i_line += 1 + else: i_line += 1 + if len(sigma) == 0: return {} - return {"nics": np.array(sigma)} + + return {"nmr_tensors": np.array(sigma)} def _parse_log_spin_exp(self, log_file_string): """Parse spin expectation values""" diff --git a/examples/example_04_nmr_nics.py b/examples/example_04_nmr_nics.py new file mode 100644 index 0000000..466f207 --- /dev/null +++ b/examples/example_04_nmr_nics.py @@ -0,0 +1,90 @@ +# pylint: disable=invalid-name +"""Run simple DFT calculation""" + + +import sys + +import ase.io +import click +from aiida.common import NotExistent +from aiida.engine import run_get_node +from aiida.orm import Code, Dict, StructureData +from aiida.plugins import CalculationFactory + +GaussianCalculation = CalculationFactory("gaussian") + + +def example_nmr_nics(gaussian_code): + """Run an NMR calculation. + + The geometry also includes fake/ghost atoms (X), + where nucleus-independent chemical shift (NICS) is calculated. + """ + + # structure + structure = StructureData(ase=ase.io.read("./napthalene_nics.xyz")) + + num_cores = 1 + memory_mb = 300 + + parameters = Dict( + { + "link0_parameters": { + "%chk": "aiida.chk", + "%mem": "%dMB" % memory_mb, + "%nprocshared": num_cores, + }, + "functional": "BLYP", + "basis_set": "6-31g", + "charge": 0, + "multiplicity": 1, + "dieze_tag": "#P", + "route_parameters": { + "nmr": None, + }, + } + ) + + # Construct process builder + + builder = GaussianCalculation.get_builder() + + builder.structure = structure + builder.parameters = parameters + builder.code = gaussian_code + + builder.metadata.options.resources = { + "num_machines": 1, + "tot_num_mpiprocs": num_cores, + } + + # The advanced parser handles the NMR output + builder.metadata.options.parser_name = "gaussian.advanced" + + # Should ask for extra +25% extra memory + builder.metadata.options.max_memory_kb = int(1.25 * memory_mb) * 1024 + builder.metadata.options.max_wallclock_seconds = 5 * 60 + + print("Running calculation...") + res, _node = run_get_node(builder) + + print("NMR tensors of each atom:") + for i, site in enumerate(structure.sites): + print(site) + print(" ", res["output_parameters"]["nmr_tensors"][i]) + + +@click.command("cli") +@click.argument("codelabel", default="gaussian@localhost") +def cli(codelabel): + """Click interface""" + try: + code = Code.get_from_string(codelabel) + except NotExistent: + print(f"The code '{codelabel}' does not exist") + sys.exit(1) + example_nmr_nics(code) + + +if __name__ == "__main__": + cli() # pylint: disable=no-value-for-parameter diff --git a/examples/napthalene_nics.xyz b/examples/napthalene_nics.xyz new file mode 100644 index 0000000..39f1f1f --- /dev/null +++ b/examples/napthalene_nics.xyz @@ -0,0 +1,22 @@ +20 +Properties=species:S:1:pos:R:3 napthalene=T pbc="F F F" +C 0.00000000 1.42000000 0.00000000 +C 0.00000000 0.00000000 0.00000000 +C 1.22975600 2.13000000 0.00000000 +C 2.45951200 1.42000000 0.00000000 +C 2.45951200 0.00000000 0.00000000 +C 1.22975600 -0.71000000 0.00000000 +C 3.68926800 -0.71000000 0.00000000 +C 4.91902400 0.00000000 0.00000000 +C 4.91902400 1.42000000 0.00000000 +H -0.98726900 1.99000000 0.00000000 +H -0.98726900 -0.57000000 0.00000000 +H 1.22975600 3.27000000 0.00000000 +H 1.22975600 -1.85000000 0.00000000 +H 3.68926800 -1.85000000 0.00000000 +H 3.68926800 3.27000000 0.00000000 +H 5.90629300 1.99000000 0.00000000 +H 5.90629300 -0.57000000 0.00000000 +C 3.68926800 2.13000000 0.00000000 +X 1.22975600 0.73000000 1.00000000 +X 3.68926800 0.73000000 1.00000000