Handle ghost atoms and parse NMR output (#41)

* added try except to_string VS to_str

* modified to allow nics calculations

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* modify a bit; add example

* bug fix nmr_parsing return None

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Kristjan Eimre <[email protected]>

aiida_gaussian/calculations/gaussian.py

@@ -171,6 +171,13 @@ def prepare_for_submission(self, folder):
             self.inputs.parameters.get_dict(), pmg_structure
         )
 
+        # 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:
             out_file.write(input_string)
 

aiida_gaussian/parsers/gaussian.py

@@ -185,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))
 
@@ -199,6 +201,60 @@ def _parse_log(self, log_file_string, inputs):
 
         return None
 
+    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 = []
+
+        def extract_values(line):
+            parts = line.split()
+            return np.array([parts[1], parts[3], parts[5]], dtype=float)
+
+        lines = log_file_string.splitlines()
+        i_line = 0
+        sigma_block = False
+        while i_line < len(lines):
+            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 {"nmr_tensors": np.array(sigma)}
+
     def _parse_log_spin_exp(self, log_file_string):
         """Parse spin expectation values"""
 

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 import orm
+from aiida.common import NotExistent
+from aiida.engine import run_get_node
+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 = orm.StructureData(ase=ase.io.read("./naphthalene_nics.xyz"))
+
+    num_cores = 1
+    memory_mb = 300
+
+    parameters = orm.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 = orm.load_code(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

examples/naphthalene_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