Spectral function implementation statred.

This commit is contained in:
Atsushi Togo 2020-09-05 22:24:45 +09:00
parent bcc1b8ced6
commit 75f3595cb5
6 changed files with 138 additions and 36 deletions

View File

@ -1421,18 +1421,19 @@ class Phono3py(object):
raise RuntimeError(msg)
if temperatures is None:
temperatures = [0.0, 300.0]
self._temperatures = [0.0, 300.0]
else:
self._temperatures = temperatures
self._grid_points = grid_points
self._temperatures = temperatures
self._scattering_event_class = scattering_event_class
vals = get_imag_self_energy(
self._interaction,
grid_points,
self._sigmas,
temperatures,
sigmas=self._sigmas,
frequency_points=frequency_points,
frequency_step=frequency_step,
num_frequency_points=num_frequency_points,
temperatures=temperatures,
scattering_event_class=scattering_event_class,
write_gamma_detail=write_gamma_detail,
return_gamma_detail=keep_gamma_detail,

View File

@ -459,6 +459,12 @@ class FrequencyShift(object):
return sum_d
def imag_to_real(im_part, frequency_points):
i2r = ImagToReal(im_part, frequency_points)
i2r.run()
return i2r.re_part, i2r.frequency_points
class ImagToReal(object):
"""Calculate real part of self-energy using Kramers-Kronig relation"""
@ -512,7 +518,7 @@ class ImagToReal(object):
def _pick_one(self):
re_part = []
fpoints = []
coef = - self._df / np.pi
coef = self._df / np.pi
i_zero = (len(self._im_part) - 1) // 2
for i, im_part_at_i in enumerate(self._im_part):
if i < i_zero:
@ -529,7 +535,7 @@ class ImagToReal(object):
def _half_shift(self):
re_part = []
fpoints = []
coef = - self._df / np.pi
coef = self._df / np.pi
i_zero = (len(self._im_part) - 1) // 2
for i, im_part_at_i in enumerate(self._im_part):
if i < i_zero:
@ -542,15 +548,6 @@ class ImagToReal(object):
return (np.array(re_part, dtype='double'),
np.array(fpoints, dtype='double'))
# def half_shift(self):
# df = gammas[1, 0] - gammas[0, 0]
# shift_kk = []
# freqs = gammas[:, 0] + df / 2
# for f in freqs:
# vals = gammas[:, 1] / (gammas[:, 0] - f)
# shift_kk.append(- vals.sum() / np.pi * df)
# return np.array([freqs, shift_kk]).T
def _expand_bubble_im_part(self, im_part, frequency_points):
if (np.abs(frequency_points[0]) > 1e-8).any():
raise RuntimeError(

View File

@ -44,11 +44,11 @@ from phono3py.file_IO import (write_gamma_detail_to_hdf5,
def get_imag_self_energy(interaction,
grid_points,
sigmas,
temperatures,
sigmas=None,
frequency_points=None,
frequency_step=None,
num_frequency_points=None,
temperatures=None,
scattering_event_class=None, # class 1 or 2
write_gamma_detail=False,
return_gamma_detail=False,
@ -65,10 +65,14 @@ def get_imag_self_energy(interaction,
grid_points : array_like
Grid-point indices where imag-self-energeis are caclculated.
dtype=int, shape=(grid_points,)
sigmas : array_like
temperatures : array_like
Temperatures where imag-self-energies are calculated.
dtype=float, shape=(temperatures,)
sigmas : array_like, optional
A set of sigmas. simgas=[None, ] means to use tetrahedron method,
otherwise smearing method with real positive value of sigma.
For example, sigmas=[None, 0.01, 0.03] is possible.
For example, sigmas=[None, 0.01, 0.03] is possible. Default is None,
which results in [None, ].
dtype=float, shape=(sigmas,)
frequency_points : array_like, optional
Frequency sampling points. Default is None. In this case,
@ -82,10 +86,6 @@ def get_imag_self_energy(interaction,
Number of sampling sampling points to be used instead of
frequency_step. This number includes end points. Default is None,
which gives 201.
temperatures : array_like, optional
Temperatures where imag-self-energies are calculated. Default is None,
which gives [0, 300].
dtype=float, shape=(temperatures,)
scattering_event_class : int, optional
Specific choice of scattering event class, 1 or 2 that is specified
1 or 2, respectively. The result is stored in gammas. Therefore
@ -105,13 +105,11 @@ def get_imag_self_energy(interaction,
(gammas, detailed_gammas, frequency_points) are returned.
"""
if temperatures is None:
temperatures = [0.0, 300.0]
if temperatures is None:
if log_level:
print("Temperatures have to be set.")
return False
if sigmas is None:
_sigmas = [None, ]
else:
_sigmas = sigmas
if (interaction.get_phonons()[2] == 0).any():
if log_level:
@ -124,7 +122,7 @@ def get_imag_self_energy(interaction,
_frequency_points = get_frequency_points(
max_phonon_freq=max_phonon_freq,
sigmas=sigmas,
sigmas=_sigmas,
frequency_points=frequency_points,
frequency_step=frequency_step,
num_frequency_points=num_frequency_points)
@ -132,7 +130,7 @@ def get_imag_self_energy(interaction,
ise = ImagSelfEnergy(
interaction, with_detail=(write_gamma_detail or return_gamma_detail))
gamma = np.zeros(
(len(grid_points), len(sigmas), len(temperatures),
(len(grid_points), len(_sigmas), len(temperatures),
len(_frequency_points), len(interaction.band_indices)),
dtype='double', order='C')
detailed_gamma = []
@ -168,11 +166,11 @@ def get_imag_self_energy(interaction,
num_band0 = len(interaction.band_indices)
num_band = frequencies.shape[1]
detailed_gamma_at_gp = np.zeros(
(len(sigmas), len(temperatures), len(_frequency_points),
(len(_sigmas), len(temperatures), len(_frequency_points),
len(weights), num_band0, num_band, num_band),
dtype='double')
for j, sigma in enumerate(sigmas):
for j, sigma in enumerate(_sigmas):
if log_level:
if sigma:
print("Sigma: %s" % sigma)

View File

@ -0,0 +1,87 @@
# Copyright (C) 2020 Atsushi Togo
# All rights reserved.
#
# This file is part of phono3py.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions
# are met:
#
# * Redistributions of source code must retain the above copyright
# notice, this list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright
# notice, this list of conditions and the following disclaimer in
# the documentation and/or other materials provided with the
# distribution.
#
# * Neither the name of the phonopy project nor the names of its
# contributors may be used to endorse or promote products derived
# from this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
import numpy as np
from phono3py.phonon3.imag_self_energy import get_imag_self_energy
from phono3py.phonon3.frequency_shift import imag_to_real
class SpectralFunction(object):
"""Calculate spectral function"""
def __init__(self,
interaction,
grid_points,
frequency_points=None,
frequency_step=None,
num_frequency_points=None,
temperatures=None):
self._interaction = interaction
self._grid_points = grid_points
self._frequency_points_in = frequency_points
self._frequency_step = frequency_step
self._num_frequency_points = num_frequency_points
self._temperatures = temperatures
def run(self):
# gammas[grid_points, sigmas, temps, freq_points, band_indices]
gammas, self._frequency_points = get_imag_self_energy(
self._interaction,
self._grid_points,
frequency_points=self._frequency_points_in,
frequency_step=self._frequency_step,
num_frequency_points=self._num_frequency_points,
temperatures=self._temperatures)
self._gammas = np.array(gammas[:, 0], dtype='double', order='C')
self._deltas = np.zeros_like(self._gammas)
for i, gp in enumerate(self._grid_points):
for j, temp in enumerate(self._temperatures):
for k, bi in enumerate(self._interaction.band_indices):
re_part, fpoints = imag_to_real(
self._gammas[i, j, :, k], self._frequency_points)
self._deltas[i, j, :, k] = -re_part
assert (np.abs(self._frequency_points - fpoints) < 1e-8).all()
@property
def shifts(self):
return self._deltas
@property
def half_linewidths(self):
return self._gammas
@property
def frequency_points(self):
return self._frequency_points

View File

@ -7,7 +7,7 @@ si_pbesol_Delta = [
[-0.02078728, -0.02102094, -0.06573269, -0.11432603,
-0.1366966, -0.14371315]]
# imag-self-energy Si-PBEsol 50x50x50 gp=5, bi=4, 101 points
# imag-self-energy Si-PBEsol 50x50x50 gp=5, bi=4, 101 points, 300K
im_part = [
[0.00000000, 0.00000000, -0.00000000, -0.17502231, 0.15346111, -0.17945834],
[0.30692220, 0.01806856, 0.30692220, -0.16920492, 0.46038331, -0.16365925],
@ -127,10 +127,10 @@ def test_RealToImag():
vals = np.array(im_part)
i2r = ImagToReal(vals[:, 1], vals[:, 0])
i2r.run()
pick_one_vals = i2r.re_part
pick_one_vals = -1 * i2r.re_part # -1 to make it freq-shift
pick_one_freqs = i2r.frequency_points
i2r.run(method='half_shift')
half_shift_vals = i2r.re_part
half_shift_vals = -1 * i2r.re_part # -1 to make it freq-shift
half_shift_freqs = i2r.frequency_points
# for f, im, f1, re1, f2, re2, in zip(
# vals[:, 0], vals[:, 1],

View File

@ -0,0 +1,19 @@
import numpy as np
from phono3py.phonon3.spectral_function import SpectralFunction
def test_SpectralFunction(si_pbesol):
si_pbesol.mesh_numbers = [9, 9, 9]
si_pbesol.init_phph_interaction()
sf = SpectralFunction(si_pbesol.phph_interaction,
[1, 103],
temperatures=[300, ],
num_frequency_points=10)
sf.run()
print(sf.half_linewidths)
# np.testing.assert_allclose(
# gammas, si_pbesol.gammas.ravel(), atol=1e-2)
# np.testing.assert_allclose(
# freq_points, si_pbesol.frequency_points.ravel(),
# atol=1e-5)