From fb7525eb14c6c93e7ba0563d5e6129d489a4ccbb Mon Sep 17 00:00:00 2001 From: Carlo Pignedoli Date: Thu, 21 Nov 2024 23:01:08 +0100 Subject: [PATCH] Nics (#166) * work in progress * not much done * working * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * added extras * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * working * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * added extras to input structure --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> --- aiida_nanotech_empa/utils/cycle_tools.py | 119 +++++++++++ aiida_nanotech_empa/utils/nmr.py | 163 +++++++++++++++ .../workflows/gaussian/__init__.py | 2 + .../gaussian/casscf_series_workchain.py | 4 +- .../workflows/gaussian/nics_workchain.py | 185 ++++++++++++++++++ .../workflows/gaussian/relax_workchain.py | 38 +++- .../workflows/gaussian/scf_workchain.py | 50 ++++- examples/workflows/benzene.xyz | 14 ++ examples/workflows/example_gaussian_nics.py | 68 +++++++ examples/workflows/example_gaussian_opt.py | 55 ++++++ examples/workflows/naphthalene.xyz | 20 ++ pyproject.toml | 1 + 12 files changed, 715 insertions(+), 4 deletions(-) create mode 100644 aiida_nanotech_empa/utils/cycle_tools.py create mode 100644 aiida_nanotech_empa/utils/nmr.py create mode 100644 aiida_nanotech_empa/workflows/gaussian/nics_workchain.py create mode 100644 examples/workflows/benzene.xyz create mode 100644 examples/workflows/example_gaussian_nics.py create mode 100644 examples/workflows/example_gaussian_opt.py create mode 100644 examples/workflows/naphthalene.xyz diff --git a/aiida_nanotech_empa/utils/cycle_tools.py b/aiida_nanotech_empa/utils/cycle_tools.py new file mode 100644 index 0000000..feaf9ca --- /dev/null +++ b/aiida_nanotech_empa/utils/cycle_tools.py @@ -0,0 +1,119 @@ +import ase +import ase.io +import ase.neighborlist +import ase.visualize +import numpy as np +import scipy as sp +import scipy.linalg + + +def convert_neighbor_list(nl): + new = {} + n_vert = np.max(nl) + 1 + + for i_v in range(n_vert): + new[i_v] = [] + + for i_v, j_v in zip(nl[0], nl[1]): + + new[i_v].append(j_v) + + return new + + +def find_cycles(i_vert, cnl, max_length, cur_path, passed_edges): + + if len(cur_path) - 1 == max_length: + return [] + + acc_cycles = [] + sort_cycles = [] + + neighbs = cnl[i_vert] + + # if we are connected to something that is not the end + # then we crossed multiple cycles + for n in neighbs: + edge = (np.min([i_vert, n]), np.max([i_vert, n])) + if edge not in passed_edges: + + if n in cur_path[1:]: + # path went too close to itself... + return [] + + # CHeck if we are at the end + for n in neighbs: + edge = (np.min([i_vert, n]), np.max([i_vert, n])) + if edge not in passed_edges: + + if n == cur_path[0]: + # found cycle + return [cur_path] + + # Continue in all possible directions + for n in neighbs: + edge = (np.min([i_vert, n]), np.max([i_vert, n])) + if edge not in passed_edges: + + cycs = find_cycles( + n, cnl, max_length, cur_path + [n], passed_edges + [edge] + ) + for cyc in cycs: + sorted_cyc = tuple(sorted(cyc)) + if sorted_cyc not in sort_cycles: + sort_cycles.append(sorted_cyc) + acc_cycles.append(cyc) + + return acc_cycles + + +def dumb_cycle_detection(ase_atoms_no_h, max_length): + + neighbor_list = ase.neighborlist.neighbor_list("ij", ase_atoms_no_h, 2.0) + + cycles = [] + sorted_cycles = [] + n_vert = np.max(neighbor_list) + 1 + + cnl = convert_neighbor_list(neighbor_list) + + for i_vert in range(n_vert): + + cycs = find_cycles(i_vert, cnl, max_length, [i_vert], []) + for cyc in cycs: + sorted_cyc = tuple(sorted(cyc)) + if sorted_cyc not in sorted_cycles: + sorted_cycles.append(sorted_cyc) + cycles.append(cyc) + + return cycles + + +def cycle_normal(cycle, h): + cycle = np.array(cycle) + centroid = np.mean(cycle, axis=0) + + points = cycle - centroid + u, s, v = np.linalg.svd(points.T) + normal = u[:, -1] + normal /= np.linalg.norm(normal) + if np.dot(normal, h * np.array([1, 1, 1])) < 0.0: + normal *= -1.0 + return normal + + +def find_cycle_centers_and_normals(ase_atoms_no_h, cycles, h=0.0): + """ + positive h means projection to z axis is positive and vice-versa + """ + if h == 0.0: + h = 1.0 + normals = [] + centers = [] + for cyc in cycles: + cyc_p = [] + for i_at in cyc: + cyc_p.append(ase_atoms_no_h[i_at].position) + normals.append(cycle_normal(cyc_p, h)) + centers.append(np.mean(cyc_p, axis=0)) + return np.array(centers), np.array(normals) diff --git a/aiida_nanotech_empa/utils/nmr.py b/aiida_nanotech_empa/utils/nmr.py new file mode 100644 index 0000000..6108d6b --- /dev/null +++ b/aiida_nanotech_empa/utils/nmr.py @@ -0,0 +1,163 @@ +import ase +import ase.io +import ase.neighborlist +import ase.visualize +import cclib +import numpy as np +import scipy as sp +import scipy.linalg + +from . import cycle_tools + +### ----------------------------------------------------------------- +### SETUP + + +def find_ref_points(ase_atoms_no_h, cycles, h=0.0): + """ + positive h means projection to z axis is positive and vice-versa + """ + + centers, normals = cycle_tools.find_cycle_centers_and_normals( + ase_atoms_no_h, cycles, h + ) + + ase_ref_p = ase.Atoms() + + for i_cyc in range(len(cycles)): + if h == 0.0: + ase_ref_p.append(ase.Atom("X", centers[i_cyc])) + else: + pos = centers[i_cyc] + np.abs(h) * normals[i_cyc] + ase_ref_p.append(ase.Atom("X", pos)) + + return ase_ref_p + + +def dist(p1, p2): + if isinstance(p1, ase.Atom): + p1 = p1.position + if isinstance(p2, ase.Atom): + p2 = p2.position + return np.linalg.norm(p2 - p1) + + +def interp_pts(p1, p2, dx): + vec = p2 - p1 + dist = np.sqrt(np.sum(vec**2)) + num_p = int(np.round(dist / dx)) + dx_real = dist / num_p + dvec = vec / dist * dx_real + + points = np.outer(np.arange(0, num_p), dvec) + p1 + return points + + +def build_path(ref_pts, dx=0.1): + + point_arr = None + + for i_rp in range(len(ref_pts) - 1): + + pt1 = ref_pts[i_rp].position + pt2 = ref_pts[i_rp + 1].position + + points = interp_pts(pt1, pt2, dx) + + if i_rp == len(ref_pts) - 2: + points = np.concatenate([points, [pt2]], axis=0) + + if point_arr is None: + point_arr = points + else: + point_arr = np.concatenate([point_arr, points], axis=0) + + ase_arr = ase.Atoms("X%d" % len(point_arr), point_arr) + return ase_arr + + +### ----------------------------------------------------------------- +### PROCESS + + +def load_nics_gaussian(nics_path): + + sigma = [] + + def extract_values(line): + parts = line.split() + return np.array([parts[1], parts[3], parts[5]], dtype=float) + + with open(nics_path) as file: + lines = file.readlines() + i_line = 0 + 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 + else: + i_line += 1 + + return np.array(sigma) + + +def is_number(x): + try: + float(x) + return True + except ValueError: + return False + + +def parse_nmr_cmo_matrix(log_file_str, property_dict): + + lines = log_file_str.splitlines() + + # build the object + n_atom = property_dict["natom"] + n_occupied_mo = property_dict["homos"][0] + 1 + + nmr_cmo_matrix = np.zeros((n_atom, n_occupied_mo, 3, 3)) + + i_line = 0 + while i_line < len(lines): + + # -------------------------------------------------------------------------- + # Full Cartesian NMR shielding tensor (ppm) for atom C( 1): + # Canonical MO contributions + # + # MO XX XY XZ YX YY YZ ZX ZY ZZ + # =============================================================================== + # 1. 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 -0.16 + # 2. 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 -0.13 + # 3. 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.24 + # 4. 0.00 0.00 0.00 0.00 -0.03 0.00 0.00 0.00 0.27 + # ... + + if "Full Cartesian NMR shielding tensor (ppm) for atom" in lines[i_line]: + + i_atom = int(lines[i_line].replace("(", ")").split(")")[-2]) - 1 + + i_line += 1 + + if "Canonical MO contributions" in lines[i_line]: + + for _i in range(2000): + i_line += 1 + if "Total" in lines[i_line]: + break + split = lines[i_line].split() + + if len(split) == 10 and is_number(split[-1]): + i_mo = int(split[0][:-1]) - 1 + + arr = np.array([float(x) for x in split[1:]]) + nmr_cmo_matrix[i_atom, i_mo, :, :] = arr.reshape(3, 3) + + i_line += 1 + + return nmr_cmo_matrix diff --git a/aiida_nanotech_empa/workflows/gaussian/__init__.py b/aiida_nanotech_empa/workflows/gaussian/__init__.py index d08b655..809a790 100644 --- a/aiida_nanotech_empa/workflows/gaussian/__init__.py +++ b/aiida_nanotech_empa/workflows/gaussian/__init__.py @@ -4,6 +4,7 @@ from .delta_scf_workchain import GaussianDeltaScfWorkChain from .hf_mp2_workchain import GaussianHfMp2WorkChain from .natorb_workchain import GaussianNatOrbWorkChain +from .nics_workchain import GaussianNicsWorkChain from .relax_workchain import GaussianRelaxWorkChain from .scf_workchain import GaussianScfWorkChain from .spin_workchain import GaussianSpinWorkChain @@ -18,4 +19,5 @@ "GaussianConstrOptChainWorkChain", "GaussianCasscfWorkChain", "GaussianCasscfSeriesWorkChain", + "GaussianNicsWorkChain", ) diff --git a/aiida_nanotech_empa/workflows/gaussian/casscf_series_workchain.py b/aiida_nanotech_empa/workflows/gaussian/casscf_series_workchain.py index 4aac31f..0baf0ec 100644 --- a/aiida_nanotech_empa/workflows/gaussian/casscf_series_workchain.py +++ b/aiida_nanotech_empa/workflows/gaussian/casscf_series_workchain.py @@ -267,5 +267,7 @@ def finalize(self): f"{var}_cube_image_folder", self.ctx[var].outputs.cube_image_folder, ) - + # Add extras. + struc = self.inputs.structure + common_utils.add_extras(struc, "surfaces", self.node.uuid) return engine.ExitCode(0) diff --git a/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py b/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py new file mode 100644 index 0000000..2486a19 --- /dev/null +++ b/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py @@ -0,0 +1,185 @@ +import ase +import numpy as np +from aiida import engine, orm + +from ...utils import common_utils, cycle_tools, nmr +from .relax_workchain import GaussianRelaxWorkChain +from .scf_workchain import GaussianScfWorkChain + + +@engine.calcfunction +def prepare_nmr_structure(structure, height): + ase_geo = structure.get_ase() + pos = ase_geo.positions + extents = np.array( + [ + np.max(pos[:, 0]) - np.min(pos[:, 0]), + np.max(pos[:, 1]) - np.min(pos[:, 1]), + np.max(pos[:, 2]) - np.min(pos[:, 2]), + ] + ) + inds = np.argsort(-extents) + new_pos = pos[:, inds] + ase_atoms = ase.Atoms(numbers=ase_geo.numbers, positions=new_pos) + ase_atoms_no_h = ase.Atoms([a for a in ase_atoms if a.symbol != "H"]) + cycles = cycle_tools.dumb_cycle_detection(ase_atoms_no_h, 8) + ref_p = nmr.find_ref_points(ase_atoms_no_h, cycles, height.value) + new_ase = ase_atoms + ref_p + return orm.StructureData(ase=new_ase) + + +class GaussianNicsWorkChain(engine.WorkChain): + @classmethod + def define(cls, spec): + super().define(spec) + + spec.input("gaussian_code", valid_type=orm.Code) + + spec.input( + "structure", + valid_type=orm.StructureData, + required=True, + help="input geometry", + ) + spec.input( + "opt", + valid_type=orm.Bool, + required=False, + help="False do not optimize structure", + ) + spec.input( + "height", + valid_type=orm.Float, + required=False, + help="Height of NICS centers", + ) + spec.input( + "functional", valid_type=orm.Str, required=True, help="xc functional" + ) + + spec.input("basis_set", valid_type=orm.Str, required=True, help="basis_set") + + spec.input( + "multiplicity", + valid_type=orm.Int, + required=False, + default=lambda: orm.Int(0), + help="spin multiplicity; 0 means RKS", + ) + spec.input( + "wfn_stable_opt", + valid_type=orm.Bool, + required=False, + default=lambda: orm.Bool(False), + help="if true, perform wfn stability optimization", + ) + spec.input( + "empirical_dispersion", + valid_type=orm.Str, + required=False, + default=lambda: orm.Str(""), + help=("Include empirical dispersion corrections" '(e.g. "GD3", "GD3BJ")'), + ) + spec.input( + "options", + valid_type=orm.Dict, + required=False, + help="Use custom metadata.options instead of the automatic ones.", + ) + + spec.outline( + cls.setup, + engine.if_(cls.should_submit_opt)(cls.submit_opt, cls.inspect_opt), + cls.submit_nics, + cls.inspect_nics, + cls.finalize, + ) + + spec.outputs.dynamic = True + + spec.exit_code( + 390, + "ERROR_TERMINATION", + message="One or more steps of the work chain failed.", + ) + + def setup(self): + self.report("Setting up...") + self.ctx.nmr_structure = self.inputs.structure + self.ctx_should_opt = getattr(self.inputs, "opt", True) + + def should_submit_opt(self): + return self.ctx_should_opt + + def submit_opt(self): + self.report("Submitting optimization...") + label = "geo_opt" + builder = GaussianRelaxWorkChain.get_builder() + builder.gaussian_code = self.inputs.gaussian_code + builder.structure = self.inputs.structure + builder.functional = self.inputs.functional + builder.empirical_dispersion = self.inputs.empirical_dispersion + builder.basis_set = self.inputs.basis_set + builder.multiplicity = self.inputs.multiplicity + builder.tight = orm.Bool(True) + builder.int = orm.Str("superfine") + builder.cdiis = orm.Bool(True) + builder.maxcycle = orm.Int(2048) + builder.conver = orm.Int(8) + if "options" in self.inputs: + builder.options = self.inputs.options + + submitted_node = self.submit(builder) + submitted_node.description = label + self.to_context(**{label: submitted_node}) + + def inspect_opt(self): + self.report("Inspecting optimization...") + label = "geo_opt" + # check if everything finished nicely + if not common_utils.check_if_calc_ok(self, self.ctx[label]): + return self.exit_codes.ERROR_TERMINATION + self.ctx.nmr_structure = self.ctx[label].outputs.output_structure + return engine.ExitCode(0) + + def submit_nics(self): + self.report("Submitting NICS calculation...") + label = "nics" + builder = GaussianScfWorkChain.get_builder() + height = getattr(self.inputs, "height", 1.0) + builder.gaussian_code = self.inputs.gaussian_code + builder.structure = prepare_nmr_structure( + self.ctx.nmr_structure, orm.Float(height) + ) + builder.functional = self.inputs.functional + builder.basis_set = self.inputs.basis_set + builder.multiplicity = self.inputs.multiplicity + builder.int = orm.Str("superfine") + builder.cdiis = orm.Bool(True) + builder.nmr = orm.Bool(True) + builder.maxcycle = orm.Int(2048) + builder.conver = orm.Int(8) + if "options" in self.inputs: + builder.options = self.inputs.options + + submitted_node = self.submit(builder) + submitted_node.description = label + self.to_context(**{label: submitted_node}) + + def inspect_nics(self): + self.report("Inspecting nics...") + label = "nics" + # check if everything finished nicely + if not common_utils.check_if_calc_ok(self, self.ctx[label]): + return self.exit_codes.ERROR_TERMINATION + + self.out("output_parameters", self.ctx[label].outputs.output_parameters) + self.out("output_structure", self.ctx.nics.inputs.structure) + return engine.ExitCode(0) + + def finalize(self): + self.report("Finalizing...") + + # Add extras. + struc = self.inputs.structure + common_utils.add_extras(struc, "surfaces", self.node.uuid) diff --git a/aiida_nanotech_empa/workflows/gaussian/relax_workchain.py b/aiida_nanotech_empa/workflows/gaussian/relax_workchain.py index 9c03573..fe6b210 100644 --- a/aiida_nanotech_empa/workflows/gaussian/relax_workchain.py +++ b/aiida_nanotech_empa/workflows/gaussian/relax_workchain.py @@ -49,7 +49,30 @@ def define(cls, spec): default=lambda: orm.Bool(False), help="Use tight optimization criteria.", ) - + spec.input( + "maxcycle", + valid_type=orm.Int, + required=False, + help="the maximum number of scf cycles", + ) + spec.input( + "cdiis", + valid_type=orm.Bool, + required=False, + help="Conjugate Direct Inversion in the Iterative Subspace", + ) + spec.input( + "conver", + valid_type=orm.Int, + required=False, + help="the scf convergece threshold", + ) + spec.input( + "int", + valid_type=orm.Str, + required=False, + help="the integral grid", + ) spec.input( "freq", valid_type=orm.Bool, @@ -290,7 +313,18 @@ def optimization(self): if len(opt_dict) != 0: parameters["route_parameters"]["opt"] = opt_dict - + conver = getattr(self.inputs, "conver", None) + maxcycle = getattr(self.inputs, "maxcycle", None) + grid = getattr(self.inputs, "int", None) + cdiis = getattr(self.inputs, "cdiis", None) + if grid: + parameters["route_parameters"]["int"] = grid + if maxcycle: + parameters["route_parameters"]["scf"]["maxcycle"] = maxcycle + if cdiis: + parameters["route_parameters"]["scf"]["cdiis"] = None + if conver: + parameters["route_parameters"]["scf"]["conver"] = conver builder.gaussian.parameters = parameters builder.gaussian.structure = self.inputs.structure builder.gaussian.code = self.inputs.gaussian_code diff --git a/aiida_nanotech_empa/workflows/gaussian/scf_workchain.py b/aiida_nanotech_empa/workflows/gaussian/scf_workchain.py index f38857a..e36adda 100644 --- a/aiida_nanotech_empa/workflows/gaussian/scf_workchain.py +++ b/aiida_nanotech_empa/workflows/gaussian/scf_workchain.py @@ -60,6 +60,36 @@ def define(cls, spec): required=False, help="the folder of a completed gaussian calculation", ) + spec.input( + "maxcycle", + valid_type=orm.Int, + required=False, + help="the maximum number of scf cycles", + ) + spec.input( + "conver", + valid_type=orm.Int, + required=False, + help="the scf convergece threshold", + ) + spec.input( + "cdiis", + valid_type=orm.Bool, + required=False, + help="Conjugate Direct Inversion in the Iterative Subspace", + ) + spec.input( + "int", + valid_type=orm.Str, + required=False, + help="the integral grid", + ) + spec.input( + "nmr", + valid_type=orm.Bool, + required=False, + help="nmr calculation", + ) # Inputs for cubes generation. spec.input("formchk_code", valid_type=orm.Code, required=False) @@ -167,7 +197,7 @@ def setup(self): # Use default convergence criterion at the start # but switch to conver=7 in case of failure. - self.ctx.conver = None + self.ctx.conver = getattr(self.inputs, "conver", None) self.ctx.scf_label = "scf" return engine.ExitCode(0) @@ -228,6 +258,24 @@ def scf(self): }, } ) + nmr = getattr(self.inputs, "nmr", None) + conver = getattr(self.inputs, "conver", None) + maxcycle = getattr(self.inputs, "maxcycle", None) + grid = getattr(self.inputs, "int", None) + cdiis = getattr(self.inputs, "cdiis", None) + if grid: + parameters["route_parameters"]["int"] = grid + if maxcycle: + parameters["route_parameters"]["scf"]["maxcycle"] = maxcycle + if cdiis: + parameters["route_parameters"]["scf"]["cdiis"] = None + if conver: + parameters["route_parameters"]["scf"]["conver"] = conver + if nmr: + if conver is None: + conver = 8 + parameters["route_parameters"]["nmr"] = None + parameters["route_parameters"]["cphf"] = {"conver": conver} if self.inputs.wfn_stable_opt: parameters["route_parameters"]["stable"] = "opt" else: diff --git a/examples/workflows/benzene.xyz b/examples/workflows/benzene.xyz new file mode 100644 index 0000000..5255bb3 --- /dev/null +++ b/examples/workflows/benzene.xyz @@ -0,0 +1,14 @@ +12 +benzene +C 8.39528138 6.66748365 5.00000291 +C 8.51925793 8.06102240 5.00000509 +C 7.12645219 6.07808103 5.00000441 +C 5.98160072 6.88221755 5.00000464 +C 6.10557756 8.27575669 5.00000369 +C 7.37440623 8.86515917 5.00000734 +H 9.50085998 8.51699940 5.00000031 +H 9.28097172 6.04538232 5.00000427 +H 7.03053864 5.00000000 5.00000036 +H 5.00000000 6.42623909 5.00000439 +H 5.21988789 8.89785938 5.00000000 +H 7.47031658 9.94324034 5.00000253 diff --git a/examples/workflows/example_gaussian_nics.py b/examples/workflows/example_gaussian_nics.py new file mode 100644 index 0000000..295741d --- /dev/null +++ b/examples/workflows/example_gaussian_nics.py @@ -0,0 +1,68 @@ +import pathlib + +import click +import numpy as np +from aiida import engine, orm +from aiida.plugins import WorkflowFactory +from ase.build import molecule +from ase.io import read + +import aiida_nanotech_empa.utils.gaussian_wcs_postprocess as pp + +GaussianNicsWorkChain = WorkflowFactory("nanotech_empa.gaussian.nics") +DATA_DIR = pathlib.Path(__file__).parent.absolute() +GEO_FILE = "naphthalene.xyz" + + +def _example_gaussian_nics(gaussian_code, opt): + # Check test geometry is already in database. + qb = orm.QueryBuilder() + qb.append(orm.Node, filters={"label": {"==": GEO_FILE}}) + structure = None + for node_tuple in qb.iterall(): + node = node_tuple[0] + structure = node + if structure is not None: + print("found existing structure: ", structure.pk) + else: + structure = orm.StructureData(ase=read(DATA_DIR / GEO_FILE)) + structure.label = GEO_FILE + structure.store() + print("created new structure: ", structure.pk) + + builder = GaussianNicsWorkChain.get_builder() + builder.gaussian_code = gaussian_code + builder.structure = structure + builder.functional = orm.Str("B3LYP") + builder.basis_set = orm.Str("6-311G(d,p)") + builder.opt = orm.Bool(opt) + builder.options = orm.Dict( + { + "resources": { + "tot_num_mpiprocs": 4, + "num_machines": 1, + }, + "max_wallclock_seconds": 1 * 60 * 60, + "max_memory_kb": 8 * 1024 * 1024, # GB + } + ) + + _, wc_node = engine.run_get_node(builder) + + assert wc_node.is_finished_ok + return wc_node.pk + + +@click.command("cli") +@click.argument("gaussian_code", default="gaussian@localhost") +def run_nics(gaussian_code): + # print("#### Running Gaussian NICS WorkChain with geo_opt ####") + # uuid = _example_gaussian_nics(orm.load_code(gaussian_code),True) + # print(f"WorkChain completed uuid: {uuid}") + print("#### Running Gaussian NICS WorkChain without geo_opt ####") + uuid = _example_gaussian_nics(orm.load_code(gaussian_code), True) + print(f"WorkChain completed uuid: {uuid}") + + +if __name__ == "__main__": + run_nics() diff --git a/examples/workflows/example_gaussian_opt.py b/examples/workflows/example_gaussian_opt.py new file mode 100644 index 0000000..7a9bb70 --- /dev/null +++ b/examples/workflows/example_gaussian_opt.py @@ -0,0 +1,55 @@ +import os + +import ase.io +import numpy as np +from aiida.engine import run_get_node +from aiida.orm import Bool, Int, Str, StructureData, load_code +from aiida.plugins import WorkflowFactory + +import aiida_nanotech_empa.utils.gaussian_wcs_postprocess as pp + +GaussianRelaxWorkChain = WorkflowFactory("nanotech_empa.gaussian.relax") + +DATA_DIR = os.path.dirname(os.path.abspath(__file__)) +OUTPUT_DIR = os.path.dirname(os.path.abspath(__file__)) + + +def _example_gaussian_spin(gaussian_code): # , formchk_code, cubegen_code): + ase_geom = ase.io.read(os.path.join(DATA_DIR, "benzene.xyz")) + ase_geom.cell = np.diag([10.0, 10.0, 10.0]) + + builder = GaussianRelaxWorkChain.get_builder() + builder.gaussian_code = gaussian_code + # builder.formchk_code = formchk_code + # builder.cubegen_code = cubegen_code + builder.structure = StructureData(ase=ase_geom) + builder.functional = Str("B3LYP") + builder.empirical_dispersion = Str("GD3") + builder.basis_set = Str("6-311+G(d,p)") + # builder.basis_set_scf = Str("STO-3G") + builder.multiplicity = Int(0) + builder.tight = Bool(True) + builder.options = Dict( + { + "resources": { + "tot_num_mpiprocs": 4, + "num_machines": 1, + }, + "max_wallclock_seconds": 1 * 60 * 60, + "max_memory_kb": 4 * 1024 * 1024, + } + ) + + _, wc_node = run_get_node(builder) + + assert wc_node.is_finished_ok + + # pp.make_report(wc_node, nb=False, save_image_loc=OUTPUT_DIR) + + +if __name__ == "__main__": + _example_gaussian_spin( + load_code("gaussian@tigu"), + # load_code("formchk@localhost"), + # load_code("cubegen@localhost"), + ) diff --git a/examples/workflows/naphthalene.xyz b/examples/workflows/naphthalene.xyz new file mode 100644 index 0000000..2f13db8 --- /dev/null +++ b/examples/workflows/naphthalene.xyz @@ -0,0 +1,20 @@ +18 +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 diff --git a/pyproject.toml b/pyproject.toml index 2229c90..5d906aa 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -70,6 +70,7 @@ dev = [ "nanotech_empa.gaussian.constr_opt_chain" = "aiida_nanotech_empa.workflows.gaussian:GaussianConstrOptChainWorkChain" "nanotech_empa.gaussian.casscf" = "aiida_nanotech_empa.workflows.gaussian:GaussianCasscfWorkChain" "nanotech_empa.gaussian.casscf_series" = "aiida_nanotech_empa.workflows.gaussian:GaussianCasscfSeriesWorkChain" +"nanotech_empa.gaussian.nics" = "aiida_nanotech_empa.workflows.gaussian:GaussianNicsWorkChain" "nanotech_empa.cp2k.geo_opt" = "aiida_nanotech_empa.workflows.cp2k:Cp2kGeoOptWorkChain" "nanotech_empa.cp2k.fragment_separation" = "aiida_nanotech_empa.workflows.cp2k:Cp2kFragmentSeparationWorkChain" "nanotech_empa.cp2k.ads_gw_ic" = "aiida_nanotech_empa.workflows.cp2k:Cp2kAdsorbedGwIcWorkChain"