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

Introduce batch for frequency points of JDOS

This commit is contained in:
Atsushi Togo 2022-03-24 17:14:40 +09:00
parent 4f89cf25a5
commit 8bd306774d
5 changed files with 199 additions and 142 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

94
phono3py/api_jointdos.py
View File

@ -38,6 +38,10 @@ from phonopy.structure.symmetry import Symmetry
from phonopy.units import VaspToTHz
from phono3py.file_IO import write_joint_dos
from phono3py.phonon3.imag_self_energy import (
get_freq_points_batches,
get_frequency_points,
)
from phono3py.phonon3.joint_dos import JointDos
from phono3py.phonon.grid import BZGrid
@ -94,9 +98,6 @@ class Phono3pyJointDos:
fc2,
nac_params=nac_params,
cutoff_frequency=cutoff_frequency,
frequency_step=frequency_step,
num_frequency_points=num_frequency_points,
temperatures=self._temperatures,
frequency_factor_to_THz=frequency_factor_to_THz,
frequency_scale_factor=frequency_scale_factor,
is_mesh_symmetry=self._is_mesh_symmetry,
@ -107,12 +108,30 @@ class Phono3pyJointDos:
)
self._joint_dos = None
self._num_frequency_points_in_batch = None
self._frequency_step = frequency_step
self._num_frequency_points = num_frequency_points
@property
def grid(self):
"""Return BZGrid class instance."""
return self._bz_grid
@property
def num_frequency_points_in_batch(self):
"""Getter and setter of num_frequency_points_in_batch.
Number of sampling frequency points per batch.
Larger value gives better concurrency in tetrahedron method,
but requires more memory.
"""
return self._num_frequency_points_in_batch
@num_frequency_points_in_batch.setter
def num_frequency_points_in_batch(self, nelems_in_batch):
self._num_frequency_points_in_batch = nelems_in_batch
def run(self, grid_points, write_jdos=False):
"""Calculate joint-density-of-states."""
if self._log_level:
@ -122,6 +141,36 @@ class Phono3pyJointDos:
)
print("Sampling mesh: [ %d %d %d ]" % tuple(self._bz_grid.D_diag))
self._jdos.run_phonon_solver(
np.arange(len(self._bz_grid.addresses), dtype="int_")
)
frequencies, _, _ = self._jdos.get_phonons()
max_phonon_freq = np.max(frequencies)
self._frequency_points = get_frequency_points(
max_phonon_freq=max_phonon_freq,
sigmas=self._sigmas,
frequency_points=None,
frequency_step=self._frequency_step,
num_frequency_points=self._num_frequency_points,
)
batches = get_freq_points_batches(
len(self._frequency_points), nelems=self._num_frequency_points_in_batch
)
if self._temperatures is None:
temperatures = [None]
else:
temperatures = self._temperatures
self._joint_dos = np.zeros(
(
len(self._sigmas),
len(temperatures),
len(self._frequency_points),
2,
),
dtype="double",
order="C",
)
for i, gp in enumerate(grid_points):
if (self._bz_grid.addresses[gp] == 0).all():
self._jdos.nac_q_direction = self._nac_q_direction
@ -147,19 +196,32 @@ class Phono3pyJointDos:
if not self._sigmas:
raise RuntimeError("sigma or tetrahedron method has to be set.")
for sigma in self._sigmas:
for i_s, sigma in enumerate(self._sigmas):
self._jdos.sigma = sigma
if self._log_level:
if sigma is None:
print("Tetrahedron method is used.")
else:
print("Smearing method with sigma=%s is used." % sigma)
self._jdos.set_sigma(sigma)
self._jdos.run()
print(
f"Calculations at {len(self._frequency_points)} "
f"frequency points are devided into {len(batches)} batches."
)
for i_t, temperature in enumerate(temperatures):
self._jdos.temperature = temperature
if write_jdos:
filename = self._write(gp, sigma=sigma)
if self._log_level:
print('JDOS is written into "%s".' % filename)
for ib, freq_indices in enumerate(batches):
print(f"{ib + 1}/{len(batches)}: {freq_indices}")
self._jdos.frequency_points = self._frequency_points[
freq_indices
]
self._jdos.run()
self._joint_dos[i_s, i_t, freq_indices] = self._jdos.joint_dos
if write_jdos:
filename = self._write(gp, i_sigma=i_s)
if self._log_level:
print('JDOS is written into "%s".' % filename)
@property
def dynamical_matrix(self):
@ -169,20 +231,20 @@ class Phono3pyJointDos:
@property
def frequency_points(self):
"""Return frequency points."""
return self._jdos.frequency_points
return self._frequency_points
@property
def joint_dos(self):
"""Return calculated joint-density-of-states."""
return self._jdos.joint_dos
return self._joint_dos
def _write(self, gp, sigma=None):
def _write(self, gp, i_sigma=0):
return write_joint_dos(
gp,
self._bz_grid.D_diag,
self._jdos.frequency_points,
self._jdos.joint_dos,
sigma=sigma,
self._frequency_points,
self._joint_dos[i_sigma],
sigma=self._sigmas[i_sigma],
temperatures=self._temperatures,
filename=self._filename,
is_mesh_symmetry=self._is_mesh_symmetry,

2
phono3py/file_IO.py
View File

@ -535,7 +535,7 @@ def write_joint_dos(
gp,
mesh,
frequencies,
jdos,
jdos[0],
sigma=sigma,
temperature=None,
filename=filename,

14
phono3py/phonon3/imag_self_energy.py
View File

@ -553,12 +553,16 @@ def run_ise_at_frequency_points_batch(
] = ise.get_detailed_imag_self_energy()
def get_freq_points_batches(tot_nelems, nelems=10):
def get_freq_points_batches(tot_nelems, nelems=None):
"""Divide frequency points into batches."""
nbatch = tot_nelems // nelems
batches = [np.arange(i * nelems, (i + 1) * nelems) for i in range(nbatch)]
if tot_nelems % nelems > 0:
batches.append(np.arange(nelems * nbatch, tot_nelems))
if nelems is None:
_nelems = 10
else:
_nelems = nelems
nbatch = tot_nelems // _nelems
batches = [np.arange(i * _nelems, (i + 1) * _nelems) for i in range(nbatch)]
if tot_nelems % _nelems > 0:
batches.append(np.arange(_nelems * nbatch, tot_nelems))
return batches

181
phono3py/phonon3/joint_dos.py
View File

@ -42,10 +42,6 @@ from phonopy.structure.cells import Primitive
from phonopy.structure.tetrahedron_method import TetrahedronMethod
from phonopy.units import VaspToTHz
from phono3py.phonon3.imag_self_energy import (
get_freq_points_batches,
get_frequency_points,
)
from phono3py.phonon3.triplets import (
get_nosym_triplets_at_q,
get_tetrahedra_vertices,
@ -70,9 +66,6 @@ class JointDos:
nac_q_direction=None,
sigma=None,
cutoff_frequency=None,
frequency_step=None,
num_frequency_points=None,
temperatures=None,
frequency_factor_to_THz=VaspToTHz,
frequency_scale_factor=1.0,
is_mesh_symmetry=True,
@ -91,15 +84,12 @@ class JointDos:
self._nac_params = nac_params
self.nac_q_direction = nac_q_direction
self._sigma = None
self.set_sigma(sigma)
self.sigma = sigma
if cutoff_frequency is None:
self._cutoff_frequency = 0
else:
self._cutoff_frequency = cutoff_frequency
self._frequency_step = frequency_step
self._num_frequency_points = num_frequency_points
self._temperatures = temperatures
self._frequency_factor_to_THz = frequency_factor_to_THz
self._frequency_scale_factor = frequency_scale_factor
self._is_mesh_symmetry = is_mesh_symmetry
@ -125,6 +115,7 @@ class JointDos:
self._g = None
self._g_zero = None
self._ones_pp_strength = None
self._temperature = None
def run(self):
"""Calculate joint-density-of-states."""
@ -207,13 +198,43 @@ class JointDos:
warnings.warn("Use attribute, nac_q_direction", DeprecationWarning)
self.nac_q_direction = nac_q_direction
def set_sigma(self, sigma):
"""Set sigma value. None means tetrahedron method."""
@property
def sigma(self):
"""Getter and setter of sigma."""
return self._sigma
@sigma.setter
def sigma(self, sigma):
if sigma is None:
self._sigma = None
else:
self._sigma = float(sigma)
def set_sigma(self, sigma):
"""Set sigma value. None means tetrahedron method."""
warnings.warn(
"Use attribute, JointDOS.sigma instead of JointDOS.set_sigma()",
DeprecationWarning,
)
self.sigma = sigma
@property
def bz_grid(self) -> BZGrid:
"""Return BZGrid."""
return self._bz_grid
@property
def temperature(self):
"""Setter and getter of temperature."""
return self._temperature
@temperature.setter
def temperature(self, temperature):
if temperature is None:
self._temperature = None
else:
self._temperature = float(temperature)
def set_grid_point(self, grid_point):
"""Set a grid point at which joint-DOS is calculated."""
self._grid_point = grid_point
@ -248,11 +269,6 @@ class JointDos:
"""Return triplets information."""
return self._triplets_at_q, self._weights_at_q
@property
def bz_grid(self) -> BZGrid:
"""Return BZGrid."""
return self._bz_grid
def run_phonon_solver(self, grid_points):
"""Calculate phonons at grid_points.
@ -260,6 +276,8 @@ class JointDos:
method name. So this name is not allowed to change.
"""
if self._phonon_done is None:
self._allocate_phonons()
run_phonon_solver_c(
self._dm,
self._frequencies,
@ -279,7 +297,7 @@ class JointDos:
self,
np.array(freq_points, dtype="double"),
self._sigma,
is_collision_matrix=(self._temperatures is None),
is_collision_matrix=(self._temperature is None),
)
def _run_c(self, lang="C"):
@ -287,7 +305,7 @@ class JointDos:
if lang == "C":
self._run_with_g()
else:
if self._temperatures is not None:
if self._temperature is not None:
print(
"JDOS with phonon occupation numbers doesn't work "
"in this option."
@ -305,36 +323,19 @@ class JointDos:
integration in JDOS. Performance benefit with lang="C" is very limited.
"""
max_phonon_freq = np.max(self._frequencies)
self._frequency_points = get_frequency_points(
max_phonon_freq=max_phonon_freq,
sigmas=[
self._sigma,
],
frequency_points=None,
frequency_step=self._frequency_step,
num_frequency_points=self._num_frequency_points,
)
if self._temperatures is None:
jdos = np.zeros((len(self._frequency_points), 2), dtype="double", order="C")
for i, freq_point in enumerate(self._frequency_points):
self.run_integration_weights([freq_point])
if self._temperatures is None:
g = self._g
jdos[i, 1] = np.sum(
np.tensordot(g[0, :, 0], self._weights_at_q, axes=(0, 0))
)
gx = g[2] - g[0]
jdos[i, 0] = np.sum(
np.tensordot(gx[:, 0], self._weights_at_q, axes=(0, 0))
)
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
jdos[i, 1] = np.sum(
np.tensordot(g[0, :, i], self._weights_at_q, axes=(0, 0))
)
gx = g[2] - g[0]
jdos[i, 0] = np.sum(
np.tensordot(gx[:, i], self._weights_at_q, axes=(0, 0))
)
else:
jdos = np.zeros(
(len(self._temperatures), len(self._frequency_points), 2),
dtype="double",
order="C",
)
if lang == "C":
num_band = len(self._primitive) * 3
self._ones_pp_strength = np.ones(
@ -342,67 +343,55 @@ class JointDos:
dtype="double",
order="C",
)
batches = get_freq_points_batches(len(self._frequency_points))
print(
f"Calculations at {len(self._frequency_points)} "
f"frequency points are devided into {len(batches)} batches."
)
for ib, freq_indices in enumerate(batches):
print(f"{ib + 1}/{len(batches)}: {freq_indices}")
fpoints = self._frequency_points[freq_indices]
self.run_integration_weights(fpoints)
for k in range(2):
g = self._g.copy()
g[k] = 0
self._run_c_with_g_at_temperature(jdos, g, k, freq_indices)
for k in range(2):
g = self._g.copy()
g[k] = 0
self._run_c_with_g_at_temperature(jdos, g, k)
else:
self._run_occupation()
for i, freq_point in enumerate(self._frequency_points):
self.run_integration_weights([freq_point])
for i, _ in enumerate(self._frequency_points):
self._run_py_with_g_at_temperature(jdos, i)
self._joint_dos = jdos / np.prod(self._bz_grid.D_diag)
def _run_c_with_g_at_temperature(self, jdos, g, k, freq_indices):
def _run_c_with_g_at_temperature(self, jdos, g, k):
import phono3py._phono3py as phono3c
jdos_elem = np.zeros(1, dtype="double")
for i, i_freq in enumerate(freq_indices):
for j, temp in enumerate(self._temperatures):
phono3c.imag_self_energy_with_g(
jdos_elem,
self._ones_pp_strength,
self._triplets_at_q,
self._weights_at_q,
self._frequencies,
float(temp),
g,
self._g_zero,
self._cutoff_frequency,
i,
)
jdos[j, i_freq, k] = jdos_elem[0]
for i, _ in enumerate(self._frequency_points):
phono3c.imag_self_energy_with_g(
jdos_elem,
self._ones_pp_strength,
self._triplets_at_q,
self._weights_at_q,
self._frequencies,
self._temperature,
g,
self._g_zero,
self._cutoff_frequency,
i,
)
jdos[i, k] = jdos_elem[0]
def _run_occupation(self):
self._occupations = []
for t in self._temperatures:
freqs = self._frequencies[self._triplets_at_q[:, 1:]]
self._occupations.append(
np.where(freqs > self._cutoff_frequency, bose_einstein(freqs, t), 0)
)
t = self._temperature
freqs = self._frequencies[self._triplets_at_q[:, 1:]]
self._occupations = np.where(
freqs > self._cutoff_frequency, bose_einstein(freqs, t), 0
)
def _run_py_with_g_at_temperature(self, jdos, i):
g = self._g
for j, n in enumerate(self._occupations):
for k, l in list(np.ndindex(g.shape[3:])):
jdos[j, i, 1] += np.dot(
(n[:, 0, k] + n[:, 1, l] + 1) * g[0, :, 0, k, l],
self._weights_at_q,
)
jdos[j, i, 0] += np.dot(
(n[:, 0, k] - n[:, 1, l]) * g[1, :, 0, k, l],
self._weights_at_q,
)
n = self._occupations
for k, l in list(np.ndindex(g.shape[3:])):
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,
)
def _run_py_tetrahedron_method(self):
thm = TetrahedronMethod(self._bz_grid.microzone_lattice)

50
test/phonon3/test_joint_dos.py
View File

@ -97,39 +97,41 @@ nacl_jdos_12 = [
0.0000000,
0.0000000,
]
nacl_freq_points_gamma = [
0.0000000,
1.6553106,
3.3106213,
4.9659319,
6.6212426,
8.2765532,
9.9318639,
11.5871745,
13.2424851,
14.8977958,
1.6322306,
3.2644613,
4.8966919,
6.5289225,
8.1611531,
9.7933838,
11.4256144,
13.0578450,
14.6900756,
]
nacl_jdos_12_gamma = [
1742452844146884.7500000,
1638874677039989.2500000,
0.0000000,
8.8165476,
0.0415488,
1.4914142,
0.3104766,
0.3679421,
1.0509976,
0.0358263,
5.8578016,
9.1146639,
0.0403982,
1.5633599,
0.2967087,
0.3866759,
0.9911241,
0.0446820,
5.2759325,
0.0000000,
7.2272898,
9.3698670,
0.0000000,
5.7740314,
5.1064113,
0.0000000,
0.6663207,
0.0000000,
0.1348658,
0.8748614,
0.0000000,
0.1474919,
0.0000000,
0.0112830,
]
nacl_freq_points_at_300K = [
@ -274,7 +276,7 @@ def test_jdos_nacl_at_300K(nacl_pbe: Phono3py):
)
# 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
nacl_jdos_12_at_300K[2:], jdos.joint_dos.ravel()[2:], atol=1e-5
)