scnc
/
phono3py
mirror of https://github.com/phonopy/phono3py.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Add test for Phono3py.symmetrize_fc3 and symmetrize_fc2

This commit is contained in:
Atsushi Togo 2025-07-08 18:12:43 +09:00
parent 8488ed0424
commit cfcd64fb9f
3 changed files with 129 additions and 36 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

50
phono3py/api_phono3py.py
View File

@ -1545,9 +1545,12 @@ class Phono3py:
if fc_calculator == "traditional" or fc_calculator is None:
if symmetrize_fc3r:
self.symmetrize_fc3(
use_symfc_projector=use_symfc_projector and fc_calculator is None
)
if use_symfc_projector and fc_calculator is None:
self.symmetrize_fc3(
use_symfc_projector=True, symfc_options=fc_calculator_options
)
else:
self.symmetrize_fc3()
elif fc_calculator == "symfc":
symfc_solver = cast(SymfcFCSolver, fc_solver.fc_solver)
fc3_nonzero_elems = symfc_solver.get_nonzero_atomic_indices_fc3()
@ -1570,13 +1573,17 @@ class Phono3py:
self._fc2 = fc2
if fc_calculator == "traditional" or fc_calculator is None:
if symmetrize_fc3r:
self.symmetrize_fc2(
use_symfc_projector=use_symfc_projector
and fc_calculator is None
)
if use_symfc_projector and fc_calculator is None:
self.symmetrize_fc2(
use_symfc_projector=True,
symfc_options=fc_calculator_options,
)
else:
self.symmetrize_fc2()
def symmetrize_fc3(
self,
level: int = 1,
use_symfc_projector: bool = False,
symfc_options: str | None = None,
):
@ -1584,6 +1591,9 @@ class Phono3py:
Parameters
----------
level : int, optional
Number of times translational and permutation symmetries are applied
for traditional symmetrization. Default is 1.
use_symfc_projector : bool, optional
If True, the force constants are symmetrized by symfc projector
instead of traditional approach.
@ -1614,12 +1624,13 @@ class Phono3py:
else:
if self._log_level:
print("Symmetrizing fc3 by traditional approach.", flush=True)
if self._fc3.shape[0] == self._fc3.shape[1]:
set_translational_invariance_fc3(self._fc3)
set_permutation_symmetry_fc3(self._fc3)
else:
set_translational_invariance_compact_fc3(self._fc3, self._primitive)
set_permutation_symmetry_compact_fc3(self._fc3, self._primitive)
for _ in range(level):
if self._fc3.shape[0] == self._fc3.shape[1]:
set_translational_invariance_fc3(self._fc3)
set_permutation_symmetry_fc3(self._fc3)
else:
set_translational_invariance_compact_fc3(self._fc3, self._primitive)
set_permutation_symmetry_compact_fc3(self._fc3, self._primitive)
def produce_fc2(
self,
@ -1686,6 +1697,7 @@ class Phono3py:
def symmetrize_fc2(
self,
level: int = 1,
use_symfc_projector: bool = False,
symfc_options: str | None = None,
):
@ -1693,6 +1705,9 @@ class Phono3py:
Parameters
----------
level : int, optional
Number of times translational and permutation symmetries are applied
for traditional symmetrization. Default is 1.
use_symfc_projector : bool, optional
If True, the force constants are symmetrized by symfc projector
instead of traditional approach.
@ -1733,10 +1748,11 @@ class Phono3py:
else:
if self._log_level:
print("Symmetrizing fc2 by traditional approach.", flush=True)
if self._fc2.shape[0] == self._fc2.shape[1]:
symmetrize_force_constants(self._fc2)
else:
symmetrize_compact_force_constants(self._fc2, primitive)
for _ in range(level):
if self._fc2.shape[0] == self._fc2.shape[1]:
symmetrize_force_constants(self._fc2)
else:
symmetrize_compact_force_constants(self._fc2, primitive)
def cutoff_fc3_by_zero(self, cutoff_distance, fc3=None):
"""Set zero to fc3 elements out of cutoff distance.

57
phono3py/phonon3/fc3.py
View File

@ -34,10 +34,13 @@
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
from __future__ import annotations
import logging
import sys
import numpy as np
from numpy.typing import NDArray
from phonopy.harmonic.force_constants import (
distribute_force_constants,
get_nsym_list_and_s2pp,
@ -573,30 +576,50 @@ def cutoff_fc3_by_zero(fc3, supercell, cutoff_distance, p2s_map=None, symprec=1e
break
def show_drift_fc3(fc3, primitive=None, name="fc3"):
"""Show drift of fc3."""
def show_drift_fc3(
fc3: NDArray, primitive: Primitive | None = None, name: str = "fc3", digit: int = 8
):
"""Show max drift of fc3."""
maxval1, maxval2, maxval3, xyz1, xyz2, xyz3 = get_drift_fc3(
fc3, primitive=primitive
)
text = f"Max drift of {name}: "
text += f"{maxval1:.{digit}f} ({'xyz'[xyz1[0]]}{'xyz'[xyz1[1]]}{'xyz'[xyz1[2]]}) "
text += f"{maxval2:.{digit}f} ({'xyz'[xyz2[0]]}{'xyz'[xyz2[1]]}{'xyz'[xyz2[2]]}) "
text += f"{maxval3:.{digit}f} ({'xyz'[xyz3[0]]}{'xyz'[xyz3[1]]}{'xyz'[xyz3[2]]})"
print(text)
def get_drift_fc3(
fc3: NDArray, primitive: Primitive | None = None
) -> tuple[float, float, float, list[int], list[int], list[int]]:
"""Return max drift of fc3."""
if fc3.shape[0] == fc3.shape[1]:
num_atom = fc3.shape[0]
maxval1 = 0
maxval2 = 0
maxval3 = 0
klm1 = [0, 0, 0]
klm2 = [0, 0, 0]
klm3 = [0, 0, 0]
xyz1 = [0, 0, 0]
xyz2 = [0, 0, 0]
xyz3 = [0, 0, 0]
for i, j, k, ll, m in list(np.ndindex((num_atom, num_atom, 3, 3, 3))):
val1 = fc3[:, i, j, k, ll, m].sum()
val2 = fc3[i, :, j, k, ll, m].sum()
val3 = fc3[i, j, :, k, ll, m].sum()
if abs(val1) > abs(maxval1):
maxval1 = val1
klm1 = [k, ll, m]
xyz1 = [k, ll, m]
if abs(val2) > abs(maxval2):
maxval2 = val2
klm2 = [k, ll, m]
xyz2 = [k, ll, m]
if abs(val3) > abs(maxval3):
maxval3 = val3
klm3 = [k, ll, m]
xyz3 = [k, ll, m]
else:
if primitive is None:
raise RuntimeError(
"Primitive cell is required to get drift of compact fc3."
)
try:
import phono3py._phono3py as phono3c
@ -614,9 +637,9 @@ def show_drift_fc3(fc3, primitive=None, name="fc3"):
maxval1 = 0
maxval2 = 0
maxval3 = 0
klm1 = [0, 0, 0]
klm2 = [0, 0, 0]
klm3 = [0, 0, 0]
xyz1 = [0, 0, 0]
xyz2 = [0, 0, 0]
xyz3 = [0, 0, 0]
phono3c.transpose_compact_fc3(
fc3, permutations, s2pp_map, p2s_map, nsym_list, 0
) # dim[0] <--> dim[1]
@ -624,7 +647,7 @@ def show_drift_fc3(fc3, primitive=None, name="fc3"):
val1 = fc3[i, :, j, k, ll, m].sum()
if abs(val1) > abs(maxval1):
maxval1 = val1
klm1 = [k, ll, m]
xyz1 = [k, ll, m]
phono3c.transpose_compact_fc3(
fc3, permutations, s2pp_map, p2s_map, nsym_list, 0
) # dim[0] <--> dim[1]
@ -633,10 +656,10 @@ def show_drift_fc3(fc3, primitive=None, name="fc3"):
val3 = fc3[i, j, :, k, ll, m].sum()
if abs(val2) > abs(maxval2):
maxval2 = val2
klm2 = [k, ll, m]
xyz2 = [k, ll, m]
if abs(val3) > abs(maxval3):
maxval3 = val3
klm3 = [k, ll, m]
xyz3 = [k, ll, m]
except ImportError as exc:
text = (
"Import error at phono3c.tranpose_compact_fc3. "
@ -644,11 +667,7 @@ def show_drift_fc3(fc3, primitive=None, name="fc3"):
)
raise RuntimeError(text) from exc
text = "Max drift of %s: " % name
text += "%f (%s%s%s) " % (maxval1, "xyz"[klm1[0]], "xyz"[klm1[1]], "xyz"[klm1[2]])
text += "%f (%s%s%s) " % (maxval2, "xyz"[klm2[0]], "xyz"[klm2[1]], "xyz"[klm2[2]])
text += "%f (%s%s%s)" % (maxval3, "xyz"[klm3[0]], "xyz"[klm3[1]], "xyz"[klm3[2]])
print(text)
return maxval1, maxval2, maxval3, xyz1, xyz2, xyz3
def _set_permutation_symmetry_fc3_elem_with_cutoff(fc3, fc3_done, a, b, c):

58
test/api/test_api_phono3py.py
View File

@ -6,9 +6,11 @@ from pathlib import Path
import numpy as np
import pytest
from phonopy.harmonic.force_constants import get_drift_force_constants
from phonopy.interface.pypolymlp import PypolymlpParams
from phono3py import Phono3py
from phono3py.phonon3.fc3 import get_drift_fc3
cwd = Path(__file__).parent
@ -207,3 +209,59 @@ def test_use_pypolymlp_mgo(mgo_222rd_444rd: Phono3py):
assert (
pytest.approx(63.0018546, abs=1e-1) == ph3.thermal_conductivity.kappa[0, 0, 0]
)
@pytest.mark.parametrize("is_compact_fc", [True, False])
def test_symmetrize_fc_traditional(si_pbesol: Phono3py, is_compact_fc: bool):
"""Test symmetrize_fc3 and symmetrize_fc2 with traditional approach."""
ph3 = Phono3py(
si_pbesol.unitcell,
supercell_matrix=si_pbesol.supercell_matrix,
primitive_matrix=si_pbesol.primitive_matrix,
log_level=2,
)
ph3.dataset = si_pbesol.dataset
ph3.produce_fc3(is_compact_fc=is_compact_fc)
assert ph3.fc3 is not None
assert ph3.fc2 is not None
v1, v2, v3, _, _, _ = get_drift_fc3(ph3.fc3, primitive=ph3.primitive)
np.testing.assert_allclose(
[v1, v2, v3], [-1.755065e-01, -1.749287e-01, 3.333333e-05], atol=1e-6
)
ph3.symmetrize_fc3(level=3)
v1_sym, v2_sym, v3_sym, _, _, _ = get_drift_fc3(ph3.fc3, primitive=ph3.primitive)
if is_compact_fc:
np.testing.assert_allclose([v1_sym, v2_sym, v3_sym], 1.217081e-05, atol=1e-6)
else:
np.testing.assert_allclose([v1_sym, v2_sym, v3_sym], 1.421085e-14, atol=1e-6)
v1_sym, v2_sym, _, _ = get_drift_force_constants(ph3.fc2, primitive=ph3.primitive)
if is_compact_fc:
np.testing.assert_allclose([v1_sym, v2_sym], -1.0e-06, atol=1e-6)
else:
np.testing.assert_allclose([v1_sym, v2_sym], -1.0e-06, atol=1e-6)
@pytest.mark.parametrize("is_compact_fc", [True, False])
def test_symmetrize_fc_symfc(si_pbesol: Phono3py, is_compact_fc: bool):
"""Test symmetrize_fc3 and symmetrize_fc2 with symfc."""
pytest.importorskip("symfc")
ph3 = Phono3py(
si_pbesol.unitcell,
supercell_matrix=si_pbesol.supercell_matrix,
primitive_matrix=si_pbesol.primitive_matrix,
log_level=2,
)
ph3.dataset = si_pbesol.dataset
ph3.produce_fc3(is_compact_fc=is_compact_fc)
assert ph3.fc3 is not None
assert ph3.fc2 is not None
ph3.symmetrize_fc3(use_symfc_projector=True)
v1_sym, v2_sym, v3_sym, _, _, _ = get_drift_fc3(ph3.fc3, primitive=ph3.primitive)
np.testing.assert_allclose([v1_sym, v2_sym, v3_sym], 0, atol=1e-6)
ph3.symmetrize_fc2(use_symfc_projector=True)
v1_sym, v2_sym, _, _ = get_drift_force_constants(ph3.fc2, primitive=ph3.primitive)
np.testing.assert_allclose([v1_sym, v2_sym], 0, atol=1e-6)