Made cutoff_force_constants accept store_in_dense_array

phonopy/harmonic/force_constants.py

@@ -129,18 +129,33 @@ def cutoff_force_constants(force_constants,
                            primitive,
                            cutoff_radius,
                            symprec=1e-5):
+    """Set zero to force constants outside of cutoff distance.
+
+    Note
+    ----
+    `force_constants` is overwritten.
+
+    """
     fc_shape = force_constants.shape
+
     if fc_shape[0] == fc_shape[1]:
-        svecs, _ = get_smallest_vectors(supercell.get_cell(),
-                                        supercell.get_scaled_positions(),
-                                        supercell.get_scaled_positions(),
-                                        symprec=symprec)
-        min_distances = np.sqrt(np.sum(
-            np.dot(svecs[:, :, 0, :], supercell.get_cell()) ** 2, axis=-1))
+        svecs, multi = get_smallest_vectors(
+            supercell.cell,
+            supercell.scaled_positions,
+            supercell.scaled_positions,
+            symprec=symprec,
+            store_in_dense_array=primitive.store_in_dense_array)
+        lattice = supercell.cell
     else:
-        svecs, _ = primitive.get_smallest_vectors()
-        min_distances = np.sqrt(np.sum(
-            np.dot(svecs[:, :, 0, :], primitive.get_cell()) ** 2, axis=-1))
+        svecs, multi = primitive.get_smallest_vectors()
+        lattice = primitive.cell
+
+    if primitive.store_in_dense_array:
+        _svecs = svecs[multi[:, :, 1]]
+    else:
+        _svecs = svecs[:, :, 0, :]
+
+    min_distances = np.sqrt(np.sum(np.dot(_svecs, lattice) ** 2, axis=-1))
 
     for i in range(fc_shape[0]):
         for j in range(fc_shape[1]):

phonopy/structure/cells.py

@@ -52,13 +52,13 @@ def get_supercell(unitcell, supercell_matrix, is_old_style=True, symprec=1e-5):
 def get_primitive(supercell,
                   primitive_frame,
                   symprec=1e-5,
-                  store_dense_vectors=False,
+                  store_in_dense_array=False,
                   positions_to_reorder=None):
     """Create primitive cell."""
     return Primitive(supercell,
                      primitive_frame,
                      symprec=symprec,
-                     store_dense_vectors=store_dense_vectors,
+                     store_in_dense_array=store_in_dense_array,
                      positions_to_reorder=positions_to_reorder)
 
 
@@ -401,6 +401,7 @@ class Primitive(PhonopyAtoms):
     s2p_map : ndarray
     p2p_map : dict
     atomic_permutations : ndarray
+    store_in_dense_array : bool
 
     """
 
@@ -408,7 +409,7 @@ class Primitive(PhonopyAtoms):
                  supercell,
                  primitive_matrix,
                  symprec=1e-5,
-                 store_dense_vectors=False,
+                 store_in_dense_array=False,
                  positions_to_reorder=None):
         """Init method.
 
@@ -423,7 +424,7 @@ class Primitive(PhonopyAtoms):
             shape=(3,3)
         symprec : float, optional, default=1e-5
             Tolerance to find overlapping atoms in primitive cell.
-        store_dense_vectors : bool, optional, default=False
+        store_in_dense_array : bool, optional, default=False
             Shortest vectors are stored in a dense array. See `ShortestPairs`.
         positions_to_reorder : array_like, optional, default=None
             If atomic positions in a created primitive cell is known and
@@ -435,7 +436,7 @@ class Primitive(PhonopyAtoms):
         self._primitive_matrix = np.array(
             primitive_matrix, dtype='double', order='C')
         self._symprec = symprec
-        self._store_dense_vectors = store_dense_vectors
+        self._store_in_dense_array = store_in_dense_array
         self._p2s_map = None
         self._s2p_map = None
         self._p2p_map = None
@@ -566,6 +567,11 @@ class Primitive(PhonopyAtoms):
                       DeprecationWarning)
         return self.atomic_permutations
 
+    @property
+    def store_in_dense_array(self):
+        """Return whether shortest vectors are stored in dense array or not."""
+        return self._store_in_dense_array
+
     def _run(self, supercell, positions_to_reorder=None):
         self._p2s_map = self._create_primitive_cell(
             supercell, positions_to_reorder=positions_to_reorder)
@@ -647,7 +653,7 @@ class Primitive(PhonopyAtoms):
             supercell_bases,
             supercell_pos,
             primitive_pos,
-            store_dense_vectors=self._store_dense_vectors,
+            store_in_dense_array=self._store_in_dense_array,
             symprec=self._symprec)
         trans_mat_float = np.dot(
             supercell_bases, np.linalg.inv(primitive_bases))
