modify a bit; add example

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:

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"""

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

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