Implemented the van der Waals radii inside the check distance function

CHANGELOG.md

@@ -7,6 +7,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased]
 ### Changed
 - vdW radii scaling parameter can now be adjusted via `mindlessgen.toml` or CLI
+- The check_distance function now checks based on the sum of the van der Waals radii and a scaling factor acessible via `mindlessgen.toml` or CLI
 - better type hints for `Callables`
 
 ### Fixed

mindlessgen.toml

@@ -31,6 +31,8 @@ increase_scaling_factor = 1.3
 dist_threshold = 1.2
 # > Scaling factor for the employed van der Waals radii. Options: <float>
 scale_vdw_radii = 1.3333
+# > Scaling factor for the minimal bondlength based on the sum of the van der Waals radii. Options: <float>
+scale_minimal_bondlength = 0.75
 # > 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

@@ -92,6 +92,12 @@ def cli_parser(argv: Sequence[str] | None = None) -> dict:
         required=False,
         help="Scaling factor for van der Waals radii.",
     )
+    parser.add_argument(
+        "--scale-minimal-bondlength",
+        type=float,
+        required=False,
+        help="Minimum bond length scaling factor.",
+    )
 
     ### Molecule generation arguments ###
     parser.add_argument(
@@ -273,6 +279,7 @@ def cli_parser(argv: Sequence[str] | None = None) -> dict:
         "element_composition": args_dict["element_composition"],
         "forbidden_elements": args_dict["forbidden_elements"],
         "scale_vdw_radii": args_dict["scale_vdw_radii"],
+        "scale_minimal_bondlength": args_dict["scale_minimal_bondlength"],
     }
     # XTB specific arguments
     rev_args_dict["xtb"] = {"xtb_path": args_dict["xtb_path"]}

src/mindlessgen/molecules/generate_molecule.py

@@ -6,6 +6,7 @@
 import numpy as np
 from ..prog import GenerateConfig
 from .molecule import Molecule
