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

Name of frequency shift is changed to real self energy

This commit is contained in:
Atsushi Togo 2020-09-06 14:59:14 +09:00
parent 41ce7fe379
commit 43fac60ae5
13 changed files with 553 additions and 455 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

72
c/_phono3py.c
View File

@ -44,7 +44,7 @@
#include <phonoc_array.h>
#include <phonoc_const.h>
#include <phonon3_h/fc3.h>
#include <phonon3_h/frequency_shift.h>
#include <phonon3_h/real_self_energy.h>
#include <phonon3_h/interaction.h>
#include <phonon3_h/imag_self_energy_with_g.h>
#include <phonon3_h/pp_collision.h>
@ -67,10 +67,10 @@ static PyObject *
py_get_imag_self_energy_with_g(PyObject *self, PyObject *args);
static PyObject *
py_get_detailed_imag_self_energy_with_g(PyObject *self, PyObject *args);
static PyObject * py_get_frequency_shift_at_bands(PyObject *self,
PyObject *args);
static PyObject * py_get_frequency_shift_at_frequency_point(PyObject *self,
PyObject *args);
static PyObject * py_get_real_self_energy_at_bands(PyObject *self,
PyObject *args);
static PyObject * py_get_real_self_energy_at_frequency_point(PyObject *self,
PyObject *args);
static PyObject * py_get_collision_matrix(PyObject *self, PyObject *args);
static PyObject * py_get_reducible_collision_matrix(PyObject *self,
PyObject *args);
@ -160,14 +160,14 @@ static PyMethodDef _phono3py_methods[] = {
(PyCFunction)py_get_detailed_imag_self_energy_with_g,
METH_VARARGS,
"Detailed contribution to imaginary part of self energy at frequency points with g"},
{"frequency_shift_at_bands",
(PyCFunction)py_get_frequency_shift_at_bands,
{"real_self_energy_at_bands",
(PyCFunction)py_get_real_self_energy_at_bands,
METH_VARARGS,
"Phonon frequency shift from third order force constants"},
{"frequency_shift_at_frequency_point",
(PyCFunction)py_get_frequency_shift_at_frequency_point,
"Real part of self energy from third order force constants"},
{"real_self_energy_at_frequency_point",
(PyCFunction)py_get_real_self_energy_at_frequency_point,
METH_VARARGS,
"Phonon frequency shift from third order force constants at a frequency point"},
"Real part of self energy from third order force constants at a frequency point"},
{"collision_matrix",
(PyCFunction)py_get_collision_matrix,
METH_VARARGS,
@ -995,8 +995,8 @@ py_get_detailed_imag_self_energy_with_g(PyObject *self, PyObject *args)
Py_RETURN_NONE;
}
static PyObject * py_get_frequency_shift_at_bands(PyObject *self,
PyObject *args)
static PyObject * py_get_real_self_energy_at_bands(PyObject *self,
PyObject *args)
{
PyArrayObject *py_shift;
PyArrayObject *py_fc3_normal_squared;
@ -1035,16 +1035,16 @@ static PyObject * py_get_frequency_shift_at_bands(PyObject *self,
triplets = (size_t(*)[3])PyArray_DATA(py_triplets);
triplet_weights = (int*)PyArray_DATA(py_triplet_weights);
get_frequency_shift_at_bands(shift,
fc3_normal_squared,
band_indices,
frequencies,
triplets,
triplet_weights,
epsilon,
temperature,
unit_conversion_factor,
cutoff_frequency);
get_real_self_energy_at_bands(shift,
fc3_normal_squared,
band_indices,
frequencies,
triplets,
triplet_weights,
epsilon,
temperature,
unit_conversion_factor,
cutoff_frequency);
free(fc3_normal_squared);
fc3_normal_squared = NULL;
@ -1052,8 +1052,8 @@ static PyObject * py_get_frequency_shift_at_bands(PyObject *self,
Py_RETURN_NONE;
}
static PyObject * py_get_frequency_shift_at_frequency_point(PyObject *self,
PyObject *args)
static PyObject * py_get_real_self_energy_at_frequency_point(PyObject *self,
PyObject *args)
{
PyArrayObject *py_shift;
PyArrayObject *py_fc3_normal_squared;
@ -1094,17 +1094,17 @@ static PyObject * py_get_frequency_shift_at_frequency_point(PyObject *self,
triplets = (size_t(*)[3])PyArray_DATA(py_triplets);
triplet_weights = (int*)PyArray_DATA(py_triplet_weights);
get_frequency_shift_at_frequency_point(shift,
frequency_point,
fc3_normal_squared,
band_indices,
frequencies,
triplets,
triplet_weights,
epsilon,
temperature,
unit_conversion_factor,
cutoff_frequency);
get_real_self_energy_at_frequency_point(shift,
frequency_point,
fc3_normal_squared,
band_indices,
frequencies,
triplets,
triplet_weights,
epsilon,
temperature,
unit_conversion_factor,
cutoff_frequency);
free(fc3_normal_squared);
fc3_normal_squared = NULL;

275
c/anharmonic/phonon3/frequency_shift.c
View File

@ -1,275 +0,0 @@
/* Copyright (C) 2015 Atsushi Togo */
/* All rights reserved. */
/* This file is part of phonopy. */
/* 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. */
#include <stdlib.h>
#include <math.h>
#include <phonoc_array.h>
#include <phonoc_utils.h>
#include <phonon3_h/frequency_shift.h>
#include <phonon3_h/real_to_reciprocal.h>
static double get_frequency_shift_at_band(const int band_index,
const Darray *fc3_normal_squared,
const double fpoint,
const double *frequencies,
const size_t (*triplets)[3],
const int *triplet_weights,
const double epsilon,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency);
static double sum_frequency_shift_at_band(const int num_band,
const double *fc3_normal_squared,
const double fpoint,
const double *freqs1,
const double *freqs2,
const double epsilon,
const double temperature,
const double cutoff_frequency);
static double sum_frequency_shift_at_band_0K(const int num_band,
const double *fc3_normal_squared,
const double fpoint,
const double *freqs1,
const double *freqs2,
const double epsilon,
const double cutoff_frequency);
void get_frequency_shift_at_bands(double *frequency_shift,
const Darray *fc3_normal_squared,
const int *band_indices,
const double *frequencies,
const size_t (*triplets)[3],
const int *triplet_weights,
const double epsilon,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency)
{
int i, num_band0, num_band, gp0;
double fpoint;
num_band0 = fc3_normal_squared->dims[1];
num_band = fc3_normal_squared->dims[2];
gp0 = triplets[0][0];
/* num_band0 and num_band_indices have to be same. */
for (i = 0; i < num_band0; i++) {
fpoint = frequencies[gp0 * num_band + band_indices[i]];
if (fpoint < cutoff_frequency) {
frequency_shift[i] = 0;
} else {
frequency_shift[i] =
get_frequency_shift_at_band(i,
fc3_normal_squared,
fpoint,
frequencies,
triplets,
triplet_weights,
epsilon,
temperature,
unit_conversion_factor,
cutoff_frequency);
}
}
}
void
get_frequency_shift_at_frequency_point(double *frequency_shift,
const double frequency_point,
const Darray *fc3_normal_squared,
const int *band_indices,
const double *frequencies,
const size_t (*triplets)[3],
const int *triplet_weights,
const double epsilon,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency)
{
int i, num_band0;
num_band0 = fc3_normal_squared->dims[1];
/* num_band0 and num_band_indices have to be same. */
for (i = 0; i < num_band0; i++) {
if (frequency_point < cutoff_frequency) {
frequency_shift[i] = 0;
} else {
frequency_shift[i] =
get_frequency_shift_at_band(i,
fc3_normal_squared,
frequency_point,
frequencies,
triplets,
triplet_weights,
epsilon,
temperature,
unit_conversion_factor,
cutoff_frequency);
}
}
}
static double get_frequency_shift_at_band(const int band_index,
const Darray *fc3_normal_squared,
const double fpoint,
const double *frequencies,
const size_t (*triplets)[3],
const int *triplet_weights,
const double epsilon,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency)
{
int i, num_triplets, num_band0, num_band, gp1, gp2;
double shift;
num_triplets = fc3_normal_squared->dims[0];
num_band0 = fc3_normal_squared->dims[1];
num_band = fc3_normal_squared->dims[2];
shift = 0;
#pragma omp parallel for private(gp1, gp2) reduction(+:shift)
for (i = 0; i < num_triplets; i++) {
gp1 = triplets[i][1];
gp2 = triplets[i][2];
if (temperature > 0) {
shift +=
sum_frequency_shift_at_band(num_band,
fc3_normal_squared->data +
i * num_band0 * num_band * num_band +
band_index * num_band * num_band,
fpoint,
frequencies + gp1 * num_band,
frequencies + gp2 * num_band,
epsilon,
temperature,
cutoff_frequency) *
triplet_weights[i] * unit_conversion_factor;
} else {
shift +=
sum_frequency_shift_at_band_0K(num_band,
fc3_normal_squared->data +
i * num_band0 * num_band * num_band +
band_index * num_band * num_band,
fpoint,
frequencies + gp1 * num_band,
frequencies + gp2 * num_band,
epsilon,
cutoff_frequency) *
triplet_weights[i] * unit_conversion_factor;
}
}
return shift;
}
static double sum_frequency_shift_at_band(const int num_band,
const double *fc3_normal_squared,
const double fpoint,
const double *freqs1,
const double *freqs2,
const double epsilon,
const double temperature,
const double cutoff_frequency)
{
int i, j;
double n1, n2, f1, f2, f3, f4, shift;
/* double sum; */
shift = 0;
for (i = 0; i < num_band; i++) {
if (freqs1[i] > cutoff_frequency) {
n1 = bose_einstein(freqs1[i], temperature);
for (j = 0; j < num_band; j++) {
if (freqs2[j] > cutoff_frequency) {
n2 = bose_einstein(freqs2[j], temperature);
f1 = fpoint + freqs1[i] + freqs2[j];
f2 = fpoint - freqs1[i] - freqs2[j];
f3 = fpoint - freqs1[i] + freqs2[j];
f4 = fpoint + freqs1[i] - freqs2[j];
/* sum = 0;
* if (fabs(f1) > epsilon) {
* sum -= (n1 + n2 + 1) / f1;
* }
* if (fabs(f2) > epsilon) {
* sum += (n1 + n2 + 1) / f2;
* }
* if (fabs(f3) > epsilon) {
* sum -= (n1 - n2) / f3;
* }
* if (fabs(f4) > epsilon) {
* sum += (n1 - n2) / f4;
* }
* shift += sum * fc3_normal_squared[i * num_band + j]; */
shift += (- (n1 + n2 + 1) * f1 / (f1 * f1 + epsilon * epsilon)
+ (n1 + n2 + 1) * f2 / (f2 * f2 + epsilon * epsilon)
- (n1 - n2) * f3 / (f3 * f3 + epsilon * epsilon)
+ (n1 - n2) * f4 / (f4 * f4 + epsilon * epsilon)) *
fc3_normal_squared[i * num_band + j];
}
}
}
}
return shift;
}
static double sum_frequency_shift_at_band_0K(const int num_band,
const double *fc3_normal_squared,
const double fpoint,
const double *freqs1,
const double *freqs2,
const double epsilon,
const double cutoff_frequency)
{
int i, j;
double f1, f2, shift;
shift = 0;
for (i = 0; i < num_band; i++) {
if (freqs1[i] > cutoff_frequency) {
for (j = 0; j < num_band; j++) {
if (freqs2[j] > cutoff_frequency) {
f1 = fpoint + freqs1[i] + freqs2[j];
f2 = fpoint - freqs1[i] - freqs2[j];
shift += (- 1 * f1 / (f1 * f1 + epsilon * epsilon)
+ 1 * f2 / (f2 * f2 + epsilon * epsilon)) *
fc3_normal_squared[i * num_band + j];
}
}
}
}
return shift;
}

275
c/anharmonic/phonon3/real_self_energy.c Normal file
View File

@ -0,0 +1,275 @@
/* Copyright (C) 2015 Atsushi Togo */
/* All rights reserved. */
/* This file is part of phonopy. */
/* 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. */
#include <stdlib.h>
#include <math.h>
#include <phonoc_array.h>
#include <phonoc_utils.h>
#include <phonon3_h/real_self_energy.h>
#include <phonon3_h/real_to_reciprocal.h>
static double get_real_self_energy_at_band(const int band_index,
const Darray *fc3_normal_squared,
const double fpoint,
const double *frequencies,
const size_t (*triplets)[3],
const int *triplet_weights,
const double epsilon,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency);
static double sum_real_self_energy_at_band(const int num_band,
const double *fc3_normal_squared,
const double fpoint,
const double *freqs1,
const double *freqs2,
const double epsilon,
const double temperature,
const double cutoff_frequency);
static double sum_real_self_energy_at_band_0K(const int num_band,
const double *fc3_normal_squared,
const double fpoint,
const double *freqs1,
const double *freqs2,
const double epsilon,
const double cutoff_frequency);
void get_real_self_energy_at_bands(double *real_self_energy,
const Darray *fc3_normal_squared,
const int *band_indices,
const double *frequencies,
const size_t (*triplets)[3],
const int *triplet_weights,
const double epsilon,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency)
{
int i, num_band0, num_band, gp0;
double fpoint;
num_band0 = fc3_normal_squared->dims[1];
num_band = fc3_normal_squared->dims[2];
gp0 = triplets[0][0];
/* num_band0 and num_band_indices have to be same. */
for (i = 0; i < num_band0; i++) {
fpoint = frequencies[gp0 * num_band + band_indices[i]];
if (fpoint < cutoff_frequency) {
real_self_energy[i] = 0;
} else {
real_self_energy[i] =
get_real_self_energy_at_band(i,
fc3_normal_squared,
fpoint,
frequencies,
triplets,
triplet_weights,
epsilon,
temperature,
unit_conversion_factor,
cutoff_frequency);
}
}
}
void
get_real_self_energy_at_frequency_point(double *real_self_energy,
const double frequency_point,
const Darray *fc3_normal_squared,
const int *band_indices,
const double *frequencies,
const size_t (*triplets)[3],
const int *triplet_weights,
const double epsilon,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency)
{
int i, num_band0;
num_band0 = fc3_normal_squared->dims[1];
/* num_band0 and num_band_indices have to be same. */
for (i = 0; i < num_band0; i++) {
if (frequency_point < cutoff_frequency) {
real_self_energy[i] = 0;
} else {
real_self_energy[i] =
get_real_self_energy_at_band(i,
fc3_normal_squared,
frequency_point,
frequencies,
triplets,
triplet_weights,
epsilon,
temperature,
unit_conversion_factor,
cutoff_frequency);
}
}
}
static double get_real_self_energy_at_band(const int band_index,
const Darray *fc3_normal_squared,
const double fpoint,
const double *frequencies,
const size_t (*triplets)[3],
const int *triplet_weights,
const double epsilon,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency)
{
int i, num_triplets, num_band0, num_band, gp1, gp2;
double shift;
num_triplets = fc3_normal_squared->dims[0];
num_band0 = fc3_normal_squared->dims[1];
num_band = fc3_normal_squared->dims[2];
shift = 0;
#pragma omp parallel for private(gp1, gp2) reduction(+:shift)
for (i = 0; i < num_triplets; i++) {
gp1 = triplets[i][1];
gp2 = triplets[i][2];
if (temperature > 0) {
shift +=
sum_real_self_energy_at_band(num_band,
fc3_normal_squared->data +
i * num_band0 * num_band * num_band +
band_index * num_band * num_band,
fpoint,
frequencies + gp1 * num_band,
frequencies + gp2 * num_band,
epsilon,
temperature,
cutoff_frequency) *
triplet_weights[i] * unit_conversion_factor;
} else {
shift +=
sum_real_self_energy_at_band_0K(num_band,
fc3_normal_squared->data +
i * num_band0 * num_band * num_band +
band_index * num_band * num_band,
fpoint,
frequencies + gp1 * num_band,
frequencies + gp2 * num_band,
epsilon,
cutoff_frequency) *
triplet_weights[i] * unit_conversion_factor;
}
}
return shift;
}
static double sum_real_self_energy_at_band(const int num_band,
const double *fc3_normal_squared,
const double fpoint,
const double *freqs1,
const double *freqs2,
const double epsilon,
const double temperature,
const double cutoff_frequency)
{
int i, j;
double n1, n2, f1, f2, f3, f4, shift;
/* double sum; */
shift = 0;
for (i = 0; i < num_band; i++) {
if (freqs1[i] > cutoff_frequency) {
n1 = bose_einstein(freqs1[i], temperature);
for (j = 0; j < num_band; j++) {
if (freqs2[j] > cutoff_frequency) {
n2 = bose_einstein(freqs2[j], temperature);
f1 = fpoint + freqs1[i] + freqs2[j];
f2 = fpoint - freqs1[i] - freqs2[j];
f3 = fpoint - freqs1[i] + freqs2[j];
f4 = fpoint + freqs1[i] - freqs2[j];
/* sum = 0;
* if (fabs(f1) > epsilon) {
* sum -= (n1 + n2 + 1) / f1;
* }
* if (fabs(f2) > epsilon) {
* sum += (n1 + n2 + 1) / f2;
* }
* if (fabs(f3) > epsilon) {
* sum -= (n1 - n2) / f3;
* }
* if (fabs(f4) > epsilon) {
* sum += (n1 - n2) / f4;
* }
* shift += sum * fc3_normal_squared[i * num_band + j]; */
shift += (- (n1 + n2 + 1) * f1 / (f1 * f1 + epsilon * epsilon)
+ (n1 + n2 + 1) * f2 / (f2 * f2 + epsilon * epsilon)
- (n1 - n2) * f3 / (f3 * f3 + epsilon * epsilon)
+ (n1 - n2) * f4 / (f4 * f4 + epsilon * epsilon)) *
fc3_normal_squared[i * num_band + j];
}
}
}
}
return shift;
}
static double sum_real_self_energy_at_band_0K(const int num_band,
const double *fc3_normal_squared,
const double fpoint,
const double *freqs1,
const double *freqs2,
const double epsilon,
const double cutoff_frequency)
{
int i, j;
double f1, f2, shift;
shift = 0;
for (i = 0; i < num_band; i++) {
if (freqs1[i] > cutoff_frequency) {
for (j = 0; j < num_band; j++) {
if (freqs2[j] > cutoff_frequency) {
f1 = fpoint + freqs1[i] + freqs2[j];
f2 = fpoint - freqs1[i] - freqs2[j];
shift += (- 1 * f1 / (f1 * f1 + epsilon * epsilon)
+ 1 * f2 / (f2 * f2 + epsilon * epsilon)) *
fc3_normal_squared[i * num_band + j];
}
}
}
}
return shift;
}

42
c/anharmonic_h/phonon3_h/frequency_shift.h → c/anharmonic_h/phonon3_h/real_self_energy.h
View File

@ -39,25 +39,25 @@
#endif
void get_frequency_shift_at_bands(double *frequency_shift,
const Darray *fc3_normal_squared,
const int *band_indices,
const double *frequencies,
const size_t (*triplets)[3],
const int *triplet_weights,
const double epsilon,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency);
void get_real_self_energy_at_bands(double *real_self_energy,
const Darray *fc3_normal_squared,
const int *band_indices,
const double *frequencies,
const size_t (*triplets)[3],
const int *triplet_weights,
const double epsilon,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency);
void
get_frequency_shift_at_frequency_point(double *frequency_shift,
const double frequency_point,
const Darray *fc3_normal_squared,
const int *band_indices,
const double *frequencies,
const size_t (*triplets)[3],
const int *triplet_weights,
const double epsilon,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency);
get_real_self_energy_at_frequency_point(double *real_self_energy,
const double frequency_point,
const Darray *fc3_normal_squared,
const int *band_indices,
const double *frequencies,
const size_t (*triplets)[3],
const int *triplet_weights,
const double epsilon,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency);

124
phono3py/api_phono3py.py
View File

@ -54,7 +54,8 @@ from phonopy.structure.grid_points import length2mesh
from phono3py.version import __version__
from phono3py.phonon3.imag_self_energy import (get_imag_self_energy,
write_imag_self_energy)
from phono3py.phonon3.frequency_shift import get_frequency_shift
from phono3py.phonon3.real_self_energy import get_real_self_energy
from phono3py.phonon3.spectral_function import SpectralFunction
from phono3py.phonon3.interaction import Interaction
from phono3py.phonon3.conductivity_RTA import get_thermal_conductivity_RTA
from phono3py.phonon3.conductivity_LBTE import get_thermal_conductivity_LBTE
@ -155,7 +156,7 @@ class Phono3py(object):
self._scattering_event_class = None
# Frequency shift (real part of bubble diagram)
self._frequency_shift = None
self._real_self_energy = None
self._grid_points = None
self._frequency_points = None
@ -1372,21 +1373,27 @@ class Phono3py(object):
def run_imag_self_energy(self,
grid_points,
temperatures,
frequency_points=None,
frequency_step=None,
num_frequency_points=None,
temperatures=None,
scattering_event_class=None,
write_gamma_detail=False,
keep_gamma_detail=False,
output_filename=None):
"""
"""Calculate imaginary part of self-energy of bubble diagram (Gamma)
Pi = Delta - i Gamma.
Parameters
----------
grid_points : array_like
Grid-point indices where imag-self-energeis are caclculated.
Grid-point indices where imaginary part of self-energies are
caclculated.
dtype=int, shape=(grid_points,)
temperatures : array_like
Temperatures where imaginary part of self-energies are calculated.
dtype=float, shape=(temperatures,)
frequency_points : array_like, optional
Frequency sampling points. Default is None. In this case,
num_frequency_points or frequency_step is used to generate uniform
@ -1399,10 +1406,6 @@ class Phono3py(object):
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
@ -1412,6 +1415,8 @@ class Phono3py(object):
Detailed gammas are written into a file in hdf5. Default is False.
keep_gamma_detail : bool, optional
With True, detailed gammas are stored. Default is False.
output_filename : str
This string is inserted in the output file names.
"""
@ -1459,25 +1464,55 @@ class Phono3py(object):
filename=filename,
is_mesh_symmetry=self._is_mesh_symmetry)
def run_frequency_shift(
def run_real_self_energy(
self,
grid_points,
temperatures,
run_on_bands=False,
frequency_points=None,
frequency_step=None,
num_frequency_points=None,
temperatures=None,
epsilons=None,
write_hdf5=True,
output_filename=None):
"""Frequency shift from lowest order diagram is calculated.
"""Calculate real-part of self-energy of bubble diagram (Delta)
Pi = Delta - i Gamma.
Parameters
----------
grid_points : array_like
Grid-point indices where real part of self-energies are
caclculated.
dtype=int, shape=(grid_points,)
temperatures : array_like
Temperatures where real part of self-energies are calculated.
dtype=float, shape=(temperatures,)
run_on_bands : bool, optional
With False, frequency shifts are calculated at frquency sampling
points. When True, they are done at the phonon frequencies.
Default is False.
frequency_points : array_like, optional
Frequency sampling points. Default is None. In this case,
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.
epsilons : array_like
The value to avoid divergence. When multiple values are given
frequency shifts for those values are returned.
Smearing widths to computer principal part. When multiple values
are given frequency shifts for those values are returned.
dtype=float, shape=(epsilons,)
write_hdf5 : bool
Results are stored in hdf5 files independently at grid points,
epsilons, and temperatures.
output_filename : str
This string is inserted in the output file names.
"""
@ -1498,7 +1533,7 @@ class Phono3py(object):
self._grid_points = grid_points
# (epsilon, grid_point, temperature, band)
self._frequency_points, self._frequency_shifts = get_frequency_shift(
self._frequency_points, self._real_self_energys = get_real_self_energy(
self._interaction,
self._grid_points,
temperatures,
@ -1511,7 +1546,64 @@ class Phono3py(object):
output_filename=output_filename,
log_level=self._log_level)
return self._frequency_points, self._frequency_shifts
return self._frequency_points, self._real_self_energys
def run_spectral_function(
self,
grid_points,
temperatures,
frequency_points=None,
frequency_step=None,
num_frequency_points=None,
write_hdf5=True,
output_filename=None):
"""Frequency shift from lowest order diagram is calculated.
Parameters
----------
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,)
frequency_points : array_like, optional
Frequency sampling points. Default is None. In this case,
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.
epsilons : array_like
The value to avoid divergence. When multiple values are given
frequency shifts for those values are returned.
dtype=float, shape=(epsilons,)
write_hdf5 : bool
Results are stored in hdf5 files independently at grid points,
epsilons, and temperatures.
output_filename : str
This string is inserted in the output file names.
"""
if self._interaction is None:
msg = ("Phono3py.init_phph_interaction has to be called "
"before running this method.")
raise RuntimeError(msg)
self._spectral_function = SpectralFunction(
self._interaction,
grid_points,
frequency_points=frequency_points,
frequency_step=frequency_step,
num_frequency_points=num_frequency_points,
temperatures=temperatures,
log_level=self._log_level)
def run_thermal_conductivity(
self,

27
phono3py/cui/phono3py_script.py
View File

@ -146,8 +146,10 @@ def get_run_mode(settings):
run_mode = "isotope"
elif settings.is_imag_self_energy:
run_mode = "imag_self_energy"
elif settings.is_frequency_shift:
run_mode = "frequency_shift"
elif settings.is_real_self_energy:
run_mode = "real_self_energy"
elif settings.is_spectral_function:
run_mode = "spectral_function"
elif settings.is_bterta:
run_mode = "conductivity-RTA"
elif settings.is_lbte:
@ -1095,9 +1097,9 @@ def main(**argparse_control):
if run_mode == "imag_self_energy":
phono3py.run_imag_self_energy(
updated_settings['grid_points'],
updated_settings['temperature_points'],
frequency_step=updated_settings['frequency_step'],
num_frequency_points=updated_settings['num_frequency_points'],
temperatures=updated_settings['temperature_points'],
scattering_event_class=settings.scattering_event_class,
write_gamma_detail=settings.write_gamma_detail,
output_filename=output_filename)
@ -1106,10 +1108,23 @@ def main(**argparse_control):
#####################################################
# Run frequency shift calculation of bubble diagram #
#####################################################
if run_mode == "frequency_shift":
phono3py.run_frequency_shift(
if run_mode == "real_self_energy":
phono3py.run_real_self_energy(
updated_settings['grid_points'],
temperatures=updated_settings['temperature_points'],
updated_settings['temperature_points'],
frequency_step=updated_settings['frequency_step'],
num_frequency_points=updated_settings['num_frequency_points'],
output_filename=output_filename)
#######################################################
# Run spectral function calculation of bubble diagram #
#######################################################
if run_mode == "spectral_function":
phono3py.run_spectral_function(
updated_settings['grid_points'],
updated_settings['temperature_points'],
frequency_step=updated_settings['frequency_step'],
num_frequency_points=updated_settings['num_frequency_points'],
output_filename=output_filename)
####################################

22
phono3py/file_IO.py
View File

@ -506,17 +506,17 @@ def write_linewidth_at_grid_point(gp,
w.write("%15.7f %20.15e\n" % (t, v))
def write_frequency_shift(gp,
band_indices,
frequency_points,
deltas,
mesh,
epsilon,
temperature,
filename=None,
is_mesh_symmetry=True):
def write_real_self_energy(gp,
band_indices,
frequency_points,
deltas,
mesh,
epsilon,
temperature,
filename=None,
is_mesh_symmetry=True):
fst_filename = "frequency_shift"
fst_filename = "real_self_energy"
suffix = _get_filename_suffix(mesh,
grid_point=gp,
band_indices=band_indices,
@ -574,7 +574,7 @@ def write_Delta_to_hdf5(grid_point,
w.create_dataset('grid_point', data=grid_point)
w.create_dataset('mesh', data=mesh)
w.create_dataset('band_index', data=_band_indices)
w.create_dataset('bubble', data=deltas)
w.create_dataset('delta', data=deltas)
w.create_dataset('temperature', data=temperatures)
w.create_dataset('epsilon', data=epsilon)
if frequency_points is not None:

110
phono3py/phonon3/frequency_shift.py → phono3py/phonon3/real_self_energy.py
View File

@ -37,21 +37,21 @@ import numpy as np
from phonopy.units import Hbar, EV, THz
from phonopy.phonon.degeneracy import degenerate_sets
from phono3py.phonon.func import bose_einstein
from phono3py.file_IO import write_frequency_shift, write_Delta_to_hdf5
from phono3py.file_IO import write_real_self_energy, write_Delta_to_hdf5
from phono3py.phonon3.imag_self_energy import get_frequency_points
def get_frequency_shift(interaction,
grid_points,
temperatures,
run_on_bands=False,
frequency_points=None,
frequency_step=None,
num_frequency_points=None,
epsilons=None,
output_filename=None,
write_hdf5=True,
log_level=0):
def get_real_self_energy(interaction,
grid_points,
temperatures,
run_on_bands=False,
frequency_points=None,
frequency_step=None,
num_frequency_points=None,
epsilons=None,
output_filename=None,
write_hdf5=True,
log_level=0):
if epsilons is None:
_epsilons = [None, ]
else:
@ -68,7 +68,7 @@ def get_frequency_shift(interaction,
print("Running harmonic phonon calculations...")
interaction.run_phonon_solver()
fst = FrequencyShift(interaction)
fst = RealSelfEnergy(interaction)
mesh = interaction.mesh_numbers
frequencies = interaction.get_phonons()[0]
max_phonon_freq = np.amax(frequencies)
@ -96,7 +96,7 @@ def get_frequency_shift(interaction,
fst.grid_point = gp
if log_level:
weights = interaction.get_triplets_at_q()[1]
print("------ Frequency shift -o- at grid point %d ------" % gp)
print("------ Re self-energy -o- at grid point %d ------" % gp)
print("Number of ir-triplets: %d / %d"
% (len(weights), weights.sum()))
@ -115,12 +115,12 @@ def get_frequency_shift(interaction,
for k, t in enumerate(_temperatures):
fst.temperature = t
fst.run()
all_deltas[i, j, k] = fst.frequency_shift
all_deltas[i, j, k] = fst.real_self_energy
if not run_on_bands:
pos = 0
for bi_set in [[bi, ] for bi in band_indices]:
filename = write_frequency_shift(
filename = write_real_self_energy(
gp,
bi_set,
_frequency_points,
@ -155,7 +155,7 @@ def get_frequency_shift(interaction,
return _frequency_points, all_deltas
class FrequencyShift(object):
class RealSelfEnergy(object):
default_epsilon = 0.05
@ -218,7 +218,7 @@ class FrequencyShift(object):
self._unit_conversion = None
self._cutoff_frequency = interaction.cutoff_frequency
self._frequency_points = None
self._frequency_shifts = None
self._real_self_energies = None
# Unit to THz of Delta
self._unit_conversion = (18 / (Hbar * EV) ** 2
@ -231,10 +231,10 @@ class FrequencyShift(object):
num_band0 = len(self._pp.band_indices)
if self._frequency_points is None:
self._frequency_shifts = np.zeros(num_band0, dtype='double')
self._real_self_energies = np.zeros(num_band0, dtype='double')
self._run_with_band_indices()
else:
self._frequency_shifts = np.zeros(
self._real_self_energies = np.zeros(
(len(self._frequency_points), num_band0), dtype='double')
self._run_with_frequency_points()
@ -248,11 +248,11 @@ class FrequencyShift(object):
self._band_indices = self._pp.band_indices
@property
def frequency_shift(self):
def real_self_energy(self):
if self._cutoff_frequency is None:
return self._frequency_shifts
return self._real_self_energies
else: # Averaging frequency shifts by degenerate bands
shifts = np.zeros_like(self._frequency_shifts)
shifts = np.zeros_like(self._real_self_energies)
freqs = self._frequencies[self._grid_point]
deg_sets = degenerate_sets(freqs) # like [[0,1], [2], [3,4,5]]
for dset in deg_sets:
@ -263,12 +263,13 @@ class FrequencyShift(object):
if len(bi_set) > 0:
for i in bi_set:
if self._frequency_points is None:
shifts[i] = (self._frequency_shifts[bi_set].sum() /
len(bi_set))
shifts[i] = (
self._real_self_energies[bi_set].sum()
/ len(bi_set))
else:
shifts[:, i] = (
self._frequency_shifts[:, bi_set].sum(axis=1) /
len(bi_set))
self._real_self_energies[:, bi_set].sum(axis=1)
/ len(bi_set))
return shifts
@property
@ -330,16 +331,16 @@ class FrequencyShift(object):
def _run_c_with_band_indices(self):
import phono3py._phono3py as phono3c
phono3c.frequency_shift_at_bands(self._frequency_shifts,
self._pp_strength,
self._triplets_at_q,
self._weights_at_q,
self._frequencies,
self._band_indices,
self._temperature,
self._epsilon,
self._unit_conversion,
self._cutoff_frequency)
phono3c.real_self_energy_at_bands(self._real_self_energies,
self._pp_strength,
self._triplets_at_q,
self._weights_at_q,
self._frequencies,
self._band_indices,
self._temperature,
self._epsilon,
self._unit_conversion,
self._cutoff_frequency)
def _run_py_with_band_indices(self):
for i, (triplet, w, interaction) in enumerate(
@ -351,21 +352,22 @@ class FrequencyShift(object):
for j, bi in enumerate(self._band_indices):
fpoint = freqs[0, bi]
if self._temperature > 0:
self._frequency_shifts[j] += (
self._frequency_shifts_at_bands(
self._real_self_energies[j] += (
self._real_self_energies_at_bands(
j, fpoint, freqs, interaction, w))
else:
self._frequency_shifts[j] += (
self._frequency_shifts_at_bands_0K(
self._real_self_energies[j] += (
self._real_self_energies_at_bands_0K(
j, fpoint, freqs, interaction, w))
self._frequency_shifts *= self._unit_conversion
self._real_self_energies *= self._unit_conversion
def _run_c_with_frequency_points(self):
import phono3py._phono3py as phono3c
for i, fpoint in enumerate(self._frequency_points):
shifts = np.zeros(self._frequency_shifts.shape[1], dtype='double')
phono3c.frequency_shift_at_frequency_point(
shifts = np.zeros(self._real_self_energies.shape[1],
dtype='double')
phono3c.real_self_energy_at_frequency_point(
shifts,
fpoint,
self._pp_strength,
@ -377,7 +379,7 @@ class FrequencyShift(object):
self._epsilon,
self._unit_conversion,
self._cutoff_frequency)
self._frequency_shifts[i][:] = shifts
self._real_self_energies[i][:] = shifts
def _run_py_with_frequency_points(self):
for k, fpoint in enumerate(self._frequency_points):
@ -389,18 +391,18 @@ class FrequencyShift(object):
freqs = self._frequencies[triplet]
for j, bi in enumerate(self._band_indices):
if self._temperature > 0:
self._frequency_shifts[k, j] += (
self._frequency_shifts_at_bands(
self._real_self_energies[k, j] += (
self._real_self_energies_at_bands(
j, fpoint, freqs, interaction, w))
else:
self._frequency_shifts[k, j] += (
self._frequency_shifts_at_bands_0K(
self._real_self_energies[k, j] += (
self._real_self_energies_at_bands_0K(
j, fpoint, freqs, interaction, w))
self._frequency_shifts *= self._unit_conversion
self._real_self_energies *= self._unit_conversion
def _frequency_shifts_at_bands(self, i, fpoint, freqs,
interaction, weight):
def _real_self_energies_at_bands(self, i, fpoint, freqs,
interaction, weight):
if fpoint < self._cutoff_frequency:
return 0
@ -435,8 +437,8 @@ class FrequencyShift(object):
sum_d += d * interaction[i, j, k] * weight
return sum_d
def _frequency_shifts_at_bands_0K(self, i, fpoint, freqs,
interaction, weight):
def _real_self_energies_at_bands_0K(self, i, fpoint, freqs,
interaction, weight):
if fpoint < self._cutoff_frequency:
return 0

37
phono3py/phonon3/spectral_function.py
View File

@ -32,29 +32,9 @@
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
import sys
import numpy as np
from phono3py.phonon3.imag_self_energy import get_imag_self_energy
from phono3py.phonon3.frequency_shift import imag_to_real
def get_spectral_function(interaction,
grid_points,
temperatures,
frequency_points=None,
frequency_step=None,
num_frequency_points=None,
output_filename=None,
write_hdf5=True,
log_level=0):
sf = SpectralFunction(interaction,
grid_points,
frequency_points=frequency_points,
frequency_step=frequency_step,
num_frequency_points=num_frequency_points,
temperatures=temperatures)
sf.run()
return sf
from phono3py.phonon3.real_self_energy import imag_to_real
class SpectralFunction(object):
@ -66,13 +46,16 @@ class SpectralFunction(object):
frequency_points=None,
frequency_step=None,
num_frequency_points=None,
temperatures=None):
temperatures=None,
log_level=0):
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
self._log_level = log_level
self._spectral_functions = None
self._gammas = None
self._deltas = None
@ -99,7 +82,14 @@ class SpectralFunction(object):
def frequency_points(self):
return self._frequency_points
@property
def grid_points(self):
return self._grid_points
def _run_gamma(self):
if self._log_level > 0:
print("Start calculation of imaginary part of self energies")
# gammas[grid_points, sigmas, temps, freq_points, band_indices]
fpoints, gammas = get_imag_self_energy(
self._interaction,
@ -107,7 +97,8 @@ class SpectralFunction(object):
frequency_points=self._frequency_points_in,
frequency_step=self._frequency_step,
num_frequency_points=self._num_frequency_points,
temperatures=self._temperatures)
temperatures=self._temperatures,
log_level=self._log_level)
self._gammas = np.array(gammas[:, 0], dtype='double', order='C')
self._frequency_points = fpoints

2
setup.py
View File

@ -41,7 +41,7 @@ sources = ['c/_phono3py.c',
'c/harmonic/phonoc_array.c',
'c/harmonic/phonoc_utils.c',
'c/anharmonic/phonon3/fc3.c',
'c/anharmonic/phonon3/frequency_shift.c',
'c/anharmonic/phonon3/real_self_energy.c',
'c/anharmonic/phonon3/interaction.c',
'c/anharmonic/phonon3/real_to_reciprocal.c',
'c/anharmonic/phonon3/reciprocal_to_normal.c',

11
test/phonon3/test_imag_self_energy.py
View File

@ -74,9 +74,8 @@ def test_imag_self_energy_npoints(si_pbesol):
si_pbesol.init_phph_interaction()
si_pbesol.run_imag_self_energy(
[1, 103],
temperatures=[300, ],
[300, ],
num_frequency_points=10)
print(np.array(si_pbesol.gammas).shape)
np.testing.assert_allclose(
gammas, si_pbesol.gammas.ravel(), atol=1e-2)
np.testing.assert_allclose(
@ -89,7 +88,7 @@ def test_imag_self_energy_freq_points(si_pbesol):
si_pbesol.init_phph_interaction()
si_pbesol.run_imag_self_energy(
[1, 103],
temperatures=[300, ],
[300, ],
frequency_points=freq_points)
np.testing.assert_allclose(
gammas, si_pbesol.gammas.ravel(), atol=1e-2)
@ -102,7 +101,7 @@ def test_imag_self_energy_detailed(si_pbesol):
si_pbesol.init_phph_interaction()
si_pbesol.run_imag_self_energy(
[1, ],
temperatures=[300, ],
[300, ],
frequency_points=freq_points,
keep_gamma_detail=True)
np.testing.assert_allclose(
@ -116,7 +115,7 @@ def test_imag_self_energy_scat_class1(si_pbesol):
si_pbesol.init_phph_interaction()
si_pbesol.run_imag_self_energy(
[1, 103],
temperatures=[300, ],
[300, ],
frequency_points=freq_points,
scattering_event_class=1)
# for line in si_pbesol.gammas.reshape(-1, 6):
@ -130,7 +129,7 @@ def test_imag_self_energy_scat_class2(si_pbesol):
si_pbesol.init_phph_interaction()
si_pbesol.run_imag_self_energy(
[1, 103],
temperatures=[300, ],
[300, ],
frequency_points=freq_points,
scattering_event_class=2)
# for line in si_pbesol.gammas.reshape(-1, 6):

8
test/phonon3/test_frequency_shift.py → test/phonon3/test_real_self_energy.py
View File

@ -1,5 +1,5 @@
import numpy as np
from phono3py.phonon3.frequency_shift import ImagToReal
from phono3py.phonon3.real_self_energy import ImagToReal
si_pbesol_Delta = [
[-0.0057666, -0.0057666, -0.01639729, -0.14809965,
@ -112,12 +112,12 @@ im_part = [
[30.69222190, 0.00000000, 30.69222190, 0.41736574, 30.84568301, 0.38214163]]
def test_frequency_shift(si_pbesol):
def test_real_self_energy(si_pbesol):
si_pbesol.mesh_numbers = [9, 9, 9]
si_pbesol.init_phph_interaction()
_, delta = si_pbesol.run_frequency_shift(
_, delta = si_pbesol.run_real_self_energy(
[1, 103],
temperatures=[300, ],
[300, ],
write_hdf5=False,
run_on_bands=True)
np.testing.assert_allclose(si_pbesol_Delta, delta[0, :, 0], atol=1e-5)

3
test/phonon3/test_spectral_function.py
View File

@ -57,7 +57,6 @@ def test_SpectralFunction(si_pbesol):
# for line in sf.spectral_functions.reshape(-1, 6):
# print(("%.7f, " * 6) % tuple(line))
# raise
print(np.where(np.abs(spec_funcs - sf.spectral_functions.ravel()) > 1e-2))
np.testing.assert_allclose(shifts, sf.shifts.ravel(), atol=1e-2)
np.testing.assert_allclose(spec_funcs, sf.spectral_functions.ravel(),
atol=1e-2)
atol=1e-4, rtol=1e-2)