Fix random number generation by removing legacy random number generation (#78)

* fix random number generation by incorporating default_rng

Signed-off-by: Marcel Müller <marcel.mueller@thch.uni-bonn.de>

* update CHANGELOG.md

Signed-off-by: Marcel Müller <marcel.mueller@thch.uni-bonn.de>

---------

Signed-off-by: Marcel Müller <marcel.mueller@thch.uni-bonn.de>

Author: marcelmbn

CHANGELOG.md

@@ -20,6 +20,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - bug leading to `UnicodeDecodeError` when reading `xtb` output files
 - bug with all atom lists being initialized with a length of 102 instead of 103
 - inconsistent default values for the `mindlessgen.toml` and the `ConfigManager` class
+- legacy pseudo random number generation removed and replaced by `np.random.default_rng()` for avoiding interference with other packages
 
 ### Added
 - support for the novel "g-xTB" method (working title: GP3-xTB)

src/mindlessgen/molecules/generate_molecule.py

@@ -65,6 +65,9 @@ def generate_atom_list(cfg: GenerateConfig, verbosity: int = 1) -> np.ndarray:
     """
     Generate a random molecule with a random number of atoms.
     """
+    # initialize a default random number generator
+    rng = np.random.default_rng()
+
     # Define a new set of all elements that can be included
     set_all_elem = set(range(0, MAX_ELEM))
     if cfg.forbidden_elements:
@@ -166,17 +169,17 @@ def add_random(min_adds: int, max_adds: int, min_nat: int, max_nat: int):
         """
         Default random atom generation.
         """
-        numatoms_all = np.random.randint(
-            min_adds, max_adds
+        numatoms_all = rng.integers(
+            low=min_adds, high=max_adds
         )  # with range(1, 7) -> mean value: 3.5
         for _ in range(numatoms_all):
             # Define the atom type to be added via a random choice from the set of valid elements
-            ati = np.random.choice(list(valid_elems))
+            ati = rng.choice(list(valid_elems))
             if verbosity > 1:
                 print(f"Adding atom type {ati}...")
             # Add a random number of atoms of the defined type
-            natoms[ati] = natoms[ati] + np.random.randint(
-                min_nat, max_nat
+            natoms[ati] = natoms[ati] + rng.integers(
+                low=min_nat, high=max_nat
             )  # with range(0, 3) -> mean value: 1
             # max value of this section with commented settings: 12
 
@@ -193,11 +196,11 @@ def add_organic(num_adds: int, min_nat: int, max_nat: int) -> None:
             return
         for _ in range(num_adds):  # with range(5) -> mean value 1.5
             # go through the elements B to F (4-9 in 0-based indexing)
-            ati = np.random.choice(valid_organic)
+            ati = rng.choice(valid_organic)
             if verbosity > 1:
                 print(f"Adding atom type {ati}...")
-            natoms[ati] = natoms[ati] + np.random.randint(
-                min_nat, max_nat
+            natoms[ati] = natoms[ati] + rng.integers(
+                low=min_nat, high=max_nat
             )  # with range(0, 3) -> mean value: 1
             # max value of this section with commented settings: 8
 
@@ -243,7 +246,7 @@ def add_hydrogen():
         # If no H is included, add H atoms
         if natoms[0] == 0:
             nat = np.sum(natoms)
-            randint = np.random.rand()
+            randint = rng.random()
             j = 1 + round(randint * nat * 1.2)
             natoms[0] = natoms[0] + j
             # Example: For 5 atoms at this point,
@@ -257,7 +260,7 @@ def check_min_max_atoms():
                     f"Minimal number of atoms: {cfg.min_num_atoms}; "
                     + f"Actual number of atoms: {np.sum(natoms)}.\nAdding atoms..."
                 )
-            ati = np.random.choice(list(valid_elems))
+            ati = rng.choice(list(valid_elems))
             max_limit = cfg.element_composition.get(ati, (None, None))[1]
             if max_limit is not None and natoms[ati] >= max_limit:
                 continue
@@ -285,7 +288,7 @@ def check_min_max_atoms():
             # randomly select an atom type from the list, thereby weighting the selection
             # for reduction by the current occurrence
             # generate a random number between 0 and the number of atoms in the list
-            random_index = np.random.randint(len(atom_list))
+            random_index = rng.integers(0, len(atom_list))
             i = atom_list[int(random_index)]
             if natoms[i] > 0:
                 min_limit = cfg.element_composition.get(i, (None, None))[0]
@@ -361,13 +364,15 @@ def generate_random_coordinates(at: np.ndarray) -> tuple[np.ndarray, np.ndarray]
     atilist: list[int] = []
     xyz = np.zeros((sum(at), 3))
     numatoms = 0
+    rng = np.random.default_rng()
     for elem, count in enumerate(at):
         for m in range(count):
             # different rules for hydrogen
             if elem == 0:
-                xyz[numatoms + m, :] = np.random.rand(3) * 3 - 1.5
+                xyz[numatoms + m, :] = rng.random(3) * 3 - 1.5
             else:
-                xyz[numatoms + m, :] = np.random.rand(3) * 2 - 1
+                xyz[numatoms + m, :] = rng.random(3) * 2 - 1
+
             atilist.append(elem)
 
         numatoms += count

src/mindlessgen/molecules/miscellaneous.py

@@ -63,7 +63,8 @@ def set_random_charge(ati: np.ndarray, verbosity: int = 1) -> tuple[int, int]:
             iseven = True
         # if the number of electrons is even, the charge is -2, 0, or 2
         # if the number of electrons is odd, the charge is -1, 1
-        randint = np.random.rand()
+        rng = np.random.default_rng()
+        randint = rng.random()
         if iseven:
             if randint < 1 / 3:
                 charge = -2

src/mindlessgen/molecules/molecule.py

@@ -156,6 +156,8 @@ def __init__(self, name: str = ""):
         self._xyz: np.ndarray = np.array([], dtype=float)
         self._ati: np.ndarray = np.array([], dtype=int)
 
+        self.rng = np.random.default_rng()
+
     def __str__(self) -> str:
         """
         Return a user-friendly string representation of the molecule.
@@ -637,5 +639,5 @@ def set_name_from_formula(self) -> None:
 
         molname = self.sum_formula()
         # add a random hash to the name
-        hashname = hashlib.sha256(np.random.bytes(32)).hexdigest()[:6]
+        hashname = hashlib.sha256(self.rng.bytes(32)).hexdigest()[:6]
         self.name = f"{molname}_{hashname}"