diff --git a/aiida_nanotech_empa/utils/cycle_tools.py b/aiida_nanotech_empa/utils/cycle_tools.py index f4bed8d..feaf9ca 100644 --- a/aiida_nanotech_empa/utils/cycle_tools.py +++ b/aiida_nanotech_empa/utils/cycle_tools.py @@ -1,38 +1,34 @@ -import numpy as np - -import scipy as sp -import scipy.linalg - import ase import ase.io -import ase.visualize import ase.neighborlist - -import os -import shutil +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 - + 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: + + 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 @@ -49,45 +45,47 @@ def find_cycles(i_vert, cnl, max_length, cur_path, passed_edges): 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]) + 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 - + 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 @@ -99,7 +97,7 @@ def cycle_normal(cycle, h): 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: + if np.dot(normal, h * np.array([1, 1, 1])) < 0.0: normal *= -1.0 return normal diff --git a/aiida_nanotech_empa/utils/nmr.py b/aiida_nanotech_empa/utils/nmr.py index 4145659..6108d6b 100644 --- a/aiida_nanotech_empa/utils/nmr.py +++ b/aiida_nanotech_empa/utils/nmr.py @@ -1,16 +1,11 @@ -import numpy as np - -import scipy as sp -import scipy.linalg - import ase import ase.io -import ase.visualize import ase.neighborlist - -import os -import shutil +import ase.visualize import cclib +import numpy as np +import scipy as sp +import scipy.linalg from . import cycle_tools @@ -18,67 +13,73 @@ ### SETUP -def find_ref_points(ase_atoms_no_h, cycles, h = 0.0): +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) - + 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])) + 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)) - + 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) + 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 - + 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): + for i_rp in range(len(ref_pts) - 1): pt1 = ref_pts[i_rp].position - pt2 = ref_pts[i_rp+1].position + pt2 = ref_pts[i_rp + 1].position points = interp_pts(pt1, pt2, dx) - if i_rp == len(ref_pts)-2: + 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) + + ase_arr = ase.Atoms("X%d" % len(point_arr), point_arr) return ase_arr ### ----------------------------------------------------------------- ### PROCESS + def load_nics_gaussian(nics_path): sigma = [] @@ -110,14 +111,15 @@ def is_number(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 + 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)) @@ -138,24 +140,24 @@ def parse_nmr_cmo_matrix(log_file_str, property_dict): if "Full Cartesian NMR shielding tensor (ppm) for atom" in lines[i_line]: - i_atom = int(lines[i_line].replace('(', ')').split(')')[-2]) - 1 + i_atom = int(lines[i_line].replace("(", ")").split(")")[-2]) - 1 i_line += 1 - if 'Canonical MO contributions' in lines[i_line]: + if "Canonical MO contributions" in lines[i_line]: for _i in range(2000): i_line += 1 - if 'Total' in lines[i_line]: + 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:]]) + 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 0b5d156..809a790 100644 --- a/aiida_nanotech_empa/workflows/gaussian/__init__.py +++ b/aiida_nanotech_empa/workflows/gaussian/__init__.py @@ -4,10 +4,10 @@ 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 -from .nics_workchain import GaussianNicsWorkChain __all__ = ( "GaussianScfWorkChain", diff --git a/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py b/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py index 966a42f..e667e83 100644 --- a/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py +++ b/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py @@ -1,29 +1,33 @@ -import numpy as np 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 -from ...utils import cycle_tools, nmr -from ...utils import common_utils + @engine.calcfunction -def prepare_nmr_structure(structure,height): +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]), - ]) + 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']) + 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): @@ -81,11 +85,11 @@ def define(cls, spec): 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), + engine.if_(cls.should_submit_opt)(cls.submit_opt, cls.inspect_opt), cls.submit_nics, cls.inspect_nics, cls.finalize, @@ -98,14 +102,15 @@ def define(cls, spec): "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" @@ -127,23 +132,25 @@ def submit_opt(self): 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 + 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.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 @@ -151,23 +158,23 @@ def submit_nics(self): builder.cdiis = orm.Bool(True) builder.nmr = orm.Bool(True) builder.maxcycle = orm.Int(2048) - builder.conver = orm.Int(8) + 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) - return engine.ExitCode(0) - + return engine.ExitCode(0) + def finalize(self): - self.report("Finalizing...") \ No newline at end of file + self.report("Finalizing...") diff --git a/aiida_nanotech_empa/workflows/gaussian/oldnics_workchain.py b/aiida_nanotech_empa/workflows/gaussian/oldnics_workchain.py index 06fe78e..7a7a002 100644 --- a/aiida_nanotech_empa/workflows/gaussian/oldnics_workchain.py +++ b/aiida_nanotech_empa/workflows/gaussian/oldnics_workchain.py @@ -1,7 +1,8 @@ import numpy as np from aiida import engine, orm, plugins -from ...utils import common_utils,nmr,cycle_tools +from ...utils import common_utils, cycle_tools, nmr + #from .delta_scf_workchain import GaussianDeltaScfWorkChain #from .natorb_workchain import GaussianNatOrbWorkChain #from .relax_workchain import GaussianRelaxWorkChain @@ -14,7 +15,7 @@ @engine.calcfunction def nics_structure(structure=None, h=1.0): - + ase_geo = structure.get_ase() # orient the geometry pos = ase_geo.positions @@ -25,10 +26,10 @@ def nics_structure(structure=None, h=1.0): ]) inds = np.argsort(-extents) ase_atoms = ase.Atoms(numbers=ase_geo.numbers, positions=pos[:, inds]) - ase_atoms_no_h = ase.Atoms([a for a in ase_atoms if a.symbol != 'H']) + 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) - ase_geom = ase_atoms + nmr.find_ref_points(ase_atoms_no_h, cycles, h.value) - return orm.StructureData(ase=ase_geom) + ase_geom = ase_atoms + nmr.find_ref_points(ase_atoms_no_h, cycles, h.value) + return orm.StructureData(ase=ase_geom) class GaussianNicsWorkChain(engine.WorkChain): @classmethod @@ -189,7 +190,7 @@ def submit_next_steps(self): builder.structure = self.ctx.gs_structure builder.functional = self.inputs.functional builder.empirical_dispersion = self.inputs.empirical_dispersion - builder.basis_set = + builder.basis_set = builder.multiplicity = self.inputs.multiplicity #if "options" in self.inputs: diff --git a/aiida_nanotech_empa/workflows/gaussian/relax_workchain.py b/aiida_nanotech_empa/workflows/gaussian/relax_workchain.py index 1ea2d76..fe6b210 100644 --- a/aiida_nanotech_empa/workflows/gaussian/relax_workchain.py +++ b/aiida_nanotech_empa/workflows/gaussian/relax_workchain.py @@ -60,13 +60,13 @@ def define(cls, spec): 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, @@ -313,10 +313,10 @@ 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) + 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: diff --git a/aiida_nanotech_empa/workflows/gaussian/scf_workchain.py b/aiida_nanotech_empa/workflows/gaussian/scf_workchain.py index db8644e..e36adda 100644 --- a/aiida_nanotech_empa/workflows/gaussian/scf_workchain.py +++ b/aiida_nanotech_empa/workflows/gaussian/scf_workchain.py @@ -71,13 +71,13 @@ def define(cls, spec): 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, @@ -254,15 +254,15 @@ def scf(self): "charge": 0, "multiplicity": self.ctx.mult, "route_parameters": { - "scf": {"maxcycle": 140 }, + "scf": {"maxcycle": 140}, }, } ) - 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) + 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: diff --git a/examples/workflows/benzene.xyz b/examples/workflows/benzene.xyz index 66fd521..5255bb3 100644 --- a/examples/workflows/benzene.xyz +++ b/examples/workflows/benzene.xyz @@ -1,14 +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 \ No newline at end of file +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 index 9c43f4b..12fa82f 100644 --- a/examples/workflows/example_gaussian_nics.py +++ b/examples/workflows/example_gaussian_nics.py @@ -1,16 +1,16 @@ -import os import click -from ase.build import molecule import numpy as np from aiida import engine, orm from aiida.plugins import WorkflowFactory +from ase.build import molecule import aiida_nanotech_empa.utils.gaussian_wcs_postprocess as pp GaussianNicsWorkChain = WorkflowFactory("nanotech_empa.gaussian.nics") -def _example_gaussian_nics(gaussian_code,opt): - ase_geom = molecule('C6H6') + +def _example_gaussian_nics(gaussian_code, opt): + ase_geom = molecule("C6H6") ase_geom.cell = np.diag([10.0, 10.0, 10.0]) builder = GaussianNicsWorkChain.get_builder() @@ -26,25 +26,26 @@ def _example_gaussian_nics(gaussian_code,opt): "num_machines": 1, }, "max_wallclock_seconds": 1 * 60 * 60, - "max_memory_kb": 8 * 1024 * 1024, #GB + "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 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),False) + uuid = _example_gaussian_nics(orm.load_code(gaussian_code), False) print(f"WorkChain completed uuid: {uuid}") if __name__ == "__main__": - run_nics() \ No newline at end of file + run_nics() diff --git a/examples/workflows/example_gaussian_opt.py b/examples/workflows/example_gaussian_opt.py index f7beade..7a9bb70 100644 --- a/examples/workflows/example_gaussian_opt.py +++ b/examples/workflows/example_gaussian_opt.py @@ -3,7 +3,7 @@ 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.orm import Bool, Int, Str, StructureData, load_code from aiida.plugins import WorkflowFactory import aiida_nanotech_empa.utils.gaussian_wcs_postprocess as pp @@ -14,19 +14,19 @@ OUTPUT_DIR = os.path.dirname(os.path.abspath(__file__)) -def _example_gaussian_spin(gaussian_code):#, formchk_code, cubegen_code): +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.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.basis_set_scf = Str("STO-3G") builder.multiplicity = Int(0) builder.tight = Bool(True) builder.options = Dict( @@ -38,18 +38,18 @@ def _example_gaussian_spin(gaussian_code):#, formchk_code, cubegen_code): "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) + # 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"), + # load_code("formchk@localhost"), + # load_code("cubegen@localhost"), )