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

Sampling frequency points for real-free-energy 

Sampling frequency points are determined by the max phonon frequency in
phonon data on grid. This commit changes the behaviour slightly.
The phonon frequencies at Gamma point is given without NAC though
previously with NAC. This change aims the uniform frequency sampling
over all grid points.
This commit is contained in:
Atsushi Togo 2022-05-09 18:56:43 +09:00
parent 4ca9daf208
commit 5b6f442ef2
2 changed files with 354 additions and 249 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

505
phono3py/phonon3/real_self_energy.py
View File

@ -48,246 +48,6 @@ from phono3py.phonon3.interaction import Interaction
from phono3py.phonon.func import bose_einstein
def get_real_self_energy(
interaction: Interaction,
grid_points,
temperatures,
epsilons=None,
frequency_points=None,
frequency_step=None,
num_frequency_points=None,
frequency_points_at_bands=False,
write_hdf5=True,
output_filename=None,
log_level=0,
):
"""Real part of self energy at frequency points.
Band indices to be calculated at are kept in Interaction instance.
Parameters
----------
interaction : Interaction
Ph-ph interaction.
grid_points : array_like
Grid-point indices where imag-self-energeis are caclculated.
dtype=int, shape=(grid_points,)
temperatures : array_like
Temperatures where imag-self-energies are calculated.
dtype=float, shape=(temperatures,)
epsilons : array_like
Smearing widths to computer principal part. When multiple values
are given frequency shifts for those values are returned.
dtype=float, shape=(epsilons,)
frequency_points : array_like, optional
Frequency sampling points. Default is None. With
frequency_points_at_bands=False and frequency_points is None,
num_frequency_points or frequency_step is used to generate uniform
frequency sampling points.
dtype=float, shape=(frequency_points,)
frequency_step : float, optional
Uniform pitch of frequency sampling points. Default is None. This
results in using num_frequency_points.
num_frequency_points: int, optional
Number of sampling sampling points to be used instead of
frequency_step. This number includes end points. Default is None,
which gives 201.
frequency_points_at_bands : bool, optional
Phonon band frequencies are used as frequency points when True.
Default is False.
num_points_in_batch: int, optional
Number of sampling points in one batch. This is for the frequency
sampling mode and the sampling points are divided into batches.
Lager number provides efficient use of multi-cores but more
memory demanding. Default is None, which give the number of 10.
log_level: int
Log level. Default is 0.
Returns
-------
tuple :
(frequency_points, all_deltas) are returned.
When frequency_points_at_bands is True,
all_deltas.shape = (epsilons, temperatures, grid_points,
band_indices)
otherwise
all_deltas.shape = (epsilons, temperatures, grid_points,
band_indices, frequency_points)
"""
if epsilons is None:
_epsilons = [
None,
]
else:
_epsilons = epsilons
_temperatures = np.array(temperatures, dtype="double")
if (interaction.get_phonons()[2] == 0).any():
if log_level:
print("Running harmonic phonon calculations...")
interaction.run_phonon_solver()
fst = RealSelfEnergy(interaction)
mesh = interaction.mesh_numbers
bz_grid = interaction.bz_grid
frequencies = interaction.get_phonons()[0]
max_phonon_freq = np.amax(frequencies)
band_indices = interaction.band_indices
if frequency_points_at_bands:
_frequency_points = None
all_deltas = np.zeros(
(len(_epsilons), len(_temperatures), len(grid_points), len(band_indices)),
dtype="double",
order="C",
)
else:
_frequency_points = get_frequency_points(
max_phonon_freq=max_phonon_freq,
sigmas=epsilons,
frequency_points=frequency_points,
frequency_step=frequency_step,
num_frequency_points=num_frequency_points,
)
all_deltas = np.zeros(
(
len(_epsilons),
len(_temperatures),
len(grid_points),
len(band_indices),
len(_frequency_points),
),
dtype="double",
order="C",
)
fst.frequency_points = _frequency_points
for j, gp in enumerate(grid_points):
fst.grid_point = gp
if log_level:
weights = interaction.get_triplets_at_q()[1]
if len(grid_points) > 1:
print(
"------------------- Real part of self energy -o- (%d/%d) "
"-------------------" % (j + 1, len(grid_points))
)
else:
print(
"----------------------- Real part of self energy -o- "
"-----------------------"
)
print("Grid point: %d" % gp)
print("Number of ir-triplets: %d / %d" % (len(weights), weights.sum()))
fst.run_interaction()
frequencies = interaction.get_phonons()[0][gp]
if log_level:
bz_grid = interaction.bz_grid
qpoint = np.dot(bz_grid.QDinv, bz_grid.addresses[gp])
print("Phonon frequencies at (%4.2f, %4.2f, %4.2f):" % tuple(qpoint))
for bi, freq in enumerate(frequencies):
print("%3d %f" % (bi + 1, freq))
sys.stdout.flush()
for i, epsilon in enumerate(_epsilons):
fst.epsilon = epsilon
for k, t in enumerate(_temperatures):
fst.temperature = t
fst.run()
all_deltas[i, k, j] = fst.real_self_energy.T
# if not frequency_points_at_bands:
# pos = 0
# for bi_set in [[bi, ] for bi in band_indices]:
# filename = write_real_self_energy(
# gp,
# bi_set,
# _frequency_points,
# all_deltas[i, k, j, pos:(pos + len(bi_set))],
# mesh,
# fst.epsilon,
# t,
# filename=output_filename)
# pos += len(bi_set)
# if log_level:
# print("Real part of self energies were stored in "
# "\"%s\"." % filename)
# sys.stdout.flush()
if write_hdf5:
filename = write_real_self_energy_to_hdf5(
gp,
band_indices,
_temperatures,
all_deltas[i, :, j],
mesh,
fst.epsilon,
bz_grid=bz_grid,
frequency_points=_frequency_points,
frequencies=frequencies,
filename=output_filename,
)
if log_level:
print('Real part of self energies were stored in "%s".' % filename)
sys.stdout.flush()
return _frequency_points, all_deltas
def write_real_self_energy(
real_self_energy,
mesh,
grid_points,
band_indices,
frequency_points,
temperatures,
epsilons,
output_filename=None,
is_mesh_symmetry=True,
log_level=0,
):
"""Write real-part of self-energies into files."""
if epsilons is None:
_epsilons = [
RealSelfEnergy.default_epsilon,
]
else:
_epsilons = epsilons
for epsilon, rse_temps in zip(_epsilons, real_self_energy):
for t, rse_gps in zip(temperatures, rse_temps):
for gp, rse in zip(grid_points, rse_gps):
for i, bi in enumerate(band_indices):
pos = 0
for j in range(i):
pos += len(band_indices[j])
filename = write_real_self_energy_at_grid_point(
gp,
bi,
frequency_points,
rse[pos : (pos + len(bi))].sum(axis=0) / len(bi),
mesh,
epsilon,
t,
filename=output_filename,
is_mesh_symmetry=is_mesh_symmetry,
)
if log_level:
print(
"Real parts of self-energies were "
'written to "%s".' % filename
)
class RealSelfEnergy:
"""Class to calculate real-part of self-energy of bubble diagram.
@ -311,7 +71,12 @@ class RealSelfEnergy:
default_epsilon = 0.05
def __init__(
self, interaction, grid_point=None, temperature=None, epsilon=None, lang="C"
self,
interaction: Interaction,
grid_point=None,
temperature=None,
epsilon=None,
lang="C",
):
"""Init method.
@ -604,13 +369,6 @@ class RealSelfEnergy:
return sum_d
def imag_to_real(im_part, frequency_points):
"""Calculate real-part of self-energy from the imaginary-part."""
i2r = ImagToReal(im_part, frequency_points)
i2r.run()
return i2r.re_part, i2r.frequency_points
class ImagToReal:
"""Calculate real part of self-energy using Kramers-Kronig relation."""
@ -710,3 +468,254 @@ class ImagToReal:
_frequency_points = np.hstack([-frequency_points[::-1], frequency_points[1:]])
return _im_part, _frequency_points, df
def get_real_self_energy(
interaction: Interaction,
grid_points,
temperatures,
epsilons=None,
frequency_points=None,
frequency_step=None,
num_frequency_points=None,
frequency_points_at_bands=False,
write_hdf5=True,
output_filename=None,
log_level=0,
):
"""Real part of self energy at frequency points.
Band indices to be calculated at are kept in Interaction instance.
Parameters
----------
interaction : Interaction
Ph-ph interaction.
grid_points : array_like
Grid-point indices where imag-self-energeis are caclculated.
dtype=int, shape=(grid_points,)
temperatures : array_like
Temperatures where imag-self-energies are calculated.
dtype=float, shape=(temperatures,)
epsilons : array_like
Smearing widths to computer principal part. When multiple values
are given frequency shifts for those values are returned.
dtype=float, shape=(epsilons,)
frequency_points : array_like, optional
Frequency sampling points. Default is None. With
frequency_points_at_bands=False and frequency_points is None,
num_frequency_points or frequency_step is used to generate uniform
frequency sampling points.
dtype=float, shape=(frequency_points,)
frequency_step : float, optional
Uniform pitch of frequency sampling points. Default is None. This
results in using num_frequency_points.
num_frequency_points: int, optional
Number of sampling sampling points to be used instead of
frequency_step. This number includes end points. Default is None,
which gives 201.
frequency_points_at_bands : bool, optional
Phonon band frequencies are used as frequency points when True.
Default is False.
num_points_in_batch: int, optional
Number of sampling points in one batch. This is for the frequency
sampling mode and the sampling points are divided into batches.
Lager number provides efficient use of multi-cores but more
memory demanding. Default is None, which give the number of 10.
log_level: int
Log level. Default is 0.
Returns
-------
tuple :
(frequency_points, all_deltas) are returned.
When frequency_points_at_bands is True,
all_deltas.shape = (epsilons, temperatures, grid_points,
band_indices)
otherwise
all_deltas.shape = (epsilons, temperatures, grid_points,
band_indices, frequency_points)
"""
if epsilons is None:
_epsilons = [
None,
]
else:
_epsilons = epsilons
_temperatures = np.array(temperatures, dtype="double")
if (interaction.get_phonons()[2] == 0).any():
if log_level:
print("Running harmonic phonon calculations...")
interaction.run_phonon_solver()
fst = RealSelfEnergy(interaction)
mesh = interaction.mesh_numbers
bz_grid = interaction.bz_grid
# Set phonon at Gamma without NAC for finding max_phonon_freq.
interaction.run_phonon_solver_at_gamma()
max_phonon_freq = np.amax(interaction.get_phonons()[0])
interaction.run_phonon_solver_at_gamma(is_nac=True)
band_indices = interaction.band_indices
if frequency_points_at_bands:
_frequency_points = None
all_deltas = np.zeros(
(len(_epsilons), len(_temperatures), len(grid_points), len(band_indices)),
dtype="double",
order="C",
)
else:
_frequency_points = get_frequency_points(
max_phonon_freq=max_phonon_freq,
sigmas=epsilons,
frequency_points=frequency_points,
frequency_step=frequency_step,
num_frequency_points=num_frequency_points,
)
all_deltas = np.zeros(
(
len(_epsilons),
len(_temperatures),
len(grid_points),
len(band_indices),
len(_frequency_points),
),
dtype="double",
order="C",
)
fst.frequency_points = _frequency_points
for j, gp in enumerate(grid_points):
fst.grid_point = gp
if log_level:
weights = interaction.get_triplets_at_q()[1]
if len(grid_points) > 1:
print(
"------------------- Real part of self energy -o- (%d/%d) "
"-------------------" % (j + 1, len(grid_points))
)
else:
print(
"----------------------- Real part of self energy -o- "
"-----------------------"
)
print("Grid point: %d" % gp)
print("Number of ir-triplets: %d / %d" % (len(weights), weights.sum()))
fst.run_interaction()
frequencies = interaction.get_phonons()[0][gp]
if log_level:
bz_grid = interaction.bz_grid
qpoint = np.dot(bz_grid.QDinv, bz_grid.addresses[gp])
print("Phonon frequencies at (%4.2f, %4.2f, %4.2f):" % tuple(qpoint))
for bi, freq in enumerate(frequencies):
print("%3d %f" % (bi + 1, freq))
sys.stdout.flush()
for i, epsilon in enumerate(_epsilons):
fst.epsilon = epsilon
for k, t in enumerate(_temperatures):
fst.temperature = t
fst.run()
all_deltas[i, k, j] = fst.real_self_energy.T
# if not frequency_points_at_bands:
# pos = 0
# for bi_set in [[bi, ] for bi in band_indices]:
# filename = write_real_self_energy(
# gp,
# bi_set,
# _frequency_points,
# all_deltas[i, k, j, pos:(pos + len(bi_set))],
# mesh,
# fst.epsilon,
# t,
# filename=output_filename)
# pos += len(bi_set)
# if log_level:
# print("Real part of self energies were stored in "
# "\"%s\"." % filename)
# sys.stdout.flush()
if write_hdf5:
filename = write_real_self_energy_to_hdf5(
gp,
band_indices,
_temperatures,
all_deltas[i, :, j],
mesh,
fst.epsilon,
bz_grid=bz_grid,
frequency_points=_frequency_points,
frequencies=frequencies,
filename=output_filename,
)
if log_level:
print('Real part of self energies were stored in "%s".' % filename)
sys.stdout.flush()
return _frequency_points, all_deltas
def write_real_self_energy(
real_self_energy,
mesh,
grid_points,
band_indices,
frequency_points,
temperatures,
epsilons,
output_filename=None,
is_mesh_symmetry=True,
log_level=0,
):
"""Write real-part of self-energies into files."""
if epsilons is None:
_epsilons = [
RealSelfEnergy.default_epsilon,
]
else:
_epsilons = epsilons
for epsilon, rse_temps in zip(_epsilons, real_self_energy):
for t, rse_gps in zip(temperatures, rse_temps):
for gp, rse in zip(grid_points, rse_gps):
for i, bi in enumerate(band_indices):
pos = 0
for j in range(i):
pos += len(band_indices[j])
filename = write_real_self_energy_at_grid_point(
gp,
bi,
frequency_points,
rse[pos : (pos + len(bi))].sum(axis=0) / len(bi),
mesh,
epsilon,
t,
filename=output_filename,
is_mesh_symmetry=is_mesh_symmetry,
)
if log_level:
print(
"Real parts of self-energies were "
'written to "%s".' % filename
)
def imag_to_real(im_part, frequency_points):
"""Calculate real-part of self-energy from the imaginary-part."""
i2r = ImagToReal(im_part, frequency_points)
i2r.run()
return i2r.re_part, i2r.frequency_points