@@ -891,7 +897,7 @@ def get_reduced_bases(lattice,
 def get_smallest_vectors(supercell_bases,
                          supercell_pos,
                          primitive_pos,
-                         store_dense_vectors=False,
+                         store_in_dense_array=False,
                          symprec=1e-5):
     """Return shortest vectors and multiplicities.
 
@@ -901,7 +907,7 @@ def get_smallest_vectors(supercell_bases,
     spairs = ShortestPairs(supercell_bases,
                            supercell_pos,
                            primitive_pos,
-                           store_dense_vectors=store_dense_vectors,
+                           store_in_dense_array=store_in_dense_array,
                            symprec=symprec)
     return spairs.shortest_vectors, spairs.multiplicities
 
@@ -939,7 +945,7 @@ class ShortestPairs(object):
                  supercell_bases,
                  supercell_pos,
                  primitive_pos,
-                 store_dense_vectors=False,
+                 store_in_dense_array=False,
                  symprec=1e-5):
         """Init method.
 
@@ -955,7 +961,7 @@ class ShortestPairs(object):
             Atomic positions in fractional coordinates of supercell. Note that
             not in fractional coodinates of primitive cell.  dtype='double',
             shape=(size_prim, 3)
-        store_dense_vectors_: bool, optional
+        store_in_dense_array_: bool, optional
             ``shortest_vectors`` are stored in the dense data structure.
             Default is False.
         symprec : float, optional
@@ -967,7 +973,7 @@ class ShortestPairs(object):
         self._primitive_pos = primitive_pos
         self._symprec = symprec
 
-        if store_dense_vectors:
+        if store_in_dense_array:
             svecs, multi = self._run_dense()
             self._smallest_vectors = svecs
             self._multiplicities = multi

test/harmonic/test_force_constants.py

@@ -1,4 +1,7 @@
+import pytest
 import numpy as np
+from phonopy.harmonic.force_constants import cutoff_force_constants
+from phonopy.structure.cells import get_primitive
 
 
 fc_1_10_ref = [-0.037549, 0.000000, 0.000000,
@@ -31,3 +34,47 @@ def test_fc_compact_fcsym(ph_nacl_compact_fcsym):
     # print("".join(["%f, " % v for v in fc_1_10.ravel()]))
     np.testing.assert_allclose(
         fc_1_10.ravel(), fc_1_10_compact_fcsym_ref, atol=1e-5)
+
+
+@pytest.mark.parametrize("is_compact", [False, True])
+def test_fc_cutoff_radius(ph_nacl, ph_nacl_compact_fcsym, is_compact):
+    if is_compact:
+        ph = ph_nacl_compact_fcsym
+    else:
+        ph = ph_nacl
+
+    # Need restore fc because fc are overwritten.
+    fc_orig = ph.force_constants.copy()
+    ph.set_force_constants_zero_with_radius(4.0)
+    changed = (np.abs(fc_orig - ph.force_constants) > 1e-8)
+    ph.force_constants = fc_orig
+
+    if is_compact:
+        assert np.sum(changed) == 534
+    else:
+        assert np.sum(changed) == 17088
+
+
+@pytest.mark.parametrize("is_compact,store_in_dense_array",
+                         [(False, False), (False, True),
+                          (True, False), (True, True)])
+def test_fc_cutoff_radius_svecs(ph_nacl, ph_nacl_compact_fcsym,
+                                is_compact, store_in_dense_array):
+    if is_compact:
+        ph = ph_nacl_compact_fcsym
+    else:
+        ph = ph_nacl
+
+    fc = ph.force_constants.copy()
+    primitive_matrix = np.dot(np.linalg.inv(ph.supercell_matrix),
+                              ph.primitive_matrix)
+    primitive = get_primitive(ph.supercell, primitive_matrix,
+                              store_in_dense_array=store_in_dense_array)
+
+    cutoff_force_constants(fc, ph.supercell, primitive, 4.0)
+    changed = (np.abs(ph.force_constants - fc) > 1e-8)
+
+    if is_compact:
+        assert np.sum(changed) == 534
+    else:
+        assert np.sum(changed) == 17088

test/structure/test_cells.py

@@ -111,15 +111,15 @@ def test_get_primitive_convcell_nacl(convcell_nacl,
     helper_methods.compare_cells_with_order(pcell, primcell_nacl)
 
 
-@pytest.mark.parametrize("store_dense_vectors", [True, False])
+@pytest.mark.parametrize("store_in_dense_array", [True, False])
 def test_get_primitive_convcell_nacl_svecs(convcell_nacl,
                                            primcell_nacl,
-                                           store_dense_vectors):
+                                           store_in_dense_array):
     pcell = get_primitive(convcell_nacl,
                           primitive_matrix_nacl,
-                          store_dense_vectors=store_dense_vectors)
+                          store_in_dense_array=store_in_dense_array)
     svecs, multi = pcell.get_smallest_vectors()
-    if store_dense_vectors:
+    if store_in_dense_array:
         assert svecs.shape == (54, 3)
         assert multi.shape == (8, 2, 2)
         assert np.sum(multi[:, :, 0]) == 54
@@ -174,7 +174,7 @@ def test_ShortestPairs_dense_nacl(ph_nacl, helper_methods):
     pcell = ph_nacl.primitive
     pos = scell.scaled_positions
     spairs = ShortestPairs(scell.cell, pos, pos[pcell.p2s_map],
-                           store_dense_vectors=True)
+                           store_in_dense_array=True)
     svecs = spairs.shortest_vectors
     multi = spairs.multiplicities
     assert multi[-1, -1, :].sum() == multi[:, :, 0].sum()