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

Isotope scattering with tetrahedron method in C

This commit is contained in:
Atsushi Togo 2014-02-12 20:18:40 +09:00
parent 92c04f2aee
commit a688c08621
8 changed files with 378 additions and 53 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

106
anharmonic/other/isotope.py
View File

@ -19,7 +19,7 @@ class Isotope:
symprec=1e-5,
cutoff_frequency=None,
lapack_zheev_uplo='L'):
self._mesh = mesh
self._mesh = np.array(mesh, dtype='intc')
self._mass_variances = np.array(mass_variances, dtype='double')
self._primitive = primitive
self._band_indices = band_indices
@ -34,6 +34,7 @@ class Isotope:
self._nac_q_direction = None
self._grid_address = None
self._bz_map = None
self._grid_points = None
self._frequencies = None
@ -47,25 +48,21 @@ class Isotope:
def set_grid_point(self, grid_point):
self._grid_point = grid_point
self._reciprocal_lattice = np.linalg.inv(self._primitive.get_cell())
num_band = self._primitive.get_number_of_atoms() * 3
if self._band_indices is None:
self._band_indices = np.arange(num_band, dtype='intc')
else:
self._band_indices = np.array(self._band_indices, dtype='intc')
num_grid = np.prod(self._mesh)
self._grid_points = np.arange(num_grid, dtype='intc')
self._grid_points = np.arange(np.prod(self._mesh), dtype='intc')
if self._phonon_done is None:
self._phonon_done = np.zeros(num_grid, dtype='byte')
self._frequencies = np.zeros((num_grid, num_band), dtype='double')
self._eigenvectors = np.zeros((num_grid, num_band, num_band),
dtype='complex128')
if self._grid_address is None:
primitive_lattice = np.linalg.inv(self._primitive.get_cell())
self._grid_address = get_bz_grid_address(self._mesh,
primitive_lattice)
self._grid_address, self._bz_map = get_bz_grid_address(
self._mesh, primitive_lattice, with_boundary=True)
if self._phonon_done is None:
self._allocate_phonon()
def set_sigma(self, sigma):
if sigma is None:
@ -114,27 +111,71 @@ class Isotope:
self._nac_q_direction = np.array(nac_q_direction, dtype='double')
def _run_c(self):
self._set_phonon_c()
self._set_phonon_c(self._grid_points)
import anharmonic._phono3py as phono3c
gamma = np.zeros(len(self._band_indices), dtype='double')
phono3c.isotope_strength(gamma,
self._grid_point,
self._mass_variances,
self._frequencies,
self._eigenvectors,
self._band_indices,
np.prod(self._mesh),
self._sigma,
self._cutoff_frequency)
if self._sigma is None:
self._set_integration_weights()
weights = np.ones(len(self._grid_points), dtype='intc')
phono3c.thm_isotope_strength(gamma,
self._grid_point,
self._grid_points,
weights,
self._mass_variances,
self._frequencies,
self._eigenvectors,
self._band_indices,
self._integration_weights,
self._cutoff_frequency)
else:
phono3c.isotope_strength(gamma,
self._grid_point,
self._mass_variances,
self._frequencies,
self._eigenvectors,
self._band_indices,
np.prod(self._mesh),
self._sigma,
self._cutoff_frequency)
self._gamma = np.pi ** 2 / np.prod(self._mesh) * gamma
def _set_integration_weights(self):
thm = TetrahedronMethod(self._reciprocal_lattice, mesh=self._mesh)
primitive_lattice = np.linalg.inv(self._primitive.get_cell())
thm = TetrahedronMethod(primitive_lattice, mesh=self._mesh)
num_grid_points = len(self._grid_points)
num_band = self._primitive.get_number_of_atoms() * 3
self._integration_weights = np.zeros(
(len(self._band_indices), num_band, num_grid_points), dtype='double')
(num_grid_points, len(self._band_indices), num_band), dtype='double')
self._set_integration_weights_c(thm)
def _set_integration_weights_c(self, thm):
import anharmonic._phono3py as phono3c
unique_vertices = thm.get_unique_tetrahedra_vertices()
neighboring_grid_points = np.zeros(
len(unique_vertices) * len(self._grid_points), dtype='intc')
phono3c.neighboring_grid_points(
neighboring_grid_points,
self._grid_points,
unique_vertices,
self._mesh,
self._grid_address,
self._bz_map)
self._set_phonon_c(np.unique(neighboring_grid_points))
freq_points = np.array(
self._frequencies[self._grid_point, self._band_indices],
dtype='double', order='C')
phono3c.integration_weights(
self._integration_weights,
freq_points,
thm.get_tetrahedra(),
self._mesh,
self._grid_points,
self._frequencies,
self._grid_address,
self._bz_map)
def _set_integration_weights_py(self, thm):
for i, gp in enumerate(self._grid_points):
tfreqs = get_tetrahedra_frequencies(
gp,
@ -149,7 +190,7 @@ class Isotope:
thm.set_tetrahedra_omegas(frequencies)
thm.run(self._frequencies[self._grid_point, self._band_indices])
iw = thm.get_integration_weight()
self._integration_weights[:, bi, i] = iw
self._integration_weights[i, :, bi] = iw
def _run_py(self):
for gp in self._grid_points:
@ -173,19 +214,19 @@ class Isotope:
ti_sum_band = np.sum(np.abs(vec * vec0) ** 2 * mass_v)
if self._sigma is None:
ti_sum += ti_sum_band * self._integration_weights[
bi, j, i]
i, bi, j]
else:
ti_sum += ti_sum_band * gaussian(f0 - f, self._sigma)
t_inv.append(np.pi ** 2 / np.prod(self._mesh) * f0 ** 2 * ti_sum)
self._gamma = np.array(t_inv, dtype='double') / 2
def _set_phonon_c(self):
def _set_phonon_c(self, grid_points):
set_phonon_c(self._dm,
self._frequencies,
self._eigenvectors,
self._phonon_done,
self._grid_points,
grid_points,
self._grid_address,
self._mesh,
self._frequency_factor_to_THz,
@ -202,3 +243,12 @@ class Isotope:
self._dm,
self._frequency_factor_to_THz,
self._lapack_zheev_uplo)
def _allocate_phonon(self):
num_band = self._primitive.get_number_of_atoms() * 3
num_grid = len(self._grid_address)
self._phonon_done = np.zeros(num_grid, dtype='byte')
self._frequencies = np.zeros((num_grid, num_band), dtype='double')
self._eigenvectors = np.zeros((num_grid, num_band, num_band),
dtype='complex128')

