Merge pull request #51 from jonathan-schoeps/dev/dist_covalent_radii

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
@@ -17,6 +18,9 @@ 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)
 
+### Breaking Changes
+- Removal of the `dist_threshold` flag and in the `-toml` file.
+
 ## [0.4.0] - 2024-09-19
 ### Changed
 - Default file name of `.xyz` file contains prefix `mlm_`

mindlessgen.toml

@@ -25,12 +25,12 @@ min_num_atoms = 5
 max_num_atoms = 10
 # > Initial coordinate scaling factor. Options: <float>
 init_scaling = 3.0
-# > Increase in the coordinate scaling factor per trial after dist_threshold was not met. Options: <float>
+# > Increase in the coordinate scaling factor per trial after check_distance was not met. Options: <float>
 increase_scaling_factor = 1.3
-# > Distance threshold for the inital, randomly generated coordinates. Options: <float>
-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(
@@ -119,12 +125,6 @@ def cli_parser(argv: Sequence[str] | None = None) -> dict:
         help="Factor with which the coordinate scaling factor is increased "
         + "after a failed attempt.",
     )
-    parser.add_argument(
-        "--dist-threshold",
-        type=float,
-        required=False,
-        help="Distance threshold for generating coordinates.",
-    )
     parser.add_argument(
         "--element-composition",
         type=str,
@@ -269,10 +269,10 @@ def cli_parser(argv: Sequence[str] | None = None) -> dict:
         "max_num_atoms": args_dict["max_num_atoms"],
         "init_coord_scaling": args_dict["init_coord_scaling"],
         "increase_scaling_factor": args_dict["increase_scaling_factor"],
-        "dist_threshold": args_dict["dist_threshold"],
         "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

@@ -175,11 +175,11 @@ def __init__(self: GenerateConfig) -> None:
         self._min_num_atoms: int = 2
         self._max_num_atoms: int = 100
         self._init_coord_scaling: float = 3.0
-        self._dist_threshold: float = 1.2
         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
 
     def get_identifier(self) -> str:
         return "generate"
@@ -238,24 +238,6 @@ def init_coord_scaling(self, init_coord_scaling: float):
             raise ValueError("Initial coordinate scaling should be greater than 0.")
         self._init_coord_scaling = init_coord_scaling
 
-    @property
-    def dist_threshold(self):
-        """
-        Get the distance threshold.
-        """
-        return self._dist_threshold
-
-    @dist_threshold.setter
-    def dist_threshold(self, dist_threshold: float):
-        """
-        Set the distance threshold.
-        """
-        if not isinstance(dist_threshold, float):
-            raise TypeError("Distance threshold should be a float.")
-        if dist_threshold <= 0:
-            raise ValueError("Distance threshold should be greater than 0.")
-        self._dist_threshold = dist_threshold
-
     @property
     def increase_scaling_factor(self):
         """
@@ -383,6 +365,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/fixtures/example_config.toml

@@ -23,10 +23,8 @@ min_num_atoms = 2
 max_num_atoms = 100
 # > Initial coordinate scaling factor. Options: <float>
 init_scaling = 3.0
-# > Increase in the coordinate scaling factor per trial after dist_threshold was not met. Options: <float>
+# > Increase in the coordinate scaling factor per trial after check_distance was not met. Options: <float>
 increase_scaling_factor = 1.3
-# > Distance threshold for the inital, randomly generated coordinates. Options: <float>
-dist_threshold = 1.2
 # > 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.

test/test_config/test_config_set_attributes.py

@@ -65,8 +65,6 @@ def test_general_config_default_values(property_name, initial_value):
         ("max_num_atoms", 80, None, TypeError),
         ("init_coord_scaling", 1.0, -0.5, ValueError),
         ("init_coord_scaling", 1.0, "1.0", TypeError),
-        ("dist_threshold", 1.5, -1.0, ValueError),
-        ("dist_threshold", 1.5, "1.5", TypeError),
         ("increase_scaling_factor", 1.1, 0.0, ValueError),
         ("increase_scaling_factor", 1.1, "1.1", TypeError),
     ],
@@ -115,7 +113,6 @@ def test_generate_config_element_composition(
         ("min_num_atoms", 2),
         ("max_num_atoms", 100),
         ("init_coord_scaling", 3.0),
-        ("dist_threshold", 1.2),
         ("increase_scaling_factor", 1.3),
         ("element_composition", {}),
         ("forbidden_elements", None),

test/test_config/test_load_from_toml.py

@@ -36,7 +36,6 @@ def test_load_generate_config(config_manager):
     assert config_manager.generate.max_num_atoms == 100
     assert config_manager.generate.init_coord_scaling == 3.0
     assert config_manager.generate.increase_scaling_factor == 1.3
-    assert config_manager.generate.dist_threshold == 1.2
     assert config_manager.generate.element_composition == {
         5: (2, 10),  # Carbon (C)
         0: (10, 20),  # Hydrogen (H)

test/test_generate/test_generate_molecule.py

@@ -218,49 +218,77 @@ 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)",
+            "Two Hydrogenes 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)",
+            "Two Hydrogenes 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",
+            "Three Hydrogenes 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",
+            "Three Hydrogenes 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",
+            "Single Hydrogene, 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",
+            "Two Hydrogenes 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 Hydrogenes with diagonal distance just above threshold (sqrt(3) ≈ 1.732, 1.7323/0.64 = 2.70625)(0.64 = sum of covalent radii for H)",
+        ),
+        (
+            np.array([[0.0, 0.0, 0.0], [2.3, 0.0, 0.0]]),
+            np.array([18, 8]),
+            0.9,
             True,
-            "Two atoms with diagonal distance just above threshold (sqrt(3) ≈ 1.732)",
+            "Potassium plus flourine with distance greater than threshold (r = 2.3, scaled_minimal_bondlength = 2.16)",
+        ),
+        (
+            np.array([[0.0, 0.0, 0.0], [2.3, 0.0, 0.0]]),
+            np.array([18, 8]),
+            2.0,
+            False,
+            "Potassium plus flourine with distance less than threshold (r = 2.3, scaled_minimal_bondlength = 4.8)",
+        ),
+        (
+            np.array([[0.0, 0.0, 0.0], [2.3, 0.0, 0.0]]),
+            np.array([18, 16]),
+            0.9,
+            False,
+            "Potassium plus chlorine with distance less than threshold (r = 2.3, scaled_minimal_bondlength = 2.61)",
         ),
     ],
     ids=[
@@ -271,10 +299,13 @@ def test_generate_coordinates() -> None:
         "single_atom",
         "two_identical",
         "diagonal_distance",
+        "different_elements_apart",
+        "different_elements_scaled_close",
+        "different_elements_close",
     ],
 )
-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):