98
test/phonon3/test_real_self_energy.py
View File

@ -1,6 +1,7 @@
"""Test for real_self_energy.py."""
import numpy as np
from phono3py import Phono3py
from phono3py.phonon3.real_self_energy import ImagToReal
si_pbesol_Delta = [
@ -167,6 +168,80 @@ im_part = [
[30.3852997, 0.0468188, 30.3852997, 0.4941608, 30.5387608, 0.4706733],
[30.6922219, 0.0000000, 30.6922219, 0.4173657, 30.8456830, 0.3821416],
]
delta_nacl_nac = [
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
-0.43835832,
-0.43835832,
-0.27175512,
0.00000000,
0.00000000,
0.00000000,
-0.24732785,
-0.24732785,
-0.15332798,
0.00000000,
0.00000000,
0.00000000,
-0.45125778,
-0.45125778,
-0.27975198,
0.00000000,
0.00000000,
0.00000000,
-0.80781739,
-0.80781739,
-0.50079699,
0.00000000,
0.00000000,
0.00000000,
-0.36217245,
-0.36217245,
-0.22452456,
0.00000000,
0.00000000,
0.00000000,
0.54848940,
0.54848940,
0.34002959,
0.00000000,
0.00000000,
0.00000000,
0.29335952,
0.29335952,
0.18186480,
0.00000000,
0.00000000,
0.00000000,
0.16329650,
0.16329650,
0.10123375,
0.00000000,
0.00000000,
0.00000000,
0.11208511,
0.11208511,
0.06948585,
]
freq_points_nacl_nac = [
0.0,
1.63223063,
3.26446125,
4.89669188,
6.5289225,
8.16115313,
9.79338375,
11.42561438,
13.057845,
14.69007563,
]
def test_real_self_energy_with_band_indices(si_pbesol):
@ -188,7 +263,7 @@ def test_real_self_energy_with_band_indices(si_pbesol):
np.testing.assert_allclose(si_pbesol_Delta, delta[0, 0, :], atol=0.01)
def test_real_self_energy_with_frequency_points(si_pbesol):
def test_real_self_energy_with_frequency_points(si_pbesol: Phono3py):
"""Real part of self energy spectrum of Si.
* specified frquency points
@ -221,6 +296,27 @@ def test_real_self_energy_with_frequency_points(si_pbesol):
# print("".join("%.8f, " % s for s in b))
def test_real_self_energy_nacl_nac_npoints(nacl_pbe: Phono3py):
"""Real part of self energy spectrum of NaCl.
* at 10 frequency points sampled uniformly.
* at q->0
"""
nacl_pbe.mesh_numbers = [9, 9, 9]
nacl_pbe.init_phph_interaction(nac_q_direction=[1, 0, 0])
fps, delta = nacl_pbe.run_real_self_energy(
[nacl_pbe.grid.gp_Gamma], [300], num_frequency_points=10
)
# for line in np.swapaxes(delta, -1, -2).ravel().reshape(-1, 6):
# print(("%10.8f, " * 6) % tuple(line))
# print(fps.ravel())
np.testing.assert_allclose(freq_points_nacl_nac, fps.ravel(), rtol=0, atol=1e-5)
np.testing.assert_allclose(
delta_nacl_nac, np.swapaxes(delta, -1, -2).ravel(), rtol=0, atol=1e-2
)
def test_ImagToReal():
"""Test ImagToReal class (Kramers–Kronig relation)."""
vals = np.array(im_part)