Merge pull request #326 from phonopy/random-disp

(WIP) Extend feature of random displacements

doc/changelog.md

@@ -2,6 +2,11 @@
 
 # Change Log
 
+##
+
+- `np.random.randn` is replaced by `np.random.Generator.standard_normal`. Note
+  that by this chang, random seed became incompatible.
+
 ## Dec-4-2023: Version 2.21.0
 
 - Maintenance release.

phonopy/api_phonopy.py

@@ -899,7 +899,7 @@ class Phonopy:
 
     def generate_displacements(
         self,
-        distance: float = 0.01,
+        distance: Optional[float] = None,
         is_plusminus: Union[str, bool] = "auto",
         is_diagonal: bool = True,
         is_trigonal: bool = False,
@@ -908,6 +908,8 @@ class Phonopy:
         temperature: Optional[float] = None,
         cutoff_frequency: Optional[float] = None,
         max_distance: Optional[float] = None,
+        is_random_distance: bool = False,
+        min_distance: Optional[float] = None,
     ):
         """Generate displacement dataset.
 
@@ -926,7 +928,9 @@ class Phonopy:
         ----------
         distance : float, optional
             Displacement distance. Unit is the same as that used for crystal
-            structure. Default is 0.01.
+            structure. Default is 0.01. For random direction and distance
+            displacements generation, this value is used when `max_distance` is
+            unspecified.
         is_plusminus : 'auto', True, or False, optional
             For each atom, displacement of one direction (False), both
             direction, i.e., one directiona and its opposite direction (True),
@@ -944,7 +948,7 @@ class Phonopy:
             'distance' parameter, i.e., all atoms in supercell are displaced
             with the same displacement distance in direct space. Default is
             None.
-        random_seed : 32bit unsigned int or None, optional
+        random_seed : int or None, optional
             Random seed for random displacements generation. Default is None.
         temperature : float or None, optional
             With given temperature, random displacements at temperature is
@@ -961,7 +965,17 @@ class Phonopy:
             In random displacements generation from canonical ensemble of
             harmonic phonons, displacements larger than max distance are
             renormalized to the max distance, i.e., a disptalcement d is shorten
-            by d -> d / |d| * max_distance if |d| > max_distance.
+            by d -> d / |d| * max_distance if |d| > max_distance. In random
+            direction and distance displacements generation, this value is is
+            specified.
+        is_random_distance : bool, optional
+            Random direction displacements are generated also with random
+            amplitudes. The maximum value is defined by `distance` and minimum
+            value is given by `min_distance`. Default is False.
+        min_distance : float or None, optional
+            In random direction displacements generation with random distance
+            (`is_random_distance=True`), the minimum distance is given by this
+            value.
 
         """
         if number_of_snapshots is not None and number_of_snapshots > 0:
@@ -972,12 +986,21 @@ class Phonopy:
                 _random_seed = None
                 displacement_dataset = {}
             if temperature is None:
+                if max_distance is None:
+                    if distance is None:
+                        _distance = 0.01
+                    else:
+                        _distance = distance
+                else:
+                    _distance = max_distance
                 d = get_random_displacements_dataset(
                     number_of_snapshots,
-                    distance,
                     len(self._supercell),
+                    _distance,
                     random_seed=_random_seed,
                     is_plusminus=(is_plusminus is True),
+                    is_random_distance=is_random_distance,
+                    min_distance=min_distance,
                 )
                 displacement_dataset["displacements"] = d
             else:
@@ -992,6 +1015,10 @@ class Phonopy:
                 )
                 displacement_dataset["displacements"] = d
         else:
+            if distance is None:
+                _distance = 0.01
+            else:
+                _distance = distance
             displacement_directions = get_least_displacements(
                 self._symmetry,
                 is_plusminus=is_plusminus,
@@ -1000,7 +1027,7 @@ class Phonopy:
                 log_level=self._log_level,
             )
             displacement_dataset = directions_to_displacement_dataset(
-                displacement_directions, distance, self._supercell
+                displacement_directions, _distance, self._supercell
             )
         self.dataset = displacement_dataset
 

