xyz file changes and more DFT options in post-processing

CHANGELOG.md

@@ -4,7 +4,18 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
-## [0.4.0] - 2024-XX-YY
+## [Unreleased]
+### Changed
+- Default file name of `.xyz` file contains prefix `mlm_`
+- Comment line of `.xyz` file contains the total charge and number of unpaired electrons
+- Default ORCA calculation changed from r2SCAN-3c to PBE/def2-SVP
+
+### Added
+- Optimization via DFT in the post-processing step
+- Detailed input of ORCA settings (functional, basis, grid size, SCF cycles, ...) possible
+- Maximum number of optimization cycles are an argument for the `QMMethod.optimize` base function
+
+## [0.3.0] - 2024-08-20
 ### Breaking Changes
 - ...
 

mindlessgen.toml

@@ -45,7 +45,11 @@ engine = "xtb"
 
 [postprocess]
 # > Engine for the post-processing part. Options: 'xtb', 'orca'
-engine = "xtb"
+engine = "orca"
+# > Optimize geometry in the post-processing part. Options: <bool>
+optimize = true
+# > Optimization cycles for the post-processing part. If not given, the program default is chosen. Options: <int>
+opt_cycles = 5
 
 [xtb]
 # > Path to the xtb executable. Options: <str | Path>
@@ -54,3 +58,11 @@ xtb_path = "/path/to/xtb"
 [orca]
 # > Path to the orca executable. Options: <str | Path>
 orca_path = "/path/to/orca"
+# > Functional/Method: Options: <str>
+functional = "PBE"
+# > Basis set: Options: <str>
+basis = "SV(P)"
+# > Gridsize for the numerical integration: Options: <int => 1> = coarse, <int = 2> => medium, <int = 3> => fine
+gridsize = 1
+# > Maximum number of SCF cycles: Options: <int>
+scf_cycles = 100

src/mindlessgen/cli/cli_parser.py

@@ -132,6 +132,19 @@ def cli_parser(argv: Sequence[str] | None = None) -> dict:
         choices=["xtb", "orca"],
         help="QM engine to use for postprocessing.",
     )
+    parser.add_argument(
+        "--postprocess-optimize",
+        action="store_true",
+        default=None,
+        required=False,
+        help="Optimize the molecule during post-processing.",
+    )
+    parser.add_argument(
+        "--postprocess-opt-cycles",
+        type=int,
+        required=False,
+        help="Number of optimization cycles in postprocessing.",
+    )
     # xTB specific arguments
     parser.add_argument(
         "--xtb-path",
@@ -146,6 +159,30 @@ def cli_parser(argv: Sequence[str] | None = None) -> dict:
         required=False,
         help="Path to the ORCA binary.",
     )
+    parser.add_argument(
+        "--orca-functional",
+        type=str,
+        required=False,
+        help="Functional to use in ORCA.",
+    )
+    parser.add_argument(
+        "--orca-basis",
+        type=str,
+        required=False,
+        help="Basis set to use in ORCA.",
+    )
+    parser.add_argument(
+        "--orca-gridsize",
+        type=int,
+        required=False,
+        help="Solvent to use in ORCA.",
+    )
+    parser.add_argument(
+        "--orca-scf-cycles",
+        type=int,
+        required=False,
+        help="Maximum number of SCF cycles in ORCA.",
+    )
     args = parser.parse_args(argv)
     args_dict = vars(args)
 
@@ -176,7 +213,17 @@ def cli_parser(argv: Sequence[str] | None = None) -> dict:
     # XTB specific arguments
     rev_args_dict["xtb"] = {"xtb_path": args_dict["xtb_path"]}
     # ORCA specific arguments
-    rev_args_dict["orca"] = {"orca_path": args_dict["orca_path"]}
-    rev_args_dict["postprocess"] = {"engine": args_dict["postprocess_engine"]}
+    rev_args_dict["orca"] = {
+        "orca_path": args_dict["orca_path"],
+        "functional": args_dict["orca_functional"],
+        "basis": args_dict["orca_basis"],
+        "gridsize": args_dict["orca_gridsize"],
+        "scf_cycles": args_dict["orca_scf_cycles"],
+    }
+    rev_args_dict["postprocess"] = {
+        "engine": args_dict["postprocess_engine"],
+        "optimize": args_dict["postprocess_optimize"],
+        "opt_cycles": args_dict["postprocess_opt_cycles"],
+    }
 
     return rev_args_dict

src/mindlessgen/cli/entrypoint.py

@@ -47,8 +47,8 @@ def console_entry_point(argv: Sequence[str] | None = None) -> int:
             try:
                 molecule.write_xyz_to_file()
                 if config.general.verbosity > 0:
