Dev/sym opt

CHANGELOG.md

@@ -16,6 +16,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - `GXTBConfig` class for the g-xTB method, supporting SCF cycles check
 - support for TURBOMOLE as QM engine.
 - updated the parallelization to work over the number of molecules
+- posibillity to generate NCI complexes with symmetrie as mindless molecules.
 
 ### Fixed
 - version string is now correctly formatted and printed

mindlessgen.toml

@@ -17,6 +17,9 @@ num_molecules = 1
 postprocess = false
 # > Switch molecule structure XYZ writing on and off. Default: true. Options: <bool>
 write_xyz = true
+# > Switch structure modification on and off. Defaul: false. Options: <bool>
+# > postprocess has to be turned on and the postprocess engine has to be turbomole.
+structure_mod = false
 
 [generate]
 # > Minimum number of atoms in the generated molecule. Options: <int>
@@ -104,3 +107,15 @@ functional = "PBE"
 basis = "def2-SVP"
 # > Maximum number of SCF cycles: Options: <int>
 scf_cycles = 100
+
+[modification]
+# > There is only the posibillity to do one modification at one molecule.
+# > define the distance of the opposit NCI fragments. Options: <float>
+distance = 3.0
+# > mirror the molecule along the y,z plane.
+mirror = false
+# > rotate the molecule around the z axis. n is the Cn rotaion. Options: rotation: <bool>, n_fold_rotation: <int> (2, 3, 4...)
+rotation = false
+n_fold_rotation = 2
+# > inverts the molecule through the origin.
+inversion = false

src/mindlessgen/generator/main.py

@@ -18,7 +18,7 @@
 from tqdm import tqdm
 
 # Internal modules