+from .refinement import get_cov_radii, COV_RADII
 from .miscellaneous import (
     set_random_charge,
     get_alkali_metals,
@@ -36,9 +37,9 @@ def generate_random_molecule(
     mol.xyz, mol.ati = generate_coordinates(
         at=mol.atlist,
         scaling=config_generate.init_coord_scaling,
-        dist_threshold=config_generate.dist_threshold,
         inc_scaling_factor=config_generate.increase_scaling_factor,
         verbosity=verbosity,
+        scale_bondlength=config_generate.scale_minimal_bondlength,
     )
     mol.charge, mol.uhf = set_random_charge(mol.ati, verbosity)
     mol.set_name_from_formula()
@@ -322,9 +323,9 @@ def check_composition():
 def generate_coordinates(
     at: np.ndarray,
     scaling: float,
-    dist_threshold: float,
     inc_scaling_factor: float = 1.3,
     verbosity: int = 1,
+    scale_bondlength: float = 0.75,
 ) -> tuple[np.ndarray, np.ndarray]:
     """
     Generate random coordinates for a molecule.
@@ -335,7 +336,7 @@ def generate_coordinates(
     xyz, ati = generate_random_coordinates(at)
     xyz = xyz * eff_scaling
     # do while check_distances is False
-    while not check_distances(xyz, dist_threshold):
+    while not check_distances(xyz, ati, scale_bondlength=scale_bondlength):
         if verbosity > 1:
             print(
                 f"Distance check failed. Increasing expansion factor by {inc_scaling_factor}..."
@@ -370,14 +371,17 @@ def generate_random_coordinates(at: np.ndarray) -> tuple[np.ndarray, np.ndarray]
     return xyz, ati
 
 
-def check_distances(xyz: np.ndarray, threshold: float) -> bool:
+def check_distances(xyz: np.ndarray, ati: np.ndarray, scale_bondlength: float) -> bool:
     """
     Check if the distances between atoms are larger than a threshold.
     """
     # go through the atoms dimension of the xyz array
     for i in range(xyz.shape[0] - 1):
         for j in range(i + 1, xyz.shape[0]):
             r = np.linalg.norm(xyz[i, :] - xyz[j, :])
-            if r < threshold:
+            sum_radii = get_cov_radii(ati[i], COV_RADII) + get_cov_radii(
+                ati[j], COV_RADII
+            )
+            if r < scale_bondlength * sum_radii:
                 return False
     return True

src/mindlessgen/prog/config.py

@@ -180,6 +180,7 @@ def __init__(self: GenerateConfig) -> None:
         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
 
     def get_identifier(self) -> str:
         return "generate"
@@ -383,6 +384,24 @@ def scale_vdw_radii(self, scale_vdw_radii: float):
             raise ValueError("Scale van der Waals radii should be greater than 0.")
         self._scale_vdw_radii = scale_vdw_radii
 
+    @property
+    def scale_minimal_bondlength(self):
+        """
+        Get the scaling factor for minimal bond length.
+        """
+        return self._scale_minimal_bondlength
+
+    @scale_minimal_bondlength.setter
+    def scale_minimal_bondlength(self, scale_minimal_bondlength: float):
+        """
+        Set the scaling factor for minimal bond length.
+        """
+        if not isinstance(scale_minimal_bondlength, float):
+            raise TypeError("Scale minimal bond length should be a float.")
+        if scale_minimal_bondlength <= 0:
+            raise ValueError("Scale minimal bond length should be greater than 0.")
+        self._scale_minimal_bondlength = scale_minimal_bondlength
+
 
 class RefineConfig(BaseConfig):
     """

test/test_generate/test_generate_molecule.py

@@ -218,49 +218,56 @@ def test_generate_coordinates() -> None:
 
 
 @pytest.mark.parametrize(
-    "xyz, threshold, expected, description",
+    "xyz, ati, scale_minimal_bondlength, expected, description",
     [
         (
             np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]]),
+            np.array([0, 0]),
             0.5,
             True,
             "Two atoms with distance greater than threshold (1.0 > 0.5)",
         ),
         (
             np.array([[0.0, 0.0, 0.0], [0.4, 0.0, 0.0]]),
-            0.5,
+            np.array([0, 0]),
+            0.75,
             False,
             "Two atoms with distance less than threshold (0.4 < 0.5)",
         ),
         (
             np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [2.0, 0.0, 0.0]]),
-            0.5,
+            np.array([0, 0, 0]),
+            1.0,
             True,
             "Three atoms in a line with distances greater than threshold",
         ),
         (
             np.array([[0.0, 0.0, 0.0], [0.4, 0.0, 0.0], [1.0, 0.0, 0.0]]),
-            0.5,
+            np.array([0, 0, 0]),
+            0.75,
             False,
             "Three atoms with one pair close together: distance between first two is less than threshold",
         ),
         (
             np.array([[0.0, 0.0, 0.0]]),
+            np.array([0]),
             0.5,
             True,
             "Single atom, no distances to compare",
         ),
         (
             np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]),
-            0.5,
+            np.array([0, 0]),
+            0.75,
             False,
             "Two atoms at identical positions: distance is zero, less than threshold",
         ),
         (
             np.array([[0.0, 0.0, 0.0], [1.0, 1.0, 1.0]]),
-            1.7320,
+            np.array([0, 0]),
+            2.70625,
             True,
-            "Two atoms with diagonal distance just above threshold (sqrt(3) ≈ 1.732)",
+            "Two atoms with diagonal distance just above threshold (sqrt(3) ≈ 1.732, 1.7323/0.64 = 2.70625)(0.64 = sum of covalent radii for H)",
         ),
     ],
     ids=[
@@ -273,8 +280,8 @@ def test_generate_coordinates() -> None:
         "diagonal_distance",
     ],
 )
-def test_check_distances(xyz, threshold, expected, description):
-    assert check_distances(xyz, threshold) == expected
+def test_check_distances(xyz, ati, scale_minimal_bondlength, expected, description):
+    assert check_distances(xyz, ati, scale_minimal_bondlength) == expected
 
 
 def test_generate_atom_list_min_larger_than_max(default_generate_config):