phonopy/harmonic/displacement.py

@@ -245,16 +245,70 @@ def _determinant(a, b, c):
 
 def get_random_displacements_dataset(
     num_supercells: int,
-    distance: float,
     num_atoms: int,
+    distance: float,
     random_seed: Optional[int] = None,
     is_plusminus: bool = False,
+    is_random_distance: bool = False,
+    min_distance: Optional[float] = None,
 ) -> np.ndarray:
-    """Return random displacements at constant displacement distance."""
-    disps = (
-        _get_random_directions(num_atoms * num_supercells, random_seed=random_seed)
-        * distance
-    )
+    """Return random displacements at constant displacement distance.
+
+    num_supercells : int
+        Number of snapshots of supercells with random displacements. Random
+        displacements are generated displacing all atoms in random directions
+        with a fixed displacement distance specified by 'distance' parameter,
+        i.e., all atoms in supercell are displaced with the same displacement
+        distance in direct space.
+    num_atoms : int
+        Number of atoms in supercell.
+    distance : float
+        Displacement distance. Unit is the same as that used for crystal
+        structure.
+    random_seed : int or None, optional
+        Random seed for random displacements generation. Default is None.
+    is_plusminus : True, or False, optional
+        In addition to sets of usual random displacements for supercell, sets
+        of the opposite displacements for supercell are concatenated.
+        Therefore, total number of sets of displacements is `2 *
+        num_supercells`. Default is False.
+    is_random_distance : bool, optional
+        Random direction displacements are generated also with random
+        amplitudes. The maximum value is defined by `distance` and minimum value
+        is given by `min_distance`. Default is False. Random distance is given
+        by `sqrt(random(distance - min_distance) + min_distance)`.
+    min_distance : float or None, optional
+        In random direction displacements generation with random distance
+        (`is_random_distance=True`), the minimum distance is given by this
+        value.
+
+    """
+    if is_random_distance:
+        if min_distance is None:
+            _min_distance = 0.0
+        else:
+            _min_distance = min_distance
+
+    if np.issubdtype(type(random_seed), np.integer):
+        rng = np.random.default_rng(seed=random_seed)
+    else:
+        rng = np.random.default_rng()
+
+    if is_random_distance:
+        if distance < _min_distance:
+            raise RuntimeError(
+                "Random displacements generation failed. min_distance is too large."
+            )
+        directions = _get_random_directions(num_atoms * num_supercells, rng)
+        rand_dists = np.array([])
+        while len(rand_dists) < num_atoms * num_supercells:
+            rd = np.sqrt(rng.random(num_atoms * num_supercells)) * distance
+            rand_dists = np.r_[rand_dists, rd[rd > _min_distance]]
+        disps = rand_dists[: num_atoms * num_supercells, None] * directions
+    else:
+        directions = _get_random_directions(num_atoms * num_supercells, rng)
+        disps = directions * distance
+
     supercell_disps = np.array(
         disps.reshape(num_supercells, num_atoms, 3), dtype="double", order="C"
     )