-from ..molecules import generate_random_molecule, Molecule
+from ..molecules import generate_random_molecule, Molecule, structure_modification_mol
 from ..qm import (
     XTB,
     get_xtb_path,
@@ -191,12 +191,12 @@ def single_molecule_generator(
     with resources.occupy_cores(ncores):
         # print a decent header for each molecule iteration
         if config.general.verbosity > 0:
-            print(f"\n{'='*80}")
+            print(f"\n{'=' * 80}")
             print(
-                f"{'='*22} Generating molecule {molcount + 1:<4} of "
-                + f"{config.general.num_molecules:<4} {'='*24}"
+                f"{'=' * 22} Generating molecule {molcount + 1:<4} of "
+                + f"{config.general.num_molecules:<4} {'=' * 24}"
             )
-            print(f"{'='*80}")
+            print(f"{'=' * 80}")
 
         with setup_managers(ncores, ncores) as (executor, manager, resources_local):
             stop_event = manager.Event()
@@ -315,6 +315,21 @@ def single_molecule_step(
         if config.refine.debug:
             stop_event.set()
 
+    if config.general.structure_mod:
+        # raise SystemExit("Structure modification is not implemented yet.")
+        try:
+            optimized_molecule = structure_modification_mol(
+                optimized_molecule, config.modification, config.general.verbosity
+            )
+        except RuntimeError as e:
+            if config.general.verbosity > 0:
+                print(f"Structure modification failed for cycle {cycle + 1}.")
+                if config.general.verbosity > 1:
+                    print(e)
+            return None
+        if config.general.verbosity > 1:
+            print("Structure modification successful.")
+
     if config.general.postprocess:
         try:
             optimized_molecule = postprocess_mol(

src/mindlessgen/molecules/__init__.py

@@ -11,6 +11,7 @@
 )
 from .refinement import iterative_optimization, detect_fragments
 from .postprocess import postprocess_mol
+from .structure_modification import structure_modification_mol
 from .miscellaneous import (
     set_random_charge,
     get_three_d_metals,
@@ -43,4 +44,5 @@
     "ati_to_atlist",
     "atlist_to_ati",
     "postprocess_mol",
+    "structure_modification_mol",
 ]

src/mindlessgen/molecules/structure_modification.py

@@ -0,0 +1,96 @@
+"""
+Modifies the generated molecules to get symmetric NCI structures.
+"""
+
+import numpy as np
+from .molecule import Molecule
+from ..prog import StructureModConfig
+
+
+def structure_modification_mol(
+    mol: Molecule, config: StructureModConfig, verbosity: int = 1
+) -> Molecule:
+    """
+    Modifies the generated molecules to get symmetric NCI structures.
+
+    Arguments:
+    mol (Molecule): Molecule to modify
+
+    Returns:
+    Molecule: Modified molecule
+    """
+
+    if verbosity > 2:
+        print("Modifying molecule structure...")
+
+    xyz = mol.xyz
+    translation_vector = xyz.min()
+    for i in range(mol.num_atoms):
+        if translation_vector < 0:
+            xyz[i, 0] += (
+                -translation_vector + config.distance / 2
+            )  # extra distance to avoid that an atom is at 0
+    mol.xyz = xyz
+
+    if config.mirror:
+        mirror_matrix = np.array([[-1, 0, 0], [0, 1, 0], [0, 0, 1]])
+        xyz_mirror = xyz.copy()
+        modified_molecule = mol.copy()
+        for i in range(mol.num_atoms):
+            xyz_mirror[i] = np.dot(mirror_matrix, xyz[i])
+
+        # Combine the original and the mirror image
+        xyz_combined = np.vstack((xyz, xyz_mirror))
+        modified_molecule = update_molecule_objekt(modified_molecule, mol, xyz_combined)
+        modified_molecule.set_name_from_formula()
+        return modified_molecule
+    elif config.rotation:
+        n = config.n_fold_rotation
+        rotation_matrix = np.array(
+            [
+                [np.cos((2 * np.pi) / n), -np.sin((2 * np.pi) / n), 0],
+                [np.sin((2 * np.pi) / n), np.cos((2 * np.pi) / n), 0],
+                [0, 0, 1],
+            ]
+        )
+        xyz_rotation = xyz.copy()
+        modified_molecule = mol.copy()
+        # For loop to get n times the molecule
+        for i in range(1, n):
+            for k in range(mol.num_atoms):
+                xyz_rotation[k] = np.dot(rotation_matrix, xyz_rotation[k])
+            xyz = np.vstack((xyz, xyz_rotation))
+            modified_molecule = update_molecule_objekt(modified_molecule, mol, xyz)
+
+        modified_molecule.set_name_from_formula()
+        return modified_molecule
+    elif config.inversion:
+        inversion_matrix = np.array([[-1, 0, 0], [0, -1, 0], [0, 0, -1]])
+        xyz_inversion = xyz.copy()
+        modified_molecule = mol.copy()
+        for i in range(mol.num_atoms):
+            xyz_inversion[i] = np.dot(inversion_matrix, xyz[i])
+
+        xyz_combined = np.vstack((xyz, xyz_inversion))
+        modified_molecule = update_molecule_objekt(modified_molecule, mol, xyz_combined)
+        modified_molecule.set_name_from_formula()
+        return modified_molecule
+
+    print("structure_modification_mol is now working as a function")
+
+    return mol
+
+
+def update_molecule_objekt(
+    modified_molecule: Molecule, mol: Molecule, xyz: np.ndarray
+) -> Molecule:
+    """
+    Updates the molecule object to the new twice as large molecule objekt.
+    """
+    modified_molecule.xyz = xyz
+    modified_molecule.num_atoms += mol.num_atoms
+    modified_molecule.ati = np.hstack((modified_molecule.ati, mol.ati))
+    modified_molecule.charge += mol.charge
+    modified_molecule.uhf += mol.uhf
+    modified_molecule.atlist += mol.atlist
+    return modified_molecule

src/mindlessgen/prog/__init__.py

@@ -13,6 +13,7 @@
     GenerateConfig,
     RefineConfig,
     PostProcessConfig,
+    StructureModConfig,
 )
 from .parallel import setup_managers, ResourceMonitor, setup_blocks
 
@@ -29,4 +30,5 @@
     "setup_managers",
     "ResourceMonitor",
     "setup_blocks",
+    "StructureModConfig",
 ]

src/mindlessgen/prog/config.py

@@ -47,6 +47,7 @@ def __init__(self: GeneralConfig) -> None:
         self._num_molecules: int = 1
         self._postprocess: bool = False
         self._write_xyz: bool = True
+        self._structure_mod: bool = False
 
     def get_identifier(self) -> str:
         return "general"
@@ -189,6 +190,22 @@ def check_config(self, verbosity: int = 1) -> None:
                 + "Set '--verbosity 0' or '-P 1' to avoid this warning, or simply ignore it."
             )
 
