Merge pull request grimme-lab#53 from jonathan-schoeps/dev/contract_c…

…oordinates

A new function to contract the coordinates after the random generation to ensure that the atoms are not to far apart.

CHANGELOG.md

@@ -17,6 +17,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 - Support for the novel "g-xTB" method (working title: GP3-xTB)
+- A function which contracts the coordinates after the initial generation.
+- A function which is able to printout the xyz coordinates to the terminal similar to the `.xyz` layout.
 
 ### Breaking Changes
 - Removal of the `dist_threshold` flag and in the `-toml` file.

mindlessgen.toml

@@ -28,9 +28,11 @@ init_scaling = 3.0
 # > Increase in the coordinate scaling factor per trial after check_distance was not met. Options: <float>
 increase_scaling_factor = 1.3
 # > Scaling factor for the employed van der Waals radii. Options: <float>
-scale_vdw_radii = 1.3333
+scale_vdw_radii = 1.25
 # > Scaling factor for the minimal bondlength based on the sum of the van der Waals radii. Options: <float>
-scale_minimal_bondlength = 0.75
+scale_minimal_bondlength = 0.8
+# > Contract the coordinates after the initial generation. Leads to more cluster-like and less extended structures. Options: <bool>
+contract_coords = false
 # > Atom types and their minimum and maximum occurrences. Format: "<element>:<min_count>-<max_count>"
 # > Elements that are not specified are only added by random selection.
 # > A star sign (*) can be used as a wildcard for integer value.

src/mindlessgen/cli/cli_parser.py

@@ -137,6 +137,12 @@ def cli_parser(argv: Sequence[str] | None = None) -> dict:
         required=False,
         help="List of forbidden elements.",
     )
+    parser.add_argument(
+        "--contract_coords",
+        type=bool,
+        required=False,
+        help="Contract the coordinates of the molecule after the coordinats generation.",
+    )
 
     ### Refinement arguments ###
     parser.add_argument(
@@ -195,19 +201,19 @@ def cli_parser(argv: Sequence[str] | None = None) -> dict:
         required=False,
         help="Path to the xTB binary.",
     )
+    parser.add_argument(
+        "--xtb-level",
+        type=int,
+        required=False,
+        help="Level of theory to use in xTB.",
+    )
     parser.add_argument(
         "--postprocess-debug",
         action="store_true",
         default=None,
         required=False,
         help="Print debug information during postprocessing.",
     )
-    parser.add_argument(
-        "--xtb-level",
-        type=int,
-        required=False,
-        help="Level of theory to use in xTB.",
-    )
 
     ### ORCA specific arguments ###
     parser.add_argument(
@@ -273,6 +279,7 @@ def cli_parser(argv: Sequence[str] | None = None) -> dict:
         "forbidden_elements": args_dict["forbidden_elements"],
         "scale_vdw_radii": args_dict["scale_vdw_radii"],
         "scale_minimal_bondlength": args_dict["scale_minimal_bondlength"],
+        "contract_coords": args_dict["contract_coords"],
     }
     # XTB specific arguments
     rev_args_dict["xtb"] = {"xtb_path": args_dict["xtb_path"]}

src/mindlessgen/molecules/generate_molecule.py

@@ -41,6 +41,12 @@ def generate_random_molecule(
         verbosity=verbosity,
         scale_bondlength=config_generate.scale_minimal_bondlength,
     )
+    if config_generate.contract_coords:
+        mol.xyz = contract_coordinates(
+            xyz=mol.xyz,
+            ati=mol.ati,
+            scale_minimal_distance=config_generate.scale_minimal_bondlength,
+        )
     mol.charge, mol.uhf = set_random_charge(mol.ati, verbosity)
     mol.set_name_from_formula()
 
@@ -371,6 +377,39 @@ def generate_random_coordinates(at: np.ndarray) -> tuple[np.ndarray, np.ndarray]
     return xyz, ati
 
 
