Equivalent treatment of imaginary modes to integrate JDOS in python implementation

This commit is contained in:
Atsushi Togo 2022-03-24 19:58:17 +09:00
parent 8bd306774d
commit 1481240578
2 changed files with 20 additions and 19 deletions

View File

@ -314,13 +314,15 @@ class JointDos:
else:
self._run_with_g()
def _run_with_g(self, lang="Py"):
def _run_with_g(self, lang="C"):
"""Calculate JDOS.
lang="Py" is the original implementation.
lang="C" calculates JDOS using C routine for imag-free-energy.
Computational efficiency is determined by tetraherdon method, but not
integration in JDOS. Performance benefit with lang="C" is very limited.
Computational efficiency is roughly determined by tetraherdon method, but not
integration in JDOS. Although performance benefit using lang="C" is limited,
using the same routine as imag-free-energy is considered a good idea.
So here, the implementation in C is used for the integration of JDOS.
"""
jdos = np.zeros((len(self._frequency_points), 2), dtype="double", order="C")
@ -377,21 +379,20 @@ class JointDos:
t = self._temperature
freqs = self._frequencies[self._triplets_at_q[:, 1:]]
self._occupations = np.where(
freqs > self._cutoff_frequency, bose_einstein(freqs, t), 0
freqs > self._cutoff_frequency, bose_einstein(freqs, t), -1
)
def _run_py_with_g_at_temperature(self, jdos, i):
g = self._g
n = self._occupations
for k, l in list(np.ndindex(g.shape[3:])):
weights = np.where(
np.logical_or(n[:, 0, k] < 0, n[:, 1, l] < 0), 0, self._weights_at_q
)
jdos[i, 1] += np.dot(
(n[:, 0, k] + n[:, 1, l] + 1) * g[0, :, i, k, l],
self._weights_at_q,
)
jdos[i, 0] += np.dot(
(n[:, 0, k] - n[:, 1, l]) * g[1, :, i, k, l],
self._weights_at_q,
(n[:, 0, k] + n[:, 1, l] + 1) * g[0, :, i, k, l], weights
)
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)

View File

@ -149,16 +149,16 @@ nacl_freq_points_at_300K = [
nacl_jdos_12_at_300K = [
0.0000000,
0.0000000,
8.4625631,
8.4768159,
0.0000000,
4.1076174,
1.5151176,
0.7992725,
6.7993659,
4.1241485,
1.4712023,
0.8016066,
6.7628440,
0.0000000,
21.2271309,
21.2134161,
0.0000000,
26.9803907,
26.9803216,
0.0000000,
14.9103483,
0.0000000,
@ -274,9 +274,9 @@ def test_jdos_nacl_at_300K(nacl_pbe: Phono3py):
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:], atol=1e-5
nacl_jdos_12_at_300K[2:], jdos.joint_dos.ravel()[2:], rtol=1e-2, atol=1e-5
)