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

Update for tetrahedron method implementing in C

This commit is contained in:
Atsushi Togo 2014-01-21 22:07:03 +09:00
parent 2ac0e9d3ff
commit d39a6381c1
4 changed files with 341 additions and 22 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

46
anharmonic/phonon3/imag_self_energy.py
View File

@ -126,7 +126,10 @@ class ImagSelfEnergy:
def _run_with_band_indices(self):
if self._sigma is None:
self._run_thm_with_band_indices()
if self._lang == 'C':
self._run_thm_c_with_band_indices()
else:
self._run_thm_py_with_band_indices()
else:
if self._lang == 'C':
self._run_c_with_band_indices()
@ -136,6 +139,10 @@ class ImagSelfEnergy:
def _run_with_frequency_points(self):
if self._sigma is None:
self._run_thm_with_frequency_points()
# if self._lang == 'C':
# self._run_c_thm_with_frequency_points()
# else:
# self._run_py_thm_with_frequency_points()
else:
if self._lang == 'C':
self._run_c_with_frequency_points()
@ -155,6 +162,19 @@ class ImagSelfEnergy:
self._unit_conversion,
self._cutoff_frequency)
def _run_thm_c_with_band_indices(self):
import anharmonic._phono3py as phono3c
phono3c.thm_imag_self_energy_at_bands(self._imag_self_energy,
self._fc3_normal_squared,
self._triplets_at_q,
self._weights_at_q,
self._frequencies,
self._band_indices,
self._temperature,
self._g,
self._unit_conversion,
self._cutoff_frequency)
def _run_c_with_frequency_points(self):
import anharmonic._phono3py as phono3c
for i, fpoint in enumerate(self._frequency_points):
@ -213,7 +233,7 @@ class ImagSelfEnergy:
return sum_g
def _run_thm_with_band_indices(self):
def _run_thm_py_with_band_indices(self):
if self._temperature > 0:
self._ise_thm_with_band_indices()
else:
@ -232,8 +252,8 @@ class ImagSelfEnergy:
f2 > self._cutoff_frequency):
n2 = n[i, 0, j]
n3 = n[i, 1, k]
g1 = self._g[0, i, j, k]
g2_g3 = self._g[1, i, j, k] # g2 - g3
g1 = self._g[i, :, j, k, 0]
g2_g3 = self._g[i, :, j, k, 1] # g2 - g3
self._imag_self_energy[:] += (
(n2 + n3 + 1) * g1 +
(n2 - n3) * (g2_g3)) * interaction[:, j, k] * w
@ -244,7 +264,7 @@ class ImagSelfEnergy:
for i, (w, interaction) in enumerate(zip(self._weights_at_q,
self._fc3_normal_squared)):
for j, k in list(np.ndindex(interaction.shape[1:])):
g1 = self._g[0, i, j, k]
g1 = self._g[i, :, j, k, 0]
self._imag_self_energy[:] += g1 * interaction[:, j, k] * w
self._imag_self_energy *= self._unit_conversion
@ -308,8 +328,8 @@ class ImagSelfEnergy:
f2 > self._cutoff_frequency):
n2 = occupation(f1, self._temperature)
n3 = occupation(f2, self._temperature)
g1 = self._g[0, i, j, k]
g2_g3 = self._g[1, i, j, k] # g2 - g3
g1 = self._g[i, :, j, k, 0]
g2_g3 = self._g[i, :, j, k, 1] # g2 - g3
for l in range(len(interaction)):
self._imag_self_energy[:, l] += (
(n2 + n3 + 1) * g1 +
@ -321,7 +341,7 @@ class ImagSelfEnergy:
for i, (w, interaction) in enumerate(zip(self._weights_at_q,
self._fc3_normal_squared)):
for j, k in list(np.ndindex(interaction.shape[1:])):
g1 = self._g[0, i, j, k]
g1 = self._g[i, :, j, k, 0]
for l in range(len(interaction)):
self._imag_self_energy[:, l] += g1 * interaction[l, j, k] * w
@ -362,7 +382,7 @@ class ImagSelfEnergy:
frequency_points = self._frequency_points
shape = self._fc3_normal_squared.shape
self._g = np.zeros(
(2, shape[0], shape[2], shape[3], len(frequency_points)),
(shape[0], len(frequency_points), shape[2], shape[3], 2),
dtype='double')
phono3c.triplets_integration_weights(
@ -392,7 +412,7 @@ class ImagSelfEnergy:
frequency_points = self._frequency_points
shape = self._fc3_normal_squared.shape
self._g = np.zeros(
(2, shape[0], shape[2], shape[3], len(frequency_points)),
(shape[0], len(frequency_points), shape[2], shape[3], 2),
dtype='double')
for i, vertices in enumerate(tetrahedra_vertices):
@ -401,10 +421,10 @@ class ImagSelfEnergy:
f2_v = self._frequencies[vertices[1], k]
thm.set_tetrahedra_omegas(f1_v + f2_v)
thm.run(frequency_points)
self._g[0, i, j, k] = thm.get_integration_weight()
self._g[i, :, j, k, 0] = thm.get_integration_weight()
thm.set_tetrahedra_omegas(-f1_v + f2_v)
thm.run(frequency_points)
self._g[1, i, j, k] = thm.get_integration_weight()
self._g[i, :, j, k, 0] = thm.get_integration_weight()
thm.set_tetrahedra_omegas(f1_v - f2_v)
thm.run(frequency_points)
self._g[1, i, j, k] -= thm.get_integration_weight()
self._g[i, :, j, k, 0] -= thm.get_integration_weight()

70
c/_phono3py.c
View File

@ -10,6 +10,7 @@
#include "phonon3_h/fc3.h"
#include "phonon3_h/interaction.h"
#include "phonon3_h/imag_self_energy.h"
#include "phonon3_h/imag_self_energy_thm.h"
#include "other_h/isotope.h"
#include "spglib_h/kpoint.h"
#include "spglib_h/tetrahedron_method.h"
@ -21,6 +22,8 @@ static PyObject * py_get_interaction(PyObject *self, PyObject *args);
static PyObject * py_get_imag_self_energy(PyObject *self, PyObject *args);
static PyObject * py_get_imag_self_energy_at_bands(PyObject *self,
PyObject *args);
static PyObject * py_get_thm_imag_self_energy_at_bands(PyObject *self,
PyObject *args);
static PyObject * py_set_phonons_at_gridpoints(PyObject *self, PyObject *args);
static PyObject * py_get_phonon(PyObject *self, PyObject *args);
static PyObject * py_distribute_fc3(PyObject *self, PyObject *args);
@ -44,6 +47,7 @@ static PyMethodDef functions[] = {
{"interaction", py_get_interaction, METH_VARARGS, "Interaction of triplets"},
{"imag_self_energy", py_get_imag_self_energy, METH_VARARGS, "Imaginary part of self energy"},
{"imag_self_energy_at_bands", py_get_imag_self_energy_at_bands, METH_VARARGS, "Imaginary part of self energy at phonon frequencies of bands"},
{"thm_imag_self_energy_at_bands", py_get_thm_imag_self_energy_at_bands, METH_VARARGS, "Imaginary part of self energy at phonon frequencies of bands for tetrahedron method"},
{"phonons_at_gridpoints", py_set_phonons_at_gridpoints, METH_VARARGS, "Set phonons at grid points"},
{"phonon", py_get_phonon, METH_VARARGS, "Get phonon"},
{"distribute_fc3", py_distribute_fc3, METH_VARARGS, "Distribute least fc3 to full fc3"},
@ -443,6 +447,57 @@ static PyObject * py_get_imag_self_energy_at_bands(PyObject *self,
Py_RETURN_NONE;
}
static PyObject * py_get_thm_imag_self_energy_at_bands(PyObject *self,
PyObject *args)
{
PyArrayObject* gamma_py;
PyArrayObject* fc3_normal_squared_py;
PyArrayObject* frequencies_py;
PyArrayObject* grid_point_triplets_py;
PyArrayObject* triplet_weights_py;
PyArrayObject* band_indices_py;
PyArrayObject* g_py;
double unit_conversion_factor, cutoff_frequency, temperature;
if (!PyArg_ParseTuple(args, "OOOOOOdOdd",
&gamma_py,
&fc3_normal_squared_py,
&grid_point_triplets_py,
&triplet_weights_py,
&frequencies_py,
&band_indices_py,
&temperature,
&g_py,
&unit_conversion_factor,
&cutoff_frequency)) {
return NULL;
}
Darray* fc3_normal_squared = convert_to_darray(fc3_normal_squared_py);
double* gamma = (double*)gamma_py->data;
double* g = (double*)g_py->data;
const double* frequencies = (double*)frequencies_py->data;
const int* band_indices = (int*)band_indices_py->data;
const int* grid_point_triplets = (int*)grid_point_triplets_py->data;
const int* triplet_weights = (int*)triplet_weights_py->data;
get_thm_imag_self_energy_at_bands(gamma,
fc3_normal_squared,
band_indices,
frequencies,
grid_point_triplets,
triplet_weights,
g,
temperature,
unit_conversion_factor,
cutoff_frequency);
free(fc3_normal_squared);
Py_RETURN_NONE;
}
static PyObject * py_get_isotope_strength(PyObject *self, PyObject *args)
{
PyArrayObject* gamma_py;
@ -698,22 +753,19 @@ py_set_triplets_integration_weights(PyObject *self, PyObject *args)
}
for (l = 0; l < num_band0; l++) {
f0 = frequency_points[l];
adrs_shift = (i * num_band * num_band * num_band0 +
b1 * num_band * num_band0 +
b2 * num_band0);
iw[adrs_shift + l] =
adrs_shift = (i * num_band0 * num_band * num_band +
l * num_band * num_band + b1 * num_band + b2) * 2;
iw[adrs_shift] =
thm_get_integration_weight(f0, freq_vertices[0], 'I');
adrs_shift += num_triplets * num_band * num_band * num_band0;
iw[adrs_shift + l] =
iw[adrs_shift + 1] =
thm_get_integration_weight(f0, freq_vertices[1], 'I');
iw[adrs_shift + l] -=
iw[adrs_shift + 1] -=
thm_get_integration_weight(f0, freq_vertices[2], 'I');
}
}
}
}
Py_RETURN_NONE;
}