+    @property
+    def structure_mod(self):
+        """
+        Get the structure_mod flag.
+        """
+        return self._structure_mod
+
+    @structure_mod.setter
+    def structure_mod(self, structure_mod: bool):
+        """
+        Set the structure_mod flag.
+        """
+        if not isinstance(structure_mod, bool):
+            raise TypeError("Structure modification should be a boolean.")
+        self._structure_mod = structure_mod
+
 
 class GenerateConfig(BaseConfig):
     """
@@ -1104,6 +1121,106 @@ def scf_cycles(self, max_scf_cycles: int):
         self._scf_cycles = max_scf_cycles
 
 
+class StructureModConfig(BaseConfig):
+    """
+    Configuration class for structure modification settings.
+    """
+
+    def __init__(self: StructureModConfig) -> None:
+        self.distance: float = 3.0
+        self._rotation: bool = False
+        self._n_fold_rotation: int = 4
+        self._mirror: bool = False
+        self._inversion: bool = False
+
+    def get_identifier(self) -> str:
+        return "modification"
+
+    @property
+    def distance(self):
+        """
+        Get the distance.
+        """
+        return self._distance
+
+    @distance.setter
+    def distance(self, distance: float):
+        """
+        Set the distance.
+        """
+        if not isinstance(distance, float):
+            raise TypeError("Distance should be a float.")
+        if distance <= 0:
+            raise ValueError("Distance should be greater than 0.")
+        self._distance = distance
+
+    @property
+    def rotation(self):
+        """
+        Get the rotations flag.
+        """
+        return self._rotation
+
+    @rotation.setter
+    def rotation(self, rotation: bool):
+        """
+        Set the rotation flag.
+        """
+        if not isinstance(rotation, bool):
+            raise TypeError("Rotation should be a boolean.")
+        self._rotation = rotation
+
+    @property
+    def n_fold_rotation(self):
+        """
+        Get the number of the rotation axis.
+        """
+        return self._n_fold_rotation
+
+    @n_fold_rotation.setter
+    def n_fold_rotation(self, n_fold_rotation: int):
+        """
+        Set the number of rotations.
+        """
+        if not isinstance(n_fold_rotation, int):
+            raise TypeError("Number of rotations should be an integer.")
+        if n_fold_rotation < 1:
+            raise ValueError("Number of rotations should be greater than 0.")
+        self._n_fold_rotation = n_fold_rotation
+
+    @property
+    def mirror(self):
+        """
+        Get the mirror flag.
+        """
+        return self._mirror
+
+    @mirror.setter
+    def mirror(self, mirror: bool):
+        """
+        Set the mirror flag.
+        """
+        if not isinstance(mirror, bool):
+            raise TypeError("Mirror should be a boolean.")
+        self._mirror = mirror
+
+    @property
+    def inversion(self):
+        """
+        Get the inversion flag.
+        """
+        return self._inversion
+
+    @inversion.setter
+    def inversion(self, inversion: bool):
+        """
+        Set the inversion flag.
+        """
+        if not isinstance(inversion, bool):
+            raise TypeError("Inversion should be a boolean.")
+        self._inversion = inversion
+
+
 class ConfigManager:
     """
     Overall configuration manager for the program.
@@ -1121,6 +1238,7 @@ def __init__(self, config_file: str | Path | None = None):
         self.postprocess = PostProcessConfig()
         self.generate = GenerateConfig()
         self.gxtb = GXTBConfig()
+        self.modification = StructureModConfig()
 
         if config_file:
             self.load_from_toml(config_file)

src/mindlessgen/qm/tm.py

@@ -67,6 +67,10 @@ def optimize(
             with open(temp_path / inputname, "w", encoding="utf8") as f:
                 f.write(tm_input)
 
+            singlepoint_logout = self.singlepoint(molecule, ncores, verbosity)
+            if verbosity > 2:
+                print(singlepoint_logout)
+
             # Setup the turbomole optimization command including the max number of optimization cycles
             arguments = [f"PARNODES={ncores} {self.jobex_path} -c {max_cycles}"]
             if verbosity > 2:
@@ -207,7 +211,10 @@ def _run_opt(self, temp_path: Path, arguments: list[str]) -> tuple[str, str, int
                 )
             return tm_log_out, tm_log_err, 0
         except sp.CalledProcessError as e:
-            tm_log_out = e.stdout.decode("utf8", errors="replace")
+            output_file = temp_path / "job.last"
+            if output_file.exists():
+                with open(output_file, encoding="utf-8") as file:
+                    tm_log_out = file.read()
             tm_log_err = e.stderr.decode("utf8", errors="replace")
             return tm_log_out, tm_log_err, e.returncode