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

Write test for JDOS integration in python

This commit is contained in:
Atsushi Togo 2022-03-24 23:00:13 +09:00
parent 1481240578
commit e8ae6ef276
2 changed files with 142 additions and 102 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

108
phono3py/phonon3/joint_dos.py
View File

@ -39,12 +39,10 @@ import warnings
import numpy as np
from phonopy.harmonic.dynamical_matrix import DynamicalMatrix, get_dynamical_matrix
from phonopy.structure.cells import Primitive
from phonopy.structure.tetrahedron_method import TetrahedronMethod
from phonopy.units import VaspToTHz
from phono3py.phonon3.triplets import (
get_nosym_triplets_at_q,
get_tetrahedra_vertices,
get_triplets_at_q,
get_triplets_integration_weights,
)
@ -65,6 +63,7 @@ class JointDos:
nac_params=None,
nac_q_direction=None,
sigma=None,
sigma_cutoff=None,
cutoff_frequency=None,
frequency_factor_to_THz=VaspToTHz,
frequency_scale_factor=1.0,
@ -85,6 +84,7 @@ class JointDos:
self.nac_q_direction = nac_q_direction
self._sigma = None
self.sigma = sigma
self._sigma_cutoff = sigma_cutoff
if cutoff_frequency is None:
self._cutoff_frequency = 0
@ -117,17 +117,6 @@ class JointDos:
self._ones_pp_strength = None
self._temperature = None
def run(self):
"""Calculate joint-density-of-states."""
self.run_phonon_solver(np.arange(len(self._bz_grid.addresses), dtype="int_"))
try:
import phono3py._phono3py as phono3c # noqa F401
self._run_c()
except ImportError:
print("Joint density of states in python is not implemented.")
return None, None
@property
def dynamical_matrix(self) -> DynamicalMatrix:
"""Return dynamical matrix class instance."""
@ -235,12 +224,15 @@ class JointDos:
else:
self._temperature = float(temperature)
def get_triplets_at_q(self):
"""Return triplets information."""
return self._triplets_at_q, self._weights_at_q
def set_grid_point(self, grid_point):
"""Set a grid point at which joint-DOS is calculated."""
self._grid_point = grid_point
self._set_triplets()
self._joint_dos = None
self._frequency_points = None
if self._phonon_done is None:
self._allocate_phonons()
@ -265,17 +257,17 @@ class JointDos:
self.run_phonon_solver(np.array([gamma_gp, grid_point], dtype="int_"))
def get_triplets_at_q(self):
"""Return triplets information."""
return self._triplets_at_q, self._weights_at_q
def run_phonon_solver(self, grid_points):
def run_phonon_solver(self, grid_points=None):
"""Calculate phonons at grid_points.
This method is used in get_triplets_integration_weights by this
method name. So this name is not allowed to change.
"""
if grid_points is None:
_grid_points = np.arange(len(self._bz_grid.addresses), dtype="int_")
else:
_grid_points = grid_points
if self._phonon_done is None:
self._allocate_phonons()
run_phonon_solver_c(
@ -283,7 +275,7 @@ class JointDos:
self._frequencies,
self._eigenvectors,
self._phonon_done,
grid_points,
_grid_points,
self._bz_grid.addresses,
self._bz_grid.QDinv,
self._frequency_factor_to_THz,
@ -291,31 +283,31 @@ class JointDos:
self._lapack_zheev_uplo,
)
def run_integration_weights(self, freq_points):
def run(self):
"""Calculate joint-density-of-states."""
self.run_phonon_solver()
try:
import phono3py._phono3py as phono3c # noqa F401
self.run_integration_weights()
self.run_jdos()
except ImportError:
print("Joint density of states in python is not implemented.")
return None, None
def run_integration_weights(self, lang="C"):
"""Compute triplets integration weights."""
self._g, self._g_zero = get_triplets_integration_weights(
self,
np.array(freq_points, dtype="double"),
self._frequency_points,
self._sigma,
sigma_cutoff=self._sigma_cutoff,
is_collision_matrix=(self._temperature is None),
lang=lang,
)
def _run_c(self, lang="C"):
if self._sigma is None:
if lang == "C":
self._run_with_g()
else:
if self._temperature is not None:
print(
"JDOS with phonon occupation numbers doesn't work "
"in this option."
)
self._run_py_tetrahedron_method()
else:
self._run_with_g()
def _run_with_g(self, lang="C"):
"""Calculate JDOS.
def run_jdos(self, lang="C"):
"""Run JDOS calculation with having integration weights.
lang="Py" is the original implementation.
lang="C" calculates JDOS using C routine for imag-free-energy.
@ -326,7 +318,6 @@ class JointDos:
"""
jdos = np.zeros((len(self._frequency_points), 2), dtype="double", order="C")
self.run_integration_weights(self._frequency_points)
if self._temperature is None:
for i, _ in enumerate(self._frequency_points):
g = self._g
@ -394,45 +385,6 @@ class JointDos:
)
jdos[i, 0] += np.dot((n[:, 0, k] - n[:, 1, l]) * g[1, :, i, k, l], weights)
def _run_py_tetrahedron_method(self):
thm = TetrahedronMethod(self._bz_grid.microzone_lattice)
self._vertices = get_tetrahedra_vertices(
np.array(
np.dot(thm.get_tetrahedra(), self._bz_grid.P.T), dtype="int_", order="C"
),
self._bz_grid.D_diag,
self._triplets_at_q,
self._bz_grid,
)
self.run_phonon_solver(self._vertices.ravel())
f_max = np.max(self._frequencies) * 2
f_max *= 1.005
f_min = 0
self._set_uniform_frequency_points(f_min, f_max)
num_freq_points = len(self._frequency_points)
jdos = np.zeros((num_freq_points, 2), dtype="double")
for vertices, w in zip(self._vertices, self._weights_at_q):
for i, j in list(np.ndindex(self._num_band, self._num_band)):
f1 = self._frequencies[vertices[0], i]
f2 = self._frequencies[vertices[1], j]
thm.set_tetrahedra_omegas(f1 + f2)
thm.run(self._frequency_points)
iw = thm.get_integration_weight()
jdos[:, 1] += iw * w
thm.set_tetrahedra_omegas(f1 - f2)
thm.run(self._frequency_points)
iw = thm.get_integration_weight()
jdos[:, 0] += iw * w
thm.set_tetrahedra_omegas(-f1 + f2)
thm.run(self._frequency_points)
iw = thm.get_integration_weight()
jdos[:, 0] += iw * w
self._joint_dos = jdos / np.prod(self._bz_grid.D_diag)
def _init_dynamical_matrix(self):
self._dm = get_dynamical_matrix(
self._fc2,

136
test/phonon3/test_joint_dos.py
View File

@ -5,6 +5,7 @@ import pytest
from phono3py import Phono3py
from phono3py.api_jointdos import Phono3pyJointDos
from phono3py.phonon3.joint_dos import JointDos
from phono3py.phonon.grid import BZGrid
si_freq_points = [
0.0000000,
@ -170,8 +171,14 @@ nacl_jdos_12_at_300K = [
]
def test_jdos_si(si_pbesol: Phono3py):
"""Test joint-DOS by Si."""
@pytest.mark.parametrize("gp,store_dense_gp_map", [(103, False), (105, True)])
def test_jdos_si(si_pbesol: Phono3py, gp: int, store_dense_gp_map: bool):
"""Test joint-DOS by Si.
store_dense_gp_map=False : 103
store_dense_gp_map=True : 105
"""
si_pbesol.mesh_numbers = [9, 9, 9]
jdos = Phono3pyJointDos(
si_pbesol.phonon_supercell,
@ -179,9 +186,11 @@ def test_jdos_si(si_pbesol: Phono3py):
si_pbesol.fc2,
mesh=si_pbesol.mesh_numbers,
num_frequency_points=10,
store_dense_gp_map=store_dense_gp_map,
log_level=1,
)
jdos.run([103])
jdos.run([gp])
# print(", ".join(["%.7f" % fp for fp in jdos.frequency_points]))
np.testing.assert_allclose(si_freq_points, jdos.frequency_points, atol=1e-5)
# print(", ".join(["%.7f" % jd for jd in jdos.joint_dos.ravel()]))
@ -190,7 +199,8 @@ def test_jdos_si(si_pbesol: Phono3py):
)
def test_jdso_si_nomeshsym(si_pbesol: Phono3py):
@pytest.mark.parametrize("gp,store_dense_gp_map", [(103, False), (105, True)])
def test_jdso_si_nomeshsym(si_pbesol: Phono3py, gp: int, store_dense_gp_map: bool):
"""Test joint-DOS without considering mesh symmetry by Si."""
si_pbesol.mesh_numbers = [9, 9, 9]
jdos = Phono3pyJointDos(
@ -200,9 +210,10 @@ def test_jdso_si_nomeshsym(si_pbesol: Phono3py):
mesh=si_pbesol.mesh_numbers,
num_frequency_points=10,
is_mesh_symmetry=False,
store_dense_gp_map=store_dense_gp_map,
log_level=1,
)
jdos.run([103])
jdos.run([gp])
# print(", ".join(["%.7f" % fp for fp in jdos.frequency_points]))
np.testing.assert_allclose(si_freq_points, jdos.frequency_points, atol=1e-5)
# print(", ".join(["%.7f" % jd for jd in jdos.joint_dos.ravel()]))
@ -211,7 +222,8 @@ def test_jdso_si_nomeshsym(si_pbesol: Phono3py):
)
def test_jdos_nacl(nacl_pbe: Phono3py):
@pytest.mark.parametrize("gp,store_dense_gp_map", [(103, False), (105, True)])
def test_jdos_nacl(nacl_pbe: Phono3py, gp: int, store_dense_gp_map: bool):
"""Test joint-DOS by NaCl."""
nacl_pbe.mesh_numbers = [9, 9, 9]
jdos = Phono3pyJointDos(
@ -221,9 +233,10 @@ def test_jdos_nacl(nacl_pbe: Phono3py):
mesh=nacl_pbe.mesh_numbers,
nac_params=nacl_pbe.nac_params,
num_frequency_points=10,
store_dense_gp_map=store_dense_gp_map,
log_level=1,
)
jdos.run([103])
jdos.run([gp])
# print(", ".join(["%.7f" % fp for fp in jdos.frequency_points]))
np.testing.assert_allclose(nacl_freq_points, jdos.frequency_points, atol=1e-5)
# print(", ".join(["%.7f" % jd for jd in jdos.joint_dos.ravel()]))
@ -232,7 +245,8 @@ def test_jdos_nacl(nacl_pbe: Phono3py):
)
def test_jdos_nacl_gamma(nacl_pbe: Phono3py):
@pytest.mark.parametrize("gp,store_dense_gp_map", [(0, False), (0, True)])
def test_jdos_nacl_gamma(nacl_pbe: Phono3py, gp: int, store_dense_gp_map: bool):
"""Test joint-DOS at Gamma-point by NaCl."""
nacl_pbe.mesh_numbers = [9, 9, 9]
jdos = Phono3pyJointDos(
@ -243,9 +257,10 @@ def test_jdos_nacl_gamma(nacl_pbe: Phono3py):
nac_params=nacl_pbe.nac_params,
nac_q_direction=[1, 0, 0],
num_frequency_points=10,
store_dense_gp_map=store_dense_gp_map,
log_level=1,
)
jdos.run([0])
jdos.run([gp])
# print(", ".join(["%.7f" % fp for fp in jdos.frequency_points]))
np.testing.assert_allclose(nacl_freq_points_gamma, jdos.frequency_points, atol=1e-5)
# print(", ".join(["%.7f" % jd for jd in jdos.joint_dos.ravel()]))
@ -254,7 +269,8 @@ def test_jdos_nacl_gamma(nacl_pbe: Phono3py):
)
def test_jdos_nacl_at_300K(nacl_pbe: Phono3py):
@pytest.mark.parametrize("gp,store_dense_gp_map", [(103, False), (105, True)])
def test_jdos_nacl_at_300K(nacl_pbe: Phono3py, gp: int, store_dense_gp_map: bool):
"""Test joint-DOS at 300K by NaCl."""
nacl_pbe.mesh_numbers = [9, 9, 9]
jdos = Phono3pyJointDos(
@ -267,44 +283,64 @@ def test_jdos_nacl_at_300K(nacl_pbe: Phono3py):
temperatures=[
300,
],
store_dense_gp_map=store_dense_gp_map,
log_level=1,
)
jdos.run([103])
jdos.run([gp])
# print(", ".join(["%.7f" % fp for fp in jdos.frequency_points]))
np.testing.assert_allclose(
nacl_freq_points_at_300K, jdos.frequency_points, atol=1e-5
)
print(", ".join(["%.7f" % jd for jd in jdos.joint_dos.ravel()]))
# print(", ".join(["%.7f" % jd for jd in jdos.joint_dos.ravel()]))
np.testing.assert_allclose(
nacl_jdos_12_at_300K[2:], jdos.joint_dos.ravel()[2:], rtol=1e-2, atol=1e-5
)
def test_jdos_nac_direction_phonon_NaCl(nacl_pbe: Phono3py):
@pytest.mark.parametrize("gp,store_dense_gp_map", [(0, False), (0, True)])
def test_jdos_nac_direction_phonon_NaCl(
nacl_pbe: Phono3py, gp: int, store_dense_gp_map: bool
):
"""Test JDOS of NaCl with nac_q_direction."""
jdos = _get_jdos(nacl_pbe, [7, 7, 7], nac_params=nacl_pbe.nac_params)
jdos = _get_jdos(
nacl_pbe,
[7, 7, 7],
nac_params=nacl_pbe.nac_params,
store_dense_gp_map=store_dense_gp_map,
)
jdos.nac_q_direction = [1, 0, 0]
jdos.set_grid_point(0)
jdos.set_grid_point(gp)
frequencies, _, _ = jdos.get_phonons()
np.testing.assert_allclose(
frequencies[0], [0, 0, 0, 4.59488262, 4.59488262, 7.41183870], rtol=0, atol=1e-6
)
def test_jdos_nac_direction_phonon_NaCl_second_error(nacl_pbe: Phono3py):
@pytest.mark.parametrize("gp,store_dense_gp_map", [(0, False), (0, True)])
def test_jdos_nac_direction_phonon_NaCl_second_error(
nacl_pbe: Phono3py, gp: int, store_dense_gp_map: bool
):
"""Test JDOS of NaCl with nac_q_direction.
Second setting non-gamma grid point must raise exception.
"""
jdos = _get_jdos(nacl_pbe, [7, 7, 7], nac_params=nacl_pbe.nac_params)
jdos = _get_jdos(
nacl_pbe,
[7, 7, 7],
nac_params=nacl_pbe.nac_params,
store_dense_gp_map=store_dense_gp_map,
)
jdos.nac_q_direction = [1, 0, 0]
jdos.set_grid_point(0)
jdos.set_grid_point(gp)
with pytest.raises(RuntimeError):
jdos.set_grid_point(1)
def test_jdos_nac_direction_phonon_NaCl_second_no_error(nacl_pbe: Phono3py):
@pytest.mark.parametrize("gp1,gp2,store_dense_gp_map", [(0, 1, False), (0, 1, True)])
def test_jdos_nac_direction_phonon_NaCl_second_no_error(
nacl_pbe: Phono3py, gp1: int, gp2: int, store_dense_gp_map: bool
):
"""Test JDOS of NaCl with nac_q_direction.
Second setting non-gamma grid point should not raise exception because
@ -312,25 +348,77 @@ def test_jdos_nac_direction_phonon_NaCl_second_no_error(nacl_pbe: Phono3py):
NAC.
"""
jdos = _get_jdos(nacl_pbe, [7, 7, 7], nac_params=nacl_pbe.nac_params)
jdos = _get_jdos(
nacl_pbe,
[7, 7, 7],
nac_params=nacl_pbe.nac_params,
store_dense_gp_map=store_dense_gp_map,
)
jdos.nac_q_direction = [1, 0, 0]
jdos.set_grid_point(0)
jdos.set_grid_point(gp1)
jdos.nac_q_direction = None
jdos.set_grid_point(1)
jdos.set_grid_point(gp2)
frequencies, _, _ = jdos.get_phonons()
np.testing.assert_allclose(
frequencies[0], [0, 0, 0, 4.59488262, 4.59488262, 4.59488262], rtol=0, atol=1e-6
)
def _get_jdos(ph3: Phono3py, mesh, nac_params=None):
@pytest.mark.parametrize("gp,store_dense_gp_map", [(103, False), (105, True)])
def test_jdos_nac_NaCl_300K_C(nacl_pbe: Phono3py, gp: int, store_dense_gp_map: bool):
"""Test running JDOS of NaCl in C mode."""
jdos = _get_jdos(
nacl_pbe,
[9, 9, 9],
nac_params=nacl_pbe.nac_params,
store_dense_gp_map=store_dense_gp_map,
)
jdos.set_grid_point(gp)
jdos.frequency_points = nacl_freq_points_at_300K
jdos.temperature = 300
jdos.run_phonon_solver()
jdos.run_integration_weights()
jdos.run_jdos()
np.testing.assert_allclose(
nacl_jdos_12_at_300K[2:], jdos.joint_dos.ravel()[2:], rtol=1e-2, atol=1e-5
)
@pytest.mark.parametrize("gp,store_dense_gp_map", [(103, False), (105, True)])
def test_jdos_nac_NaCl_300K_Py(nacl_pbe: Phono3py, gp: int, store_dense_gp_map: bool):
"""Test running JDOS of NaCl in Py mode."""
jdos = _get_jdos(
nacl_pbe,
[9, 9, 9],
nac_params=nacl_pbe.nac_params,
store_dense_gp_map=store_dense_gp_map,
)
jdos.set_grid_point(gp)
jdos.frequency_points = nacl_freq_points_at_300K
jdos.temperature = 300
jdos.run_phonon_solver()
jdos.run_integration_weights()
jdos.run_jdos(lang="Py")
np.testing.assert_allclose(
nacl_jdos_12_at_300K[2:], jdos.joint_dos.ravel()[2:], rtol=1e-2, atol=1e-5
)
def _get_jdos(ph3: Phono3py, mesh, nac_params=None, store_dense_gp_map=False):
bz_grid = BZGrid(
mesh,
lattice=ph3.primitive.cell,
symmetry_dataset=ph3.primitive_symmetry.dataset,
store_dense_gp_map=store_dense_gp_map,
)
ph3.mesh_numbers = mesh
jdos = JointDos(
ph3.primitive,
ph3.supercell,
ph3.grid,
bz_grid,
ph3.fc2,
nac_params=nac_params,
store_dense_gp_map=store_dense_gp_map,
cutoff_frequency=1e-4,
)
return jdos