221
c/anharmonic/phonon3/imag_self_energy_thm.c Normal file
View File

@ -0,0 +1,221 @@
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "phonoc_array.h"
#include "phonoc_utils.h"
#include "phonon3_h/imag_self_energy_thm.h"
static double
get_thm_imag_self_energy_at_band(double *imag_self_energy,
const int band_index,
const Darray *fc3_normal_sqared,
const double fpoint,
const double *frequencies,
const int *grid_point_triplets,
const int *triplet_weights,
const double *g,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency);
static double
sum_thm_imag_self_energy_at_band(const int num_band,
const double *fc3_normal_sqared,
const double fpoint,
const double *freqs0,
const double *freqs1,
const double *g,
const double temperature,
const double cutoff_frequency);
static double
sum_thm_imag_self_energy_at_band_0K(const int num_band,
const double *fc3_normal_sqared,
const double fpoint,
const double *freqs0,
const double *freqs1,
const double *g,
const double cutoff_frequency);
/* imag_self_energy[num_band0] */
/* fc3_normal_sqared[num_triplets, num_band0, num_band, num_band] */
void get_thm_imag_self_energy(double *imag_self_energy,
const Darray *fc3_normal_sqared,
const double fpoint,
const double *frequencies,
const int *grid_point_triplets,
const int *triplet_weights,
const double *g,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency)
{
int i, num_band0;
num_band0 = fc3_normal_sqared->dims[1];
for (i = 0; i < num_band0; i++) {
imag_self_energy[i] =
get_thm_imag_self_energy_at_band(imag_self_energy,
i,
fc3_normal_sqared,
fpoint,
frequencies,
grid_point_triplets,
triplet_weights,
g,
temperature,
unit_conversion_factor,
cutoff_frequency);
}
}
void get_thm_imag_self_energy_at_bands(double *imag_self_energy,
const Darray *fc3_normal_sqared,
const int *band_indices,
const double *frequencies,
const int *grid_point_triplets,
const int *triplet_weights,
const double *g,
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_sqared->dims[1];
num_band = fc3_normal_sqared->dims[2];
gp0 = grid_point_triplets[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]];
imag_self_energy[i] =
get_thm_imag_self_energy_at_band(imag_self_energy,
i,
fc3_normal_sqared,
fpoint,
frequencies,
grid_point_triplets,
triplet_weights,
g,
temperature,
unit_conversion_factor,
cutoff_frequency);
}
}
static double
get_thm_imag_self_energy_at_band(double *imag_self_energy,
const int band_index,
const Darray *fc3_normal_sqared,
const double fpoint,
const double *frequencies,
const int *grid_point_triplets,
const int *triplet_weights,
const double *g,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency)
{
int i, num_triplets, num_band0, num_band, gp1, gp2;
double sum_g;
num_triplets = fc3_normal_sqared->dims[0];
num_band0 = fc3_normal_sqared->dims[1];
num_band = fc3_normal_sqared->dims[2];
sum_g = 0;
#pragma omp parallel for private(gp1, gp2) reduction(+:sum_g)
for (i = 0; i < num_triplets; i++) {
gp1 = grid_point_triplets[i * 3 + 1];
gp2 = grid_point_triplets[i * 3 + 2];
if (temperature > 0) {
sum_g +=
sum_thm_imag_self_energy_at_band
(num_band,
fc3_normal_sqared->data +
i * num_band0 * num_band * num_band +
band_index * num_band * num_band,
fpoint,
frequencies + gp1 * num_band,
frequencies + gp2 * num_band,
g +
(i * num_band0 * num_band * num_band +
band_index * num_band * num_band) * 2,
temperature,
cutoff_frequency) * triplet_weights[i] * unit_conversion_factor;
} else {
sum_g +=
sum_thm_imag_self_energy_at_band_0K
(num_band,
fc3_normal_sqared->data +
i * num_band0 * num_band * num_band +
band_index * num_band * num_band,
fpoint,
frequencies + gp1 * num_band,
frequencies + gp2 * num_band,
g +
(i * num_band0 * num_band * num_band +
band_index * num_band * num_band) * 2,
cutoff_frequency) * triplet_weights[i] * unit_conversion_factor;
}
}
return sum_g;
}
static double
sum_thm_imag_self_energy_at_band(const int num_band,
const double *fc3_normal_sqared,
const double fpoint,
const double *freqs0,
const double *freqs1,
const double *g,
const double temperature,
const double cutoff_frequency)
{
int i, j;
double n2, n3, g1, g2_3, sum_g;
sum_g = 0;
for (i = 0; i < num_band; i++) {
if (freqs0[i] > cutoff_frequency) {
n2 = bose_einstein(freqs0[i], temperature);
for (j = 0; j < num_band; j++) {
if (freqs1[j] > cutoff_frequency) {
n3 = bose_einstein(freqs1[j], temperature);
g1 = g[i * num_band * 2 + j * 2];
g2_3 = g[i * num_band * 2 + j * 2 + 1];
sum_g += ((n2 + n3 + 1) * g1 + (n2 - n3) * (g2_3)) *
fc3_normal_sqared[i * num_band + j];
}
}
}
}
return sum_g;
}
static double
sum_thm_imag_self_energy_at_band_0K(const int num_band,
const double *fc3_normal_sqared,
const double fpoint,
const double *freqs0,
const double *freqs1,
const double *g,
const double cutoff_frequency)
{
int i, j;
double g1, sum_g;
sum_g = 0;
for (i = 0; i < num_band; i++) {
if (freqs0[i] > cutoff_frequency) {
for (j = 0; j < num_band; j++) {
if (freqs1[j] > cutoff_frequency) {
g1 = g[i * num_band * 2 + j * 2];
sum_g += g1 * fc3_normal_sqared[i * num_band + j];
}
}
}
}
return sum_g;
}

26
c/anharmonic_h/phonon3_h/imag_self_energy_thm.h Normal file
View File

@ -0,0 +1,26 @@
#ifndef __imag_self_energy_thm_H__
#define __imag_self_energy_thm_H__
#include "phonoc_array.h"
void get_thm_imag_self_energy(double *gamma,
const Darray *fc3_normal_sqared,
const double fpoint,
const double *frequencies,
const int *grid_point_triplets,
const int *triplet_weights,
const double *g,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency);
void get_thm_imag_self_energy_at_bands(double *imag_self_energy,
const Darray *fc3_normal_sqared,
const int *band_indices,
const double *frequencies,
const int *grid_point_triplets,
const int *triplet_weights,
const double *g,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency);
#endif