-                    print(f"Written molecule file '{molecule.name}.xyz'.")
-                molecules_written.append(molecule.name)
+                    print(f"Written molecule file 'mlm_{molecule.name}.xyz'.")
+                molecules_written.append("mlm_" + molecule.name)
             except Exception as e:
                 warnings.warn(f"Failed to write molecule file: {e}")
         if exitcode != 0:

src/mindlessgen/generator/main.py

@@ -185,11 +185,11 @@ def single_molecule_generator(
                 if config.general.verbosity > 1:
                     print(e)
             return None
+        if config.general.verbosity > 1:
+            print("Postprocessing successful.")
 
     if not stop_event.is_set():
         stop_event.set()  # Signal other processes to stop
-        if config.general.verbosity > 1:
-            print("Postprocessing successful. Optimized molecule:")
         return optimized_molecule
     else:
         return None
@@ -223,7 +223,7 @@ def header(version: str) -> str:
 # Define a utility function to set up the required engine
 def setup_engines(
     engine_type: str,
-    engine_config: ConfigManager,
+    cfg: ConfigManager,
     xtb_path_func,
     orca_path_func,
 ):
@@ -232,19 +232,19 @@ def setup_engines(
     """
     if engine_type == "xtb":
         try:
-            path = xtb_path_func(engine_config.xtb.xtb_path)
+            path = xtb_path_func(cfg.xtb.xtb_path)
             if not path:
                 raise ImportError("xtb not found.")
         except ImportError as e:
             raise ImportError("xtb not found.") from e
-        return XTB(path)
+        return XTB(path, cfg.xtb)
     elif engine_type == "orca":
         try:
-            path = orca_path_func(engine_config.orca.orca_path)
+            path = orca_path_func(cfg.orca.orca_path)
             if not path:
                 raise ImportError("orca not found.")
         except ImportError as e:
             raise ImportError("orca not found.") from e
-        return ORCA(path)
+        return ORCA(path, cfg.orca)
     else:
         raise NotImplementedError("Engine not implemented.")

src/mindlessgen/molecules/molecule.py

@@ -458,11 +458,20 @@ def write_xyz_to_file(self, filename: str | Path | None = None):
             if not isinstance(filename, Path):
                 filename = Path(filename).resolve()
         else:
-            filename = Path(self.name + ".xyz").resolve()
+            filename = Path("mlm_" + self.name + ".xyz").resolve()
 
         with open(filename, "w", encoding="utf8") as f:
             f.write(f"{self.num_atoms}\n")
-            f.write(f"Generated by mindlessgen-v{__version__}\n")
+            try:
+                commentline = f"Total charge: {self.charge} ; "
+            except ValueError:
+                commentline = ""
+            try:
+                commentline = commentline + f"Unpaired electrons: {self.uhf} ; "
+            except ValueError:
+                pass
+            commentline = commentline + f"Generated by mindlessgen-v{__version__}\n"
+            f.write(commentline)
             for i in range(self.num_atoms):
                 f.write(
                     f"{PSE[self.ati[i]+1]:<5} "

src/mindlessgen/molecules/postprocess.py

@@ -22,8 +22,20 @@ def postprocess_mol(
     """
 
     # Run a singlepoint calculation with the postprocess engine to see if it is stable
-    try:
-        engine.singlepoint(mol, verbosity=verbosity)
-    except RuntimeError as e:
-        raise RuntimeError("Single point calculation in postprocessing failed.") from e
-    return mol
+    if verbosity > 2:
+        print("Postprocessing molecule...")
+    if config.optimize:
+        try:
+            postprocmol = engine.optimize(
+                mol, max_cycles=config.opt_cycles, verbosity=verbosity
+            )
+        except RuntimeError as e:
+            raise RuntimeError("Optimization in postprocessing failed.") from e
+    else:
+        try:
+            engine.singlepoint(mol, verbosity=verbosity)
+        except RuntimeError as e:
+            raise RuntimeError(
+                "Single point calculation in postprocessing failed."
+            ) from e
+    return postprocmol

src/mindlessgen/molecules/refinement.py

@@ -22,12 +22,16 @@ def iterative_optimization(
     for cycle in range(config_refine.max_frag_cycles):
         # Optimize the current molecule
         try:
-            rev_mol = engine.optimize(rev_mol, verbosity)
+            rev_mol = engine.optimize(rev_mol, None, verbosity)
         except RuntimeError as e:
             raise RuntimeError(
                 f"Optimization failed at fragmentation cycle {cycle}: {e}"
             ) from e
 
+        if verbosity > 2:
+            # Print coordinates of optimized molecule
+            print(f"Optimized molecule in cycle {cycle + 1}:\n{rev_mol.xyz}")
+
         # Detect fragments from the optimized molecule
         fragmols = detect_fragments(rev_mol, verbosity)
 

src/mindlessgen/prog/config.py

@@ -387,6 +387,8 @@ class PostProcessConfig(BaseConfig):
 
     def __init__(self):
         self._engine: str = "orca"
+        self._opt_cycles: int | None = None
+        self._optimize: bool = False
 
     def get_identifier(self) -> str:
         return "postprocess"
@@ -409,6 +411,40 @@ def engine(self, engine: str):
             raise ValueError("Postprocess engine can only be xtb or orca.")
         self._engine = engine
 
+    @property
+    def optimize(self):
+        """
+        Get the optimization flag for post-processing.
+        """
+        return self._optimize
+
+    @optimize.setter
+    def optimize(self, optimize: bool):
+        """
+        Set the optimization flag for post-processing.
+        """
+        if not isinstance(optimize, bool):
+            raise TypeError("Optimize should be a boolean.")
+        self._optimize = optimize
+
+    @property
+    def opt_cycles(self):
+        """
+        Get the optimization cycles for post-processing.
+        """
+        return self._opt_cycles
+
+    @opt_cycles.setter
+    def opt_cycles(self, opt_cycles: int):
+        """
+        Set the optimization cycles for post-processing.
+        """
+        if not isinstance(opt_cycles, int):
+            raise TypeError("Optimization cycles should be an integer.")
+        if opt_cycles < 0:
+            raise ValueError("Optimization cycles should be 0 or greater.")
+        self._opt_cycles = opt_cycles
+
 
 class XTBConfig(BaseConfig):
     """
@@ -445,6 +481,10 @@ class ORCAConfig(BaseConfig):
 
     def __init__(self):
         self._orca_path: str | Path = "orca"
+        self._functional: str = "PBE"
+        self._basis: str = ""
+        self._gridsize: int = 1
+        self._scf_cycles: int = 100
 
     def get_identifier(self) -> str:
         return "orca"
@@ -465,6 +505,76 @@ def orca_path(self, orca_path: str | Path):
             raise TypeError("orca_path should be a string or Path.")
         self._orca_path = orca_path
 
+    @property
+    def functional(self):
+        """
+        Get the ORCA functional/method.
+        """
+        return self._functional
+
+    @functional.setter
+    def functional(self, functional: str):
+        """
+        Set the ORCA functional/method.
+        """
+        if not isinstance(functional, str):
+            raise TypeError("Functional should be a string.")
+        self._functional = functional
+
+    @property
+    def basis(self):
+        """
+        Get the ORCA basis set.
+        """
+        return self._basis
+
+    @basis.setter
+    def basis(self, basis: str):
+        """
+        Set the ORCA basis set.
+        """
+        if not isinstance(basis, str):
+            raise TypeError("Basis should be a string.")
+        self._basis = basis
+
+    @property
+    def gridsize(self):
+        """
+        Get the ORCA gridsize for numerical integration.
+        """
+        return self._gridsize
+
+    @gridsize.setter
+    def gridsize(self, gridsize: int):
+        """
+        Set the ORCA gridsize for numerical integration.
+        """
+        if not isinstance(gridsize, int):
+            raise TypeError("Gridsize should be an integer.")
+        if gridsize not in [1, 2, 3]:
+            raise ValueError(
+                "Gridsize should be 1, 2, or 3. (-> DEFGRID1, DEFGRID2, DEFGRID3)"
+            )
+        self._gridsize = gridsize
+
+    @property
+    def scf_cycles(self):
+        """
+        Get the maximum number of SCF cycles.
+        """
+        return self._scf_cycles
+
+    @scf_cycles.setter
+    def scf_cycles(self, max_scf_cycles: int):
+        """
+        Set the maximum number of SCF cycles.
+        """
+        if not isinstance(max_scf_cycles, int):
+            raise TypeError("Max SCF cycles should be an integer.")
+        if max_scf_cycles < 1:
+            raise ValueError("Max SCF cycles should be greater than 0.")
+        self._scf_cycles = max_scf_cycles
+
 
 class ConfigManager:
     """

src/mindlessgen/qm/base.py

@@ -20,7 +20,9 @@ def __init__(self, path: str | Path, verbosity: int = 1):
         self.verbosity = verbosity
 
     @abstractmethod
-    def optimize(self, molecule: Molecule, verbosity: int = 1) -> Molecule:
+    def optimize(
+        self, molecule: Molecule, max_cycles: int | None = None, verbosity: int = 1
+    ) -> Molecule:
         """
         Define the optimization process.