6
anharmonic/phonon3/conductivity_RTA.py
View File

@ -421,7 +421,11 @@ class conductivity_RTA:
def _set_gamma_isotope_at_sigmas(self, i):
for j, sigma in enumerate(self._sigmas):
if self._log_level:
print "Calculating Gamma of ph-isotope with sigma=%s" % sigma
print "Calculating Gamma of ph-isotope with",
if sigma is None:
print "tetrahedron method"
else:
print "sigma=%s" % sigma
pp_freqs, pp_eigvecs, pp_phonon_done = self._pp.get_phonons()
self._isotope.set_sigma(sigma)
self._isotope.set_phonons(pp_freqs,

2
anharmonic/phonon3/interaction.py
View File

@ -335,7 +335,7 @@ class Interaction:
def _allocate_phonon(self):
primitive_lattice = np.linalg.inv(self._primitive.get_cell())
self._grid_address = get_bz_grid_address(
self._grid_address, self._bz_map = get_bz_grid_address(
self._mesh, primitive_lattice, with_boundary=True)
num_band = self._primitive.get_number_of_atoms() * 3
num_grid = len(self._grid_address)

2
anharmonic/phonon3/triplets.py
View File

@ -69,7 +69,7 @@ def get_bz_grid_address(mesh, primitive_lattice, with_boundary=False):
mesh,
primitive_lattice)
if with_boundary:
return bz_grid_address
return bz_grid_address, bz_map
else:
return bz_grid_address[:np.prod(mesh)]

208
c/_phono3py.c
View File

@ -29,9 +29,11 @@ static PyObject * py_get_phonon(PyObject *self, PyObject *args);
static PyObject * py_distribute_fc3(PyObject *self, PyObject *args);
static PyObject * py_phonopy_zheev(PyObject *self, PyObject *args);
static PyObject * py_get_isotope_strength(PyObject *self, PyObject *args);
static PyObject * py_get_thm_isotope_strength(PyObject *self, PyObject *args);
static PyObject * py_set_permutation_symmetry_fc3(PyObject *self,
PyObject *args);
static PyObject *py_get_neighboring_gird_points(PyObject *self, PyObject *args);
static PyObject * py_set_integration_weights(PyObject *self, PyObject *args);
static PyObject *
py_set_triplets_integration_weights(PyObject *self, PyObject *args);
static void get_triplet_tetrahedra_vertices
@ -53,8 +55,10 @@ static PyMethodDef functions[] = {
{"distribute_fc3", py_distribute_fc3, METH_VARARGS, "Distribute least fc3 to full fc3"},
{"zheev", py_phonopy_zheev, METH_VARARGS, "Lapack zheev wrapper"},
{"isotope_strength", py_get_isotope_strength, METH_VARARGS, "Isotope scattering strength"},
{"thm_isotope_strength", py_get_thm_isotope_strength, METH_VARARGS, "Isotope scattering strength for tetrahedron_method"},
{"permutation_symmetry_fc3", py_set_permutation_symmetry_fc3, METH_VARARGS, "Set permutation symmetry for fc3"},
{"neighboring_grid_points", py_get_neighboring_gird_points, METH_VARARGS, "Neighboring grid points by relative grid addresses"},
{"integration_weights", py_set_integration_weights, METH_VARARGS, "Integration weights of tetrahedron method"},
{"triplets_integration_weights", py_set_triplets_integration_weights, METH_VARARGS, "Integration weights of tetrahedron metod for triplets"},
{NULL, NULL, 0, NULL}
};
@ -499,7 +503,7 @@ static PyObject * py_get_isotope_strength(PyObject *self, PyObject *args)
PyArrayObject* gamma_py;
PyArrayObject* frequencies_py;
PyArrayObject* eigenvectors_py;
PyArrayObject* band_indicies_py;
PyArrayObject* band_indices_py;
PyArrayObject* mass_variances_py;
int grid_point;
int num_grid_points;
@ -512,7 +516,7 @@ static PyObject * py_get_isotope_strength(PyObject *self, PyObject *args)
&mass_variances_py,
&frequencies_py,
&eigenvectors_py,
&band_indicies_py,
&band_indices_py,
&num_grid_points,
&sigma,
&cutoff_frequency)) {
@ -524,22 +528,122 @@ static PyObject * py_get_isotope_strength(PyObject *self, PyObject *args)
const double* frequencies = (double*)frequencies_py->data;
const lapack_complex_double* eigenvectors =
(lapack_complex_double*)eigenvectors_py->data;
const int* band_indicies = (int*)band_indicies_py->data;
const int* band_indices = (int*)band_indices_py->data;
const double* mass_variances = (double*)mass_variances_py->data;
const int num_band = (int)frequencies_py->dimensions[1];
const int num_band0 = (int)band_indicies_py->dimensions[0];
const int num_band0 = (int)band_indices_py->dimensions[0];
get_isotope_scattering_strength(gamma,
grid_point,
mass_variances,
frequencies,
eigenvectors,
num_grid_points,
band_indicies,
num_band,
num_band0,
sigma,
cutoff_frequency);
int i, j, k;
double f, f0;
int *weights, *ir_grid_points;
double *integration_weights;
ir_grid_points = (int*)malloc(sizeof(int) * num_grid_points);
weights = (int*)malloc(sizeof(int) * num_grid_points);
integration_weights = (double*)malloc(sizeof(double) *
num_grid_points * num_band0 * num_band);
for (i = 0; i < num_grid_points; i++) {
ir_grid_points[i] = i;
weights[i] = 1;
for (j = 0; j < num_band0; j++) {
f0 = frequencies[grid_point * num_band + band_indices[j]];
for (k = 0; k < num_band; k++) {
f = frequencies[i * num_band + k];
integration_weights[i * num_band0 * num_band +
j * num_band + k] = gaussian(f - f0, sigma);
}
}
}
get_thm_isotope_scattering_strength(gamma,
grid_point,
ir_grid_points,
weights,
mass_variances,
frequencies,
eigenvectors,
num_grid_points,
band_indices,
num_band,
num_band0,
integration_weights,
cutoff_frequency);
free(ir_grid_points);
free(weights);
free(integration_weights);
/* get_isotope_scattering_strength(gamma, */
/* grid_point, */
/* mass_variances, */
/* frequencies, */
/* eigenvectors, */
/* num_grid_points, */
/* band_indices, */
/* num_band, */
/* num_band0, */
/* sigma, */
/* cutoff_frequency); */
Py_RETURN_NONE;
}
static PyObject * py_get_thm_isotope_strength(PyObject *self, PyObject *args)
{
PyArrayObject* gamma_py;
PyArrayObject* frequencies_py;
PyArrayObject* eigenvectors_py;
PyArrayObject* band_indices_py;
PyArrayObject* mass_variances_py;
PyArrayObject* ir_grid_points_py;
PyArrayObject* weights_py;
int grid_point;
double cutoff_frequency;
PyArrayObject* integration_weights_py;
if (!PyArg_ParseTuple(args, "OiOOOOOOOd",
&gamma_py,
&grid_point,
&ir_grid_points_py,
&weights_py,
&mass_variances_py,
&frequencies_py,
&eigenvectors_py,
&band_indices_py,
&integration_weights_py,
&cutoff_frequency)) {
return NULL;
}
double* gamma = (double*)gamma_py->data;
const double* frequencies = (double*)frequencies_py->data;
const int* ir_grid_points = (int*)ir_grid_points_py->data;
const int* weights = (int*)weights_py->data;
const lapack_complex_double* eigenvectors =
(lapack_complex_double*)eigenvectors_py->data;
const int* band_indices = (int*)band_indices_py->data;
const double* mass_variances = (double*)mass_variances_py->data;
const int num_band = (int)frequencies_py->dimensions[1];
const int num_band0 = (int)band_indices_py->dimensions[0];
const double* integration_weights = (double*)integration_weights_py->data;
const int num_ir_grid_points = (int)ir_grid_points_py->dimensions[0];
get_thm_isotope_scattering_strength(gamma,
grid_point,
ir_grid_points,
weights,
mass_variances,
frequencies,
eigenvectors,
num_ir_grid_points,
band_indices,
num_band,
num_band0,
integration_weights,
cutoff_frequency);
Py_RETURN_NONE;
}
@ -673,6 +777,72 @@ static PyObject * py_get_neighboring_gird_points(PyObject *self, PyObject *args)
Py_RETURN_NONE;
}
static PyObject * py_set_integration_weights(PyObject *self, PyObject *args)
{
PyArrayObject* iw_py;
PyArrayObject* frequency_points_py;
PyArrayObject* relative_grid_address_py;
PyArrayObject* mesh_py;
PyArrayObject* grid_points_py;
PyArrayObject* frequencies_py;
PyArrayObject* bz_grid_address_py;
PyArrayObject* bz_map_py;
if (!PyArg_ParseTuple(args, "OOOOOOOO",
&iw_py,
&frequency_points_py,
&relative_grid_address_py,
&mesh_py,
&grid_points_py,
&frequencies_py,
&bz_grid_address_py,
&bz_map_py)) {
return NULL;
}
double *iw = (double*)iw_py->data;
const double *frequency_points = (double*)frequency_points_py->data;
const int num_band0 = frequency_points_py->dimensions[0];
SPGCONST int (*relative_grid_address)[4][3] =
(int(*)[4][3])relative_grid_address_py->data;
const int *mesh = (int*)mesh_py->data;
SPGCONST int *grid_points = (int*)grid_points_py->data;
const int num_gp = (int)grid_points_py->dimensions[0];
SPGCONST int (*bz_grid_address)[3] = (int(*)[3])bz_grid_address_py->data;
const int *bz_map = (int*)bz_map_py->data;
const double *frequencies = (double*)frequencies_py->data;
const int num_band = (int)frequencies_py->dimensions[1];
int i, j, k, bi;
int vertices[24][4];
double freq_vertices[24][4];
#pragma omp parallel for private(j, k, bi, vertices, freq_vertices)
for (i = 0; i < num_gp; i++) {
for (j = 0; j < 24; j++) {
kpt_get_neighboring_grid_points(vertices[j],
grid_points[i],
relative_grid_address[j],
4,
mesh,
bz_grid_address,
bz_map);
}
for (bi = 0; bi < num_band; bi++) {
for (j = 0; j < 24; j++) {
for (k = 0; k < 4; k++) {
freq_vertices[j][k] = frequencies[vertices[j][k] * num_band + bi];
}
}
for (j = 0; j < num_band0; j++) {
iw[i * num_band0 * num_band + j * num_band + bi] =
thm_get_integration_weight(frequency_points[j], freq_vertices, 'I');
}
}
}
Py_RETURN_NONE;
}
static PyObject *
py_set_triplets_integration_weights(PyObject *self, PyObject *args)
{
@ -728,7 +898,7 @@ py_set_triplets_integration_weights(PyObject *self, PyObject *args)
}
}
#pragma omp parallel for private(j, k, l, b1, b2, sign, vertices, adrs_shift, f0, f1, f2, freq_vertices)
#pragma omp parallel for private(j, k, b1, b2, sign, vertices, adrs_shift, f0, f1, f2, freq_vertices)
for (i = 0; i < num_triplets; i++) {
get_triplet_tetrahedra_vertices(vertices,
tp_relative_grid_address,
@ -747,10 +917,10 @@ py_set_triplets_integration_weights(PyObject *self, PyObject *args)
freq_vertices[2][j][k] = f1 - f2;
}
}
for (l = 0; l < num_band0; l++) {
f0 = frequency_points[l];
for (j = 0; j < num_band0; j++) {
f0 = frequency_points[j];
adrs_shift = (i * num_band0 * num_band * num_band +
l * num_band * num_band + b1 * num_band + b2) * 2;
j * num_band * num_band + b1 * num_band + b2) * 2;
iw[adrs_shift] =
thm_get_integration_weight(f0, freq_vertices[0], 'I');
iw[adrs_shift + 1] =

86
c/anharmonic/other/isotope.c
View File

@ -75,4 +75,90 @@ void get_isotope_scattering_strength(double *gamma,
free(e0_i);
}
void
get_thm_isotope_scattering_strength(double *gamma,
const int grid_point,
const int *ir_grid_points,
const int *weights,
const double *mass_variances,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const int num_grid_points,
const int *band_indices,
const int num_band,
const int num_band0,
const double *integration_weights,
const double cutoff_frequency)
{
int i, j, k, l, gp;
double *e0_r, *e0_i, *f0, *gamma_ij;
double e1_r, e1_i, a, b, f, dist, sum_g_k;
e0_r = (double*)malloc(sizeof(double) * num_band * num_band0);
e0_i = (double*)malloc(sizeof(double) * num_band * num_band0);
f0 = (double*)malloc(sizeof(double) * num_band0);
for (i = 0; i < num_band0; i++) {
f0[i] = frequencies[grid_point * num_band + band_indices[i]];
for (j = 0; j < num_band; j++) {
e0_r[i * num_band + j] = lapack_complex_double_real
(eigenvectors[grid_point * num_band * num_band +
j * num_band + band_indices[i]]);
e0_i[i * num_band + j] = lapack_complex_double_imag
(eigenvectors[grid_point * num_band * num_band +
j * num_band + band_indices[i]]);
}
}
gamma_ij = (double*)malloc(sizeof(double) * num_grid_points * num_band0);
#pragma omp parallel for private(j, k, l, f, gp, e1_r, e1_i, a, b, dist, sum_g_k)
for (i = 0; i < num_grid_points; i++) {
gp = ir_grid_points[i];
for (j = 0; j < num_band0; j++) { /* band index0 */
if (f0[j] < cutoff_frequency) {
continue;
}
sum_g_k = 0;
for (k = 0; k < num_band; k++) { /* band index */
f = frequencies[gp * num_band + k];
if (f < cutoff_frequency) {
continue;
}
dist = integration_weights[gp * num_band0 * num_band + j * num_band + k];
for (l = 0; l < num_band; l++) { /* elements */
e1_r = lapack_complex_double_real
(eigenvectors[gp * num_band * num_band + l * num_band + k]);
e1_i = lapack_complex_double_imag
(eigenvectors[gp * num_band * num_band + l * num_band + k]);
a = e0_r[j * num_band + l] * e1_r + e0_i[j * num_band + l] * e1_i;
b = e0_i[j * num_band + l] * e1_r - e0_r[j * num_band + l] * e1_i;
sum_g_k += (a * a + b * b) * mass_variances[l / 3] * dist;
}
}
gamma_ij[gp * num_band0 + j] = sum_g_k * weights[gp];
}
}
for (i = 0; i < num_band0; i++) {
gamma[i] = 0;
}
for (i = 0; i < num_grid_points; i++) {
gp = ir_grid_points[i];
for (j = 0; j < num_band0; j++) {
gamma[j] += gamma_ij[gp * num_band0 + j];
}
}
for (i = 0; i < num_band0; i++) {
gamma[i] *= f0[i] * f0[i] / 2; /* gamma = 1/2t */
}
free(gamma_ij);
free(f0);
free(e0_r);
free(e0_i);
}

5
c/anharmonic/phonon3/imag_self_energy_thm.c
View File

@ -29,7 +29,7 @@ void get_thm_imag_self_energy_at_bands(double *imag_self_energy,
const double cutoff_frequency)
{
int i, j, num_triplets, num_band0, num_band, gp1, gp2;
double sum_g, f1, f2;
double f1, f2;
double *n1, *n2, *ise;
num_triplets = fc3_normal_sqared->dims[0];
@ -38,7 +38,7 @@ void get_thm_imag_self_energy_at_bands(double *imag_self_energy,
ise = (double*)malloc(sizeof(double) * num_triplets * num_band0);
#pragma omp parallel for private(j, gp1, gp2, sum_g, n1, n2, f1, f2)
#pragma omp parallel for private(j, gp1, gp2, n1, n2, f1, f2)
for (i = 0; i < num_triplets; i++) {
gp1 = triplets[i * 3 + 1];
gp2 = triplets[i * 3 + 2];
@ -140,4 +140,3 @@ sum_thm_imag_self_energy_at_band_0K(const int num_band,
}
return sum_g;
}

16
c/anharmonic_h/other_h/isotope.h
View File

@ -12,4 +12,20 @@ void get_isotope_scattering_strength(double *gamma,
const int num_band0,
const double sigma,
const double cutoff_frequency);
void
get_thm_isotope_scattering_strength(double *gamma,
const int grid_point,
const int *ir_grid_points,
const int *weights,
const double *mass_variances,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const int num_grid_points,
const int *band_indices,
const int num_band,
const int num_band0,
const double *integration_weights,
const double cutoff_frequency);
#endif