+def contract_coordinates(
+    xyz: np.ndarray, ati: np.ndarray, scale_minimal_distance: float
+) -> np.ndarray:
+    """
+    Pull the atoms towards the origin.
+    """
+    # Initialize the old coordinates as an array of zeros
+    xyz_old: np.ndarray = np.zeros_like(xyz)
+    cycle = 0
+    # Break if the coordinates do not change
+    while not np.array_equal(xyz_old, xyz):
+        cycle += 1
+        if cycle > 2500:
+            raise RuntimeError(
+                "Could not contract the coordinates in a reasonable amount of cycles."
+            )
+        xyz_old = xyz.copy()
+        # Go through the atoms dimension of the xyz array in a reversed order
+        # Justification: First atoms are most likely hydrogen atoms, which should be moved last
+        for i in range(len(xyz) - 1, -1, -1):
+            atom_xyz = xyz[i]
+            atom_xyz_norm = np.linalg.norm(atom_xyz)
+            normalized_atom_xyz = atom_xyz / atom_xyz_norm
+            # Shift the atom only if it is closer to the origin after the shift
+            if atom_xyz_norm > 0.1:
+                # Pull the atom towards the origin
+                xyz[i] -= normalized_atom_xyz * 0.2
+                # When the check_distances function returns False, reset the atom coordinates
+                if not check_distances(xyz, ati, scale_minimal_distance):
+                    xyz[i] = xyz_old[i]
+    return xyz
+
+
 def check_distances(xyz: np.ndarray, ati: np.ndarray, scale_bondlength: float) -> bool:
     """
     Check if the distances between atoms are larger than a threshold.

src/mindlessgen/molecules/molecule.py

@@ -459,12 +459,29 @@ def atlist(self, value: np.ndarray):
 
         self._atlist = value
 
-    def print_xyz(self):
+    def get_xyz_str(self) -> str:
         """
-        Print the XYZ coordinates of the molecule.
+        Obtain a string with the full XYZ file information of the molecule.
         """
-        # TODO: Implement a better way to print the coordinates
-        print(self._xyz)
+        xyz_str = f"{self.num_atoms}\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"
+        xyz_str += commentline
+        for i in range(self.num_atoms):
+            xyz_str += (
+                f"{PSE[self.ati[i]+1]:<5} "
+                + f"{self.xyz[i, 0]:>12.7f} "
+                + f"{self.xyz[i, 1]:>12.7f} "
+                + f"{self.xyz[i, 2]:>12.7f}\n"
+            )
+        return xyz_str
 
     def write_xyz_to_file(self, filename: str | Path | None = None):
         """
@@ -496,24 +513,7 @@ def write_xyz_to_file(self, filename: str | Path | None = None):
             filename = Path("mlm_" + self.name + ".xyz").resolve()
 
         with open(filename, "w", encoding="utf8") as f:
-            f.write(f"{self.num_atoms}\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} "
-                    + f"{self.xyz[i, 0]:>12.7f} "
-                    + f"{self.xyz[i, 1]:>12.7f} "
-                    + f"{self.xyz[i, 2]:>12.7f}\n"
-                )
+            f.write(self.get_xyz_str())
         # if the charge is set, write it to a '.CHRG' file
         if self._charge is not None:
             with open(filename.with_suffix(".CHRG"), "w", encoding="utf8") as f:

src/mindlessgen/prog/config.py

@@ -178,8 +178,9 @@ def __init__(self: GenerateConfig) -> None:
         self._increase_scaling_factor: float = 1.3
         self._element_composition: dict[int, tuple[int | None, int | None]] = {}
         self._forbidden_elements: list[int] | None = None
-        self._scale_vdw_radii: float = 4.0 / 3.0
-        self._scale_minimal_bondlength: float = 0.75
+        self._scale_vdw_radii: float = 1.25
+        self._scale_minimal_bondlength: float = 0.8
+        self._contract_coords: bool = False
 
     def get_identifier(self) -> str:
         return "generate"
@@ -383,6 +384,22 @@ def scale_minimal_bondlength(self, scale_minimal_bondlength: float):
             raise ValueError("Scale minimal bond length should be greater than 0.")
         self._scale_minimal_bondlength = scale_minimal_bondlength
 
+    @property
+    def contract_coords(self):
+        """
+        Get the contract_coords flag.
+        """
+        return self._contract_coords
+
+    @contract_coords.setter
+    def contract_coords(self, contract_coords: bool):
+        """
+        Set the contract_coords flag.
+        """
+        if not isinstance(contract_coords, bool):
+            raise TypeError("Contract coords should be a boolean.")
+        self._contract_coords = contract_coords
+
 
 class RefineConfig(BaseConfig):
     """