@@ -267,15 +321,9 @@ def get_random_displacements_dataset(
     return supercell_disps
 
 
-def _get_random_directions(num_atoms, random_seed=None):
+def _get_random_directions(num_atoms: int, rng: np.random.Generator) -> np.ndarray:
     """Return random directions in sphere with radius 1."""
-    if (
-        np.issubdtype(type(random_seed), np.integer)
-        and random_seed >= 0
-        and random_seed < 2**32
-    ):
-        np.random.seed(random_seed)
-
-    xy = np.random.randn(3, num_atoms)
-    r = np.sqrt((xy**2).sum(axis=0))
-    return (xy / r).T
+    xy = rng.standard_normal(size=(3, num_atoms * 2))
+    r = np.linalg.norm(xy, axis=0)
+    condition = r > 1e-10
+    return (xy[:, condition][:, :num_atoms] / r[condition][:num_atoms]).T

phonopy/phonon/random_displacements.py

@@ -242,33 +242,42 @@ class RandomDisplacements:
         number_of_snapshots : int
             Number of snapshots to be generated.
         random_seed : int or None, optional
-            Random seed passed to np.random.seed. Default is None. Integer
-            number has to be positive.
+            Random seed passed to np.random.default_rng(seed). Default is None.
         randn : tuple
-            (randn_ii, randn_ij).
-            Used for testing purpose for the fixed random numbers of
-            np.random.normal that can depends on system.
+            (randn_ii, randn_ij). Used for testing purpose for the fixed random
+            numbers of random.Generator.standard_normal that can depends on
+            system.
 
         """
-        np.random.seed(seed=random_seed)
+        if np.issubdtype(type(random_seed), np.integer):
+            rng = np.random.default_rng(seed=random_seed)
+        else:
+            rng = np.random.default_rng()
 
         N = len(self._comm_points)
 
         # This randn is used only for testing purpose.
         if randn is None:
-            randn_ii = None
-            randn_ij = None
+            shape = (
+                len(self._eigvals_ii),
+                number_of_snapshots,
+                len(self._eigvals_ii[0]),
+            )
+            randn_ii = rng.standard_normal(size=shape)
+            shape = (
+                len(self._eigvals_ij),
+                2,
+                number_of_snapshots,
+                len(self._eigvals_ij[0]),
+            )
+            randn_ij = rng.standard_normal(size=shape)
         else:
             randn_ii = randn[0]
             randn_ij = randn[1]
 
-        u_ii, self._conditions_ii = self._solve_ii(
-            T, number_of_snapshots, randn=randn_ii
-        )
+        u_ii, self._conditions_ii = self._solve_ii(T, number_of_snapshots, randn_ii)
         if self._ij:
-            u_ij, self._conditions_ij = self._solve_ij(
-                T, number_of_snapshots, randn=randn_ij
-            )
+            u_ij, self._conditions_ij = self._solve_ij(T, number_of_snapshots, randn_ij)
         else:
             u_ij = 0
             self._conditions_ij = None
@@ -481,7 +490,12 @@ class RandomDisplacements:
         self._comm_points = get_commensurate_points_in_integers(smat)
         self._ii, self._ij = categorize_commensurate_points(self._comm_points)
 
-    def _solve_ii(self, T, number_of_snapshots, randn=None):
+    def _solve_ii(
+        self,
+        T: float,
+        number_of_snapshots: int,
+        randn: np.ndarray,
+    ):
         """Solve ii terms.
 
         randn parameter is used for the test.
@@ -490,14 +504,9 @@ class RandomDisplacements:
         natom = len(self._dynmat.supercell)
         u = np.zeros((number_of_snapshots, natom, 3), dtype="double")
 
-        shape = (len(self._eigvals_ii), number_of_snapshots, len(self._eigvals_ii[0]))
-        if randn is None:
-            _randn = np.random.normal(size=shape)
-        else:
-            _randn = randn
         sigmas, conditions = self._get_sigma(self._eigvals_ii, T)
         for norm_dist, sigma, eigvecs, phase in zip(
-            _randn, sigmas, self._eigvecs_ii, self._phase_ii
+            randn, sigmas, self._eigvecs_ii, self._phase_ii
         ):
             u_red = np.dot(norm_dist * sigma, eigvecs.T).reshape(
                 number_of_snapshots, -1, 3
@@ -508,7 +517,12 @@ class RandomDisplacements:
 
         return u, conditions
 
-    def _solve_ij(self, T, number_of_snapshots, randn=None):
+    def _solve_ij(
+        self,
+        T: float,
+        number_of_snapshots: int,
+        randn: Optional[np.ndarray] = None,
+    ):
         """Solve ij terms.
 
         randn parameter is used for the test.
@@ -516,19 +530,9 @@ class RandomDisplacements:
         """
         natom = len(self._dynmat.supercell)
         u = np.zeros((number_of_snapshots, natom, 3), dtype="double")
-        shape = (
-            len(self._eigvals_ij),
-            2,
-            number_of_snapshots,
-            len(self._eigvals_ij[0]),
-        )
-        if randn is None:
-            _randn = np.random.normal(size=shape)
-        else:
-            _randn = randn
         sigmas, conditions = self._get_sigma(self._eigvals_ij, T)
         for norm_dist, sigma, eigvecs, phase in zip(
-            _randn, sigmas, self._eigvecs_ij, self._phase_ij
+            randn, sigmas, self._eigvecs_ij, self._phase_ij
         ):
             u_red = np.dot(norm_dist * sigma, eigvecs.T).reshape(
                 2, number_of_snapshots, -1, 3

test/harmonic/test_displacements.py

@@ -1,7 +1,10 @@
 """Tests for displacements."""
+import itertools
 from copy import deepcopy
+from typing import Optional
 
 import numpy as np
+import pytest
 
 from phonopy import Phonopy
 
@@ -50,15 +53,6 @@ def test_sno2(ph_sno2: Phonopy):
 def test_tio2(ph_tio2: Phonopy):
     """Test displacements of TiO2."""
     dataset = deepcopy(ph_tio2.dataset)
-    disp_ref = [
-        [0, 0.01, 0.0, 0.0],
-        [0, 0.0, 0.01, 0.0],
-        [0, 0.0, 0.0, 0.01],
-        [0, 0.0, 0.0, -0.01],
-        [72, 0.01, 0.0, 0.0],
-        [72, 0.0, 0.0, 0.01],
-    ]
-    np.testing.assert_allclose(ph_tio2.displacements, disp_ref, atol=1e-8)
     ph_tio2.generate_displacements()
     disp_gen = [
         [0, 0.0060687317141537135, 0.0060687317141537135, 0.0051323474905008],
@@ -70,21 +64,113 @@ def test_tio2(ph_tio2: Phonopy):
     ph_tio2.dataset = dataset
 
 
-def test_tio2_random_disp(ph_tio2: Phonopy):
-    """Test random displacements of TiO2."""
+@pytest.mark.parametrize(
+    "is_plusminus,distance", itertools.product([False, True], [None, 0.03])
+)
+def test_tio2_random_disp(
+    ph_tio2: Phonopy, is_plusminus: bool, distance: Optional[float]
+):
+    """Test random displacements of TiO2.
+
+    Currently default displacement distance = 0.01.
+
+    """
     dataset = deepcopy(ph_tio2.dataset)
-    disp_ref = [
-        [0, 0.01, 0.0, 0.0],
-        [0, 0.0, 0.01, 0.0],
-        [0, 0.0, 0.0, 0.01],
-        [0, 0.0, 0.0, -0.01],
-        [72, 0.01, 0.0, 0.0],
-        [72, 0.0, 0.0, 0.01],
-    ]
-    np.testing.assert_allclose(ph_tio2.displacements, disp_ref, atol=1e-8)
-    ph_tio2.generate_displacements(number_of_snapshots=4, distance=0.03)
+    ph_tio2.generate_displacements(
+        number_of_snapshots=4, distance=distance, is_plusminus=is_plusminus
+    )
     d = ph_tio2.displacements
-    np.testing.assert_allclose(np.linalg.norm(d, axis=2).ravel(), 0.03, atol=1e-8)
+    assert len(d) == 4 * (is_plusminus + 1)
+    if distance is None:
+        np.testing.assert_allclose(np.linalg.norm(d, axis=2).ravel(), 0.01, atol=1e-8)
+    else:
+        np.testing.assert_allclose(np.linalg.norm(d, axis=2).ravel(), 0.03, atol=1e-8)
+    if is_plusminus:
+        np.testing.assert_allclose(d[:4], -d[4:], atol=1e-8)
+    ph_tio2.dataset = dataset
+
+
+@pytest.mark.parametrize("min_distance", [None, 0.05, 0.2])
+def test_tio2_random_disp_with_random_dist(
+    ph_tio2: Phonopy, min_distance: Optional[float]
+):
+    """Test random displacements with random distance of TiO2."""
+    dataset = deepcopy(ph_tio2.dataset)
+
+    if min_distance is not None and min_distance > 0.1:
+        with pytest.raises(RuntimeError):
+            ph_tio2.generate_displacements(
+                number_of_snapshots=1,
+                distance=0.1,
+                is_random_distance=True,
+                min_distance=min_distance,
+            )
+    else:
+        n_snapshots = 100
+        ph_tio2.generate_displacements(
+            number_of_snapshots=n_snapshots,
+            distance=0.1,
+            is_random_distance=True,
+            min_distance=min_distance,
+        )
+        d = ph_tio2.displacements
+        assert len(d) == n_snapshots
+        dists = np.linalg.norm(d, axis=2).ravel()
+        assert (dists < 0.1 + 1e-8).all()
+        if min_distance is not None:
+            assert (dists > min_distance - 1e-8).all()
+
+    ph_tio2.dataset = dataset
+
+
+@pytest.mark.parametrize("max_distance", [None, 0.1])
+def test_tio2_random_disp_with_random_dist_defualt(
+    ph_tio2: Phonopy, max_distance: Optional[float]
+):
+    """Test random displacements with random distance of TiO2.
+
+    Combination of default distance and max_distance.
+
+    """
+    dataset = deepcopy(ph_tio2.dataset)
+
+    n_snapshots = 100
+    ph_tio2.generate_displacements(
+        number_of_snapshots=n_snapshots,
+        is_random_distance=True,
+        max_distance=max_distance,
+    )
+    d = ph_tio2.displacements
+    assert len(d) == n_snapshots
+    dists = np.linalg.norm(d, axis=2).ravel()
+    if max_distance is None:
+        assert (dists < 0.01 + 1e-8).all()
+    else:
+        assert (dists < max_distance + 1e-8).all()
+
+    ph_tio2.dataset = dataset
+
+
+def test_tio2_random_disp_with_random_max_distance(ph_tio2: Phonopy):
+    """Test random displacements with random distance of TiO2.
+
+    Combination of distance and max_distance parameters.
+
+    """
+    dataset = deepcopy(ph_tio2.dataset)
+
+    n_snapshots = 100
+    ph_tio2.generate_displacements(
+        number_of_snapshots=n_snapshots,
+        distance=0.01,
+        max_distance=0.1,
+        is_random_distance=True,
+    )
+    d = ph_tio2.displacements
+    assert len(d) == n_snapshots
+    dists = np.linalg.norm(d, axis=2).ravel()
+    assert (dists < 0.1 + 1e-8).all()
+
     ph_tio2.dataset = dataset
 
 

test/phonon/test_random_displacements.py

@@ -736,9 +736,9 @@ def _test_random_displacements(ph, answer, ntimes=30, d_max=1, nbins=51):
         hist[0] += h_1
         hist[1] += h_2
         argmaxs = (np.argmax(hist[0]), np.argmax(hist[1]))
-        # print(hist[0])
-        # print(hist[1])
-        # print(i, argmaxs)
+        print(hist[0])
+        print(hist[1])
+        print(i, argmaxs)
         if i >= ntimes and argmaxs == answer:
             break
     return argmaxs
@@ -863,7 +863,10 @@ def test_tio2_random_disp_plusminus(ph_tio2: Phonopy, is_plusminus: bool):
     ]
     np.testing.assert_allclose(ph_tio2.displacements, disp_ref, atol=1e-8)
     ph_tio2.generate_displacements(
-        number_of_snapshots=4, distance=0.03, is_plusminus=is_plusminus, temperature=300
+        number_of_snapshots=4,
+        distance=0.03,
+        is_plusminus=is_plusminus,
+        temperature=300,
     )
     d = ph_tio2.displacements
     if is_plusminus: