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

Formatted c files by VScode

This commit is contained in:
Atsushi Togo 2021-10-12 15:56:14 +09:00
parent e4c2d1396a
commit 2a474567d8
32 changed files with 3938 additions and 2755 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

90
c/_lapackepy.c
View File

@ -38,57 +38,59 @@
#include <numpy/arrayobject.h>
#include "lapack_wrapper.h"
static PyObject * py_phonopy_pinv(PyObject *self, PyObject *args);
static PyObject * py_phonopy_zheev(PyObject *self, PyObject *args);
static PyObject *py_phonopy_pinv(PyObject *self, PyObject *args);
static PyObject *py_phonopy_zheev(PyObject *self, PyObject *args);
struct module_state {
struct module_state
{
PyObject *error;
};
#if PY_MAJOR_VERSION >= 3
#define GETSTATE(m) ((struct module_state*)PyModule_GetState(m))
#define GETSTATE(m) ((struct module_state *)PyModule_GetState(m))
#else
#define GETSTATE(m) (&_state)
static struct module_state _state;
#endif
static PyObject *
error_out(PyObject *m) {
error_out(PyObject *m)
{
struct module_state *st = GETSTATE(m);
PyErr_SetString(st->error, "something bad happened");
return NULL;
}
static PyMethodDef _lapackepy_methods[] = {
{"error_out", (PyCFunction)error_out, METH_NOARGS, NULL},
{"pinv", py_phonopy_pinv, METH_VARARGS, "Pseudo-inverse using Lapack dgesvd"},
{"zheev", py_phonopy_zheev, METH_VARARGS, "Lapack zheev wrapper"},
{NULL, NULL, 0, NULL}
};
{"error_out", (PyCFunction)error_out, METH_NOARGS, NULL},
{"pinv", py_phonopy_pinv, METH_VARARGS, "Pseudo-inverse using Lapack dgesvd"},
{"zheev", py_phonopy_zheev, METH_VARARGS, "Lapack zheev wrapper"},
{NULL, NULL, 0, NULL}};
#if PY_MAJOR_VERSION >= 3
static int _lapackepy_traverse(PyObject *m, visitproc visit, void *arg) {
static int _lapackepy_traverse(PyObject *m, visitproc visit, void *arg)
{
Py_VISIT(GETSTATE(m)->error);
return 0;
}
static int _lapackepy_clear(PyObject *m) {
static int _lapackepy_clear(PyObject *m)
{
Py_CLEAR(GETSTATE(m)->error);
return 0;
}
static struct PyModuleDef moduledef = {
PyModuleDef_HEAD_INIT,
"_lapackepy",
NULL,
sizeof(struct module_state),
_lapackepy_methods,
NULL,
_lapackepy_traverse,
_lapackepy_clear,
NULL
};
PyModuleDef_HEAD_INIT,
"_lapackepy",
NULL,
sizeof(struct module_state),
_lapackepy_methods,
NULL,
_lapackepy_traverse,
_lapackepy_clear,
NULL};
#define INITERROR return NULL
@ -98,8 +100,7 @@ PyInit__lapackepy(void)
#else
#define INITERROR return
void
init_lapackepy(void)
void init_lapackepy(void)
#endif
{
#if PY_MAJOR_VERSION >= 3
@ -114,7 +115,8 @@ PyInit__lapackepy(void)
st = GETSTATE(module);
st->error = PyErr_NewException("_lapackepy.Error", NULL, NULL);
if (st->error == NULL) {
if (st->error == NULL)
{
Py_DECREF(module);
INITERROR;
}
@ -124,10 +126,10 @@ PyInit__lapackepy(void)
#endif
}
static PyObject * py_phonopy_zheev(PyObject *self, PyObject *args)
static PyObject *py_phonopy_zheev(PyObject *self, PyObject *args)
{
PyArrayObject* dynamical_matrix;
PyArrayObject* eigenvalues;
PyArrayObject *dynamical_matrix;
PyArrayObject *eigenvalues;
int dimension;
npy_cdouble *dynmat;
@ -137,36 +139,39 @@ static PyObject * py_phonopy_zheev(PyObject *self, PyObject *args)
if (!PyArg_ParseTuple(args, "OO",
&dynamical_matrix,
&eigenvalues)) {
&eigenvalues))
{
return NULL;
}
dimension = (int)PyArray_DIMS(dynamical_matrix)[0];
dynmat = (npy_cdouble*)PyArray_DATA(dynamical_matrix);
eigvals = (double*)PyArray_DATA(eigenvalues);
dynmat = (npy_cdouble *)PyArray_DATA(dynamical_matrix);
eigvals = (double *)PyArray_DATA(eigenvalues);
a = (lapack_complex_double*) malloc(sizeof(lapack_complex_double) *
a = (lapack_complex_double *)malloc(sizeof(lapack_complex_double) *
dimension * dimension);
for (i = 0; i < dimension * dimension; i++) {
for (i = 0; i < dimension * dimension; i++)
{
a[i] = lapack_make_complex_double(dynmat[i].real, dynmat[i].imag);
}
info = phonopy_zheev(eigvals, a, dimension, 'L');
for (i = 0; i < dimension * dimension; i++) {
for (i = 0; i < dimension * dimension; i++)
{
dynmat[i].real = lapack_complex_double_real(a[i]);
dynmat[i].imag = lapack_complex_double_imag(a[i]);
}
free(a);
return PyLong_FromLong((long) info);
return PyLong_FromLong((long)info);
}
static PyObject * py_phonopy_pinv(PyObject *self, PyObject *args)
static PyObject *py_phonopy_pinv(PyObject *self, PyObject *args)
{
PyArrayObject* data_in_py;
PyArrayObject* data_out_py;
PyArrayObject *data_in_py;
PyArrayObject *data_out_py;
double cutoff;
int m;
@ -178,16 +183,17 @@ static PyObject * py_phonopy_pinv(PyObject *self, PyObject *args)
if (!PyArg_ParseTuple(args, "OOd",
&data_out_py,
&data_in_py,
&cutoff)) {
&cutoff))
{
return NULL;
}
m = (int)PyArray_DIMS(data_in_py)[0];
n = (int)PyArray_DIMS(data_in_py)[1];
data_in = (double*)PyArray_DATA(data_in_py);
data_out = (double*)PyArray_DATA(data_out_py);
data_in = (double *)PyArray_DATA(data_in_py);
data_out = (double *)PyArray_DATA(data_out_py);
info = phonopy_pinv(data_out, data_in, m, n, cutoff);
return PyLong_FromLong((long) info);
return PyLong_FromLong((long)info);
}

1205
c/_phono3py.c
View File

File diff suppressed because it is too large Load Diff

199
c/_phononmod.c
View File

@ -37,59 +37,60 @@
#include "lapack_wrapper.h"
#include "phononmod.h"
static PyObject *py_get_phonons_at_gridpoints(PyObject *self, PyObject *args);
static PyObject * py_get_phonons_at_gridpoints(PyObject *self, PyObject *args);
struct module_state {
struct module_state
{
PyObject *error;
};
#if PY_MAJOR_VERSION >= 3
#define GETSTATE(m) ((struct module_state*)PyModule_GetState(m))
#define GETSTATE(m) ((struct module_state *)PyModule_GetState(m))
#else
#define GETSTATE(m) (&_state)
static struct module_state _state;
#endif
static PyObject *
error_out(PyObject *m) {
error_out(PyObject *m)
{
struct module_state *st = GETSTATE(m);
PyErr_SetString(st->error, "something bad happened");
return NULL;
}
static PyMethodDef _phononmod_methods[] = {
{"error_out", (PyCFunction)error_out, METH_NOARGS, NULL},
{"phonons_at_gridpoints",
py_get_phonons_at_gridpoints,
METH_VARARGS,
"Set phonons at grid points"},
{NULL, NULL, 0, NULL}
};
{"error_out", (PyCFunction)error_out, METH_NOARGS, NULL},
{"phonons_at_gridpoints",
py_get_phonons_at_gridpoints,
METH_VARARGS,
"Set phonons at grid points"},
{NULL, NULL, 0, NULL}};
#if PY_MAJOR_VERSION >= 3
static int _phononmod_traverse(PyObject *m, visitproc visit, void *arg) {
static int _phononmod_traverse(PyObject *m, visitproc visit, void *arg)
{
Py_VISIT(GETSTATE(m)->error);
return 0;
}
static int _phononmod_clear(PyObject *m) {
static int _phononmod_clear(PyObject *m)
{
Py_CLEAR(GETSTATE(m)->error);
return 0;
}
static struct PyModuleDef moduledef = {
PyModuleDef_HEAD_INIT,
"_phononmod",
NULL,
sizeof(struct module_state),
_phononmod_methods,
NULL,
_phononmod_traverse,
_phononmod_clear,
NULL
};
PyModuleDef_HEAD_INIT,
"_phononmod",
NULL,
sizeof(struct module_state),
_phononmod_methods,
NULL,
_phononmod_traverse,
_phononmod_clear,
NULL};
#define INITERROR return NULL
@ -99,8 +100,7 @@ PyInit__phononmod(void)
#else
#define INITERROR return
void
init_phononmod(void)
void init_phononmod(void)
#endif
{
#if PY_MAJOR_VERSION >= 3
@ -115,7 +115,8 @@ PyInit__phononmod(void)
st = GETSTATE(module);
st->error = PyErr_NewException("_phononmod.Error", NULL, NULL);
if (st->error == NULL) {
if (st->error == NULL)
{
Py_DECREF(module);
INITERROR;
}
@ -125,51 +126,51 @@ PyInit__phononmod(void)
#endif
}
static PyObject * py_get_phonons_at_gridpoints(PyObject *self, PyObject *args)
static PyObject *py_get_phonons_at_gridpoints(PyObject *self, PyObject *args)
{
PyArrayObject* py_frequencies;
PyArrayObject* py_eigenvectors;
PyArrayObject* py_phonon_done;
PyArrayObject* py_grid_points;
PyArrayObject* py_grid_address;
PyArrayObject* py_QDinv;
PyArrayObject* py_shortest_vectors_fc2;
PyArrayObject* py_multiplicity_fc2;
PyArrayObject* py_positions_fc2;
PyArrayObject* py_fc2;
PyArrayObject* py_masses_fc2;
PyArrayObject* py_p2s_map_fc2;
PyArrayObject* py_s2p_map_fc2;
PyArrayObject* py_reciprocal_lattice;
PyArrayObject* py_born_effective_charge;
PyArrayObject* py_q_direction;
PyArrayObject* py_dielectric_constant;
PyArrayObject* py_dd_q0;
PyArrayObject* py_G_list;
PyArrayObject *py_frequencies;
PyArrayObject *py_eigenvectors;
PyArrayObject *py_phonon_done;
PyArrayObject *py_grid_points;
PyArrayObject *py_grid_address;
PyArrayObject *py_QDinv;
PyArrayObject *py_shortest_vectors_fc2;
PyArrayObject *py_multiplicity_fc2;
PyArrayObject *py_positions_fc2;
PyArrayObject *py_fc2;
PyArrayObject *py_masses_fc2;
PyArrayObject *py_p2s_map_fc2;
PyArrayObject *py_s2p_map_fc2;
PyArrayObject *py_reciprocal_lattice;
PyArrayObject *py_born_effective_charge;
PyArrayObject *py_q_direction;
PyArrayObject *py_dielectric_constant;
PyArrayObject *py_dd_q0;
PyArrayObject *py_G_list;
double nac_factor;
double unit_conversion_factor;
double lambda;
char* uplo;
char *uplo;
double (*born)[3][3];
double (*dielectric)[3];
double(*born)[3][3];
double(*dielectric)[3];
double *q_dir;
double* freqs;
lapack_complex_double* eigvecs;
char* phonon_done;
long* grid_points;
long (*grid_address)[3];
double (*QDinv)[3];
double* fc2;
double *freqs;
lapack_complex_double *eigvecs;
char *phonon_done;
long *grid_points;
long(*grid_address)[3];
double(*QDinv)[3];
double *fc2;
double(*svecs_fc2)[3];
long (*multi_fc2)[2];
double (*positions_fc2)[3];
double* masses_fc2;
long* p2s_fc2;
long* s2p_fc2;
double (*rec_lat)[3];
double * dd_q0;
double (*G_list)[3];
long(*multi_fc2)[2];
double(*positions_fc2)[3];
double *masses_fc2;
long *p2s_fc2;
long *s2p_fc2;
double(*rec_lat)[3];
double *dd_q0;
double(*G_list)[3];
long num_patom, num_satom, num_phonons, num_grid_points, num_G_points;
if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOdOOOOdOOds",
@ -195,62 +196,82 @@ static PyObject * py_get_phonons_at_gridpoints(PyObject *self, PyObject *args)
&py_dd_q0,
&py_G_list,
&lambda,
&uplo)) {
&uplo))
{
return NULL;
}
freqs = (double*)PyArray_DATA(py_frequencies);
eigvecs = (lapack_complex_double*)PyArray_DATA(py_eigenvectors);
phonon_done = (char*)PyArray_DATA(py_phonon_done);
grid_points = (long*)PyArray_DATA(py_grid_points);
freqs = (double *)PyArray_DATA(py_frequencies);
eigvecs = (lapack_complex_double *)PyArray_DATA(py_eigenvectors);
phonon_done = (char *)PyArray_DATA(py_phonon_done);
grid_points = (long *)PyArray_DATA(py_grid_points);
grid_address = (long(*)[3])PyArray_DATA(py_grid_address);
QDinv = (double(*)[3])PyArray_DATA(py_QDinv);
fc2 = (double*)PyArray_DATA(py_fc2);
fc2 = (double *)PyArray_DATA(py_fc2);
svecs_fc2 = (double(*)[3])PyArray_DATA(py_shortest_vectors_fc2);
multi_fc2 = (long(*)[2])PyArray_DATA(py_multiplicity_fc2);
masses_fc2 = (double*)PyArray_DATA(py_masses_fc2);
p2s_fc2 = (long*)PyArray_DATA(py_p2s_map_fc2);
s2p_fc2 = (long*)PyArray_DATA(py_s2p_map_fc2);
masses_fc2 = (double *)PyArray_DATA(py_masses_fc2);
p2s_fc2 = (long *)PyArray_DATA(py_p2s_map_fc2);
s2p_fc2 = (long *)PyArray_DATA(py_s2p_map_fc2);
rec_lat = (double(*)[3])PyArray_DATA(py_reciprocal_lattice);
num_patom = (long)PyArray_DIMS(py_multiplicity_fc2)[1];
num_satom = (long)PyArray_DIMS(py_multiplicity_fc2)[0];
num_phonons = (long)PyArray_DIMS(py_frequencies)[0];
num_grid_points = (long)PyArray_DIMS(py_grid_points)[0];
if ((PyObject*)py_born_effective_charge == Py_None) {
if ((PyObject *)py_born_effective_charge == Py_None)
{
born = NULL;
} else {
}
else
{
born = (double(*)[3][3])PyArray_DATA(py_born_effective_charge);
}
if ((PyObject*)py_dielectric_constant == Py_None) {
if ((PyObject *)py_dielectric_constant == Py_None)
{
dielectric = NULL;
} else {
}
else
{
dielectric = (double(*)[3])PyArray_DATA(py_dielectric_constant);
}
if ((PyObject*)py_q_direction == Py_None) {
if ((PyObject *)py_q_direction == Py_None)
{
q_dir = NULL;
} else {
q_dir = (double*)PyArray_DATA(py_q_direction);
}
else
{
q_dir = (double *)PyArray_DATA(py_q_direction);
if (fabs(q_dir[0]) < 1e-10 &&
fabs(q_dir[1]) < 1e-10 &&
fabs(q_dir[2]) < 1e-10) {
fabs(q_dir[2]) < 1e-10)
{
q_dir = NULL;
}
}
if ((PyObject*)py_dd_q0 == Py_None) {
if ((PyObject *)py_dd_q0 == Py_None)
{
dd_q0 = NULL;
} else {
dd_q0 = (double*)PyArray_DATA(py_dd_q0);
}
if ((PyObject*)py_G_list == Py_None) {
else
{
dd_q0 = (double *)PyArray_DATA(py_dd_q0);
}
if ((PyObject *)py_G_list == Py_None)
{
G_list = NULL;
num_G_points = 0;
} else {
}
else
{
G_list = (double(*)[3])PyArray_DATA(py_G_list);
num_G_points = (long)PyArray_DIMS(py_G_list)[0];
}
if ((PyObject*)py_positions_fc2 == Py_None) {
if ((PyObject *)py_positions_fc2 == Py_None)
{
positions_fc2 = NULL;
} else {
}
else
{
positions_fc2 = (double(*)[3])PyArray_DATA(py_positions_fc2);
}

415
c/bzgrid.c
View File

@ -42,133 +42,131 @@
#define BZG_NUM_BZ_SEARCH_SPACE 125
#define GRID_TOLERANCE_FACTOR 0.01
static long bz_search_space[BZG_NUM_BZ_SEARCH_SPACE][3] = {
{ 0, 0, 0},
{ 0, 0, 1},
{ 0, 0, 2},
{ 0, 0, -2},
{ 0, 0, -1},
{ 0, 1, 0},
{ 0, 1, 1},
{ 0, 1, 2},
{ 0, 1, -2},
{ 0, 1, -1},
{ 0, 2, 0},
{ 0, 2, 1},
{ 0, 2, 2},
{ 0, 2, -2},
{ 0, 2, -1},
{ 0, -2, 0},
{ 0, -2, 1},
{ 0, -2, 2},
{ 0, -2, -2},
{ 0, -2, -1},
{ 0, -1, 0},
{ 0, -1, 1},
{ 0, -1, 2},
{ 0, -1, -2},
{ 0, -1, -1},
{ 1, 0, 0},
{ 1, 0, 1},
{ 1, 0, 2},
{ 1, 0, -2},
{ 1, 0, -1},
{ 1, 1, 0},
{ 1, 1, 1},
{ 1, 1, 2},
{ 1, 1, -2},
{ 1, 1, -1},
{ 1, 2, 0},
{ 1, 2, 1},
{ 1, 2, 2},
{ 1, 2, -2},
{ 1, 2, -1},
{ 1, -2, 0},
{ 1, -2, 1},
{ 1, -2, 2},
{ 1, -2, -2},
{ 1, -2, -1},
{ 1, -1, 0},
{ 1, -1, 1},
{ 1, -1, 2},
{ 1, -1, -2},
{ 1, -1, -1},
{ 2, 0, 0},
{ 2, 0, 1},
{ 2, 0, 2},
{ 2, 0, -2},
{ 2, 0, -1},
{ 2, 1, 0},
{ 2, 1, 1},
{ 2, 1, 2},
{ 2, 1, -2},
{ 2, 1, -1},
{ 2, 2, 0},
{ 2, 2, 1},
{ 2, 2, 2},
{ 2, 2, -2},
{ 2, 2, -1},
{ 2, -2, 0},
{ 2, -2, 1},
{ 2, -2, 2},
{ 2, -2, -2},
{ 2, -2, -1},
{ 2, -1, 0},
{ 2, -1, 1},
{ 2, -1, 2},
{ 2, -1, -2},
{ 2, -1, -1},
{-2, 0, 0},
{-2, 0, 1},
{-2, 0, 2},
{-2, 0, -2},
{-2, 0, -1},
{-2, 1, 0},
{-2, 1, 1},
{-2, 1, 2},
{-2, 1, -2},
{-2, 1, -1},
{-2, 2, 0},
{-2, 2, 1},
{-2, 2, 2},
{-2, 2, -2},
{-2, 2, -1},
{-2, -2, 0},
{-2, -2, 1},
{-2, -2, 2},
{-2, -2, -2},
{-2, -2, -1},
{-2, -1, 0},
{-2, -1, 1},
{-2, -1, 2},
{-2, -1, -2},
{-2, -1, -1},
{-1, 0, 0},
{-1, 0, 1},
{-1, 0, 2},
{-1, 0, -2},
{-1, 0, -1},
{-1, 1, 0},
{-1, 1, 1},
{-1, 1, 2},
{-1, 1, -2},
{-1, 1, -1},
{-1, 2, 0},
{-1, 2, 1},
{-1, 2, 2},
{-1, 2, -2},
{-1, 2, -1},
{-1, -2, 0},
{-1, -2, 1},
{-1, -2, 2},
{-1, -2, -2},
{-1, -2, -1},
{-1, -1, 0},
{-1, -1, 1},
{-1, -1, 2},
{-1, -1, -2},
{-1, -1, -1}
};
{0, 0, 0},
{0, 0, 1},
{0, 0, 2},
{0, 0, -2},
{0, 0, -1},
{0, 1, 0},
{0, 1, 1},
{0, 1, 2},
{0, 1, -2},
{0, 1, -1},
{0, 2, 0},
{0, 2, 1},
{0, 2, 2},
{0, 2, -2},
{0, 2, -1},
{0, -2, 0},
{0, -2, 1},
{0, -2, 2},
{0, -2, -2},
{0, -2, -1},
{0, -1, 0},
{0, -1, 1},
{0, -1, 2},
{0, -1, -2},
{0, -1, -1},
{1, 0, 0},
{1, 0, 1},
{1, 0, 2},
{1, 0, -2},
{1, 0, -1},
{1, 1, 0},
{1, 1, 1},
{1, 1, 2},
{1, 1, -2},
{1, 1, -1},
{1, 2, 0},
{1, 2, 1},
{1, 2, 2},
{1, 2, -2},
{1, 2, -1},
{1, -2, 0},
{1, -2, 1},
{1, -2, 2},
{1, -2, -2},
{1, -2, -1},
{1, -1, 0},
{1, -1, 1},
{1, -1, 2},
{1, -1, -2},
{1, -1, -1},
{2, 0, 0},
{2, 0, 1},
{2, 0, 2},
{2, 0, -2},
{2, 0, -1},
{2, 1, 0},
{2, 1, 1},
{2, 1, 2},
{2, 1, -2},
{2, 1, -1},
{2, 2, 0},
{2, 2, 1},
{2, 2, 2},
{2, 2, -2},
{2, 2, -1},
{2, -2, 0},
{2, -2, 1},
{2, -2, 2},
{2, -2, -2},
{2, -2, -1},
{2, -1, 0},
{2, -1, 1},
{2, -1, 2},
{2, -1, -2},
{2, -1, -1},
{-2, 0, 0},
{-2, 0, 1},
{-2, 0, 2},
{-2, 0, -2},
{-2, 0, -1},
{-2, 1, 0},
{-2, 1, 1},
{-2, 1, 2},
{-2, 1, -2},
{-2, 1, -1},
{-2, 2, 0},
{-2, 2, 1},
{-2, 2, 2},
{-2, 2, -2},
{-2, 2, -1},
{-2, -2, 0},
{-2, -2, 1},
{-2, -2, 2},
{-2, -2, -2},
{-2, -2, -1},
{-2, -1, 0},
{-2, -1, 1},
{-2, -1, 2},
{-2, -1, -2},
{-2, -1, -1},
{-1, 0, 0},
{-1, 0, 1},
{-1, 0, 2},
{-1, 0, -2},
{-1, 0, -1},
{-1, 1, 0},
{-1, 1, 1},
{-1, 1, 2},
{-1, 1, -2},
{-1, 1, -1},
{-1, 2, 0},
{-1, 2, 1},
{-1, 2, 2},
{-1, 2, -2},
{-1, 2, -1},
{-1, -2, 0},
{-1, -2, 1},
{-1, -2, 2},
{-1, -2, -2},
{-1, -2, -1},
{-1, -1, 0},
{-1, -1, 1},
{-1, -1, 2},
{-1, -1, -2},
{-1, -1, -1}};
static void get_bz_grid_addresses_type1(BZGrid *bzgrid,
const long Qinv[3][3]);
@ -192,7 +190,6 @@ static long get_inverse_unimodular_matrix_l3(long m[3][3],
const long a[3][3]);
static double norm_squared_d3(const double a[3]);
long bzg_rotate_grid_index(const long bz_grid_index,
const long rotation[3][3],
const ConstBZGrid *bzgrid)
@ -205,27 +202,35 @@ long bzg_rotate_grid_index(const long bz_grid_index,
grg_get_grid_address(adrs_rot, dadrs_rot, bzgrid->PS);
gp = grg_get_grid_index(adrs_rot, bzgrid->D_diag);
if (bzgrid->type == 1) {
if (bzgrid->addresses[gp][0] == adrs_rot[0] &&
bzgrid->addresses[gp][1] == adrs_rot[1] &&
bzgrid->addresses[gp][2] == adrs_rot[2]) {
if (bzgrid->type == 1)
{
if (bzgrid->addresses[gp][0] == adrs_rot[0] &&
bzgrid->addresses[gp][1] == adrs_rot[1] &&
bzgrid->addresses[gp][2] == adrs_rot[2])
{
return gp;
}
num_grgp = bzgrid->D_diag[0] * bzgrid->D_diag[1] * bzgrid->D_diag[2];
num_bzgp = num_grgp * 8;
for (i = bzgrid->gp_map[num_bzgp + gp] + num_grgp;
i < bzgrid->gp_map[num_bzgp + gp + 1] + num_grgp; i++) {
if (bzgrid->addresses[i][0] == adrs_rot[0] &&
bzgrid->addresses[i][1] == adrs_rot[1] &&
bzgrid->addresses[i][2] == adrs_rot[2]) {
i < bzgrid->gp_map[num_bzgp + gp + 1] + num_grgp; i++)
{
if (bzgrid->addresses[i][0] == adrs_rot[0] &&
bzgrid->addresses[i][1] == adrs_rot[1] &&
bzgrid->addresses[i][2] == adrs_rot[2])
{
return i;
}
}
} else {
for (i = bzgrid->gp_map[gp]; i < bzgrid->gp_map[gp + 1]; i++) {
if (bzgrid->addresses[i][0] == adrs_rot[0] &&
bzgrid->addresses[i][1] == adrs_rot[1] &&
bzgrid->addresses[i][2] == adrs_rot[2]) {
}
else
{
for (i = bzgrid->gp_map[gp]; i < bzgrid->gp_map[gp + 1]; i++)
{
if (bzgrid->addresses[i][0] == adrs_rot[0] &&
bzgrid->addresses[i][1] == adrs_rot[1] &&
bzgrid->addresses[i][2] == adrs_rot[2])
{
return i;
}
}
@ -235,27 +240,29 @@ long bzg_rotate_grid_index(const long bz_grid_index,
return bzgrid->gp_map[gp];
}
long bzg_get_bz_grid_addresses(BZGrid *bzgrid)
{
long det;
long Qinv[3][3];
det = get_inverse_unimodular_matrix_l3(Qinv, bzgrid->Q);
if (det == 0) {
if (det == 0)
{
return 0;
}
if (bzgrid->type == 1) {
if (bzgrid->type == 1)
{
get_bz_grid_addresses_type1(bzgrid, Qinv);
} else {
}
else
{
get_bz_grid_addresses_type2(bzgrid, Qinv);
}
return 1;
}
/* Note: Tolerance in squared distance. */
double bzg_get_tolerance_for_BZ_reduction(const BZGrid *bzgrid)
{
@ -264,25 +271,29 @@ double bzg_get_tolerance_for_BZ_reduction(const BZGrid *bzgrid)
double length[3];
double reclatQ[3][3];
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
for (i = 0; i < 3; i++)
{
for (j = 0; j < 3; j++)
{
reclatQ[i][j] =
bzgrid->reclat[i][0] * bzgrid->Q[0][j]
+ bzgrid->reclat[i][1] * bzgrid->Q[1][j]
+ bzgrid->reclat[i][2] * bzgrid->Q[2][j];
bzgrid->reclat[i][0] * bzgrid->Q[0][j] + bzgrid->reclat[i][1] * bzgrid->Q[1][j] + bzgrid->reclat[i][2] * bzgrid->Q[2][j];
}
}
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
length[i] = 0;
for (j = 0; j < 3; j++) {
for (j = 0; j < 3; j++)
{
length[i] += reclatQ[j][i] * reclatQ[j][i];
}
length[i] /= bzgrid->D_diag[i] * bzgrid->D_diag[i];
}
tolerance = length[0];
for (i = 1; i < 3; i++) {
if (tolerance < length[i]) {
for (i = 1; i < 3; i++)
{
if (tolerance < length[i])
{
tolerance = length[i];
}
}
@ -291,21 +302,23 @@ double bzg_get_tolerance_for_BZ_reduction(const BZGrid *bzgrid)
return tolerance;
}
RotMats * bzg_alloc_RotMats(const long size)
RotMats *bzg_alloc_RotMats(const long size)
{
RotMats *rotmats;
rotmats = NULL;
if ((rotmats = (RotMats*) malloc(sizeof(RotMats))) == NULL) {
if ((rotmats = (RotMats *)malloc(sizeof(RotMats))) == NULL)
{
warning_print("Memory could not be allocated.");
return NULL;
}
rotmats->size = size;
if (size > 0) {
if ((rotmats->mat = (long (*)[3][3]) malloc(sizeof(long[3][3]) * size))
== NULL) {
if (size > 0)
{
if ((rotmats->mat = (long(*)[3][3])malloc(sizeof(long[3][3]) * size)) == NULL)
{
warning_print("Memory could not be allocated ");
warning_print("(RotMats, line %d, %s).\n", __LINE__, __FILE__);
free(rotmats);
@ -316,9 +329,10 @@ RotMats * bzg_alloc_RotMats(const long size)
return rotmats;
}
void bzg_free_RotMats(RotMats * rotmats)
void bzg_free_RotMats(RotMats *rotmats)
{
if (rotmats->size > 0) {
if (rotmats->size > 0)
{
free(rotmats->mat);
rotmats->mat = NULL;
}
@ -331,10 +345,12 @@ void bzg_multiply_matrix_vector_ld3(double v[3],
{
long i;
double c[3];
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
c[i] = a[i][0] * b[0] + a[i][1] * b[1] + a[i][2] * b[2];
}
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
v[i] = c[i];
}
}
@ -349,12 +365,14 @@ static void get_bz_grid_addresses_type1(BZGrid *bzgrid,
long count, id_shift;
tolerance = bzg_get_tolerance_for_BZ_reduction(bzgrid);
for (j = 0; j < 3; j++) {
for (j = 0; j < 3; j++)
{
bzmesh[j] = bzgrid->D_diag[j] * 2;
}
num_bzmesh = bzmesh[0] * bzmesh[1] * bzmesh[2];
for (i = 0; i < num_bzmesh; i++) {
for (i = 0; i < num_bzmesh; i++)
{
bzgrid->gp_map[i] = num_bzmesh;
}
@ -365,16 +383,22 @@ static void get_bz_grid_addresses_type1(BZGrid *bzgrid,
/* with boundary_num_gp is unstable to store bz_grid_address. */
bzgrid->gp_map[num_bzmesh] = 0;
id_shift = 0;
for (i = 0; i < total_num_gp; i++) {
for (i = 0; i < total_num_gp; i++)
{
grg_get_grid_address_from_index(gr_adrs, i, bzgrid->D_diag);
min_distance = get_bz_distances(nint, distances, bzgrid,
gr_adrs, tolerance);
count = 0;
for (j = 0; j < BZG_NUM_BZ_SEARCH_SPACE; j++) {
if (distances[j] < min_distance + tolerance) {
if (count == 0) {
for (j = 0; j < BZG_NUM_BZ_SEARCH_SPACE; j++)
{
if (distances[j] < min_distance + tolerance)
{
if (count == 0)
{
gp = i;
} else {
}
else
{
gp = boundary_num_gp + total_num_gp;
boundary_num_gp++;
}
@ -385,11 +409,12 @@ static void get_bz_grid_addresses_type1(BZGrid *bzgrid,
bzgrid->D_diag,
nint,
Qinv);
for (k = 0; k < 3; k++) {
for (k = 0; k < 3; k++)
{
bz_address_double[k] = bzgrid->addresses[gp][k] * 2 + bzgrid->PS[k];
}
bzgp = grg_get_double_grid_index(
bz_address_double, bzmesh, bzgrid->PS);
bz_address_double, bzmesh, bzgrid->PS);
bzgrid->gp_map[bzgp] = gp;
bzgrid->bzg2grg[gp] = i;
}
@ -417,12 +442,15 @@ static void get_bz_grid_addresses_type2(BZGrid *bzgrid,
bzgrid->gp_map[0] = 0;
for (i = 0;
i < bzgrid->D_diag[0] * bzgrid->D_diag[1] * bzgrid->D_diag[2]; i++) {
i < bzgrid->D_diag[0] * bzgrid->D_diag[1] * bzgrid->D_diag[2]; i++)
{
grg_get_grid_address_from_index(gr_adrs, i, bzgrid->D_diag);
min_distance = get_bz_distances(nint, distances, bzgrid,
gr_adrs, tolerance);
for (j = 0; j < BZG_NUM_BZ_SEARCH_SPACE; j++) {
if (distances[j] < min_distance + tolerance) {
for (j = 0; j < BZG_NUM_BZ_SEARCH_SPACE; j++)
{
if (distances[j] < min_distance + tolerance)
{
set_bz_address(bzgrid->addresses[num_gp],
j,
gr_adrs,
@ -433,7 +461,7 @@ static void get_bz_grid_addresses_type2(BZGrid *bzgrid,
num_gp++;
}
}
bzgrid->gp_map[i + 1] = num_gp;
bzgrid->gp_map[i + 1] = num_gp;
}
bzgrid->size = num_gp;
@ -449,11 +477,13 @@ static void set_bz_address(long address[3],
long i;
long deltaG[3];
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
deltaG[i] = bz_search_space[bz_index][i] - nint[i];
}
lagmat_multiply_matrix_vector_l3(deltaG, Qinv, deltaG);
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
address[i] = grid_address[i] + deltaG[i] * D_diag[i];
}
}
@ -471,17 +501,21 @@ static double get_bz_distances(long nint[3],
grg_get_double_grid_address(dadrs, grid_address, bzgrid->PS);
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
q_red[i] = dadrs[i] / (2.0 * bzgrid->D_diag[i]);
}
bzg_multiply_matrix_vector_ld3(q_red, bzgrid->Q, q_red);
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
nint[i] = lagmat_Nint(q_red[i]);
q_red[i] -= nint[i];
}
for (i = 0; i < BZG_NUM_BZ_SEARCH_SPACE; i++) {
for (j = 0; j < 3; j++) {
for (i = 0; i < BZG_NUM_BZ_SEARCH_SPACE; i++)
{
for (j = 0; j < 3; j++)
{
q_vec[j] = q_red[j] + bz_search_space[i][j];
}
multiply_matrix_vector_d3(q_vec, bzgrid->reclat, q_vec);
@ -493,8 +527,10 @@ static double get_bz_distances(long nint[3],
* those translationally equivalent can change by very tiny numerical
* fluctuation. */
min_distance = distances[0];
for (i = 1; i < BZG_NUM_BZ_SEARCH_SPACE; i++) {
if (distances[i] < min_distance - tolerance) {
for (i = 1; i < BZG_NUM_BZ_SEARCH_SPACE; i++)
{
if (distances[i] < min_distance - tolerance)
{
min_distance = distances[i];
}
}
@ -508,10 +544,12 @@ static void multiply_matrix_vector_d3(double v[3],
{
long i;
double c[3];
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
c[i] = a[i][0] * b[0] + a[i][1] * b[1] + a[i][2] * b[2];
}
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
v[i] = c[i];
}
}
@ -523,7 +561,8 @@ static long get_inverse_unimodular_matrix_l3(long m[3][3],
long c[3][3];
det = lagmat_get_determinant_l3(a);
if (labs(det) != 1) {
if (labs(det) != 1)
{
return 0;
}

14
c/bzgrid.h
View File

@ -35,7 +35,8 @@
#ifndef __bzgrid_H__
#define __bzgrid_H__
typedef struct {
typedef struct
{
long size;
long (*mat)[3][3];
} RotMats;
@ -67,7 +68,8 @@ typedef struct {
* shape=(3, 3)
* type : long
* 1 or 2. */
typedef struct {
typedef struct
{
long size;
long D_diag[3];
long Q[3][3];
@ -79,7 +81,8 @@ typedef struct {
long type;
} BZGrid;
typedef struct {
typedef struct
{
long size;
long D_diag[3];
long Q[3][3];
@ -91,14 +94,13 @@ typedef struct {
long type;
} ConstBZGrid;
long bzg_rotate_grid_index(const long grid_index,
const long rotation[3][3],
const ConstBZGrid *bzgrid);
long bzg_get_bz_grid_addresses(BZGrid *bzgrid);
double bzg_get_tolerance_for_BZ_reduction(const BZGrid *bzgrid);
RotMats * bzg_alloc_RotMats(const long size);
void bzg_free_RotMats(RotMats * rotmats);
RotMats *bzg_alloc_RotMats(const long size);
void bzg_free_RotMats(RotMats *rotmats);
void bzg_multiply_matrix_vector_ld3(double v[3],
const long a[3][3],
const double b[3]);

156
c/collision_matrix.c
View File

@ -105,23 +105,23 @@ void col_get_collision_matrix(double *collision_matrix,
num_band = fc3_normal_squared->dims[2];
get_collision_matrix(
collision_matrix,
fc3_normal_squared->data,
num_band0,
num_band,
frequencies,
triplets,
triplets_map,
num_gp,
map_q,
rot_grid_points,
num_ir_gp,
num_rot,
rotations_cartesian,
g + 2 * num_triplets * num_band0 * num_band * num_band,
temperature,
unit_conversion_factor,
cutoff_frequency);
collision_matrix,
fc3_normal_squared->data,
num_band0,
num_band,
frequencies,
triplets,
triplets_map,
num_gp,
map_q,
rot_grid_points,
num_ir_gp,
num_rot,
rotations_cartesian,
g + 2 * num_triplets * num_band0 * num_band * num_band,
temperature,
unit_conversion_factor,
cutoff_frequency);
}
void col_get_reducible_collision_matrix(double *collision_matrix,
@ -143,19 +143,19 @@ void col_get_reducible_collision_matrix(double *collision_matrix,
num_band = fc3_normal_squared->dims[2];
get_reducible_collision_matrix(
collision_matrix,
fc3_normal_squared->data,
num_band0,
num_band,
frequencies,
triplets,
triplets_map,
num_gp,
map_q,
g + 2 * num_triplets * num_band0 * num_band * num_band,
temperature,
unit_conversion_factor,
cutoff_frequency);
collision_matrix,
fc3_normal_squared->data,
num_band0,
num_band,
frequencies,
triplets,
triplets_map,
num_gp,
map_q,
g + 2 * num_triplets * num_band0 * num_band * num_band,
temperature,
unit_conversion_factor,
cutoff_frequency);
}
static void get_collision_matrix(double *collision_matrix,
@ -186,9 +186,11 @@ static void get_collision_matrix(double *collision_matrix,
#ifdef PHPYOPENMP
#pragma omp parallel for private(j, k, l, m, n, ti, r_gp, collision, inv_sinh)
#endif
for (i = 0; i < num_ir_gp; i++) {
inv_sinh = (double*)malloc(sizeof(double) * num_band);
for (j = 0; j < num_rot; j++) {
for (i = 0; i < num_ir_gp; i++)
{
inv_sinh = (double *)malloc(sizeof(double) * num_band);
for (j = 0; j < num_rot; j++)
{
r_gp = rot_grid_points[i * num_rot + j];
ti = gp2tp_map[triplets_map[r_gp]];
get_inv_sinh(inv_sinh,
@ -201,25 +203,30 @@ static void get_collision_matrix(double *collision_matrix,
num_band,
cutoff_frequency);
for (k = 0; k < num_band0; k++) {
for (l = 0; l < num_band; l++) {
for (k = 0; k < num_band0; k++)
{
for (l = 0; l < num_band; l++)
{
collision = 0;
for (m = 0; m < num_band; m++) {
for (m = 0; m < num_band; m++)
{
collision +=
fc3_normal_squared[ti * num_band0 * num_band * num_band +
k * num_band * num_band +
l * num_band + m] *
g[ti * num_band0 * num_band * num_band +
k * num_band * num_band +
l * num_band + m] *
inv_sinh[m] * unit_conversion_factor;
fc3_normal_squared[ti * num_band0 * num_band * num_band +
k * num_band * num_band +
l * num_band + m] *
g[ti * num_band0 * num_band * num_band +
k * num_band * num_band +
l * num_band + m] *
inv_sinh[m] * unit_conversion_factor;
}
for (m = 0; m < 3; m++) {
for (n = 0; n < 3; n++) {
for (m = 0; m < 3; m++)
{
for (n = 0; n < 3; n++)
{
collision_matrix[k * 3 * num_ir_gp * num_band * 3 +
m * num_ir_gp * num_band * 3 +
i * num_band * 3 + l * 3 + n] +=
collision * rotations_cartesian[j * 9 + m * 3 + n];
collision * rotations_cartesian[j * 9 + m * 3 + n];
}
}
}
@ -258,8 +265,9 @@ get_reducible_collision_matrix(double *collision_matrix,
#ifdef PHPYOPENMP
#pragma omp parallel for private(j, k, l, ti, collision, inv_sinh)
#endif
for (i = 0; i < num_gp; i++) {
inv_sinh = (double*)malloc(sizeof(double) * num_band);
for (i = 0; i < num_gp; i++)
{
inv_sinh = (double *)malloc(sizeof(double) * num_band);
ti = gp2tp_map[triplets_map[i]];
get_inv_sinh(inv_sinh,
i,
@ -271,18 +279,21 @@ get_reducible_collision_matrix(double *collision_matrix,
num_band,
cutoff_frequency);
for (j = 0; j < num_band0; j++) {
for (k = 0; k < num_band; k++) {
for (j = 0; j < num_band0; j++)
{
for (k = 0; k < num_band; k++)
{
collision = 0;
for (l = 0; l < num_band; l++) {
for (l = 0; l < num_band; l++)
{
collision +=
fc3_normal_squared[ti * num_band0 * num_band * num_band +
j * num_band * num_band +
k * num_band + l] *
g[ti * num_band0 * num_band * num_band +
j * num_band * num_band +
k * num_band + l] *
inv_sinh[l] * unit_conversion_factor;
fc3_normal_squared[ti * num_band0 * num_band * num_band +
j * num_band * num_band +
k * num_band + l] *
g[ti * num_band0 * num_band * num_band +
j * num_band * num_band +
k * num_band + l] *
inv_sinh[l] * unit_conversion_factor;
}
collision_matrix[j * num_gp * num_band + i * num_band + k] += collision;
}
@ -312,17 +323,24 @@ static void get_inv_sinh(double *inv_sinh,
/* This assumes the algorithm of get_ir_triplets_at_q_perm_q1q2, */
/* where defined triplets_map[gp] == triplets_map[map_q[gp]]. */
/* If triplets_map[map_q[gp]] != map_q[gp], q1 and q2 are permuted. */
if (triplets_map[gp] == map_q[gp]) {
if (triplets_map[gp] == map_q[gp])
{
gp2 = triplet[2];
} else {
}
else
{
gp2 = triplet[1];
}
for (i = 0; i < num_band; i++) {
for (i = 0; i < num_band; i++)
{
f = frequencies[gp2 * num_band + i];
if (f > cutoff_frequency) {
if (f > cutoff_frequency)
{
inv_sinh[i] = phonoc_inv_sinh_occupation(f, temperature);
} else {
}
else
{
inv_sinh[i] = 0;
}
}
@ -337,13 +355,17 @@ static long *create_gp2tp_map(const long *triplets_map,
long i, num_ir;
long *gp2tp_map;
gp2tp_map = (long*)malloc(sizeof(long) * num_gp);
gp2tp_map = (long *)malloc(sizeof(long) * num_gp);
num_ir = 0;
for (i = 0; i < num_gp; i++) {
if (triplets_map[i] == i) {
for (i = 0; i < num_gp; i++)
{
if (triplets_map[i] == i)
{
gp2tp_map[i] = num_ir;
num_ir++;
} else { /* This should not be used. */
}
else
{ /* This should not be used. */
gp2tp_map[i] = -1;
}
}

267
c/dynmat.c
View File

@ -65,7 +65,7 @@ static void get_dm(double dm_real[3][3],
const long k);
static double get_dielectric_part(const double q_cart[3],
const double dielectric[3][3]);
static void get_KK(double *dd_part, /* [natom, 3, natom, 3, (real,imag)] */
static void get_KK(double *dd_part, /* [natom, 3, natom, 3, (real,imag)] */
const double (*G_list)[3], /* [num_G, 3] */
const long num_G,
const long num_patom,
@ -96,11 +96,13 @@ long dym_get_dynamical_matrix_at_q(double *dynamical_matrix,
{
long i, j, ij;
if (with_openmp) {
if (with_openmp)
{
#ifdef PHPYOPENMP
#pragma omp parallel for
#endif
for (ij = 0; ij < num_patom * num_patom ; ij++) {
for (ij = 0; ij < num_patom * num_patom; ij++)
{
get_dynmat_ij(dynamical_matrix,
num_patom,
num_satom,
@ -115,9 +117,13 @@ long dym_get_dynamical_matrix_at_q(double *dynamical_matrix,
ij / num_patom, /* i */
ij % num_patom); /* j */
}
} else {
for (i = 0; i < num_patom; i++) {
for (j = 0; j < num_patom; j++) {
}
else
{
for (i = 0; i < num_patom; i++)
{
for (j = 0; j < num_patom; j++)
{
get_dynmat_ij(dynamical_matrix,
num_patom,
num_satom,
@ -140,8 +146,8 @@ long dym_get_dynamical_matrix_at_q(double *dynamical_matrix,
return 0;
}
void dym_get_recip_dipole_dipole(double *dd, /* [natom, 3, natom, 3, (real,imag)] */
const double *dd_q0, /* [natom, 3, 3, (real,imag)] */
void dym_get_recip_dipole_dipole(double *dd, /* [natom, 3, natom, 3, (real,imag)] */
const double *dd_q0, /* [natom, 3, 3, (real,imag)] */
const double (*G_list)[3], /* [num_G, 3] */
const long num_G,
const long num_patom,
@ -150,7 +156,7 @@ void dym_get_recip_dipole_dipole(double *dd, /* [natom, 3, natom, 3, (real,imag)
const double (*born)[3][3],
const double dielectric[3][3],
const double (*pos)[3], /* [num_patom, 3] */
const double factor, /* 4pi/V*unit-conv */
const double factor, /* 4pi/V*unit-conv */
const double lambda,
const double tolerance)
{
@ -158,9 +164,10 @@ void dym_get_recip_dipole_dipole(double *dd, /* [natom, 3, natom, 3, (real,imag)
double *dd_tmp;
dd_tmp = NULL;
dd_tmp = (double*) malloc(sizeof(double) * num_patom * num_patom * 18);
dd_tmp = (double *)malloc(sizeof(double) * num_patom * num_patom * 18);
for (i = 0; i < num_patom * num_patom * 18; i++) {
for (i = 0; i < num_patom * num_patom * 18; i++)
{
dd[i] = 0;
dd_tmp[i] = 0;
}
@ -178,9 +185,12 @@ void dym_get_recip_dipole_dipole(double *dd, /* [natom, 3, natom, 3, (real,imag)
multiply_borns(dd, dd_tmp, num_patom, born);
for (i = 0; i < num_patom; i++) {
for (k = 0; k < 3; k++) { /* alpha */
for (l = 0; l < 3; l++) { /* beta */
for (i = 0; i < num_patom; i++)
{
for (k = 0; k < 3; k++)
{ /* alpha */
for (l = 0; l < 3; l++)
{ /* beta */
adrs = i * num_patom * 9 + k * num_patom * 3 + i * 3 + l;
adrs_sum = i * 9 + k * 3 + l;
dd[adrs * 2] -= dd_q0[adrs_sum * 2];
@ -189,7 +199,8 @@ void dym_get_recip_dipole_dipole(double *dd, /* [natom, 3, natom, 3, (real,imag)
}
}
for (i = 0; i < num_patom * num_patom * 18; i++) {
for (i = 0; i < num_patom * num_patom * 18; i++)
{
dd[i] *= factor;
}
@ -200,7 +211,7 @@ void dym_get_recip_dipole_dipole(double *dd, /* [natom, 3, natom, 3, (real,imag)
dd_tmp = NULL;
}
void dym_get_recip_dipole_dipole_q0(double *dd_q0, /* [natom, 3, 3, (real,imag)] */
void dym_get_recip_dipole_dipole_q0(double *dd_q0, /* [natom, 3, 3, (real,imag)] */
const double (*G_list)[3], /* [num_G, 3] */
const long num_G,
const long num_patom,
@ -215,11 +226,12 @@ void dym_get_recip_dipole_dipole_q0(double *dd_q0, /* [natom, 3, 3, (real,imag)]
double *dd_tmp1, *dd_tmp2;
dd_tmp1 = NULL;
dd_tmp1 = (double*) malloc(sizeof(double) * num_patom * num_patom * 18);
dd_tmp1 = (double *)malloc(sizeof(double) * num_patom * num_patom * 18);
dd_tmp2 = NULL;
dd_tmp2 = (double*) malloc(sizeof(double) * num_patom * num_patom * 18);
dd_tmp2 = (double *)malloc(sizeof(double) * num_patom * num_patom * 18);
for (i = 0; i < num_patom * num_patom * 18; i++) {
for (i = 0; i < num_patom * num_patom * 18; i++)
{
dd_tmp1[i] = 0;
dd_tmp2[i] = 0;
}
@ -241,16 +253,21 @@ void dym_get_recip_dipole_dipole_q0(double *dd_q0, /* [natom, 3, 3, (real,imag)]
multiply_borns(dd_tmp2, dd_tmp1, num_patom, born);
for (i = 0; i < num_patom * 18; i++) {
for (i = 0; i < num_patom * 18; i++)
{
dd_q0[i] = 0;
}
for (i = 0; i < num_patom; i++) {
for (k = 0; k < 3; k++) { /* alpha */
for (l = 0; l < 3; l++) { /* beta */
for (i = 0; i < num_patom; i++)
{
for (k = 0; k < 3; k++)
{ /* alpha */
for (l = 0; l < 3; l++)
{ /* beta */
adrs = i * 9 + k * 3 + l;
for (j = 0; j < num_patom; j++) {
adrs_tmp = i * num_patom * 9 + k * num_patom * 3 + j * 3 + l ;
for (j = 0; j < num_patom; j++)
{
adrs_tmp = i * num_patom * 9 + k * num_patom * 3 + j * 3 + l;
dd_q0[adrs * 2] += dd_tmp2[adrs_tmp * 2];
dd_q0[adrs * 2 + 1] += dd_tmp2[adrs_tmp * 2 + 1];
}
@ -272,9 +289,12 @@ void dym_get_recip_dipole_dipole_q0(double *dd_q0, /* [natom, 3, 3, (real,imag)]
/* } */
/* } */
for (i = 0; i < num_patom; i++) {
for (k = 0; k < 3; k++) { /* alpha */
for (l = 0; l < 3; l++) { /* beta */
for (i = 0; i < num_patom; i++)
{
for (k = 0; k < 3; k++)
{ /* alpha */
for (l = 0; l < 3; l++)
{ /* beta */
adrs = i * 9 + k * 3 + l;
adrsT = i * 9 + l * 3 + k;
dd_q0[adrs * 2] += dd_q0[adrsT * 2];
@ -300,29 +320,38 @@ void dym_get_charge_sum(double (*charge_sum)[3][3],
const double (*born)[3][3])
{
long i, j, k, a, b;
double (*q_born)[3];
double(*q_born)[3];
q_born = (double (*)[3]) malloc(sizeof(double[3]) * num_patom);
for (i = 0; i < num_patom; i++) {
for (j = 0; j < 3; j++) {
q_born = (double(*)[3])malloc(sizeof(double[3]) * num_patom);
for (i = 0; i < num_patom; i++)
{
for (j = 0; j < 3; j++)
{
q_born[i][j] = 0;
}
}
for (i = 0; i < num_patom; i++) {
for (j = 0; j < 3; j++) {
for (k = 0; k < 3; k++) {
for (i = 0; i < num_patom; i++)
{
for (j = 0; j < 3; j++)
{
for (k = 0; k < 3; k++)
{
q_born[i][j] += q_cart[k] * born[i][k][j];
}
}
}
for (i = 0; i < num_patom; i++) {
for (j = 0; j < num_patom; j++) {
for (a = 0; a < 3; a++) {
for (b = 0; b < 3; b++) {
for (i = 0; i < num_patom; i++)
{
for (j = 0; j < num_patom; j++)
{
for (a = 0; a < 3; a++)
{
for (b = 0; b < 3; b++)
{
charge_sum[i * num_patom + j][a][b] =
q_born[i][a] * q_born[j][b] * factor;
q_born[i][a] * q_born[j][b] * factor;
}
}
}
@ -352,35 +381,43 @@ void dym_transform_dynmat_to_fc(double *fc,
double coef, phase, cos_phase, sin_phase;
N = num_satom / num_patom;
for (i = 0; i < num_patom * num_satom * 9; i++) {
for (i = 0; i < num_patom * num_satom * 9; i++)
{
fc[i] = 0;
}
for (i = 0; i < num_patom; i++) {
for (j = 0; j < num_satom; j++) {
for (i = 0; i < num_patom; i++)
{
for (j = 0; j < num_satom; j++)
{
i_pair = j * num_patom + i;
m_pair = multi[i_pair][0];
svecs_adrs = multi[i_pair][1];
coef = sqrt(masses[i] * masses[s2pp_map[j]]) / N;
for (k = 0; k < N; k++) {
for (k = 0; k < N; k++)
{
cos_phase = 0;
sin_phase = 0;
for (l = 0; l < m_pair; l++) {
for (l = 0; l < m_pair; l++)
{
phase = 0;
for (m = 0; m < 3; m++) {
for (m = 0; m < 3; m++)
{
phase -= comm_points[k][m] * svecs[svecs_adrs + l][m];
}
cos_phase += cos(phase * 2 * PI);
sin_phase += sin(phase * 2 * PI);
}
cos_phase /= m_pair;
sin_phase /= m_pair;
for (l = 0; l < 3; l++) {
for (m = 0; m < 3; m++) {
cos_phase /= m_pair;
sin_phase /= m_pair;
for (l = 0; l < 3; l++)
{
for (m = 0; m < 3; m++)
{
adrs = k * num_patom * num_patom * 18 + i * num_patom * 18 +
l * num_patom * 6 + s2pp_map[j] * 6 + m * 2;
l * num_patom * 6 + s2pp_map[j] * 6 + m * 2;
fc[fc_index_map[i] * num_satom * 9 + j * 9 + l * 3 + m] +=
(dm[adrs] * cos_phase - dm[adrs + 1] * sin_phase) * coef;
(dm[adrs] * cos_phase - dm[adrs + 1] * sin_phase) * coef;
}
}
}
@ -388,7 +425,6 @@ void dym_transform_dynmat_to_fc(double *fc,
}
}
static void get_dynmat_ij(double *dynamical_matrix,
const long num_patom,
const long num_satom,
@ -409,15 +445,19 @@ static void get_dynmat_ij(double *dynamical_matrix,
mass_sqrt = sqrt(mass[i] * mass[j]);
for (k = 0; k < 3; k++) {
for (l = 0; l < 3; l++) {
for (k = 0; k < 3; k++)
{
for (l = 0; l < 3; l++)
{
dm_real[k][l] = 0;
dm_imag[k][l] = 0;
}
}
for (k = 0; k < num_satom; k++) { /* Lattice points of right index of fc */
if (s2p_map[k] != p2s_map[j]) {
for (k = 0; k < num_satom; k++)
{ /* Lattice points of right index of fc */
if (s2p_map[k] != p2s_map[j])
{
continue;
}
get_dm(dm_real,
@ -435,8 +475,10 @@ static void get_dynmat_ij(double *dynamical_matrix,
k);
}
for (k = 0; k < 3; k++) {
for (l = 0; l < 3; l++) {
for (k = 0; k < 3; k++)
{
for (l = 0; l < 3; l++)
{
adrs = (i * 3 + k) * num_patom * 3 + j * 3 + l;
dynamical_matrix[adrs * 2] = dm_real[k][l] / mass_sqrt;
dynamical_matrix[adrs * 2 + 1] = dm_imag[k][l] / mass_sqrt;
@ -468,21 +510,28 @@ static void get_dm(double dm_real[3][3],
m_pair = multi[i_pair][0];
adrs = multi[i_pair][1];
for (l = 0; l < m_pair; l++) {
for (l = 0; l < m_pair; l++)
{
phase = 0;
for (m = 0; m < 3; m++) {
for (m = 0; m < 3; m++)
{
phase += q[m] * svecs[adrs + l][m];
}
cos_phase += cos(phase * 2 * PI) / m_pair;
sin_phase += sin(phase * 2 * PI) / m_pair;
}
for (l = 0; l < 3; l++) {
for (m = 0; m < 3; m++) {
if (charge_sum) {
for (l = 0; l < 3; l++)
{
for (m = 0; m < 3; m++)
{
if (charge_sum)
{
fc_elem = (fc[p2s_map[i] * num_satom * 9 + k * 9 + l * 3 + m] +
charge_sum[i * num_patom + j][l][m]);
} else {
}
else
{
fc_elem = fc[p2s_map[i] * num_satom * 9 + k * 9 + l * 3 + m];
}
dm_real[l][m] += fc_elem * cos_phase;
@ -498,22 +547,25 @@ static double get_dielectric_part(const double q_cart[3],
double x[3];
double sum;
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
x[i] = 0;
for (j = 0; j < 3; j++) {
for (j = 0; j < 3; j++)
{
x[i] += dielectric[i][j] * q_cart[j];
}
}
sum = 0;
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
sum += q_cart[i] * x[i];
}
return sum;
}
static void get_KK(double *dd_part, /* [natom, 3, natom, 3, (real,imag)] */
static void get_KK(double *dd_part, /* [natom, 3, natom, 3, (real,imag)] */
const double (*G_list)[3], /* [num_G, 3] */
const long num_G,
const long num_patom,
@ -534,39 +586,54 @@ static void get_KK(double *dd_part, /* [natom, 3, natom, 3, (real,imag)] */
/* sum over K = G + q and over G (i.e. q=0) */
/* q_direction has values for summation over K at Gamma point. */
/* q_direction is NULL for summation over G */
for (g = 0; g < num_G; g++) {
for (g = 0; g < num_G; g++)
{
norm = 0;
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
q_K[i] = G_list[g][i] + q_cart[i];
norm += q_K[i] * q_K[i];
}
if (sqrt(norm) < tolerance) {
if (!q_direction_cart) {
if (sqrt(norm) < tolerance)
{
if (!q_direction_cart)
{
continue;
} else {
}
else
{
dielectric_part = get_dielectric_part(q_direction_cart, dielectric);
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
for (i = 0; i < 3; i++)
{
for (j = 0; j < 3; j++)
{
KK[i][j] =
q_direction_cart[i] * q_direction_cart[j] / dielectric_part;
q_direction_cart[i] * q_direction_cart[j] / dielectric_part;
}
}
}
} else {
}
else
{
dielectric_part = get_dielectric_part(q_K, dielectric);
exp_damp = exp(-dielectric_part / L2);
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
for (i = 0; i < 3; i++)
{
for (j = 0; j < 3; j++)
{
KK[i][j] = q_K[i] * q_K[j] / dielectric_part * exp_damp;
}
}
}
for (i = 0; i < num_patom; i++) {
for (j = 0; j < num_patom; j++) {
for (i = 0; i < num_patom; i++)
{
for (j = 0; j < num_patom; j++)
{
phase = 0;
for (k = 0; k < 3; k++) {
for (k = 0; k < 3; k++)
{
/* For D-type dynamical matrix */
/* phase += (pos[i][k] - pos[j][k]) * q_K[k]; */
/* For C-type dynamical matrix */
@ -575,8 +642,10 @@ static void get_KK(double *dd_part, /* [natom, 3, natom, 3, (real,imag)] */
phase *= 2 * PI;
cos_phase = cos(phase);
sin_phase = sin(phase);
for (k = 0; k < 3; k++) {
for (l = 0; l < 3; l++) {
for (k = 0; k < 3; k++)
{
for (l = 0; l < 3; l++)
{
adrs = i * num_patom * 9 + k * num_patom * 3 + j * 3 + l;
dd_part[adrs * 2] += KK[k][l] * cos_phase;
dd_part[adrs * 2 + 1] += KK[k][l] * sin_phase;
@ -591,8 +660,10 @@ static void make_Hermitian(double *mat, const long num_band)
{
long i, j, adrs, adrsT;
for (i = 0; i < num_band; i++) {
for (j = i; j < num_band; j++) {
for (i = 0; i < num_band; i++)
{
for (j = i; j < num_band; j++)
{
adrs = i * num_band + j * 1;
adrs *= 2;
adrsT = j * num_band + i * 1;
@ -601,7 +672,7 @@ static void make_Hermitian(double *mat, const long num_band)
mat[adrs] += mat[adrsT];
mat[adrs] /= 2;
/* imaginary part */
mat[adrs + 1] -= mat[adrsT+ 1];
mat[adrs + 1] -= mat[adrsT + 1];
mat[adrs + 1] /= 2;
/* store */
mat[adrsT] = mat[adrs];
@ -618,14 +689,20 @@ static void multiply_borns(double *dd,
long i, j, k, l, m, n, adrs, adrs_in;
double zz;
for (i = 0; i < num_patom; i++) {
for (j = 0; j < num_patom; j++) {
for (k = 0; k < 3; k++) { /* alpha */
for (l = 0; l < 3; l++) { /* beta */
for (i = 0; i < num_patom; i++)
{
for (j = 0; j < num_patom; j++)
{
for (k = 0; k < 3; k++)
{ /* alpha */
for (l = 0; l < 3; l++)
{ /* beta */
adrs = i * num_patom * 9 + k * num_patom * 3 + j * 3 + l;
for (m = 0; m < 3; m++) { /* alpha' */
for (n = 0; n < 3; n++) { /* beta' */
adrs_in = i * num_patom * 9 + m * num_patom * 3 + j * 3 + n ;
for (m = 0; m < 3; m++)
{ /* alpha' */
for (n = 0; n < 3; n++)
{ /* beta' */
adrs_in = i * num_patom * 9 + m * num_patom * 3 + j * 3 + n;
zz = born[i][m][k] * born[j][n][l];
dd[adrs * 2] += dd_in[adrs_in * 2] * zz;
dd[adrs * 2 + 1] += dd_in[adrs_in * 2 + 1] * zz;

8
c/dynmat.h
View File

@ -47,8 +47,8 @@ long dym_get_dynamical_matrix_at_q(double *dynamical_matrix,
const long *p2s_map,
const double (*charge_sum)[3][3],
const long with_openmp);
void dym_get_recip_dipole_dipole(double *dd, /* [natom, 3, natom, 3, (real,imag)] */
const double *dd_q0, /* [natom, 3, 3, (real,imag)] */
void dym_get_recip_dipole_dipole(double *dd, /* [natom, 3, natom, 3, (real,imag)] */
const double *dd_q0, /* [natom, 3, 3, (real,imag)] */
const double (*G_list)[3], /* [num_G, 3] */
const long num_G,
const long num_patom,
@ -57,10 +57,10 @@ void dym_get_recip_dipole_dipole(double *dd, /* [natom, 3, natom, 3, (real,imag)
const double (*born)[3][3],
const double dielectric[3][3],
const double (*pos)[3], /* [num_patom, 3] */
const double factor, /* 4pi/V*unit-conv */
const double factor, /* 4pi/V*unit-conv */
const double lambda,
const double tolerance);
void dym_get_recip_dipole_dipole_q0(double *dd_q0, /* [natom, 3, 3, (real,imag)] */
void dym_get_recip_dipole_dipole_q0(double *dd_q0, /* [natom, 3, 3, (real,imag)] */
const double (*G_list)[3], /* [num_G, 3] */
const long num_G,
const long num_patom,

343
c/fc3.c
View File

@ -65,14 +65,14 @@ static void set_permutation_symmetry_fc3_elem(double *fc3_elem,
const long b,
const long c,
const long num_atom);
static void set_permutation_symmetry_compact_fc3(double * fc3,
static void set_permutation_symmetry_compact_fc3(double *fc3,
const long p2s[],
const long s2pp[],
const long nsym_list[],
const long perms[],
const long n_satom,
const long n_patom);
static void transpose_compact_fc3_type01(double * fc3,
static void transpose_compact_fc3_type01(double *fc3,
const long p2s[],
const long s2pp[],
const long nsym_list[],
@ -80,7 +80,7 @@ static void transpose_compact_fc3_type01(double * fc3,
const long n_satom,
const long n_patom,
const long t_type);
static void transpose_compact_fc3_type2(double * fc3,
static void transpose_compact_fc3_type2(double *fc3,
const long p2s[],
const long s2pp[],
const long nsym_list[],
@ -97,12 +97,12 @@ void fc3_distribute_fc3(double *fc3,
{
long i, j, adrs_out, adrs_in;
for (i = 0; i < num_atom; i++) {
for (j = 0; j < num_atom; j++) {
for (i = 0; i < num_atom; i++)
{
for (j = 0; j < num_atom; j++)
{
adrs_out = (num_atom * num_atom * target + num_atom * i + j) * 27;
adrs_in = (num_atom * num_atom * source
+ num_atom * atom_mapping[i]
+ atom_mapping[j]) * 27;
adrs_in = (num_atom * num_atom * source + num_atom * atom_mapping[i] + atom_mapping[j]) * 27;
tensor3_rotation(fc3 + adrs_out, fc3 + adrs_in, rot_cart);
}
}
@ -118,11 +118,11 @@ void fc3_rotate_delta_fc2(double (*fc3)[3][3][3],
const long num_disp)
{
long i_atoms, i, j, k, l, m, n;
double (*rot_delta_fc2s)[3][3];
double(*rot_delta_fc2s)[3][3];
rot_delta_fc2s = (double(*)[3][3]) malloc(sizeof(double[3][3])
* num_site_sym * num_disp);
for (i_atoms = 0; i_atoms < num_atom * num_atom; i_atoms++) {
rot_delta_fc2s = (double(*)[3][3])malloc(sizeof(double[3][3]) * num_site_sym * num_disp);
for (i_atoms = 0; i_atoms < num_atom * num_atom; i_atoms++)
{
i = i_atoms / num_atom;
j = i_atoms % num_atom;
rotate_delta_fc2s(rot_delta_fc2s,
@ -134,14 +134,17 @@ void fc3_rotate_delta_fc2(double (*fc3)[3][3][3],
num_atom,
num_site_sym,
num_disp);
for (k = 0; k < 3; k++) {
for (l = 0; l < 3; l++) {
for (m = 0; m < 3; m++) {
for (k = 0; k < 3; k++)
{
for (l = 0; l < 3; l++)
{
for (m = 0; m < 3; m++)
{
fc3[i_atoms][k][l][m] = 0;
for (n = 0; n < num_site_sym * num_disp; n++) {
for (n = 0; n < num_site_sym * num_disp; n++)
{
fc3[i_atoms][k][l][m] +=
inv_U[k * num_site_sym * num_disp + n]
* rot_delta_fc2s[n][l][m];
inv_U[k * num_site_sym * num_disp + n] * rot_delta_fc2s[n][l][m];
}
}
}
@ -160,9 +163,12 @@ void fc3_set_permutation_symmetry_fc3(double *fc3, const long num_atom)
#ifdef PHPYOPENMP
#pragma omp parallel for private(j, k, fc3_elem)
#endif
for (i = 0; i < num_atom; i++) {
for (j = i; j < num_atom; j++) {
for (k = j; k < num_atom; k++) {
for (i = 0; i < num_atom; i++)
{
for (j = i; j < num_atom; j++)
{
for (k = j; k < num_atom; k++)
{
set_permutation_symmetry_fc3_elem(fc3_elem, fc3, i, j, k, num_atom);
copy_permutation_symmetry_fc3_elem(fc3, fc3_elem,
i, j, k, num_atom);
@ -171,7 +177,7 @@ void fc3_set_permutation_symmetry_fc3(double *fc3, const long num_atom)
}
}
void fc3_set_permutation_symmetry_compact_fc3(double * fc3,
void fc3_set_permutation_symmetry_compact_fc3(double *fc3,
const long p2s[],
const long s2pp[],
const long nsym_list[],
@ -188,7 +194,7 @@ void fc3_set_permutation_symmetry_compact_fc3(double * fc3,
n_patom);
}
void fc3_transpose_compact_fc3(double * fc3,
void fc3_transpose_compact_fc3(double *fc3,
const long p2s[],
const long s2pp[],
const long nsym_list[],
@ -201,7 +207,8 @@ void fc3_transpose_compact_fc3(double * fc3,
/* t_type=0: dim[0] <-> dim[1] */
/* t_type=1: dim[0] <-> dim[2] */
/* t_type=2: dim[1] <-> dim[2] */
if (t_type == 0 || t_type == 1) {
if (t_type == 0 || t_type == 1)
{
transpose_compact_fc3_type01(fc3,
p2s,
s2pp,
@ -210,8 +217,11 @@ void fc3_transpose_compact_fc3(double * fc3,
n_satom,
n_patom,
t_type);
} else {
if (t_type == 2) {
}
else
{
if (t_type == 2)
{
transpose_compact_fc3_type2(fc3,
p2s,
s2pp,
@ -235,12 +245,12 @@ static void rotate_delta_fc2s(double (*rot_delta_fc2s)[3][3],
{
long i, j;
for (i = 0; i < num_disp; i++) {
for (j = 0; j < num_site_sym; j++) {
for (i = 0; i < num_disp; i++)
{
for (j = 0; j < num_site_sym; j++)
{
tensor2_rotation(rot_delta_fc2s[i * num_site_sym + j],
delta_fc2s[i * num_atom * num_atom
+ rot_map_sym[j * num_atom + i_atom] * num_atom
+ rot_map_sym[j * num_atom + j_atom]],
delta_fc2s[i * num_atom * num_atom + rot_map_sym[j * num_atom + i_atom] * num_atom + rot_map_sym[j * num_atom + j_atom]],
site_sym_cart[j]);
}
}
@ -252,16 +262,22 @@ static void tensor2_rotation(double rot_tensor[3][3],
{
long i, j, k, l;
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
for (i = 0; i < 3; i++)
{
for (j = 0; j < 3; j++)
{
rot_tensor[i][j] = 0;
}
}
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
for (k = 0; k < 3; k++) {
for (l = 0; l < 3; l++) {
for (i = 0; i < 3; i++)
{
for (j = 0; j < 3; j++)
{
for (k = 0; k < 3; k++)
{
for (l = 0; l < 3; l++)
{
rot_tensor[i][j] += r[i][k] * r[j][l] * tensor[k][l];
}
}
@ -275,7 +291,8 @@ static void tensor3_rotation(double *rot_tensor,
{
long l;
for (l = 0; l < 27; l++) {
for (l = 0; l < 27; l++)
{
rot_tensor[l] = tensor3_rotation_elem(tensor, rot_cartesian, l);
}
}
@ -292,11 +309,14 @@ static double tensor3_rotation_elem(const double *tensor,
n = pos % 3;
sum = 0.0;
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
for (k = 0; k < 3; k++) {
for (i = 0; i < 3; i++)
{
for (j = 0; j < 3; j++)
{
for (k = 0; k < 3; k++)
{
sum += r[l * 3 + i] * r[m * 3 + j] * r[n * 3 + k] *
tensor[i * 9 + j * 3 + k];
tensor[i * 9 + j * 3 + k];
}
}
}
@ -312,33 +332,36 @@ static void copy_permutation_symmetry_fc3_elem(double *fc3,
{
long i, j, k;
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
for (k = 0; k < 3; k++) {
for (i = 0; i < 3; i++)
{
for (j = 0; j < 3; j++)
{
for (k = 0; k < 3; k++)
{
fc3[a * num_atom * num_atom * 27 +
b * num_atom * 27 +
c * 27 + i * 9 + j * 3 + k] =
fc3_elem[i * 9 + j * 3 + k];
fc3_elem[i * 9 + j * 3 + k];
fc3[a * num_atom * num_atom * 27 +
c * num_atom * 27 +
b * 27 + i * 9 + k * 3 + j] =
fc3_elem[i * 9 + j * 3 + k];
fc3_elem[i * 9 + j * 3 + k];
fc3[b * num_atom * num_atom * 27 +
a * num_atom * 27 +
c * 27 + j * 9 + i * 3 + k] =
fc3_elem[i * 9 + j * 3 + k];
fc3_elem[i * 9 + j * 3 + k];
fc3[b * num_atom * num_atom * 27 +
c * num_atom * 27 +
a * 27 + j * 9 + k * 3 + i] =
fc3_elem[i * 9 + j * 3 + k];
fc3_elem[i * 9 + j * 3 + k];
fc3[c * num_atom * num_atom * 27 +
a * num_atom * 27 +
b * 27 + k * 9 + i * 3 + j] =
fc3_elem[i * 9 + j * 3 + k];
fc3_elem[i * 9 + j * 3 + k];
fc3[c * num_atom * num_atom * 27 +
b * num_atom * 27 +
a * 27 + k * 9 + j * 3 + i] =
fc3_elem[i * 9 + j * 3 + k];
fc3_elem[i * 9 + j * 3 + k];
}
}
}
@ -353,34 +376,38 @@ static void set_permutation_symmetry_fc3_elem(double *fc3_elem,
{
long i, j, k;
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
for (k = 0; k < 3; k++) {
for (i = 0; i < 3; i++)
{
for (j = 0; j < 3; j++)
{
for (k = 0; k < 3; k++)
{
fc3_elem[i * 9 + j * 3 + k] =
(fc3[a * num_atom * num_atom * 27 +
b * num_atom * 27 +
c * 27 + i * 9 + j * 3 + k] +
fc3[a * num_atom * num_atom * 27 +
c * num_atom * 27 +
b * 27 + i * 9 + k * 3 + j] +
fc3[b * num_atom * num_atom * 27 +
a * num_atom * 27 +
c * 27 + j * 9 + i * 3 + k] +
fc3[b * num_atom * num_atom * 27 +
c * num_atom * 27 +
a * 27 + j * 9 + k * 3 + i] +
fc3[c * num_atom * num_atom * 27 +
a * num_atom * 27 +
b * 27 + k * 9 + i * 3 + j] +
fc3[c * num_atom * num_atom * 27 +
b * num_atom * 27 +
a * 27 + k * 9 + j * 3 + i]) / 6;
(fc3[a * num_atom * num_atom * 27 +
b * num_atom * 27 +
c * 27 + i * 9 + j * 3 + k] +
fc3[a * num_atom * num_atom * 27 +
c * num_atom * 27 +
b * 27 + i * 9 + k * 3 + j] +
fc3[b * num_atom * num_atom * 27 +
a * num_atom * 27 +
c * 27 + j * 9 + i * 3 + k] +
fc3[b * num_atom * num_atom * 27 +
c * num_atom * 27 +
a * 27 + j * 9 + k * 3 + i] +
fc3[c * num_atom * num_atom * 27 +
a * num_atom * 27 +
b * 27 + k * 9 + i * 3 + j] +
fc3[c * num_atom * num_atom * 27 +
b * num_atom * 27 +
a * 27 + k * 9 + j * 3 + i]) /
6;
}
}
}
}
static void set_permutation_symmetry_compact_fc3(double * fc3,
static void set_permutation_symmetry_compact_fc3(double *fc3,
const long p2s[],
const long s2pp[],
const long nsym_list[],
@ -400,17 +427,21 @@ static void set_permutation_symmetry_compact_fc3(double * fc3,
char *done;
done = NULL;
done = (char*)malloc(sizeof(char) * n_patom * n_satom * n_satom);
for (i = 0; i < n_patom * n_satom * n_satom; i++) {
done = (char *)malloc(sizeof(char) * n_patom * n_satom * n_satom);
for (i = 0; i < n_patom * n_satom * n_satom; i++)
{
done[i] = 0;
}
for (i_p = 0; i_p < n_patom; i_p++) {
for (i_p = 0; i_p < n_patom; i_p++)
{
i = p2s[i_p];
for (j = 0; j < n_satom; j++) {
for (j = 0; j < n_satom; j++)
{
j_p = s2pp[j];
i_trans_j = perms[nsym_list[j] * n_satom + i];
for (k = 0; k < n_satom; k++) {
for (k = 0; k < n_satom; k++)
{
k_p = s2pp[k];
k_trans_j = perms[nsym_list[j] * n_satom + k];
i_trans_k = perms[nsym_list[k] * n_satom + i];
@ -425,24 +456,31 @@ static void set_permutation_symmetry_compact_fc3(double * fc3,
adrs[5] = k_p * n_satom * n_satom + j_trans_k * n_satom + i_trans_k;
done_any = 0;
for (l = 0; l < 6; l++) {
if (done[adrs[l]]) {
for (l = 0; l < 6; l++)
{
if (done[adrs[l]])
{
done_any = 1;
break;
}
}
if (done_any) {
if (done_any)
{
continue;
}
for (l = 0; l < 6; l++) {
for (l = 0; l < 6; l++)
{
done[adrs[l]] = 1;
adrs[l] *= 27;
}
for (l = 0; l < 3; l++) {
for (m = 0; m < 3; m++) {
for (n = 0; n < 3; n++) {
for (l = 0; l < 3; l++)
{
for (m = 0; m < 3; m++)
{
for (n = 0; n < 3; n++)
{
fc3_elem[l][m][n] = fc3[adrs[0] + l * 9 + m * 3 + n];
fc3_elem[l][m][n] += fc3[adrs[1] + l * 9 + n * 3 + m];
fc3_elem[l][m][n] += fc3[adrs[2] + m * 9 + l * 3 + n];
@ -453,9 +491,12 @@ static void set_permutation_symmetry_compact_fc3(double * fc3,
}
}
}
for (l = 0; l < 3; l++) {
for (m = 0; m < 3; m++) {
for (n = 0; n < 3; n++) {
for (l = 0; l < 3; l++)
{
for (m = 0; m < 3; m++)
{
for (n = 0; n < 3; n++)
{
fc3[adrs[0] + l * 9 + m * 3 + n] = fc3_elem[l][m][n];
fc3[adrs[1] + l * 9 + n * 3 + m] = fc3_elem[l][m][n];
fc3[adrs[2] + m * 9 + l * 3 + n] = fc3_elem[l][m][n];
@ -471,11 +512,9 @@ static void set_permutation_symmetry_compact_fc3(double * fc3,
free(done);
done = NULL;
}
static void transpose_compact_fc3_type01(double * fc3,
static void transpose_compact_fc3_type01(double *fc3,
const long p2s[],
const long s2pp[],
const long nsym_list[],
@ -494,16 +533,20 @@ static void transpose_compact_fc3_type01(double * fc3,
char *done;
done = NULL;
done = (char*)malloc(sizeof(char) * n_satom * n_patom);
for (i = 0; i < n_satom * n_patom; i++) {
done = (char *)malloc(sizeof(char) * n_satom * n_patom);
for (i = 0; i < n_satom * n_patom; i++)
{
done[i] = 0;
}
for (i_p = 0; i_p < n_patom; i_p++) {
for (i_p = 0; i_p < n_patom; i_p++)
{
i = p2s[i_p];
for (j = 0; j < n_satom; j++) {
for (j = 0; j < n_satom; j++)
{
j_p = s2pp[j];
if (!done[i_p * n_satom + j]) {
if (!done[i_p * n_satom + j])
{
/* (j, i) -- nsym_list[j] --> (j', i') */
/* nsym_list[j] translates j to j' where j' is in */
/* primitive cell. The same translation sends i to i' */
@ -512,32 +555,44 @@ static void transpose_compact_fc3_type01(double * fc3,
i_trans = perms[nsym_list[j] * n_satom + i];
done[i_p * n_satom + j] = 1;
done[j_p * n_satom + i_trans] = 1;
for (k = 0; k < n_satom; k++) {
for (k = 0; k < n_satom; k++)
{
k_trans = perms[nsym_list[j] * n_satom + k];
switch (t_type) {
switch (t_type)
{
case 0:
adrs = (i_p * n_satom * n_satom + j * n_satom + k) * 27;
adrs_t = (j_p * n_satom * n_satom + i_trans * n_satom + k_trans) * 27;
for (l = 0; l < 3; l++) {
for (m = 0; m < 3; m++) {
for (n = 0; n < 3; n++) {
for (l = 0; l < 3; l++)
{
for (m = 0; m < 3; m++)
{
for (n = 0; n < 3; n++)
{
fc3_elem[l][m][n] = fc3[adrs + l * 9 + m * 3 + n];
}
}
}
if (adrs != adrs_t) {
for (l = 0; l < 3; l++) {
for (m = 0; m < 3; m++) {
for (n = 0; n < 3; n++) {
if (adrs != adrs_t)
{
for (l = 0; l < 3; l++)
{
for (m = 0; m < 3; m++)
{
for (n = 0; n < 3; n++)
{
fc3[adrs + l * 9 + m * 3 + n] = fc3[adrs_t + m * 9 + l * 3 + n];
}
}
}
}
for (l = 0; l < 3; l++) {
for (m = 0; m < 3; m++) {
for (n = 0; n < 3; n++) {
for (l = 0; l < 3; l++)
{
for (m = 0; m < 3; m++)
{
for (n = 0; n < 3; n++)
{
fc3[adrs_t + m * 9 + l * 3 + n] = fc3_elem[l][m][n];
}
}
@ -546,32 +601,41 @@ static void transpose_compact_fc3_type01(double * fc3,
case 1:
adrs = (i_p * n_satom * n_satom + k * n_satom + j) * 27;
adrs_t = (j_p * n_satom * n_satom + k_trans * n_satom + i_trans) * 27;
for (l = 0; l < 3; l++) {
for (m = 0; m < 3; m++) {
for (n = 0; n < 3; n++) {
for (l = 0; l < 3; l++)
{
for (m = 0; m < 3; m++)
{
for (n = 0; n < 3; n++)
{
fc3_elem[l][m][n] = fc3[adrs + l * 9 + m * 3 + n];
}
}
}
if (adrs != adrs_t) {
for (l = 0; l < 3; l++) {
for (m = 0; m < 3; m++) {
for (n = 0; n < 3; n++) {
if (adrs != adrs_t)
{
for (l = 0; l < 3; l++)
{
for (m = 0; m < 3; m++)
{
for (n = 0; n < 3; n++)
{
fc3[adrs + l * 9 + m * 3 + n] = fc3[adrs_t + n * 9 + m * 3 + l];
}
}
}
}
for (l = 0; l < 3; l++) {
for (m = 0; m < 3; m++) {
for (n = 0; n < 3; n++) {
for (l = 0; l < 3; l++)
{
for (m = 0; m < 3; m++)
{
for (n = 0; n < 3; n++)
{
fc3[adrs_t + n * 9 + m * 3 + l] = fc3_elem[l][m][n];
}
}
}
break;
} /* end switch */
}
}
}
@ -581,7 +645,7 @@ static void transpose_compact_fc3_type01(double * fc3,
done = NULL;
}
static void transpose_compact_fc3_type2(double * fc3,
static void transpose_compact_fc3_type2(double *fc3,
const long p2s[],
const long s2pp[],
const long nsym_list[],
@ -593,30 +657,43 @@ static void transpose_compact_fc3_type2(double * fc3,
long adrs, adrs_t;
double fc3_elem[3][3][3];
for (i_p = 0; i_p < n_patom; i_p++) {
for (j = 0; j < n_satom; j++) {
for (k = j; k < n_satom; k++) { /* k >= j */
for (i_p = 0; i_p < n_patom; i_p++)
{
for (j = 0; j < n_satom; j++)
{
for (k = j; k < n_satom; k++)
{ /* k >= j */
adrs = (i_p * n_satom * n_satom + j * n_satom + k) * 27;
adrs_t = (i_p * n_satom * n_satom + k * n_satom + j) * 27;
for (l = 0; l < 3; l++) {
for (m = 0; m < 3; m++) {
for (n = 0; n < 3; n++) {
for (l = 0; l < 3; l++)
{
for (m = 0; m < 3; m++)
{
for (n = 0; n < 3; n++)
{
fc3_elem[l][m][n] = fc3[adrs + l * 9 + m * 3 + n];
}
}
}
if (k != j) {
for (l = 0; l < 3; l++) {
for (m = 0; m < 3; m++) {
for (n = 0; n < 3; n++) {
if (k != j)
{
for (l = 0; l < 3; l++)
{
for (m = 0; m < 3; m++)
{
for (n = 0; n < 3; n++)
{
fc3[adrs + l * 9 + m * 3 + n] = fc3[adrs_t + l * 9 + n * 3 + m];
}
}
}
}
for (l = 0; l < 3; l++) {
for (m = 0; m < 3; m++) {
for (n = 0; n < 3; n++) {
for (l = 0; l < 3; l++)
{
for (m = 0; m < 3; m++)
{
for (n = 0; n < 3; n++)
{
fc3[adrs_t + l * 9 + n * 3 + m] = fc3_elem[l][m][n];
}
}

4
c/fc3.h
View File

@ -50,14 +50,14 @@ void fc3_rotate_delta_fc2(double (*fc3)[3][3][3],
const long num_site_sym,
const long num_disp);
void fc3_set_permutation_symmetry_fc3(double *fc3, const long num_atom);
void fc3_set_permutation_symmetry_compact_fc3(double * fc3,
void fc3_set_permutation_symmetry_compact_fc3(double *fc3,
const long p2s[],
const long s2pp[],
const long nsym_list[],
const long perms[],
const long n_satom,
const long n_patom);
void fc3_transpose_compact_fc3(double * fc3,
void fc3_transpose_compact_fc3(double *fc3,
const long p2s[],
const long s2pp[],
const long nsym_list[],

114
c/grgrid.c
View File

@ -68,7 +68,6 @@ static void get_ir_grid_map(long *ir_grid_map,
const long D_diag[3],
const long PS[3]);
long grg_get_snf3x3(long D_diag[3],
long P[3][3],
long Q[3][3],
@ -79,18 +78,22 @@ long grg_get_snf3x3(long D_diag[3],
succeeded = 0;
if (lagmat_get_determinant_l3(A) == 0) {
if (lagmat_get_determinant_l3(A) == 0)
{
goto err;
}
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
for (i = 0; i < 3; i++)
{
for (j = 0; j < 3; j++)
{
D[i][j] = A[i][j];
}
}
succeeded = snf3x3(D, P, Q);
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
D_diag[i] = D[i][i];
}
@ -122,16 +125,20 @@ long grg_transform_rotations(long (*transformed_rots)[3][3],
/* 2. Compute D(Q^-1)RQ */
/* 3. Compute D(Q^-1)RQ(D^-1) */
lagmat_cast_matrix_3l_to_3d(Q_double, Q);
for (i = 0; i < num_rot; i++) {
for (i = 0; i < num_rot; i++)
{
lagmat_get_similar_matrix_ld3(r, rotations[i], Q_double, 0);
for (j = 0; j < 3; j++) {
for (k = 0; k < 3; k++) {
for (j = 0; j < 3; j++)
{
for (k = 0; k < 3; k++)
{
r[j][k] *= D_diag[j];
r[j][k] /= D_diag[k];
}
}
lagmat_cast_matrix_3d_to_3l(transformed_rots[i], r);
if (!lagmat_check_identity_matrix_ld3(transformed_rots[i], r, IDENTITY_TOL)) {
if (!lagmat_check_identity_matrix_ld3(transformed_rots[i], r, IDENTITY_TOL))
{
return 0;
}
}
@ -218,7 +225,6 @@ void grg_get_grid_address_from_index(long address[3],
get_grid_address_from_index(address, grid_index, D_diag);
}
/* ---------------------------*/
/* Rotate grid point by index */
/* ---------------------------*/
@ -263,49 +269,59 @@ long grg_get_reciprocal_point_group(long rec_rotations[48][3][3],
{
long i, j, num_rot_ret, inv_exist;
const long inversion[3][3] = {
{-1, 0, 0 },
{ 0,-1, 0 },
{ 0, 0,-1 }
};
{-1, 0, 0},
{0, -1, 0},
{0, 0, -1}};
/* Collect unique rotations */
num_rot_ret = 0;
for (i = 0; i < num_rot; i++) {
for (j = 0; j < num_rot_ret; j++) {
if (lagmat_check_identity_matrix_l3(rotations[i], rec_rotations[j])) {
for (i = 0; i < num_rot; i++)
{
for (j = 0; j < num_rot_ret; j++)
{
if (lagmat_check_identity_matrix_l3(rotations[i], rec_rotations[j]))
{
goto escape;
}
}
if (num_rot_ret == 48) {
if (num_rot_ret == 48)
{
goto err;
}
lagmat_copy_matrix_l3(rec_rotations[num_rot_ret], rotations[i]);
num_rot_ret++;
escape:
;
escape:;
}
if (is_transpose) {
for (i = 0; i < num_rot_ret; i++) {
if (is_transpose)
{
for (i = 0; i < num_rot_ret; i++)
{
lagmat_transpose_matrix_l3(rec_rotations[i], rec_rotations[i]);
}
}
if (is_time_reversal) {
if (is_time_reversal)
{
inv_exist = 0;
for (i = 0; i < num_rot_ret; i++) {
if (lagmat_check_identity_matrix_l3(inversion, rec_rotations[i])) {
for (i = 0; i < num_rot_ret; i++)
{
if (lagmat_check_identity_matrix_l3(inversion, rec_rotations[i]))
{
inv_exist = 1;
break;
}
}
if (!inv_exist) {
if (num_rot_ret > 24) {
if (!inv_exist)
{
if (num_rot_ret > 24)
{
goto err;
}
for (i = 0; i < num_rot_ret; i++) {
for (i = 0; i < num_rot_ret; i++)
{
lagmat_multiply_matrix_l3(rec_rotations[num_rot_ret + i],
inversion, rec_rotations[i]);
}
@ -318,17 +334,16 @@ err:
return 0;
}
static void reduce_grid_address(long address[3], const long D_diag[3])
{
long i;
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
address[i] = lagmat_modulo_l(address[i], D_diag[i]);
}
}
static long get_double_grid_index(const long address_double[3],
const long D_diag[3],
const long PS[3])
@ -348,11 +363,9 @@ static long get_grid_index_from_address(const long address[3],
const long D_diag[3])
{
#ifndef GRID_ORDER_XYZ
return (address[2] * D_diag[0] * D_diag[1]
+ address[1] * D_diag[0] + address[0]);
return (address[2] * D_diag[0] * D_diag[1] + address[1] * D_diag[0] + address[0]);
#else
return (address[0] * D_diag[1] * D_diag[2]
+ address[1] * D_diag[2] + address[2]);
return (address[0] * D_diag[1] * D_diag[2] + address[1] * D_diag[2] + address[2]);
#endif
}
@ -362,11 +375,14 @@ static void get_all_grid_addresses(long grid_address[][3],
long i, j, k, grid_index;
long address[3];
for (i = 0; i < D_diag[0]; i++) {
for (i = 0; i < D_diag[0]; i++)
{
address[0] = i;
for (j = 0; j < D_diag[1]; j++) {
for (j = 0; j < D_diag[1]; j++)
{
address[1] = j;
for (k = 0; k < D_diag[2]; k++) {
for (k = 0; k < D_diag[2]; k++)
{
address[2] = k;
grid_index = get_grid_index_from_address(address, D_diag);
lagmat_copy_vector_l3(grid_address[grid_index], address);
@ -403,7 +419,8 @@ static void get_grid_address(long address[3],
{
long i;
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
address[i] = (address_double[i] - PS[i]) / 2;
}
}
@ -416,7 +433,8 @@ static void get_double_grid_address(long address_double[3],
{
long i;
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
address_double[i] = address[i] * 2 + PS[i];
}
}
@ -448,23 +466,27 @@ static void get_ir_grid_map(long *ir_grid_map,
num_gp = D_diag[0] * D_diag[1] * D_diag[2];
for (gp = 0; gp < num_gp; gp++) {
for (gp = 0; gp < num_gp; gp++)
{
ir_grid_map[gp] = num_gp;
}
/* Do not simply multithreaded this for-loop. */
/* This algorithm contains race condition in different gp's. */
for (gp = 0; gp < num_gp; gp++) {
for (i = 0; i < num_rot; i++) {
for (gp = 0; gp < num_gp; gp++)
{
for (i = 0; i < num_rot; i++)
{
r_gp = rotate_grid_index(gp, rotations[i], D_diag, PS);
if (r_gp < gp) {
if (r_gp < gp)
{
ir_grid_map[gp] = ir_grid_map[r_gp];
break;
}
}
if (ir_grid_map[gp] == num_gp) {
if (ir_grid_map[gp] == num_gp)
{
ir_grid_map[gp] = gp;
}
}
}

44
c/gridsys.c
View File

@ -50,7 +50,6 @@ void gridsys_get_all_grid_addresses(long (*gr_grid_addresses)[3],
grg_get_all_grid_addresses(gr_grid_addresses, D_diag);
}
void gridsys_get_double_grid_address(long address_double[3],
const long address[3],
const long PS[3])
@ -58,7 +57,6 @@ void gridsys_get_double_grid_address(long address_double[3],
grg_get_double_grid_address(address_double, address, PS);
}
void gridsys_get_grid_address_from_index(long address[3],
const long grid_index,
const long D_diag[3])
@ -66,7 +64,6 @@ void gridsys_get_grid_address_from_index(long address[3],
grg_get_grid_address_from_index(address, grid_index, D_diag);
}
long gridsys_get_double_grid_index(const long address_double[3],
const long D_diag[3],
const long PS[3])
@ -74,7 +71,6 @@ long gridsys_get_double_grid_index(const long address_double[3],
return grg_get_double_grid_index(address_double, D_diag, PS);
}
/* From single address to grid index */
long gridsys_get_grid_index_from_address(const long address[3],
const long D_diag[3])
@ -82,7 +78,6 @@ long gridsys_get_grid_index_from_address(const long address[3],
return grg_get_grid_index(address, D_diag);
}
long gridsys_rotate_grid_index(const long grid_index,
const long rotation[3][3],
const long D_diag[3],
@ -103,7 +98,6 @@ long gridsys_get_reciprocal_point_group(long rec_rotations[48][3][3],
1);
}
long gridsys_get_snf3x3(long D_diag[3],
long P[3][3],
long Q[3][3],
@ -112,7 +106,6 @@ long gridsys_get_snf3x3(long D_diag[3],
return grg_get_snf3x3(D_diag, P, Q, A);
}
/* Rotation matrices with respect to reciprocal basis vectors are
* transformed to those for GRGrid. This set of the rotations are
* used always in GRGrid handling. */
@ -131,7 +124,6 @@ long gridsys_transform_rotations(long (*transformed_rots)[3][3],
return succeeded;
}
double gridsys_get_thm_integration_weight(const double omega,
const double tetrahedra_omegas[24][4],
const char function)
@ -139,7 +131,6 @@ double gridsys_get_thm_integration_weight(const double omega,
return thm_get_integration_weight(omega, tetrahedra_omegas, function);
}
/* Get one dataset of relative grid address used for tetrahedron */
/* method. rec_lattice is used to choose the one. */
/* rec_lattice : microzone basis vectors in column vectors */
@ -149,7 +140,6 @@ void gridsys_get_thm_relative_grid_address(long relative_grid_addresses[24][4][3
thm_get_relative_grid_address(relative_grid_addresses, rec_lattice);
}
/* The rotations are those after proper transformation in GRGrid. */
void gridsys_get_ir_grid_map(long *ir_grid_map,
const long (*rotations)[3][3],
@ -160,7 +150,6 @@ void gridsys_get_ir_grid_map(long *ir_grid_map,
grg_get_ir_grid_map(ir_grid_map, rotations, num_rot, D_diag, PS);
}
/* Find shortest grid points from Gamma considering periodicity of */
/* reciprocal lattice. See the details in docstring of BZGrid. */
/* */
@ -196,7 +185,8 @@ long gridsys_get_bz_grid_addresses(long (*bz_grid_addresses)[3],
BZGrid *bzgrid;
long i, j, size;
if ((bzgrid = (BZGrid*) malloc(sizeof(BZGrid))) == NULL) {
if ((bzgrid = (BZGrid *)malloc(sizeof(BZGrid))) == NULL)
{
warning_print("Memory could not be allocated.");
return 0;
}
@ -205,18 +195,23 @@ long gridsys_get_bz_grid_addresses(long (*bz_grid_addresses)[3],
bzgrid->gp_map = bz_map;
bzgrid->bzg2grg = bzg2grg;
bzgrid->type = type;
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
bzgrid->D_diag[i] = D_diag[i];
bzgrid->PS[i] = PS[i];
for (j = 0; j < 3; j++) {
for (j = 0; j < 3; j++)
{
bzgrid->Q[i][j] = Q[i][j];
bzgrid->reclat[i][j] = rec_lattice[i][j];
}
}
if (bzg_get_bz_grid_addresses(bzgrid)) {
if (bzg_get_bz_grid_addresses(bzgrid))
{
size = bzgrid->size;
} else {
}
else
{
size = 0;
}
@ -226,7 +221,6 @@ long gridsys_get_bz_grid_addresses(long (*bz_grid_addresses)[3],
return size;
}
long gridsys_get_triplets_at_q(long *map_triplets,
long *map_q,
const long grid_point,
@ -246,7 +240,6 @@ long gridsys_get_triplets_at_q(long *map_triplets,
swappable);
}
long gridsys_get_BZ_triplets_at_q(long (*triplets)[3],
const long grid_point,
const long (*bz_grid_addresses)[3],
@ -260,7 +253,8 @@ long gridsys_get_BZ_triplets_at_q(long (*triplets)[3],
ConstBZGrid *bzgrid;
long i, j, num_ir;
if ((bzgrid = (ConstBZGrid*) malloc(sizeof(ConstBZGrid))) == NULL) {
if ((bzgrid = (ConstBZGrid *)malloc(sizeof(ConstBZGrid))) == NULL)
{
warning_print("Memory could not be allocated.");
return 0;
}
@ -268,10 +262,12 @@ long gridsys_get_BZ_triplets_at_q(long (*triplets)[3],
bzgrid->addresses = bz_grid_addresses;
bzgrid->gp_map = bz_map;
bzgrid->type = bz_grid_type;
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
bzgrid->D_diag[i] = D_diag[i];
bzgrid->PS[i] = 0;
for (j = 0; j < 3; j++) {
for (j = 0; j < 3; j++)
{
bzgrid->Q[i][j] = Q[i][j];
}
}
@ -312,7 +308,8 @@ long gridsys_get_integration_weight(double *iw,
ConstBZGrid *bzgrid;
long i;
if ((bzgrid = (ConstBZGrid*) malloc(sizeof(ConstBZGrid))) == NULL) {
if ((bzgrid = (ConstBZGrid *)malloc(sizeof(ConstBZGrid))) == NULL)
{
warning_print("Memory could not be allocated.");
return 0;
}
@ -320,7 +317,8 @@ long gridsys_get_integration_weight(double *iw,
bzgrid->addresses = bz_grid_addresses;
bzgrid->gp_map = bz_map;
bzgrid->type = bz_grid_type;
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
bzgrid->D_diag[i] = D_diag[i];
}

180
c/gridsys.h
View File

@ -36,104 +36,104 @@
#define __gridsys_H__
#ifdef __cplusplus
extern "C" {
extern "C"
{
#endif
void gridsys_get_all_grid_addresses(long (*gr_grid_addresses)[3],
const long D_diag[3]);
void gridsys_get_double_grid_address(long address_double[3],
const long address[3],
const long PS[3]);
void gridsys_get_grid_address_from_index(long address[3],
const long grid_index,
const long D_diag[3]);
long gridsys_get_double_grid_index(const long address_double[3],
void gridsys_get_all_grid_addresses(long (*gr_grid_addresses)[3],
const long D_diag[3]);
void gridsys_get_double_grid_address(long address_double[3],
const long address[3],
const long PS[3]);
void gridsys_get_grid_address_from_index(long address[3],
const long grid_index,
const long D_diag[3]);
long gridsys_get_double_grid_index(const long address_double[3],
const long D_diag[3],
const long PS[3]);
long gridsys_get_grid_index_from_address(const long address[3],
const long D_diag[3]);
long gridsys_rotate_grid_index(const long grid_index,
const long rotation[3][3],
const long D_diag[3],
const long PS[3]);
long gridsys_get_grid_index_from_address(const long address[3],
const long D_diag[3]);
long gridsys_rotate_grid_index(const long grid_index,
const long rotation[3][3],
const long D_diag[3],
const long PS[3]);
long gridsys_get_reciprocal_point_group(long rec_rotations[48][3][3],
const long (*rotations)[3][3],
const long num_rot,
const long is_time_reversal);
long gridsys_get_snf3x3(long D_diag[3],
long P[3][3],
long Q[3][3],
const long A[3][3]);
long gridsys_transform_rotations(long (*transformed_rots)[3][3],
long gridsys_get_reciprocal_point_group(long rec_rotations[48][3][3],
const long (*rotations)[3][3],
const long num_rot,
const long is_time_reversal);
long gridsys_get_snf3x3(long D_diag[3],
long P[3][3],
long Q[3][3],
const long A[3][3]);
long gridsys_transform_rotations(long (*transformed_rots)[3][3],
const long (*rotations)[3][3],
const long num_rot,
const long D_diag[3],
const long Q[3][3]);
double gridsys_get_thm_integration_weight(const double omega,
const double tetrahedra_omegas[24][4],
const char function);
void gridsys_get_thm_relative_grid_address(long relative_grid_addresses[24][4][3],
const double rec_lattice[3][3]);
void gridsys_get_ir_grid_map(long *ir_grid_map,
const long (*rotations)[3][3],
const long num_rot,
const long D_diag[3],
const long Q[3][3]);
double gridsys_get_thm_integration_weight(const double omega,
const double tetrahedra_omegas[24][4],
const char function);
void gridsys_get_thm_relative_grid_address(long relative_grid_addresses[24][4][3],
const double rec_lattice[3][3]);
void gridsys_get_ir_grid_map(long *ir_grid_map,
const long (*rotations)[3][3],
const long num_rot,
const long D_diag[3],
const long PS[3]);
long gridsys_get_bz_grid_addresses(long (*bz_grid_addresses)[3],
long *bz_map,
long *bzg2grg,
const long PS[3]);
long gridsys_get_bz_grid_addresses(long (*bz_grid_addresses)[3],
long *bz_map,
long *bzg2grg,
const long D_diag[3],
const long Q[3][3],
const long PS[3],
const double rec_lattice[3][3],
const long type);
long gridsys_get_triplets_at_q(long *map_triplets,
long *map_q,
const long grid_point,
const long D_diag[3],
const long Q[3][3],
const long PS[3],
const double rec_lattice[3][3],
const long type);
long gridsys_get_triplets_at_q(long *map_triplets,
long *map_q,
const long grid_point,
const long D_diag[3],
const long is_time_reversal,
const long num_rot,
const long (*rec_rotations)[3][3],
const long swappable);
long gridsys_get_BZ_triplets_at_q(long (*triplets)[3],
const long grid_point,
const long (*bz_grid_addresses)[3],
const long *bz_map,
const long *map_triplets,
const long num_map_triplets,
const long D_diag[3],
const long Q[3][3],
const long bz_grid_type);
long gridsys_get_integration_weight(double *iw,
char *iw_zero,
const double *frequency_points,
const long num_band0,
const long relative_grid_address[24][4][3],
const long D_diag[3],
const long (*triplets)[3],
const long num_triplets,
const long (*bz_grid_addresses)[3],
const long *bz_map,
const long bz_grid_type,
const double *frequencies1,
const long num_band1,
const double *frequencies2,
const long num_band2,
const long tp_type,
const long openmp_per_triplets,
const long openmp_per_bands);
void gridsys_get_integration_weight_with_sigma(double *iw,
char *iw_zero,
const double sigma,
const double sigma_cutoff,
const double *frequency_points,
const long num_band0,
const long (*triplets)[3],
const long num_triplets,
const double *frequencies,
const long num_band,
const long tp_type);
const long is_time_reversal,
const long num_rot,
const long (*rec_rotations)[3][3],
const long swappable);
long gridsys_get_BZ_triplets_at_q(long (*triplets)[3],
const long grid_point,
const long (*bz_grid_addresses)[3],
const long *bz_map,
const long *map_triplets,
const long num_map_triplets,
const long D_diag[3],
const long Q[3][3],
const long bz_grid_type);
long gridsys_get_integration_weight(double *iw,
char *iw_zero,
const double *frequency_points,
const long num_band0,
const long relative_grid_address[24][4][3],
const long D_diag[3],
const long (*triplets)[3],
const long num_triplets,
const long (*bz_grid_addresses)[3],
const long *bz_map,
const long bz_grid_type,
const double *frequencies1,
const long num_band1,
const double *frequencies2,
const long num_band2,
const long tp_type,
const long openmp_per_triplets,
const long openmp_per_bands);
void gridsys_get_integration_weight_with_sigma(double *iw,
char *iw_zero,
const double sigma,
const double sigma_cutoff,
const double *frequency_points,
const long num_band0,
const long (*triplets)[3],
const long num_triplets,
const double *frequencies,
const long num_band,
const long tp_type);
#ifdef __cplusplus
}

324
c/imag_self_energy_with_g.c
View File

@ -98,7 +98,7 @@ void ise_get_imag_self_energy_at_bands_with_g(double *imag_self_energy,
{
long i, j, num_triplets, num_band0, num_band, num_band_prod;
long num_g_pos, g_index_dims, g_index_shift;
long (*g_pos)[4];
long(*g_pos)[4];
double *ise;
long at_a_frequency_point;
@ -109,14 +109,17 @@ void ise_get_imag_self_energy_at_bands_with_g(double *imag_self_energy,
num_band0 = fc3_normal_squared->dims[1];
num_band = fc3_normal_squared->dims[2];
num_band_prod = num_band0 * num_band * num_band;
ise = (double*)malloc(sizeof(double) * num_triplets * num_band0);
ise = (double *)malloc(sizeof(double) * num_triplets * num_band0);
if (frequency_point_index < 0) {
if (frequency_point_index < 0)
{
/* frequency_points == frequencies at bands */
at_a_frequency_point = 0;
g_index_dims = num_band_prod;
g_index_shift = 0;
} else {
}
else
{
/* At an arbitrary frequency point. */
at_a_frequency_point = 1;
g_index_dims = num_frequency_points * num_band * num_band;
@ -126,19 +129,23 @@ void ise_get_imag_self_energy_at_bands_with_g(double *imag_self_energy,
#ifdef PHPYOPENMP
#pragma omp parallel for private(num_g_pos, j, g_pos)
#endif
for (i = 0; i < num_triplets; i++) {
for (i = 0; i < num_triplets; i++)
{
g_pos = (long(*)[4])malloc(sizeof(long[4]) * num_band_prod);
/* ise_set_g_pos only works for the case of frquency points at */
/* bands. For frequency sampling mode, g_zero is assumed all */
/* with the array shape of (num_triplets, num_band0, num_band, */
/* num_band). */
if (at_a_frequency_point) {
if (at_a_frequency_point)
{
num_g_pos = ise_set_g_pos_frequency_point(
g_pos,
num_band0,
num_band,
g_zero + i * g_index_dims + g_index_shift);
} else {
g_pos,
num_band0,
num_band,
g_zero + i * g_index_dims + g_index_shift);
}
else
{
num_g_pos = ise_set_g_pos(g_pos,
num_band0,
num_band,
@ -146,33 +153,36 @@ void ise_get_imag_self_energy_at_bands_with_g(double *imag_self_energy,
}
ise_imag_self_energy_at_triplet(
ise + i * num_band0,
num_band0,
num_band,
fc3_normal_squared->data + i * num_band_prod,
frequencies,
triplets[i],
triplet_weights[i],
g + i * g_index_dims + g_index_shift,
g + (i + num_triplets) * g_index_dims + g_index_shift,
g_pos,
num_g_pos,
&temperature,
1,
cutoff_frequency,
0,
at_a_frequency_point);
ise + i * num_band0,
num_band0,
num_band,
fc3_normal_squared->data + i * num_band_prod,
frequencies,
triplets[i],
triplet_weights[i],
g + i * g_index_dims + g_index_shift,
g + (i + num_triplets) * g_index_dims + g_index_shift,
g_pos,
num_g_pos,
&temperature,
1,
cutoff_frequency,
0,
at_a_frequency_point);
free(g_pos);
g_pos = NULL;
}
for (i = 0; i < num_band0; i++) {
for (i = 0; i < num_band0; i++)
{
imag_self_energy[i] = 0;
}
for (i = 0; i < num_triplets; i++) {
for (j = 0; j < num_band0; j++) {
for (i = 0; i < num_triplets; i++)
{
for (j = 0; j < num_band0; j++)
{
imag_self_energy[j] += ise[i * num_band0 + j];
}
}
@ -181,19 +191,18 @@ void ise_get_imag_self_energy_at_bands_with_g(double *imag_self_energy,
ise = NULL;
}
void ise_get_detailed_imag_self_energy_at_bands_with_g
(double *detailed_imag_self_energy,
double *imag_self_energy_N,
double *imag_self_energy_U,
const Darray *fc3_normal_squared,
const double *frequencies,
const long (*triplets)[3],
const long *triplet_weights,
const long (*bz_grid_addresses)[3],
const double *g,
const char *g_zero,
const double temperature,
const double cutoff_frequency)
void ise_get_detailed_imag_self_energy_at_bands_with_g(double *detailed_imag_self_energy,
double *imag_self_energy_N,
double *imag_self_energy_U,
const Darray *fc3_normal_squared,
const double *frequencies,
const long (*triplets)[3],
const long *triplet_weights,
const long (*bz_grid_addresses)[3],
const double *g,
const char *g_zero,
const double temperature,
const double cutoff_frequency)
{
double *ise;
long i, j, num_triplets, num_band0, num_band, num_band_prod;
@ -207,46 +216,52 @@ void ise_get_detailed_imag_self_energy_at_bands_with_g
num_band0 = fc3_normal_squared->dims[1];
num_band = fc3_normal_squared->dims[2];
num_band_prod = num_band0 * num_band * num_band;
ise = (double*)malloc(sizeof(double) * num_triplets * num_band0);
ise = (double *)malloc(sizeof(double) * num_triplets * num_band0);
/* detailed_imag_self_energy has the same shape as fc3_normal_squared. */
#ifdef PHPYOPENMP
#pragma omp parallel for
#endif
for (i = 0; i < num_triplets; i++) {
detailed_imag_self_energy_at_triplet
(detailed_imag_self_energy + i * num_band_prod,
ise + i * num_band0,
num_band0,
num_band,
fc3_normal_squared->data + i * num_band_prod,
frequencies,
triplets[i],
g + i * num_band_prod,
g + (i + num_triplets) * num_band_prod,
g_zero + i * num_band_prod,
&temperature,
1,
cutoff_frequency);
for (i = 0; i < num_triplets; i++)
{
detailed_imag_self_energy_at_triplet(detailed_imag_self_energy + i * num_band_prod,
ise + i * num_band0,
num_band0,
num_band,
fc3_normal_squared->data + i * num_band_prod,
frequencies,
triplets[i],
g + i * num_band_prod,
g + (i + num_triplets) * num_band_prod,
g_zero + i * num_band_prod,
&temperature,
1,
cutoff_frequency);
}
is_N = (long*)malloc(sizeof(long) * num_triplets);
for (i = 0; i < num_triplets; i++) {
is_N = (long *)malloc(sizeof(long) * num_triplets);
for (i = 0; i < num_triplets; i++)
{
is_N[i] = tpl_is_N(triplets[i], bz_grid_addresses);
}
for (i = 0; i < num_band0; i++) {
for (i = 0; i < num_band0; i++)
{
N = 0;
U = 0;
/* #ifdef PHPYOPENMP */
/* #pragma omp parallel for private(ise_tmp) reduction(+:N,U) */
/* #endif */
for (j = 0; j < num_triplets; j++) {
/* #ifdef PHPYOPENMP */
/* #pragma omp parallel for private(ise_tmp) reduction(+:N,U) */
/* #endif */
for (j = 0; j < num_triplets; j++)
{
ise_tmp = ise[j * num_band0 + i] * triplet_weights[j];
if (is_N[j]) {
if (is_N[j])
{
N += ise_tmp;
} else {
}
else
{
U += ise_tmp;
}
}
@ -281,9 +296,10 @@ void ise_imag_self_energy_at_triplet(double *imag_self_energy,
double *n1, *n2;
long g_pos_3;
n1 = (double*)malloc(sizeof(double) * num_temps * num_band);
n2 = (double*)malloc(sizeof(double) * num_temps * num_band);
for (i = 0; i < num_temps; i++) {
n1 = (double *)malloc(sizeof(double) * num_temps * num_band);
n2 = (double *)malloc(sizeof(double) * num_temps * num_band);
for (i = 0; i < num_temps; i++)
{
set_occupations(n1 + i * num_band,
n2 + i * num_band,
num_band,
@ -293,36 +309,50 @@ void ise_imag_self_energy_at_triplet(double *imag_self_energy,
cutoff_frequency);
}
for (i = 0; i < num_band0 * num_temps; i++) {
for (i = 0; i < num_band0 * num_temps; i++)
{
imag_self_energy[i] = 0;
}
/* Do not use OpenMP here!! */
/* g_pos[i][0] takes value 0 <= x < num_band0 only, */
/* which causes race condition. */
for (i = 0; i < num_g_pos; i++) {
if (at_a_frequency_point) {
/* Do not use OpenMP here!! */
/* g_pos[i][0] takes value 0 <= x < num_band0 only, */
/* which causes race condition. */
for (i = 0; i < num_g_pos; i++)
{
if (at_a_frequency_point)
{
/* At an arbitrary frequency point */
g_pos_3 = g_pos[i][3] % (num_band * num_band);
} else {
}
else
{
/* frequency_points == frequencies at bands */
g_pos_3 = g_pos[i][3];
}
for (j = 0; j < num_temps; j++) {
for (j = 0; j < num_temps; j++)
{
if (n1[j * num_band + g_pos[i][1]] < 0 ||
n2[j * num_band + g_pos[i][2]] < 0) {
n2[j * num_band + g_pos[i][2]] < 0)
{
;
} else {
if (temperatures[j] > 0) {
}
else
{
if (temperatures[j] > 0)
{
imag_self_energy[j * num_band0 + g_pos[i][0]] +=
((n1[j * num_band + g_pos[i][1]] +
n2[j * num_band + g_pos[i][2]] + 1) * g1[g_pos_3] +
(n1[j * num_band + g_pos[i][1]] -
n2[j * num_band + g_pos[i][2]]) * g2_3[g_pos_3]) *
fc3_normal_squared[g_pos[i][3]] * triplet_weight;
} else {
((n1[j * num_band + g_pos[i][1]] +
n2[j * num_band + g_pos[i][2]] + 1) *
g1[g_pos_3] +
(n1[j * num_band + g_pos[i][1]] -
n2[j * num_band + g_pos[i][2]]) *
g2_3[g_pos_3]) *
fc3_normal_squared[g_pos[i][3]] * triplet_weight;
}
else
{
imag_self_energy[j * num_band0 + g_pos[i][0]] +=
g1[g_pos_3] * fc3_normal_squared[g_pos[i][3]] * triplet_weight;
g1[g_pos_3] * fc3_normal_squared[g_pos[i][3]] * triplet_weight;
}
}
}
@ -343,10 +373,14 @@ long ise_set_g_pos(long (*g_pos)[4],
num_g_pos = 0;
jkl = 0;
for (j = 0; j < num_band0; j++) {
for (k = 0; k < num_band; k++) {
for (l = 0; l < num_band; l++) {
if (!g_zero[jkl]) {
for (j = 0; j < num_band0; j++)
{
for (k = 0; k < num_band; k++)
{
for (l = 0; l < num_band; l++)
{
if (!g_zero[jkl])
{
g_pos[num_g_pos][0] = j;
g_pos[num_g_pos][1] = k;
g_pos[num_g_pos][2] = l;
@ -369,11 +403,15 @@ static long ise_set_g_pos_frequency_point(long (*g_pos)[4],
num_g_pos = 0;
jkl = 0;
for (j = 0; j < num_band0; j++) {
for (j = 0; j < num_band0; j++)
{
kl = 0;
for (k = 0; k < num_band; k++) {
for (l = 0; l < num_band; l++) {
if (!g_zero[kl]) {
for (k = 0; k < num_band; k++)
{
for (l = 0; l < num_band; l++)
{
if (!g_zero[kl])
{
g_pos[num_g_pos][0] = j;
g_pos[num_g_pos][1] = k;
g_pos[num_g_pos][2] = l;
@ -409,10 +447,11 @@ detailed_imag_self_energy_at_triplet(double *detailed_imag_self_energy,
n1 = NULL;
n2 = NULL;
n1 = (double*)malloc(sizeof(double) * num_band);
n2 = (double*)malloc(sizeof(double) * num_band);
n1 = (double *)malloc(sizeof(double) * num_band);
n2 = (double *)malloc(sizeof(double) * num_band);
for (i = 0; i < num_temps; i++) {
for (i = 0; i < num_temps; i++)
{
set_occupations(n1,
n2,
num_band,
@ -421,29 +460,31 @@ detailed_imag_self_energy_at_triplet(double *detailed_imag_self_energy,
frequencies,
cutoff_frequency);
for (j = 0; j < num_band0; j++) {
for (j = 0; j < num_band0; j++)
{
adrs_shift = j * num_band * num_band;
if (temperatures[i] > 0) {
if (temperatures[i] > 0)
{
imag_self_energy[i * num_band0 + j] =
collect_detailed_imag_self_energy
(detailed_imag_self_energy + adrs_shift,
num_band,
fc3_normal_squared + adrs_shift,
n1,
n2,
g1 + adrs_shift,
g2_3 + adrs_shift,
g_zero + adrs_shift);
} else {
collect_detailed_imag_self_energy(detailed_imag_self_energy + adrs_shift,
num_band,
fc3_normal_squared + adrs_shift,
n1,
n2,
g1 + adrs_shift,
g2_3 + adrs_shift,
g_zero + adrs_shift);
}
else
{
imag_self_energy[i * num_band0 + j] =
collect_detailed_imag_self_energy_0K
(detailed_imag_self_energy + adrs_shift,
num_band,
fc3_normal_squared + adrs_shift,
n1,
n2,
g1 + adrs_shift,
g_zero + adrs_shift);
collect_detailed_imag_self_energy_0K(detailed_imag_self_energy + adrs_shift,
num_band,
fc3_normal_squared + adrs_shift,
n1,
n2,
g1 + adrs_shift,
g_zero + adrs_shift);
}
}
}
@ -468,12 +509,19 @@ collect_detailed_imag_self_energy(double *imag_self_energy,
double sum_g;
sum_g = 0;
for (ij = 0; ij < num_band * num_band; ij++) {
for (ij = 0; ij < num_band * num_band; ij++)
{
imag_self_energy[ij] = 0;
if (g_zero[ij]) {continue;}
if (g_zero[ij])
{
continue;
}
i = ij / num_band;
j = ij % num_band;
if (n1[i] < 0 || n2[j] < 0) {continue;}
if (n1[i] < 0 || n2[j] < 0)
{
continue;
}
imag_self_energy[ij] = (((n1[i] + n2[j] + 1) * g1[ij] +
(n1[i] - n2[j]) * g2_3[ij]) *
fc3_normal_squared[ij]);
@ -496,12 +544,19 @@ collect_detailed_imag_self_energy_0K(double *imag_self_energy,
double sum_g;
sum_g = 0;
for (ij = 0; ij < num_band * num_band; ij++) {
for (ij = 0; ij < num_band * num_band; ij++)
{
imag_self_energy[ij] = 0;
if (g_zero[ij]) {continue;}
if (g_zero[ij])
{
continue;
}
i = ij / num_band;
j = ij % num_band;
if (n1[i] < 0 || n2[j] < 0) {continue;}
if (n1[i] < 0 || n2[j] < 0)
{
continue;
}
imag_self_energy[ij] = g1[ij] * fc3_normal_squared[ij];
sum_g += imag_self_energy[ij];
}
@ -520,17 +575,24 @@ static void set_occupations(double *n1,
long j;
double f1, f2;
for (j = 0; j < num_band; j++) {
for (j = 0; j < num_band; j++)
{
f1 = frequencies[triplet[1] * num_band + j];
f2 = frequencies[triplet[2] * num_band + j];
if (f1 > cutoff_frequency) {
if (f1 > cutoff_frequency)
{
n1[j] = phonoc_bose_einstein(f1, temperature);
} else {
}
else
{
n1[j] = -1;
}
if (f2 > cutoff_frequency) {
if (f2 > cutoff_frequency)
{
n2[j] = phonoc_bose_einstein(f2, temperature);
} else {
}
else
{
n2[j] = -1;
}
}

25
c/imag_self_energy_with_g.h
View File

@ -49,19 +49,18 @@ void ise_get_imag_self_energy_at_bands_with_g(double *imag_self_energy,
const double cutoff_frequency,
const long num_frequency_points,
const long frequency_point_index);
void ise_get_detailed_imag_self_energy_at_bands_with_g
(double *detailed_imag_self_energy,
double *imag_self_energy_N,
double *imag_self_energy_U,
const Darray *fc3_normal_squared,
const double *frequencies,
const long (*triplets)[3],
const long *triplet_weights,
const long (*bz_grid_addresses)[3],
const double *g,
const char *g_zero,
const double temperature,
const double cutoff_frequency);
void ise_get_detailed_imag_self_energy_at_bands_with_g(double *detailed_imag_self_energy,
double *imag_self_energy_N,
double *imag_self_energy_U,
const Darray *fc3_normal_squared,
const double *frequencies,
const long (*triplets)[3],
const long *triplet_weights,
const long (*bz_grid_addresses)[3],
const double *g,
const char *g_zero,
const double temperature,
const double cutoff_frequency);
void ise_imag_self_energy_at_triplet(double *imag_self_energy,
const long num_band0,
const long num_band,

182
c/isotope.c
View File

@ -38,68 +38,74 @@
#include "isotope.h"
#include "lapack_wrapper.h"
void
iso_get_isotope_scattering_strength(double *gamma,
const long grid_point,
const double *mass_variances,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const long num_grid_points,
const long *band_indices,
const long num_band,
const long num_band0,
const double sigma,
const double cutoff_frequency)
void iso_get_isotope_scattering_strength(double *gamma,
const long grid_point,
const double *mass_variances,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const long num_grid_points,
const long *band_indices,
const long num_band,
const long num_band0,
const double sigma,
const double cutoff_frequency)
{
long i, j, k, l, m;
double *e0_r, *e0_i, e1_r, e1_i, a, b, f, *f0, dist, sum_g, 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);
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++) {
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]]);
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]]);
}
}
for (i = 0; i < num_band0; i++) {
for (i = 0; i < num_band0; i++)
{
gamma[i] = 0;
}
for (i = 0; i < num_band0; i++) { /* band index0 */
if (f0[i] < cutoff_frequency) {
for (i = 0; i < num_band0; i++)
{ /* band index0 */
if (f0[i] < cutoff_frequency)
{
continue;
}
sum_g = 0;
#ifdef PHPYOPENMP
#pragma omp parallel for private(k, l, m, f, e1_r, e1_i, a, b, dist, sum_g_k) reduction(+:sum_g)
#pragma omp parallel for private(k, l, m, f, e1_r, e1_i, a, b, dist, sum_g_k) reduction(+ \
: sum_g)
#endif
for (j = 0; j < num_grid_points; j++) {
for (j = 0; j < num_grid_points; j++)
{
sum_g_k = 0;
for (k = 0; k < num_band; k++) { /* band index */
for (k = 0; k < num_band; k++)
{ /* band index */
f = frequencies[j * num_band + k];
if (f < cutoff_frequency) {
if (f < cutoff_frequency)
{
continue;
}
dist = phonoc_gaussian(f - f0[i], sigma);
for (l = 0; l < num_band / 3; l++) { /* elements */
for (l = 0; l < num_band / 3; l++)
{ /* elements */
a = 0;
b = 0;
for (m = 0; m < 3; m++) {
e1_r = lapack_complex_double_real
(eigenvectors[j * num_band * num_band +
(l * 3 + m) * num_band + k]);
e1_i = lapack_complex_double_imag
(eigenvectors[j * num_band * num_band +
(l * 3 + m) * num_band + k]);
for (m = 0; m < 3; m++)
{
e1_r = lapack_complex_double_real(eigenvectors[j * num_band * num_band +
(l * 3 + m) * num_band + k]);
e1_i = lapack_complex_double_imag(eigenvectors[j * num_band * num_band +
(l * 3 + m) * num_band + k]);
a += (e0_r[i * num_band + l * 3 + m] * e1_r +
e0_i[i * num_band + l * 3 + m] * e1_i);
b += (e0_i[i * num_band + l * 3 + m] * e1_r -
@ -113,7 +119,8 @@ iso_get_isotope_scattering_strength(double *gamma,
gamma[i] = sum_g;
}
for (i = 0; i < num_band0; i++) {
for (i = 0; i < num_band0; i++)
{
/* Frequency unit to ang-freq: *(2pi)**2/(2pi) */
/* Ang-freq to freq unit (for lifetime): /2pi */
/* gamma = 1/2t */
@ -128,76 +135,81 @@ iso_get_isotope_scattering_strength(double *gamma,
e0_i = NULL;
}
void iso_get_thm_isotope_scattering_strength
(double *gamma,
const long grid_point,
const long *ir_grid_points,
const long *weights,
const double *mass_variances,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const long num_grid_points,
const long *band_indices,
const long num_band,
const long num_band0,
const double *integration_weights,
const double cutoff_frequency)
void iso_get_thm_isotope_scattering_strength(double *gamma,
const long grid_point,
const long *ir_grid_points,
const long *weights,
const double *mass_variances,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const long num_grid_points,
const long *band_indices,
const long num_band,
const long num_band0,
const double *integration_weights,
const double cutoff_frequency)
{
long i, j, k, l, m, 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);
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++) {
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]]);
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);
gamma_ij = (double *)malloc(sizeof(double) * num_grid_points * num_band0);
#ifdef PHPYOPENMP
#pragma omp parallel for
#endif
for (i = 0; i < num_grid_points * num_band0; i++) {
for (i = 0; i < num_grid_points * num_band0; i++)
{
gamma_ij[i] = 0;
}
#ifdef PHPYOPENMP
#pragma omp parallel for private(j, k, l, m, f, gp, e1_r, e1_i, a, b, dist, sum_g_k)
#endif
for (i = 0; i < num_grid_points; i++) {
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) {
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 */
for (k = 0; k < num_band; k++)
{ /* band index */
f = frequencies[gp * num_band + k];
if (f < cutoff_frequency) {
if (f < cutoff_frequency)
{
continue;
}
dist = integration_weights[gp * num_band0 * num_band +
j * num_band + k];
for (l = 0; l < num_band / 3; l++) { /* elements */
for (l = 0; l < num_band / 3; l++)
{ /* elements */
a = 0;
b = 0;
for (m = 0; m < 3; m++) {
e1_r = lapack_complex_double_real
(eigenvectors
[gp * num_band * num_band + (l * 3 + m) * num_band + k]);
e1_i = lapack_complex_double_imag
(eigenvectors
[gp * num_band * num_band + (l * 3 + m) * num_band + k]);
for (m = 0; m < 3; m++)
{
e1_r = lapack_complex_double_real(eigenvectors
[gp * num_band * num_band + (l * 3 + m) * num_band + k]);
e1_i = lapack_complex_double_imag(eigenvectors
[gp * num_band * num_band + (l * 3 + m) * num_band + k]);
a += (e0_r[j * num_band + l * 3 + m] * e1_r +
e0_i[j * num_band + l * 3 + m] * e1_i);
b += (e0_i[j * num_band + l * 3 + m] * e1_r -
@ -210,18 +222,22 @@ void iso_get_thm_isotope_scattering_strength
}
}
for (i = 0; i < num_band0; i++) {
for (i = 0; i < num_band0; i++)
{
gamma[i] = 0;
}
for (i = 0; i < num_grid_points; i++) {
for (i = 0; i < num_grid_points; i++)
{
gp = ir_grid_points[i];
for (j = 0; j < num_band0; j++) {
for (j = 0; j < num_band0; j++)
{
gamma[j] += gamma_ij[gp * num_band0 + j];
}
}
for (i = 0; i < num_band0; i++) {
for (i = 0; i < num_band0; i++)
{
/* Frequency unit to ang-freq: *(2pi)**2/(2pi) */
/* Ang-freq to freq unit (for lifetime): /2pi */
/* gamma = 1/2t */

50
c/isotope.h
View File

@ -37,30 +37,28 @@
#include "lapack_wrapper.h"
void
iso_get_isotope_scattering_strength(double *gamma,
const long grid_point,
const double *mass_variances,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const long num_grid_points,
const long *band_indices,
const long num_band,
const long num_band0,
const double sigma,
const double cutoff_frequency);
void iso_get_thm_isotope_scattering_strength
(double *gamma,
const long grid_point,
const long *ir_grid_points,
const long *weights,
const double *mass_variances,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const long num_grid_points,
const long *band_indices,
const long num_band,
const long num_band0,
const double *integration_weights,
const double cutoff_frequency);
void iso_get_isotope_scattering_strength(double *gamma,
const long grid_point,
const double *mass_variances,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const long num_grid_points,
const long *band_indices,
const long num_band,
const long num_band0,
const double sigma,
const double cutoff_frequency);
void iso_get_thm_isotope_scattering_strength(double *gamma,
const long grid_point,
const long *ir_grid_points,
const long *weights,
const double *mass_variances,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const long num_grid_points,
const long *band_indices,
const long num_band,
const long num_band0,
const double *integration_weights,
const double cutoff_frequency);
#endif

69
c/lagrid.c
View File

@ -36,16 +36,12 @@
long lagmat_get_determinant_l3(const long a[3][3])
{
return a[0][0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1])
+ a[0][1] * (a[1][2] * a[2][0] - a[1][0] * a[2][2])
+ a[0][2] * (a[1][0] * a[2][1] - a[1][1] * a[2][0]);
return a[0][0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1]) + a[0][1] * (a[1][2] * a[2][0] - a[1][0] * a[2][2]) + a[0][2] * (a[1][0] * a[2][1] - a[1][1] * a[2][0]);
}
double lagmat_get_determinant_d3(const double a[3][3])
{
return a[0][0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1])
+ a[0][1] * (a[1][2] * a[2][0] - a[1][0] * a[2][2])
+ a[0][2] * (a[1][0] * a[2][1] - a[1][1] * a[2][0]);
return a[0][0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1]) + a[0][1] * (a[1][2] * a[2][0] - a[1][0] * a[2][2]) + a[0][2] * (a[1][0] * a[2][1] - a[1][1] * a[2][0]);
}
void lagmat_cast_matrix_3l_to_3d(double m[3][3], const long a[3][3])
@ -80,7 +76,8 @@ long lagmat_get_similar_matrix_ld3(double m[3][3],
const double precision)
{
double c[3][3];
if (!lagmat_inverse_matrix_d3(c, b, precision)) {
if (!lagmat_inverse_matrix_d3(c, b, precision))
{
warning_print("No similar matrix due to 0 determinant.\n");
return 0;
}
@ -100,10 +97,12 @@ long lagmat_check_identity_matrix_l3(const long a[3][3],
a[1][2] - b[1][2] ||
a[2][0] - b[2][0] ||
a[2][1] - b[2][1] ||
a[2][2] - b[2][2]) {
a[2][2] - b[2][2])
{
return 0;
}
else {
else
{
return 1;
}
}
@ -120,10 +119,12 @@ long lagmat_check_identity_matrix_ld3(const long a[3][3],
lagmat_Dabs(a[1][2] - b[1][2]) > symprec ||
lagmat_Dabs(a[2][0] - b[2][0]) > symprec ||
lagmat_Dabs(a[2][1] - b[2][1]) > symprec ||
lagmat_Dabs(a[2][2] - b[2][2]) > symprec) {
lagmat_Dabs(a[2][2] - b[2][2]) > symprec)
{
return 0;
}
else {
else
{
return 1;
}
}
@ -135,7 +136,8 @@ long lagmat_inverse_matrix_d3(double m[3][3],
double det;
double c[3][3];
det = lagmat_get_determinant_d3(a);
if (lagmat_Dabs(det) < precision) {
if (lagmat_Dabs(det) < precision)
{
warning_print("No inverse matrix (det=%f)\n", det);
return 0;
}
@ -174,10 +176,12 @@ void lagmat_multiply_matrix_vector_l3(long v[3],
{
long i;
long c[3];
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
c[i] = a[i][0] * b[0] + a[i][1] * b[1] + a[i][2] * b[2];
}
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
v[i] = c[i];
}
}
@ -186,12 +190,14 @@ void lagmat_multiply_matrix_l3(long m[3][3],
const long a[3][3],
const long b[3][3])
{
long i, j; /* a_ij */
long i, j; /* a_ij */
long c[3][3];
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
for (i = 0; i < 3; i++)
{
for (j = 0; j < 3; j++)
{
c[i][j] =
a[i][0] * b[0][j] + a[i][1] * b[1][j] + a[i][2] * b[2][j];
a[i][0] * b[0][j] + a[i][1] * b[1][j] + a[i][2] * b[2][j];
}
}
lagmat_copy_matrix_l3(m, c);
@ -201,12 +207,14 @@ void lagmat_multiply_matrix_ld3(double m[3][3],
const long a[3][3],
const double b[3][3])
{
long i, j; /* a_ij */
long i, j; /* a_ij */
double c[3][3];
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
for (i = 0; i < 3; i++)
{
for (j = 0; j < 3; j++)
{
c[i][j] =
a[i][0] * b[0][j] + a[i][1] * b[1][j] + a[i][2] * b[2][j];
a[i][0] * b[0][j] + a[i][1] * b[1][j] + a[i][2] * b[2][j];
}
}
lagmat_copy_matrix_d3(m, c);
@ -216,12 +224,14 @@ void lagmat_multiply_matrix_d3(double m[3][3],
const double a[3][3],
const double b[3][3])
{
long i, j; /* a_ij */
long i, j; /* a_ij */
double c[3][3];
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
for (i = 0; i < 3; i++)
{
for (j = 0; j < 3; j++)
{
c[i][j] =
a[i][0] * b[0][j] + a[i][1] * b[1][j] + a[i][2] * b[2][j];
a[i][0] * b[0][j] + a[i][1] * b[1][j] + a[i][2] * b[2][j];
}
}
lagmat_copy_matrix_d3(m, c);
@ -264,7 +274,8 @@ long lagmat_modulo_l(const long a, const long b)
{
long c;
c = a % b;
if (c < 0) {
if (c < 0)
{
c += b;
}
return c;
@ -273,9 +284,9 @@ long lagmat_modulo_l(const long a, const long b)
long lagmat_Nint(const double a)
{
if (a < 0.0)
return (long) (a - 0.5);
return (long)(a - 0.5);
else
return (long) (a + 0.5);
return (long)(a + 0.5);
}
double lagmat_Dabs(const double a)

72
c/lapack_wrapper.c
View File

@ -34,20 +34,21 @@
#include "lapack_wrapper.h"
#define min(a,b) ((a)>(b)?(b):(a))
#define min(a, b) ((a) > (b) ? (b) : (a))
#ifdef MKL_LAPACKE
MKL_Complex16 lapack_make_complex_double( double re, double im ) {
MKL_Complex16 lapack_make_complex_double(double re, double im)
{
MKL_Complex16 z;
z.real = re;
z.imag = im;
return z;
}
#ifndef LAPACKE_malloc
#define LAPACKE_malloc( size ) malloc( size )
#define LAPACKE_malloc(size) malloc(size)
#endif
#ifndef LAPACKE_free
#define LAPACKE_free( p ) free( p )
#define LAPACKE_free(p) free(p)
#endif
#endif
@ -57,7 +58,7 @@ int phonopy_zheev(double *w,
const char uplo)
{
lapack_int info;
info = LAPACKE_zheev(LAPACK_ROW_MAJOR,'V', uplo,
info = LAPACKE_zheev(LAPACK_ROW_MAJOR, 'V', uplo,
(lapack_int)n, a, (lapack_int)n, w);
return (int)info;
}
@ -72,13 +73,14 @@ int phonopy_pinv(double *data_out,
lapack_int info;
double *s, *a, *u, *vt, *superb;
a = (double*)malloc(sizeof(double) * m * n);
s = (double*)malloc(sizeof(double) * min(m,n));
u = (double*)malloc(sizeof(double) * m * m);
vt = (double*)malloc(sizeof(double) * n * n);
superb = (double*)malloc(sizeof(double) * (min(m,n) - 1));
a = (double *)malloc(sizeof(double) * m * n);
s = (double *)malloc(sizeof(double) * min(m, n));
u = (double *)malloc(sizeof(double) * m * m);
vt = (double *)malloc(sizeof(double) * n * n);
superb = (double *)malloc(sizeof(double) * (min(m, n) - 1));
for (i = 0; i < m * n; i++) {
for (i = 0; i < m * n; i++)
{
a[i] = data_in[i];
}
@ -96,14 +98,19 @@ int phonopy_pinv(double *data_out,
(lapack_int)n,
superb);
for (i = 0; i < n * m; i++) {
for (i = 0; i < n * m; i++)
{
data_out[i] = 0;
}
for (i = 0; i < m; i++) {
for (j = 0; j < n; j++) {
for (k = 0; k < min(m,n); k++) {
if (s[k] > cutoff) {
for (i = 0; i < m; i++)
{
for (j = 0; j < n; j++)
{
for (k = 0; k < min(m, n); k++)
{
if (s[k] > cutoff)
{
data_out[j * m + i] += u[i * m + k] / s[k] * vt[k * n + j];
}
}
@ -133,7 +140,8 @@ void phonopy_pinv_mt(double *data_out,
#ifdef PHPYOPENMP
#pragma omp parallel for
#endif
for (i = 0; i < num_thread; i++) {
for (i = 0; i < num_thread; i++)
{
info_out[i] = phonopy_pinv(data_out + i * max_row_num * column_num,
data_in + i * max_row_num * column_num,
row_nums[i],
@ -151,8 +159,8 @@ int phonopy_dsyev(double *data,
lapack_int liwork;
long lwork;
lapack_int* iwork;
double* work;
lapack_int *iwork;
double *work;
lapack_int iwork_query;
double work_query;
@ -162,7 +170,8 @@ int phonopy_dsyev(double *data,
iwork = NULL;
work = NULL;
switch (algorithm) {
switch (algorithm)
{
case 0: /* dsyev */
info = LAPACKE_dsyev(LAPACK_COL_MAJOR,
'V',
@ -187,11 +196,12 @@ int phonopy_dsyev(double *data,
liwork = iwork_query;
lwork = (long)work_query;
/* printf("liwork %d, lwork %ld\n", liwork, lwork); */
if ((iwork = (lapack_int*)LAPACKE_malloc(sizeof(lapack_int) * liwork))
== NULL) {
if ((iwork = (lapack_int *)LAPACKE_malloc(sizeof(lapack_int) * liwork)) == NULL)
{
goto end;
};
if ((work = (double*)LAPACKE_malloc(sizeof(double) * lwork)) == NULL) {
if ((work = (double *)LAPACKE_malloc(sizeof(double) * lwork)) == NULL)
{
goto end;
}
@ -208,11 +218,13 @@ int phonopy_dsyev(double *data,
liwork);
end:
if (iwork) {
if (iwork)
{
LAPACKE_free(iwork);
iwork = NULL;
}
if (work) {
if (work)
{
LAPACKE_free(work);
work = NULL;
}
@ -230,16 +242,14 @@ int phonopy_dsyev(double *data,
return (int)info;
}
lapack_complex_double
phonoc_complex_prod(const lapack_complex_double a,
const lapack_complex_double b)
{
lapack_complex_double c;
c = lapack_make_complex_double
(lapack_complex_double_real(a) * lapack_complex_double_real(b) -
lapack_complex_double_imag(a) * lapack_complex_double_imag(b),
lapack_complex_double_imag(a) * lapack_complex_double_real(b) +
lapack_complex_double_real(a) * lapack_complex_double_imag(b));
c = lapack_make_complex_double(lapack_complex_double_real(a) * lapack_complex_double_real(b) -
lapack_complex_double_imag(a) * lapack_complex_double_imag(b),
lapack_complex_double_imag(a) * lapack_complex_double_real(b) +
lapack_complex_double_real(a) * lapack_complex_double_imag(b));
return c;
}

399
c/phono3py.c
View File

@ -53,7 +53,6 @@
#include <stdio.h>
#include <stdlib.h>
long ph3py_get_interaction(Darray *fc3_normal_squared,
const char *g_zero,
const Darray *frequencies,
@ -78,16 +77,19 @@ long ph3py_get_interaction(Darray *fc3_normal_squared,
ConstBZGrid *bzgrid;
long i, j;
if ((bzgrid = (ConstBZGrid*) malloc(sizeof(ConstBZGrid))) == NULL) {
if ((bzgrid = (ConstBZGrid *)malloc(sizeof(ConstBZGrid))) == NULL)
{
warning_print("Memory could not be allocated.");
return 0;
}
bzgrid->addresses = bz_grid_addresses;
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
bzgrid->D_diag[i] = D_diag[i];
bzgrid->PS[i] = 0;
for (j = 0; j < 3; j++) {
for (j = 0; j < 3; j++)
{
bzgrid->Q[i][j] = Q[i][j];
}
}
@ -116,7 +118,6 @@ long ph3py_get_interaction(Darray *fc3_normal_squared,
return 1;
}
long ph3py_get_pp_collision(double *imag_self_energy,
const long relative_grid_address[24][4][3], /* thm */
const double *frequencies,
@ -125,7 +126,7 @@ long ph3py_get_pp_collision(double *imag_self_energy,
const long num_triplets,
const long *triplet_weights,
const long (*bz_grid_addresses)[3], /* thm */
const long *bz_map, /* thm */
const long *bz_map, /* thm */
const long bz_grid_type,
const long D_diag[3],
const long Q[3][3],
@ -146,7 +147,8 @@ long ph3py_get_pp_collision(double *imag_self_energy,
ConstBZGrid *bzgrid;
long i, j;
if ((bzgrid = (ConstBZGrid*) malloc(sizeof(ConstBZGrid))) == NULL) {
if ((bzgrid = (ConstBZGrid *)malloc(sizeof(ConstBZGrid))) == NULL)
{
warning_print("Memory could not be allocated.");
return 0;
}
@ -154,10 +156,12 @@ long ph3py_get_pp_collision(double *imag_self_energy,
bzgrid->addresses = bz_grid_addresses;
bzgrid->gp_map = bz_map;
bzgrid->type = bz_grid_type;
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
bzgrid->D_diag[i] = D_diag[i];
bzgrid->PS[i] = 0;
for (j = 0; j < 3; j++) {
for (j = 0; j < 3; j++)
{
bzgrid->Q[i][j] = Q[i][j];
}
}
@ -190,46 +194,48 @@ long ph3py_get_pp_collision(double *imag_self_energy,
return 1;
}
long ph3py_get_pp_collision_with_sigma(
double *imag_self_energy,
const double sigma,
const double sigma_cutoff,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const long (*triplets)[3],
const long num_triplets,
const long *triplet_weights,
const long (*bz_grid_addresses)[3],
const long D_diag[3],
const long Q[3][3],
const double *fc3,
const long is_compact_fc3,
const double (*svecs)[3],
const long multi_dims[2],
const long (*multiplicity)[2],
const double *masses,
const long *p2s_map,
const long *s2p_map,
const Larray *band_indices,
const Darray *temperatures,
const long is_NU,
const long symmetrize_fc3_q,
const double cutoff_frequency)
double *imag_self_energy,
const double sigma,
const double sigma_cutoff,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const long (*triplets)[3],
const long num_triplets,
const long *triplet_weights,
const long (*bz_grid_addresses)[3],
const long D_diag[3],
const long Q[3][3],
const double *fc3,
const long is_compact_fc3,
const double (*svecs)[3],
const long multi_dims[2],
const long (*multiplicity)[2],
const double *masses,
const long *p2s_map,
const long *s2p_map,
const Larray *band_indices,
const Darray *temperatures,
const long is_NU,
const long symmetrize_fc3_q,
const double cutoff_frequency)
{
ConstBZGrid *bzgrid;
long i, j;
if ((bzgrid = (ConstBZGrid*) malloc(sizeof(ConstBZGrid))) == NULL) {
if ((bzgrid = (ConstBZGrid *)malloc(sizeof(ConstBZGrid))) == NULL)
{
warning_print("Memory could not be allocated.");
return 0;
}
bzgrid->addresses = bz_grid_addresses;
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
bzgrid->D_diag[i] = D_diag[i];
bzgrid->PS[i] = 0;
for (j = 0; j < 3; j++) {
for (j = 0; j < 3; j++)
{
bzgrid->Q[i][j] = Q[i][j];
}
}
@ -263,19 +269,18 @@ long ph3py_get_pp_collision_with_sigma(
return 1;
}
void ph3py_get_imag_self_energy_at_bands_with_g(
double *imag_self_energy,
const Darray *fc3_normal_squared,
const double *frequencies,
const long (*triplets)[3],
const long *triplet_weights,
const double *g,
const char *g_zero,
const double temperature,
const double cutoff_frequency,
const long num_frequency_points,
const long frequency_point_index)
double *imag_self_energy,
const Darray *fc3_normal_squared,
const double *frequencies,
const long (*triplets)[3],
const long *triplet_weights,
const double *g,
const char *g_zero,
const double temperature,
const double cutoff_frequency,
const long num_frequency_points,
const long frequency_point_index)
{
ise_get_imag_self_energy_at_bands_with_g(imag_self_energy,
fc3_normal_squared,
@ -290,20 +295,19 @@ void ph3py_get_imag_self_energy_at_bands_with_g(
frequency_point_index);
}
void ph3py_get_detailed_imag_self_energy_at_bands_with_g(
double *detailed_imag_self_energy,
double *imag_self_energy_N,
double *imag_self_energy_U,
const Darray *fc3_normal_squared,
const double *frequencies,
const long (*triplets)[3],
const long *triplet_weights,
const long (*bz_grid_addresses)[3],
const double *g,
const char *g_zero,
const double temperature,
const double cutoff_frequency)
double *detailed_imag_self_energy,
double *imag_self_energy_N,
double *imag_self_energy_U,
const Darray *fc3_normal_squared,
const double *frequencies,
const long (*triplets)[3],
const long *triplet_weights,
const long (*bz_grid_addresses)[3],
const double *g,
const char *g_zero,
const double temperature,
const double cutoff_frequency)
{
ise_get_detailed_imag_self_energy_at_bands_with_g(detailed_imag_self_energy,
imag_self_energy_N,
@ -319,7 +323,6 @@ void ph3py_get_detailed_imag_self_energy_at_bands_with_g(
cutoff_frequency);
}
void ph3py_get_real_self_energy_at_bands(double *real_self_energy,
const Darray *fc3_normal_squared,
const long *band_indices,
@ -343,19 +346,18 @@ void ph3py_get_real_self_energy_at_bands(double *real_self_energy,
cutoff_frequency);
}
void ph3py_get_real_self_energy_at_frequency_point(
double *real_self_energy,
const double frequency_point,
const Darray *fc3_normal_squared,
const long *band_indices,
const double *frequencies,
const long (*triplets)[3],
const long *triplet_weights,
const double epsilon,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency)
double *real_self_energy,
const double frequency_point,
const Darray *fc3_normal_squared,
const long *band_indices,
const double *frequencies,
const long (*triplets)[3],
const long *triplet_weights,
const double epsilon,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency)
{
rse_get_real_self_energy_at_frequency_point(real_self_energy,
frequency_point,
@ -370,7 +372,6 @@ void ph3py_get_real_self_energy_at_frequency_point(
cutoff_frequency);
}
void ph3py_get_collision_matrix(double *collision_matrix,
const Darray *fc3_normal_squared,
const double *frequencies,
@ -404,7 +405,6 @@ void ph3py_get_collision_matrix(double *collision_matrix,
cutoff_frequency);
}
void ph3py_get_reducible_collision_matrix(double *collision_matrix,
const Darray *fc3_normal_squared,
const double *frequencies,
@ -430,19 +430,18 @@ void ph3py_get_reducible_collision_matrix(double *collision_matrix,
cutoff_frequency);
}
void ph3py_get_isotope_scattering_strength(
double *gamma,
const long grid_point,
const double *mass_variances,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const long num_grid_points,
const long *band_indices,
const long num_band,
const long num_band0,
const double sigma,
const double cutoff_frequency)
double *gamma,
const long grid_point,
const double *mass_variances,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const long num_grid_points,
const long *band_indices,
const long num_band,
const long num_band0,
const double sigma,
const double cutoff_frequency)
{
iso_get_isotope_scattering_strength(gamma,
grid_point,
@ -457,21 +456,19 @@ void ph3py_get_isotope_scattering_strength(
cutoff_frequency);
}
void ph3py_get_thm_isotope_scattering_strength
(double *gamma,
const long grid_point,
const long *ir_grid_points,
const long *weights,
const double *mass_variances,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const long num_ir_grid_points,
const long *band_indices,
const long num_band,
const long num_band0,
const double *integration_weights,
const double cutoff_frequency)
void ph3py_get_thm_isotope_scattering_strength(double *gamma,
const long grid_point,
const long *ir_grid_points,
const long *weights,
const double *mass_variances,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const long num_ir_grid_points,
const long *band_indices,
const long num_band,
const long num_band0,
const double *integration_weights,
const double cutoff_frequency)
{
iso_get_thm_isotope_scattering_strength(gamma,
grid_point,
@ -503,7 +500,6 @@ void ph3py_distribute_fc3(double *fc3,
rot_cart);
}
void ph3py_rotate_delta_fc2(double (*fc3)[3][3][3],
const double (*delta_fc2s)[3][3],
const double *inv_U,
@ -523,14 +519,12 @@ void ph3py_rotate_delta_fc2(double (*fc3)[3][3][3],
num_disp);
}
void ph3py_get_permutation_symmetry_fc3(double *fc3, const long num_atom)
{
fc3_set_permutation_symmetry_fc3(fc3, num_atom);
}
void ph3py_get_permutation_symmetry_compact_fc3(double * fc3,
void ph3py_get_permutation_symmetry_compact_fc3(double *fc3,
const long p2s[],
const long s2pp[],
const long nsym_list[],
@ -547,7 +541,7 @@ void ph3py_get_permutation_symmetry_compact_fc3(double * fc3,
n_patom);
}
void ph3py_transpose_compact_fc3(double * fc3,
void ph3py_transpose_compact_fc3(double *fc3,
const long p2s[],
const long s2pp[],
const long nsym_list[],
@ -566,7 +560,6 @@ void ph3py_transpose_compact_fc3(double * fc3,
t_type);
}
long ph3py_get_triplets_reciprocal_mesh_at_q(long *map_triplets,
long *map_q,
const long grid_point,
@ -586,7 +579,6 @@ long ph3py_get_triplets_reciprocal_mesh_at_q(long *map_triplets,
swappable);
}
long ph3py_get_BZ_triplets_at_q(long (*triplets)[3],
const long grid_point,
const long (*bz_grid_addresses)[3],
@ -600,7 +592,8 @@ long ph3py_get_BZ_triplets_at_q(long (*triplets)[3],
ConstBZGrid *bzgrid;
long i, j, num_ir;
if ((bzgrid = (ConstBZGrid*) malloc(sizeof(ConstBZGrid))) == NULL) {
if ((bzgrid = (ConstBZGrid *)malloc(sizeof(ConstBZGrid))) == NULL)
{
warning_print("Memory could not be allocated.");
return 0;
}
@ -608,10 +601,12 @@ long ph3py_get_BZ_triplets_at_q(long (*triplets)[3],
bzgrid->addresses = bz_grid_addresses;
bzgrid->gp_map = bz_map;
bzgrid->type = bz_grid_type;
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
bzgrid->D_diag[i] = D_diag[i];
bzgrid->PS[i] = 0;
for (j = 0; j < 3; j++) {
for (j = 0; j < 3; j++)
{
bzgrid->Q[i][j] = Q[i][j];
}
}
@ -652,7 +647,8 @@ long ph3py_get_integration_weight(double *iw,
ConstBZGrid *bzgrid;
long i;
if ((bzgrid = (ConstBZGrid*) malloc(sizeof(ConstBZGrid))) == NULL) {
if ((bzgrid = (ConstBZGrid *)malloc(sizeof(ConstBZGrid))) == NULL)
{
warning_print("Memory could not be allocated.");
return 0;
}
@ -660,7 +656,8 @@ long ph3py_get_integration_weight(double *iw,
bzgrid->addresses = bz_grid_addresses;
bzgrid->gp_map = bz_map;
bzgrid->type = bz_grid_type;
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
bzgrid->D_diag[i] = D_diag[i];
}
@ -710,7 +707,6 @@ void ph3py_get_integration_weight_with_sigma(double *iw,
tp_type);
}
/* From single address to grid index */
long ph3py_get_grid_index_from_address(const long address[3],
const long D_diag[3])
@ -718,14 +714,12 @@ long ph3py_get_grid_index_from_address(const long address[3],
return grg_get_grid_index(address, D_diag);
}
void ph3py_get_gr_grid_addresses(long gr_grid_addresses[][3],
const long D_diag[3])
{
grg_get_all_grid_addresses(gr_grid_addresses, D_diag);
}
long ph3py_get_reciprocal_rotations(long rec_rotations[48][3][3],
const long (*rotations)[3][3],
const long num_rot,
@ -738,7 +732,6 @@ long ph3py_get_reciprocal_rotations(long rec_rotations[48][3][3],
1);
}
/* Rotation matrices with respect to reciprocal basis vectors are
* transformed to those for GRGrid. This set of the rotations are
* used always in GRGrid handling. */
@ -775,8 +768,10 @@ long ph3py_get_ir_grid_map(long *ir_grid_map,
grg_get_ir_grid_map(ir_grid_map, grg_rotations, num_rot, D_diag, PS);
num_ir = 0;
for (i = 0; i < D_diag[0] * D_diag[1] * D_diag[2]; i ++) {
if (ir_grid_map[i] == i) {
for (i = 0; i < D_diag[0] * D_diag[1] * D_diag[2]; i++)
{
if (ir_grid_map[i] == i)
{
num_ir++;
}
}
@ -796,7 +791,8 @@ long ph3py_get_bz_grid_addresses(long (*bz_grid_addresses)[3],
BZGrid *bzgrid;
long i, j, size;
if ((bzgrid = (BZGrid*) malloc(sizeof(BZGrid))) == NULL) {
if ((bzgrid = (BZGrid *)malloc(sizeof(BZGrid))) == NULL)
{
warning_print("Memory could not be allocated.");
return 0;
}
@ -805,18 +801,23 @@ long ph3py_get_bz_grid_addresses(long (*bz_grid_addresses)[3],
bzgrid->gp_map = bz_map;
bzgrid->bzg2grg = bzg2grg;
bzgrid->type = type;
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
bzgrid->D_diag[i] = D_diag[i];
bzgrid->PS[i] = PS[i];
for (j = 0; j < 3; j++) {
for (j = 0; j < 3; j++)
{
bzgrid->Q[i][j] = Q[i][j];
bzgrid->reclat[i][j] = rec_lattice[i][j];
}
}
if (bzg_get_bz_grid_addresses(bzgrid)) {
if (bzg_get_bz_grid_addresses(bzgrid))
{
size = bzgrid->size;
} else {
}
else
{
size = 0;
}
@ -837,7 +838,8 @@ long ph3py_rotate_bz_grid_index(const long bz_grid_index,
ConstBZGrid *bzgrid;
long i, rot_bz_gp;
if ((bzgrid = (ConstBZGrid*) malloc(sizeof(ConstBZGrid))) == NULL) {
if ((bzgrid = (ConstBZGrid *)malloc(sizeof(ConstBZGrid))) == NULL)
{
warning_print("Memory could not be allocated.");
return 0;
}
@ -845,7 +847,8 @@ long ph3py_rotate_bz_grid_index(const long bz_grid_index,
bzgrid->addresses = bz_grid_addresses;
bzgrid->gp_map = bz_map;
bzgrid->type = bz_grid_type;
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
bzgrid->D_diag[i] = D_diag[i];
bzgrid->PS[i] = 0;
}
@ -858,7 +861,6 @@ long ph3py_rotate_bz_grid_index(const long bz_grid_index,
return rot_bz_gp;
}
void ph3py_symmetrize_collision_matrix(double *collision_matrix,
const long num_column,
const long num_temp,
@ -867,18 +869,23 @@ void ph3py_symmetrize_collision_matrix(double *collision_matrix,
double val;
long i, j, k, l, adrs_shift;
for (i = 0; i < num_sigma; i++) {
for (j = 0; j < num_temp; j++) {
for (i = 0; i < num_sigma; i++)
{
for (j = 0; j < num_temp; j++)
{
adrs_shift = (i * num_column * num_column * num_temp +
j * num_column * num_column);
/* show_colmat_info(py_collision_matrix, i, j, adrs_shift); */
#ifdef PHPYOPENMP
#pragma omp parallel for schedule(guided) private(l, val)
#endif
for (k = 0; k < num_column; k++) {
for (l = k + 1; l < num_column; l++) {
for (k = 0; k < num_column; k++)
{
for (l = k + 1; l < num_column; l++)
{
val = (collision_matrix[adrs_shift + k * num_column + l] +
collision_matrix[adrs_shift + l * num_column + k]) / 2;
collision_matrix[adrs_shift + l * num_column + k]) /
2;
collision_matrix[adrs_shift + k * num_column + l] = val;
collision_matrix[adrs_shift + l * num_column + k] = val;
}
@ -887,7 +894,6 @@ void ph3py_symmetrize_collision_matrix(double *collision_matrix,
}
}
void ph3py_expand_collision_matrix(double *collision_matrix,
const long *rot_grid_points,
const long *ir_grid_points,
@ -904,7 +910,7 @@ void ph3py_expand_collision_matrix(double *collision_matrix,
long *multi;
double *colmat_copy;
multi = (long*)malloc(sizeof(long) * num_ir_gp);
multi = (long *)malloc(sizeof(long) * num_ir_gp);
colmat_copy = NULL;
num_column = num_grid_points * num_band;
@ -913,40 +919,49 @@ void ph3py_expand_collision_matrix(double *collision_matrix,
#ifdef PHPYOPENMP
#pragma omp parallel for schedule(guided) private(j, ir_gp)
#endif
for (i = 0; i < num_ir_gp; i++) {
for (i = 0; i < num_ir_gp; i++)
{
ir_gp = ir_grid_points[i];
multi[i] = 0;
for (j = 0; j < num_rot; j++) {
if (rot_grid_points[j * num_grid_points + ir_gp] == ir_gp) {
for (j = 0; j < num_rot; j++)
{
if (rot_grid_points[j * num_grid_points + ir_gp] == ir_gp)
{
multi[i]++;
}
}
}
for (i = 0; i < num_sigma; i++) {
for (j = 0; j < num_temp; j++) {
for (i = 0; i < num_sigma; i++)
{
for (j = 0; j < num_temp; j++)
{
adrs_shift = (i * num_column * num_column * num_temp +
j * num_column * num_column);
#ifdef PHPYOPENMP
#pragma omp parallel for private(ir_gp, adrs_shift_plus, colmat_copy, l, gp_r, m, n, p)
#endif
for (k = 0; k < num_ir_gp; k++) {
for (k = 0; k < num_ir_gp; k++)
{
ir_gp = ir_grid_points[k];
adrs_shift_plus = adrs_shift + ir_gp * num_bgb;
colmat_copy = (double*)malloc(sizeof(double) * num_bgb);
for (l = 0; l < num_bgb; l++) {
colmat_copy = (double *)malloc(sizeof(double) * num_bgb);
for (l = 0; l < num_bgb; l++)
{
colmat_copy[l] = collision_matrix[adrs_shift_plus + l] / multi[k];
collision_matrix[adrs_shift_plus + l] = 0;
}
for (l = 0; l < num_rot; l++) {
for (l = 0; l < num_rot; l++)
{
gp_r = rot_grid_points[l * num_grid_points + ir_gp];
for (m = 0; m < num_band; m++) {
for (n = 0; n < num_grid_points; n++) {
for (p = 0; p < num_band; p++) {
collision_matrix[
adrs_shift + gp_r * num_bgb + m * num_grid_points * num_band
+ rot_grid_points[l * num_grid_points + n] * num_band + p] +=
colmat_copy[m * num_grid_points * num_band + n * num_band + p];
for (m = 0; m < num_band; m++)
{
for (n = 0; n < num_grid_points; n++)
{
for (p = 0; p < num_band; p++)
{
collision_matrix[adrs_shift + gp_r * num_bgb + m * num_grid_points * num_band + rot_grid_points[l * num_grid_points + n] * num_band + p] +=
colmat_copy[m * num_grid_points * num_band + n * num_band + p];
}
}
}
@ -979,7 +994,8 @@ long ph3py_get_neighboring_gird_points(long *relative_grid_points,
long i;
ConstBZGrid *bzgrid;
if ((bzgrid = (ConstBZGrid*) malloc(sizeof(ConstBZGrid))) == NULL) {
if ((bzgrid = (ConstBZGrid *)malloc(sizeof(ConstBZGrid))) == NULL)
{
warning_print("Memory could not be allocated.");
return 0;
}
@ -987,20 +1003,21 @@ long ph3py_get_neighboring_gird_points(long *relative_grid_points,
bzgrid->addresses = bz_grid_addresses;
bzgrid->gp_map = bz_map;
bzgrid->type = bz_grid_type;
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
bzgrid->D_diag[i] = D_diag[i];
}
#ifdef PHPYOPENMP
#pragma omp parallel for
#endif
for (i = 0; i < num_grid_points; i++) {
tpi_get_neighboring_grid_points
(relative_grid_points + i * num_relative_grid_address,
grid_points[i],
relative_grid_address,
num_relative_grid_address,
bzgrid);
for (i = 0; i < num_grid_points; i++)
{
tpi_get_neighboring_grid_points(relative_grid_points + i * num_relative_grid_address,
grid_points[i],
relative_grid_address,
num_relative_grid_address,
bzgrid);
}
free(bzgrid);
@ -1009,34 +1026,34 @@ long ph3py_get_neighboring_gird_points(long *relative_grid_points,
return 1;
}
/* thm_get_integration_weight at multiple grid points for using openmp
*
* relative_grid_addresses are given as P multipled with those from dataset,
* i.e.,
* np.dot(relative_grid_addresses, P.T) */
long ph3py_get_thm_integration_weights_at_grid_points(
double *iw,
const double *frequency_points,
const long num_frequency_points,
const long num_band,
const long num_gp,
const long (*relative_grid_address)[4][3],
const long D_diag[3],
const long *grid_points,
const long (*bz_grid_addresses)[3],
const long *bz_map,
const long bz_grid_type,
const double *frequencies,
const long *gp2irgp_map,
const char function)
double *iw,
const double *frequency_points,
const long num_frequency_points,
const long num_band,
const long num_gp,
const long (*relative_grid_address)[4][3],
const long D_diag[3],
const long *grid_points,
const long (*bz_grid_addresses)[3],
const long *bz_map,
const long bz_grid_type,
const double *frequencies,
const long *gp2irgp_map,
const char function)
{
long i, j, k, bi;
long vertices[24][4];
double freq_vertices[24][4];
ConstBZGrid *bzgrid;
if ((bzgrid = (ConstBZGrid*) malloc(sizeof(ConstBZGrid))) == NULL) {
if ((bzgrid = (ConstBZGrid *)malloc(sizeof(ConstBZGrid))) == NULL)
{
warning_print("Memory could not be allocated.");
return 0;
}
@ -1044,31 +1061,37 @@ long ph3py_get_thm_integration_weights_at_grid_points(
bzgrid->addresses = bz_grid_addresses;
bzgrid->gp_map = bz_map;
bzgrid->type = bz_grid_type;
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
bzgrid->D_diag[i] = D_diag[i];
}
#ifdef PHPYOPENMP
#pragma omp parallel for private(j, k, bi, vertices, freq_vertices)
#endif
for (i = 0; i < num_gp; i++) {
for (j = 0; j < 24; j++) {
for (i = 0; i < num_gp; i++)
{
for (j = 0; j < 24; j++)
{
tpi_get_neighboring_grid_points(vertices[j],
grid_points[i],
relative_grid_address[j],
4,
bzgrid);
}
for (bi = 0; bi < num_band; bi++) {
for (j = 0; j < 24; j++) {
for (k = 0; k < 4; k++) {
freq_vertices[j][k] = frequencies[
gp2irgp_map[vertices[j][k]] * num_band + bi];
for (bi = 0; bi < num_band; bi++)
{
for (j = 0; j < 24; j++)
{
for (k = 0; k < 4; k++)
{
freq_vertices[j][k] = frequencies[gp2irgp_map[vertices[j][k]] * num_band + bi];
}
}
for (j = 0; j < num_frequency_points; j++) {
for (j = 0; j < num_frequency_points; j++)
{
iw[i * num_frequency_points * num_band + j * num_band + bi] =
thm_get_integration_weight(frequency_points[j], freq_vertices, function);
thm_get_integration_weight(frequency_points[j], freq_vertices, function);
}
}
}

198
c/phono3py.h
View File

@ -70,7 +70,7 @@ long ph3py_get_pp_collision(double *imag_self_energy,
const long num_triplets,
const long *triplet_weights,
const long (*bz_grid_addresses)[3], /* thm */
const long *bz_map, /* thm */
const long *bz_map, /* thm */
const long bz_grid_type,
const long D_diag[3],
const long Q[3][3],
@ -88,55 +88,55 @@ long ph3py_get_pp_collision(double *imag_self_energy,
const long symmetrize_fc3_q,
const double cutoff_frequency);
long ph3py_get_pp_collision_with_sigma(
double *imag_self_energy,
const double sigma,
const double sigma_cutoff,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const long (*triplets)[3],
const long num_triplets,
const long *triplet_weights,
const long (*bz_grid_addresses)[3],
const long D_diag[3],
const long Q[3][3],
const double *fc3,
const long is_compact_fc3,
const double (*svecs)[3],
const long multi_dims[2],
const long (*multi)[2],
const double *masses,
const long *p2s_map,
const long *s2p_map,
const Larray *band_indices,
const Darray *temperatures,
const long is_NU,
const long symmetrize_fc3_q,
const double cutoff_frequency);
double *imag_self_energy,
const double sigma,
const double sigma_cutoff,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const long (*triplets)[3],
const long num_triplets,
const long *triplet_weights,
const long (*bz_grid_addresses)[3],
const long D_diag[3],
const long Q[3][3],
const double *fc3,
const long is_compact_fc3,
const double (*svecs)[3],
const long multi_dims[2],
const long (*multi)[2],
const double *masses,
const long *p2s_map,
const long *s2p_map,
const Larray *band_indices,
const Darray *temperatures,
const long is_NU,
const long symmetrize_fc3_q,
const double cutoff_frequency);
void ph3py_get_imag_self_energy_at_bands_with_g(
double *imag_self_energy,
const Darray *fc3_normal_squared,
const double *frequencies,
const long (*triplets)[3],
const long *triplet_weights,
const double *g,
const char *g_zero,
const double temperature,
const double cutoff_frequency,
const long num_frequency_points,
const long frequency_point_index);
double *imag_self_energy,
const Darray *fc3_normal_squared,
const double *frequencies,
const long (*triplets)[3],
const long *triplet_weights,
const double *g,
const char *g_zero,
const double temperature,
const double cutoff_frequency,
const long num_frequency_points,
const long frequency_point_index);
void ph3py_get_detailed_imag_self_energy_at_bands_with_g(
double *detailed_imag_self_energy,
double *imag_self_energy_N,
double *imag_self_energy_U,
const Darray *fc3_normal_squared,
const double *frequencies,
const long (*triplets)[3],
const long *triplet_weights,
const long (*bz_grid_addresses)[3],
const double *g,
const char *g_zero,
const double temperature,
const double cutoff_frequency);
double *detailed_imag_self_energy,
double *imag_self_energy_N,
double *imag_self_energy_U,
const Darray *fc3_normal_squared,
const double *frequencies,
const long (*triplets)[3],
const long *triplet_weights,
const long (*bz_grid_addresses)[3],
const double *g,
const char *g_zero,
const double temperature,
const double cutoff_frequency);
void ph3py_get_real_self_energy_at_bands(double *real_self_energy,
const Darray *fc3_normal_squared,
const long *band_indices,
@ -148,17 +148,17 @@ void ph3py_get_real_self_energy_at_bands(double *real_self_energy,
const double unit_conversion_factor,
const double cutoff_frequency);
void ph3py_get_real_self_energy_at_frequency_point(
double *real_self_energy,
const double frequency_point,
const Darray *fc3_normal_squared,
const long *band_indices,
const double *frequencies,
const long (*triplets)[3],
const long *triplet_weights,
const double epsilon,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency);
double *real_self_energy,
const double frequency_point,
const Darray *fc3_normal_squared,
const long *band_indices,
const double *frequencies,
const long (*triplets)[3],
const long *triplet_weights,
const double epsilon,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency);
void ph3py_get_collision_matrix(double *collision_matrix,
const Darray *fc3_normal_squared,
const double *frequencies,
@ -186,31 +186,31 @@ void ph3py_get_reducible_collision_matrix(double *collision_matrix,
const double unit_conversion_factor,
const double cutoff_frequency);
void ph3py_get_isotope_scattering_strength(
double *gamma,
const long grid_point,
const double *mass_variances,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const long num_grid_points,
const long *band_indices,
const long num_band,
const long num_band0,
const double sigma,
const double cutoff_frequency);
double *gamma,
const long grid_point,
const double *mass_variances,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const long num_grid_points,
const long *band_indices,
const long num_band,
const long num_band0,
const double sigma,
const double cutoff_frequency);
void ph3py_get_thm_isotope_scattering_strength(
double *gamma,
const long grid_point,
const long *ir_grid_points,
const long *weights,
const double *mass_variances,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const long num_ir_grid_points,
const long *band_indices,
const long num_band,
const long num_band0,
const double *integration_weights,
const double cutoff_frequency);
double *gamma,
const long grid_point,
const long *ir_grid_points,
const long *weights,
const double *mass_variances,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const long num_ir_grid_points,
const long *band_indices,
const long num_band,
const long num_band0,
const double *integration_weights,
const double cutoff_frequency);
void ph3py_distribute_fc3(double *fc3,
const long target,
const long source,
@ -226,14 +226,14 @@ void ph3py_rotate_delta_fc2(double (*fc3)[3][3][3],
const long num_site_sym,
const long num_disp);
void ph3py_get_permutation_symmetry_fc3(double *fc3, const long num_atom);
void ph3py_get_permutation_symmetry_compact_fc3(double * fc3,
void ph3py_get_permutation_symmetry_compact_fc3(double *fc3,
const long p2s[],
const long s2pp[],
const long nsym_list[],
const long perms[],
const long n_satom,
const long n_patom);
void ph3py_transpose_compact_fc3(double * fc3,
void ph3py_transpose_compact_fc3(double *fc3,
const long p2s[],
const long s2pp[],
const long nsym_list[],
@ -352,19 +352,19 @@ long ph3py_get_neighboring_gird_points(long *relative_grid_points,
const long num_grid_points,
const long num_relative_grid_address);
long ph3py_get_thm_integration_weights_at_grid_points(
double *iw,
const double *frequency_points,
const long num_band0,
const long num_band,
const long num_gp,
const long (*relative_grid_address)[4][3],
const long D_diag[3],
const long *grid_points,
const long (*bz_grid_addresses)[3],
const long *bz_map,
const long bz_grid_type,
const double *frequencies,
const long *gp2irgp_map,
const char function);
double *iw,
const double *frequency_points,
const long num_band0,
const long num_band,
const long num_gp,
const long (*relative_grid_address)[4][3],
const long D_diag[3],
const long *grid_points,
const long (*bz_grid_addresses)[3],
const long *bz_map,
const long bz_grid_type,
const double *frequencies,
const long *gp2irgp_map,
const char function);
#endif

280
c/phonon.c
View File

@ -50,7 +50,7 @@ static void get_undone_phonons(double *frequencies,
const long (*grid_address)[3],
const double QDinv[3][3],
const double *fc2,
const double(*svecs_fc2)[3],
const double (*svecs_fc2)[3],
const long (*multi_fc2)[2],
const long num_patom,
const long num_satom,
@ -71,7 +71,7 @@ static void get_gonze_undone_phonons(double *frequencies,
const long (*grid_address)[3],
const double QDinv[3][3],
const double *fc2,
const double(*svecs_fc2)[3],
const double (*svecs_fc2)[3],
const long (*multi_fc2)[2],
const double (*positions)[3],
const long num_patom,
@ -86,7 +86,7 @@ static void get_gonze_undone_phonons(double *frequencies,
const double *q_direction,
const double nac_factor,
const double *dd_q0,
const double(*G_list)[3],
const double (*G_list)[3],
const long num_G_points,
const double lambda,
const char uplo);
@ -99,7 +99,7 @@ static void get_phonons(lapack_complex_double *eigvecs,
const long (*multi)[2],
const long num_patom,
const long num_satom,
const double(*svecs)[3],
const double (*svecs)[3],
const long is_nac,
const double (*born)[3][3],
const double dielectric[3][3],
@ -117,7 +117,7 @@ static void get_gonze_phonons(lapack_complex_double *eigvecs,
const double (*positions)[3],
const long num_patom,
const long num_satom,
const double(*svecs)[3],
const double (*svecs)[3],
const long is_nac,
const double (*born)[3][3],
const double dielectric[3][3],
@ -125,7 +125,7 @@ static void get_gonze_phonons(lapack_complex_double *eigvecs,
const double *q_direction,
const double nac_factor,
const double *dd_q0,
const double(*G_list)[3],
const double (*G_list)[3],
const long num_G_points,
const double lambda);
static void
@ -138,7 +138,7 @@ get_dynamical_matrix(lapack_complex_double *dynmat,
const long (*multi)[2],
const long num_patom,
const long num_satom,
const double(*svecs)[3],
const double (*svecs)[3],
const long is_nac,
const double (*born)[3][3], /* Wang NAC unless NULL */
const double dielectric[3][3],
@ -159,35 +159,34 @@ static long needs_nac(const double (*born)[3][3],
const long gp,
const double *q_direction);
void
phn_get_phonons_at_gridpoints(double *frequencies,
lapack_complex_double *eigenvectors,
char *phonon_done,
const long num_phonons,
const long *grid_points,
const long num_grid_points,
const long (*grid_address)[3],
const double QDinv[3][3],
const double *fc2,
const double(*svecs_fc2)[3],
const long (*multi_fc2)[2],
const long num_patom,
const long num_satom,
const double *masses_fc2,
const long *p2s_fc2,
const long *s2p_fc2,
const double unit_conversion_factor,
const double (*born)[3][3],
const double dielectric[3][3],
const double reciprocal_lattice[3][3],
const double *q_direction, /* must be pointer */
const double nac_factor,
const char uplo)
void phn_get_phonons_at_gridpoints(double *frequencies,
lapack_complex_double *eigenvectors,
char *phonon_done,
const long num_phonons,
const long *grid_points,
const long num_grid_points,
const long (*grid_address)[3],
const double QDinv[3][3],
const double *fc2,
const double (*svecs_fc2)[3],
const long (*multi_fc2)[2],
const long num_patom,
const long num_satom,
const double *masses_fc2,
const long *p2s_fc2,
const long *s2p_fc2,
const double unit_conversion_factor,
const double (*born)[3][3],
const double dielectric[3][3],
const double reciprocal_lattice[3][3],
const double *q_direction, /* must be pointer */
const double nac_factor,
const char uplo)
{
long num_undone;
long *undone;
undone = (long*)malloc(sizeof(long) * num_phonons);
undone = (long *)malloc(sizeof(long) * num_phonons);
num_undone = collect_undone_grid_points(undone,
phonon_done,
num_grid_points,
@ -219,40 +218,39 @@ phn_get_phonons_at_gridpoints(double *frequencies,
undone = NULL;
}
void
phn_get_gonze_phonons_at_gridpoints(double *frequencies,
lapack_complex_double *eigenvectors,
char *phonon_done,
const long num_phonons,
const long *grid_points,
const long num_grid_points,
const long (*grid_address)[3],
const double QDinv[3][3],
const double *fc2,
const double(*svecs_fc2)[3],
const long (*multi_fc2)[2],
const double (*positions)[3],
const long num_patom,
const long num_satom,
const double *masses_fc2,
const long *p2s_fc2,
const long *s2p_fc2,
const double unit_conversion_factor,
const double (*born)[3][3],
const double dielectric[3][3],
const double reciprocal_lattice[3][3],
const double *q_direction, /* pointer */
const double nac_factor,
const double *dd_q0,
const double(*G_list)[3],
const long num_G_points,
const double lambda,
const char uplo)
void phn_get_gonze_phonons_at_gridpoints(double *frequencies,
lapack_complex_double *eigenvectors,
char *phonon_done,
const long num_phonons,
const long *grid_points,
const long num_grid_points,
const long (*grid_address)[3],
const double QDinv[3][3],
const double *fc2,
const double (*svecs_fc2)[3],
const long (*multi_fc2)[2],
const double (*positions)[3],
const long num_patom,
const long num_satom,
const double *masses_fc2,
const long *p2s_fc2,
const long *s2p_fc2,
const double unit_conversion_factor,
const double (*born)[3][3],
const double dielectric[3][3],
const double reciprocal_lattice[3][3],
const double *q_direction, /* pointer */
const double nac_factor,
const double *dd_q0,
const double (*G_list)[3],
const long num_G_points,
const double lambda,
const char uplo)
{
long num_undone;
long *undone;
undone = (long*)malloc(sizeof(long) * num_phonons);
undone = (long *)malloc(sizeof(long) * num_phonons);
num_undone = collect_undone_grid_points(undone,
phonon_done,
num_grid_points,
@ -297,9 +295,11 @@ static long collect_undone_grid_points(long *undone,
long i, gp, num_undone;
num_undone = 0;
for (i = 0; i < num_grid_points; i++) {
for (i = 0; i < num_grid_points; i++)
{
gp = grid_points[i];
if (phonon_done[gp] == 0) {
if (phonon_done[gp] == 0)
{
undone[num_undone] = gp;
num_undone++;
phonon_done[gp] = 1;
@ -316,7 +316,7 @@ static void get_undone_phonons(double *frequencies,
const long (*grid_address)[3],
const double QDinv[3][3],
const double *fc2,
const double(*svecs_fc2)[3],
const double (*svecs_fc2)[3],
const long (*multi_fc2)[2],
const long num_patom,
const long num_satom,
@ -341,12 +341,12 @@ static void get_undone_phonons(double *frequencies,
#ifdef PHPYOPENMP
#pragma omp parallel for private(j, q, gp, is_nac)
#endif
for (i = 0; i < num_undone_grid_points; i++) {
for (i = 0; i < num_undone_grid_points; i++)
{
gp = undone_grid_points[i];
for (j = 0; j < 3; j++) {
q[j] = QDinv[j][0] * grid_address[gp][0]
+ QDinv[j][1] * grid_address[gp][1]
+ QDinv[j][2] * grid_address[gp][2];
for (j = 0; j < 3; j++)
{
q[j] = QDinv[j][0] * grid_address[gp][0] + QDinv[j][1] * grid_address[gp][1] + QDinv[j][2] * grid_address[gp][2];
}
is_nac = needs_nac(born, grid_address, gp, q_direction);
@ -375,7 +375,8 @@ static void get_undone_phonons(double *frequencies,
#pragma omp parallel for private(j, gp, freqs_tmp, info)
#endif
#endif
for (i = 0; i < num_undone_grid_points; i++) {
for (i = 0; i < num_undone_grid_points; i++)
{
gp = undone_grid_points[i];
freqs_tmp = frequencies + num_band * gp;
/* Store eigenvalues in freqs array. */
@ -386,9 +387,10 @@ static void get_undone_phonons(double *frequencies,
uplo);
/* Sqrt of eigenvalues are re-stored in freqs array.*/
for (j = 0; j < num_band; j++) {
for (j = 0; j < num_band; j++)
{
freqs_tmp[j] = sqrt(fabs(freqs_tmp[j])) *
((freqs_tmp[j] > 0) - (freqs_tmp[j] < 0)) * unit_conversion_factor;
((freqs_tmp[j] > 0) - (freqs_tmp[j] < 0)) * unit_conversion_factor;
}
}
}
@ -400,7 +402,7 @@ static void get_gonze_undone_phonons(double *frequencies,
const long (*grid_address)[3],
const double QDinv[3][3],
const double *fc2,
const double(*svecs_fc2)[3],
const double (*svecs_fc2)[3],
const long (*multi_fc2)[2],
const double (*positions)[3],
const long num_patom,
@ -415,7 +417,7 @@ static void get_gonze_undone_phonons(double *frequencies,
const double *q_direction,
const double nac_factor,
const double *dd_q0,
const double(*G_list)[3],
const double (*G_list)[3],
const long num_G_points,
const double lambda,
const char uplo)
@ -430,12 +432,12 @@ static void get_gonze_undone_phonons(double *frequencies,
#ifdef PHPYOPENMP
#pragma omp parallel for private(j, q, gp, is_nac)
#endif
for (i = 0; i < num_undone_grid_points; i++) {
for (i = 0; i < num_undone_grid_points; i++)
{
gp = undone_grid_points[i];
for (j = 0; j < 3; j++) {
q[j] = QDinv[j][0] * grid_address[gp][0]
+ QDinv[j][1] * grid_address[gp][1]
+ QDinv[j][2] * grid_address[gp][2];
for (j = 0; j < 3; j++)
{
q[j] = QDinv[j][0] * grid_address[gp][0] + QDinv[j][1] * grid_address[gp][1] + QDinv[j][2] * grid_address[gp][2];
}
is_nac = needs_nac(born, grid_address, gp, q_direction);
get_gonze_phonons(eigenvectors + num_band * num_band * gp,
@ -461,14 +463,14 @@ static void get_gonze_undone_phonons(double *frequencies,
lambda);
}
/* To avoid multithreaded BLAS in OpenMP loop */
#ifdef PHPYOPENMP
#ifndef MULTITHREADED_BLAS
#pragma omp parallel for private(j, gp, freqs_tmp, info)
#endif
#endif
for (i = 0; i < num_undone_grid_points; i++) {
for (i = 0; i < num_undone_grid_points; i++)
{
gp = undone_grid_points[i];
/* Store eigenvalues in freqs array. */
/* Eigenvectors are overwritten on eigvecs array. */
@ -479,9 +481,10 @@ static void get_gonze_undone_phonons(double *frequencies,
uplo);
/* Sqrt of eigenvalues are re-stored in freqs array.*/
for (j = 0; j < num_band; j++) {
for (j = 0; j < num_band; j++)
{
freqs_tmp[j] = sqrt(fabs(freqs_tmp[j])) *
((freqs_tmp[j] > 0) - (freqs_tmp[j] < 0)) * unit_conversion_factor;
((freqs_tmp[j] > 0) - (freqs_tmp[j] < 0)) * unit_conversion_factor;
}
}
}
@ -495,7 +498,7 @@ static void get_phonons(lapack_complex_double *eigvecs,
const long (*multi)[2],
const long num_patom,
const long num_satom,
const double(*svecs)[3],
const double (*svecs)[3],
const long is_nac,
const double (*born)[3][3],
const double dielectric[3][3],
@ -533,7 +536,7 @@ static void get_gonze_phonons(lapack_complex_double *eigvecs,
const double (*positions)[3],
const long num_patom,
const long num_satom,
const double(*svecs)[3],
const double (*svecs)[3],
const long is_nac,
const double (*born)[3][3],
const double dielectric[3][3],
@ -541,7 +544,7 @@ static void get_gonze_phonons(lapack_complex_double *eigvecs,
const double *q_direction,
const double nac_factor,
const double *dd_q0,
const double(*G_list)[3],
const double (*G_list)[3],
const long num_G_points,
const double lambda)
{
@ -555,7 +558,7 @@ static void get_gonze_phonons(lapack_complex_double *eigvecs,
q_dir_cart = NULL;
num_band = num_patom * 3;
dym_get_dynamical_matrix_at_q((double*)eigvecs,
dym_get_dynamical_matrix_at_q((double *)eigvecs,
num_patom,
num_satom,
fc2,
@ -568,26 +571,31 @@ static void get_gonze_phonons(lapack_complex_double *eigvecs,
NULL,
0);
dd = (lapack_complex_double*)
malloc(sizeof(lapack_complex_double) * num_band * num_band);
for (i = 0; i < 3; i++) {
dd = (lapack_complex_double *)
malloc(sizeof(lapack_complex_double) * num_band * num_band);
for (i = 0; i < 3; i++)
{
q_cart[i] = 0;
for (j = 0; j < 3; j++) {
for (j = 0; j < 3; j++)
{
q_cart[i] += reciprocal_lattice[i][j] * q[j];
}
}
if (q_direction) {
q_dir_cart = (double*)malloc(sizeof(double) * 3);
for (i = 0; i < 3; i++) {
if (q_direction)
{
q_dir_cart = (double *)malloc(sizeof(double) * 3);
for (i = 0; i < 3; i++)
{
q_dir_cart[i] = 0;
for (j = 0; j < 3; j++) {
for (j = 0; j < 3; j++)
{
q_dir_cart[i] += reciprocal_lattice[i][j] * q_direction[j];
}
}
}
dym_get_recip_dipole_dipole((double*)dd,
dym_get_recip_dipole_dipole((double *)dd,
dd_q0,
G_list,
num_G_points,
@ -601,22 +609,27 @@ static void get_gonze_phonons(lapack_complex_double *eigvecs,
lambda,
1e-5);
if (q_direction) {
if (q_direction)
{
free(q_dir_cart);
q_dir_cart = NULL;
}
for (i = 0; i < num_patom; i++) {
for (j = 0; j < num_patom; j++) {
for (i = 0; i < num_patom; i++)
{
for (j = 0; j < num_patom; j++)
{
mm = sqrt(masses[i] * masses[j]);
for (k = 0; k < 3; k++) {
for (l = 0; l < 3; l++) {
for (k = 0; k < 3; k++)
{
for (l = 0; l < 3; l++)
{
adrs = i * num_patom * 9 + k * num_patom * 3 + j * 3 + l;
eigvecs[adrs] = lapack_make_complex_double(
lapack_complex_double_real(eigvecs[adrs]) +
lapack_complex_double_real(dd[adrs]) / mm,
lapack_complex_double_imag(eigvecs[adrs]) +
lapack_complex_double_imag(dd[adrs]) / mm);
lapack_complex_double_real(eigvecs[adrs]) +
lapack_complex_double_real(dd[adrs]) / mm,
lapack_complex_double_imag(eigvecs[adrs]) +
lapack_complex_double_imag(dd[adrs]) / mm);
}
}
}
@ -636,7 +649,7 @@ get_dynamical_matrix(lapack_complex_double *dynmat,
const long (*multi)[2],
const long num_patom,
const long num_satom,
const double(*svecs)[3],
const double (*svecs)[3],
const long is_nac,
const double (*born)[3][3], /* Wang NAC unless NULL */
const double dielectric[3][3],
@ -644,13 +657,14 @@ get_dynamical_matrix(lapack_complex_double *dynmat,
const double *q_direction,
const double nac_factor)
{
double (*charge_sum)[3][3];
double(*charge_sum)[3][3];
charge_sum = NULL;
if (is_nac) {
if (is_nac)
{
charge_sum = (double(*)[3][3])
malloc(sizeof(double[3][3]) * num_patom * num_patom * 9);
malloc(sizeof(double[3][3]) * num_patom * num_patom * 9);
get_charge_sum(charge_sum,
num_patom,
num_satom,
@ -662,7 +676,7 @@ get_dynamical_matrix(lapack_complex_double *dynmat,
nac_factor);
}
dym_get_dynamical_matrix_at_q((double*)dynmat,
dym_get_dynamical_matrix_at_q((double *)dynmat,
num_patom,
num_satom,
fc2,
@ -674,7 +688,8 @@ get_dynamical_matrix(lapack_complex_double *dynmat,
p2s,
charge_sum,
0);
if (is_nac) {
if (is_nac)
{
free(charge_sum);
charge_sum = NULL;
}
@ -694,26 +709,35 @@ static void get_charge_sum(double (*charge_sum)[3][3],
double inv_dielectric_factor, dielectric_factor, tmp_val;
double q_cart[3];
if (q_direction) {
for (i = 0; i < 3; i++) {
if (q_direction)
{
for (i = 0; i < 3; i++)
{
q_cart[i] = 0.0;
for (j = 0; j < 3; j++) {
for (j = 0; j < 3; j++)
{
q_cart[i] += reciprocal_lattice[i][j] * q_direction[j];
}
}
} else {
for (i = 0; i < 3; i++) {
}
else
{
for (i = 0; i < 3; i++)
{
q_cart[i] = 0.0;
for (j = 0; j < 3; j++) {
for (j = 0; j < 3; j++)
{
q_cart[i] += reciprocal_lattice[i][j] * q[j];
}
}
}
inv_dielectric_factor = 0.0;
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
tmp_val = 0.0;
for (j = 0; j < 3; j++) {
for (j = 0; j < 3; j++)
{
tmp_val += dielectric[i][j] * q_cart[j];
}
inv_dielectric_factor += tmp_val * q_cart[i];
@ -721,7 +745,7 @@ static void get_charge_sum(double (*charge_sum)[3][3],
/* N = num_satom / num_patom = number of prim-cell in supercell */
/* N is used for Wang's method. */
dielectric_factor = nac_factor /
inv_dielectric_factor / num_satom * num_patom;
inv_dielectric_factor / num_satom * num_patom;
dym_get_charge_sum(charge_sum,
num_patom,
dielectric_factor,
@ -736,16 +760,22 @@ static long needs_nac(const double (*born)[3][3],
{
long is_nac;
if (born) {
if (born)
{
if (grid_address[gp][0] == 0 &&
grid_address[gp][1] == 0 &&
grid_address[gp][2] == 0 &&
q_direction == NULL) {
q_direction == NULL)
{
is_nac = 0;
} else {
}
else
{
is_nac = 1;
}
} else {
}
else
{
is_nac = 0;
}

104
c/phonon.h
View File

@ -38,57 +38,55 @@
#include "dynmat.h"
#include "lapack_wrapper.h"
void
phn_get_phonons_at_gridpoints(double *frequencies,
lapack_complex_double *eigenvectors,
char *phonon_done,
const long num_phonons,
const long *grid_points,
const long num_grid_points,
const long (*grid_address)[3],
const double QDinv[3][3],
const double *fc2,
const double(*svecs_fc2)[3],
const long (*multi_fc2)[2],
const long num_patom,
const long num_satom,
const double *masses_fc2,
const long *p2s_fc2,
const long *s2p_fc2,
const double unit_conversion_factor,
const double (*born)[3][3],
const double dielectric[3][3],
const double reciprocal_lattice[3][3],
const double *q_direction, /* must be pointer */
const double nac_factor,
const char uplo);
void
phn_get_gonze_phonons_at_gridpoints(double *frequencies,
lapack_complex_double *eigenvectors,
char *phonon_done,
const long num_phonons,
const long *grid_points,
const long num_grid_points,
const long (*grid_address)[3],
const double QDinv[3][3],
const double *fc2,
const double(*svecs_fc2)[3],
const long (*multi_fc2)[2],
const double (*positions)[3],
const long num_patom,
const long num_satom,
const double *masses_fc2,
const long *p2s_fc2,
const long *s2p_fc2,
const double unit_conversion_factor,
const double (*born)[3][3],
const double dielectric[3][3],
const double reciprocal_lattice[3][3],
const double *q_direction, /* pointer */
const double nac_factor,
const double *dd_q0,
const double(*G_list)[3],
const long num_G_points,
const double lambda,
const char uplo);
void phn_get_phonons_at_gridpoints(double *frequencies,
lapack_complex_double *eigenvectors,
char *phonon_done,
const long num_phonons,
const long *grid_points,
const long num_grid_points,
const long (*grid_address)[3],
const double QDinv[3][3],
const double *fc2,
const double (*svecs_fc2)[3],
const long (*multi_fc2)[2],
const long num_patom,
const long num_satom,
const double *masses_fc2,
const long *p2s_fc2,
const long *s2p_fc2,
const double unit_conversion_factor,
const double (*born)[3][3],
const double dielectric[3][3],
const double reciprocal_lattice[3][3],
const double *q_direction, /* must be pointer */
const double nac_factor,
const char uplo);
void phn_get_gonze_phonons_at_gridpoints(double *frequencies,
lapack_complex_double *eigenvectors,
char *phonon_done,
const long num_phonons,
const long *grid_points,
const long num_grid_points,
const long (*grid_address)[3],
const double QDinv[3][3],
const double *fc2,
const double (*svecs_fc2)[3],
const long (*multi_fc2)[2],
const double (*positions)[3],
const long num_patom,
const long num_satom,
const double *masses_fc2,
const long *p2s_fc2,
const long *s2p_fc2,
const double unit_conversion_factor,
const double (*born)[3][3],
const double dielectric[3][3],
const double reciprocal_lattice[3][3],
const double *q_direction, /* pointer */
const double nac_factor,
const double *dd_q0,
const double (*G_list)[3],
const long num_G_points,
const double lambda,
const char uplo);
#endif

11
c/phononmod.c
View File

@ -45,7 +45,7 @@ void phmod_get_phonons_at_gridpoints(double *frequencies,
const long (*grid_address)[3],
const double QDinv[3][3],
const double *fc2,
const double(*svecs_fc2)[3],
const double (*svecs_fc2)[3],
const long (*multi_fc2)[2],
const double (*positions_fc2)[3],
const long num_patom,
@ -60,12 +60,13 @@ void phmod_get_phonons_at_gridpoints(double *frequencies,
const double *q_direction, /* pointer */
const double nac_factor,
const double *dd_q0,
const double(*G_list)[3],
const double (*G_list)[3],
const long num_G_points,
const double lambda,
const char uplo)
{
if (!dd_q0) {
if (!dd_q0)
{
phn_get_phonons_at_gridpoints(frequencies,
eigenvectors,
phonon_done,
@ -89,7 +90,9 @@ void phmod_get_phonons_at_gridpoints(double *frequencies,
q_direction,
nac_factor,
uplo);
} else {
}
else
{
phn_get_gonze_phonons_at_gridpoints(frequencies,
eigenvectors,
phonon_done,

4
c/phononmod.h
View File

@ -46,7 +46,7 @@ void phmod_get_phonons_at_gridpoints(double *frequencies,
const long (*grid_address)[3],
const double QDinv[3][3],
const double *fc2,
const double(*svecs_fc2)[3],
const double (*svecs_fc2)[3],
const long (*multi_fc2)[2],
const double (*positions_fc2)[3],
const long num_patom,
@ -61,7 +61,7 @@ void phmod_get_phonons_at_gridpoints(double *frequencies,
const double *q_direction, /* pointer */
const double nac_factor,
const double *dd_q0,
const double(*G_list)[3],
const double (*G_list)[3],
const long num_G_points,
const double lambda,
const char uplo);

173
c/real_self_energy.c
View File

@ -84,59 +84,67 @@ void rse_get_real_self_energy_at_bands(double *real_self_energy,
gp0 = triplets[0][0];
/* num_band0 and num_band_indices have to be same. */
for (i = 0; i < num_band0; i++) {
for (i = 0; i < num_band0; i++)
{
fpoint = frequencies[gp0 * num_band + band_indices[i]];
if (fpoint < cutoff_frequency) {
if (fpoint < cutoff_frequency)
{
real_self_energy[i] = 0;
} else {
}
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);
get_real_self_energy_at_band(i,
fc3_normal_squared,
fpoint,
frequencies,
triplets,
triplet_weights,
epsilon,
temperature,
unit_conversion_factor,
cutoff_frequency);
}
}
}
void rse_get_real_self_energy_at_frequency_point(
double *real_self_energy,
const double frequency_point,
const Darray *fc3_normal_squared,
const long *band_indices,
const double *frequencies,
const long (*triplets)[3],
const long *triplet_weights,
const double epsilon,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency)
double *real_self_energy,
const double frequency_point,
const Darray *fc3_normal_squared,
const long *band_indices,
const double *frequencies,
const long (*triplets)[3],
const long *triplet_weights,
const double epsilon,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency)
{
long 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) {
for (i = 0; i < num_band0; i++)
{
if (frequency_point < cutoff_frequency)
{
real_self_energy[i] = 0;
} else {
}
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);
get_real_self_energy_at_band(i,
fc3_normal_squared,
frequency_point,
frequencies,
triplets,
triplet_weights,
epsilon,
temperature,
unit_conversion_factor,
cutoff_frequency);
}
}
}
@ -161,36 +169,41 @@ static double get_real_self_energy_at_band(const long band_index,
shift = 0;
#ifdef PHPYOPENMP
#pragma omp parallel for private(gp1, gp2) reduction(+:shift)
#pragma omp parallel for private(gp1, gp2) reduction(+ \
: shift)
#endif
for (i = 0; i < num_triplets; i++) {
for (i = 0; i < num_triplets; i++)
{
gp1 = triplets[i][1];
gp2 = triplets[i][2];
if (temperature > 0) {
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 {
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;
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;
@ -210,11 +223,15 @@ static double sum_real_self_energy_at_band(const long num_band,
/* double sum; */
shift = 0;
for (i = 0; i < num_band; i++) {
if (freqs1[i] > cutoff_frequency) {
for (i = 0; i < num_band; i++)
{
if (freqs1[i] > cutoff_frequency)
{
n1 = phonoc_bose_einstein(freqs1[i], temperature);
for (j = 0; j < num_band; j++) {
if (freqs2[j] > cutoff_frequency) {
for (j = 0; j < num_band; j++)
{
if (freqs2[j] > cutoff_frequency)
{
n2 = phonoc_bose_einstein(freqs2[j], temperature);
f1 = fpoint + freqs1[i] + freqs2[j];
f2 = fpoint - freqs1[i] - freqs2[j];
@ -236,11 +253,8 @@ static double sum_real_self_energy_at_band(const long num_band,
* }
* 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];
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];
}
}
}
@ -260,15 +274,18 @@ static double sum_real_self_energy_at_band_0K(const long num_band,
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) {
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];
shift += (-1 * f1 / (f1 * f1 + epsilon * epsilon) + 1 * f2 / (f2 * f2 + epsilon * epsilon)) *
fc3_normal_squared[i * num_band + j];
}
}
}

22
c/real_self_energy.h
View File

@ -50,14 +50,14 @@ void rse_get_real_self_energy_at_bands(double *real_self_energy,
const double unit_conversion_factor,
const double cutoff_frequency);
void rse_get_real_self_energy_at_frequency_point(
double *real_self_energy,
const double frequency_point,
const Darray *fc3_normal_squared,
const long *band_indices,
const double *frequencies,
const long (*triplets)[3],
const long *triplet_weights,
const double epsilon,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency);
double *real_self_energy,
const double frequency_point,
const Darray *fc3_normal_squared,
const long *band_indices,
const double *frequencies,
const long (*triplets)[3],
const long *triplet_weights,
const double epsilon,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency);

139
c/snf3x3.c
View File

@ -38,7 +38,6 @@ static void matmul(long m[3][3],
SNF3x3CONST long b[3][3]);
static long det(SNF3x3CONST long m[3][3]);
/* static void test_set_A(long A[3][3]);
* static void test_show_A(SNF3x3CONST long A[3][3]);
* static void test_extended_gcd(void);
@ -57,9 +56,12 @@ int snf3x3(long A[3][3], long P[3][3], long Q[3][3])
int i;
initialize_PQ(P, Q);
for (i = 0; i < 100; i++) {
if (first(A, P, Q)) {
if (second(A, P, Q)) {
for (i = 0; i < 100; i++)
{
if (first(A, P, Q))
{
if (second(A, P, Q))
{
finalize(A, P, Q);
transpose(Q);
goto succeeded;
@ -76,12 +78,17 @@ static void initialize_PQ(long P[3][3], long Q[3][3])
{
int i, j;
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
if (i == j) {
for (i = 0; i < 3; i++)
{
for (j = 0; j < 3; j++)
{
if (i == j)
{
P[i][j] = 1;
Q[i][j] = 1;
} else {
}
else
{
P[i][j] = 0;
Q[i][j] = 0;
}
@ -96,14 +103,16 @@ static int first(long A[3][3], long P[3][3], long Q[3][3])
first_one_loop(A, P, Q);
/* rows and columns are all zero except for the pivot */
if ((A[1][0] == 0) && (A[2][0] == 0)) {
if ((A[1][0] == 0) && (A[2][0] == 0))
{
return 1;
}
/* columns of the pivot are assumed zero because of first_one_loop. */
/* rows of the pivot are non-zero, but divisible by the pivot. */
/* first_finalize makes the rows be zero. */
if ((A[1][0] % A[0][0] == 0) && (A[2][0] % A[0][0] == 0)) {
if ((A[1][0] % A[0][0] == 0) && (A[2][0] % A[0][0] == 0))
{
first_finalize(L, A);
matmul(A, L, A);
matmul(P, L, P);
@ -126,28 +135,31 @@ static void first_column(long A[3][3], long P[3][3])
long L[3][3];
i = search_first_pivot(A);
if (i > 0) {
if (i > 0)
{
swap_rows(L, 0, i);
matmul(A, L, A);
matmul(P, L, P);
}
if (i < 0) {
if (i < 0)
{
goto err;
}
if (A[1][0] != 0) {
if (A[1][0] != 0)
{
zero_first_column(L, 1, A);
matmul(A, L, A);
matmul(P, L, P);
}
if (A[2][0] != 0) {
if (A[2][0] != 0)
{
zero_first_column(L, 2, A);
matmul(A, L, A);
matmul(P, L, P);
}
err:
;
err:;
}
static void zero_first_column(long L[3][3], const int j,
@ -163,8 +175,10 @@ static int search_first_pivot(SNF3x3CONST long A[3][3])
{
int i;
for (i = 0; i < 3; i++) {
if (A[i][0] != 0) {
for (i = 0; i < 3; i++)
{
if (A[i][0] != 0)
{
return i;
}
}
@ -190,11 +204,13 @@ static int second(long A[3][3], long P[3][3], long Q[3][3])
second_one_loop(A, P, Q);
if (A[2][1] == 0) {
if (A[2][1] == 0)
{
return 1;
}
if (A[2][1] % A[1][1] == 0) {
if (A[2][1] % A[1][1] == 0)
{
second_finalize(L, A);
matmul(A, L, A);
matmul(P, L, P);
@ -216,13 +232,15 @@ static void second_column(long A[3][3], long P[3][3])
{
long L[3][3];
if ((A[1][1] == 0) && (A[2][1] != 0)) {
if ((A[1][1] == 0) && (A[2][1] != 0))
{
swap_rows(L, 1, 2);
matmul(A, L, A);
matmul(P, L, P);
}
if (A[2][1] != 0) {
if (A[2][1] != 0)
{
zero_second_column(L, A);
matmul(A, L, A);
matmul(P, L, P);
@ -265,13 +283,16 @@ static void finalize(long A[3][3], long P[3][3], long Q[3][3])
static void finalize_sort(long A[3][3], long P[3][3], long Q[3][3])
{
if (A[0][0] > A[1][1]) {
if (A[0][0] > A[1][1])
{
swap_diag_elems(A, P, Q, 0, 1);
}
if (A[1][1] > A[2][2]) {
if (A[1][1] > A[2][2])
{
swap_diag_elems(A, P, Q, 1, 2);
}
if (A[0][0] > A[1][1]) {
if (A[0][0] > A[1][1])
{
swap_diag_elems(A, P, Q, 0, 1);
}
}
@ -280,7 +301,8 @@ static void finalize_disturb(long A[3][3], long Q[3][3], const int i, const int
{
long L[3][3];
if (A[j][j] % A[i][i] != 0) {
if (A[j][j] % A[i][i] != 0)
{
transpose(A);
disturb_rows(L, i, j);
matmul(A, L, A);
@ -325,8 +347,10 @@ static void make_diagA_positive(long A[3][3], long P[3][3])
int i;
long L[3][3];
for (i = 0; i < 3; i++) {
if (A[i][i] < 0) {
for (i = 0; i < 3; i++)
{
if (A[i][i] < 0)
{
flip_sign_row(L, i);
matmul(A, L, A);
matmul(P, L, P);
@ -338,9 +362,12 @@ static void flip_PQ(long P[3][3], long Q[3][3])
{
int i, j;
if (det(P) < 0) {
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
if (det(P) < 0)
{
for (i = 0; i < 3; i++)
{
for (j = 0; j < 3; j++)
{
P[i][j] *= -1;
Q[i][j] *= -1;
}
@ -393,16 +420,18 @@ static void extended_gcd(long retvals[3], const long a, const long b)
int i;
long vals[6];
vals[0] = a; /* r0 */
vals[1] = b; /* r1 */
vals[2] = 1; /* s0 */
vals[3] = 0; /* s1 */
vals[4] = 0; /* t0 */
vals[5] = 1; /* t1 */
vals[0] = a; /* r0 */
vals[1] = b; /* r1 */
vals[2] = 1; /* s0 */
vals[3] = 0; /* s1 */
vals[4] = 0; /* t0 */
vals[5] = 1; /* t1 */
for (i = 0; i < 1000; i++) {
for (i = 0; i < 1000; i++)
{
extended_gcd_step(vals);
if (vals[1] == 0) {
if (vals[1] == 0)
{
break;
}
}
@ -420,12 +449,15 @@ static void extended_gcd_step(long vals[6])
q = vals[0] / vals[1];
r2 = vals[0] % vals[1];
if (r2 < 0) {
if (vals[1] > 0) {
if (r2 < 0)
{
if (vals[1] > 0)
{
r2 += vals[1];
q -= 1;
}
if (vals[1] < 0) {
if (vals[1] < 0)
{
r2 -= vals[1];
q += 1;
}
@ -462,8 +494,10 @@ static void transpose(long m[3][3])
long tmp;
int i, j;
for (i = 0; i < 3; i++) {
for (j = i; j < 3; j++) {
for (i = 0; i < 3; i++)
{
for (j = i; j < 3; j++)
{
tmp = m[i][j];
m[i][j] = m[j][i];
m[j][i] = tmp;
@ -478,14 +512,18 @@ static void matmul(long m[3][3],
int i, j;
long c[3][3];
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
for (i = 0; i < 3; i++)
{
for (j = 0; j < 3; j++)
{
c[i][j] = a[i][0] * b[0][j] + a[i][1] * b[1][j] + a[i][2] * b[2][j];
}
}
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
for (i = 0; i < 3; i++)
{
for (j = 0; j < 3; j++)
{
m[i][j] = c[i][j];
}
}
@ -493,12 +531,9 @@ static void matmul(long m[3][3],
static long det(SNF3x3CONST long m[3][3])
{
return m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1])
+ m[0][1] * (m[1][2] * m[2][0] - m[1][0] * m[2][2])
+ m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0]);
return m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1]) + m[0][1] * (m[1][2] * m[2][0] - m[1][0] * m[2][2]) + m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0]);
}
/* int main()
* {
* test_extended_gcd();

1080
c/tetrahedron_method.c
View File

File diff suppressed because it is too large Load Diff

70
c/triplet.c
View File

@ -99,25 +99,25 @@ void tpl_get_integration_weight(double *iw,
#ifdef PHPYOPENMP
#pragma omp parallel for if (openmp_per_triplets)
#endif
for (i = 0; i < num_triplets; i++) {
for (i = 0; i < num_triplets; i++)
{
tpi_get_integration_weight(iw + i * num_band_prod,
iw_zero + i * num_band_prod,
frequency_points, /* f0 */
frequency_points, /* f0 */
num_band0,
tp_relative_grid_address,
triplets[i],
num_triplets,
bzgrid,
frequencies1, /* f1 */
frequencies1, /* f1 */
num_band1,
frequencies2, /* f2 */
frequencies2, /* f2 */
num_band2,
tp_type,
openmp_per_bands);
}
}
void tpl_get_integration_weight_with_sigma(double *iw,
char *iw_zero,
const double sigma,
@ -140,35 +140,38 @@ void tpl_get_integration_weight_with_sigma(double *iw,
#ifdef PHPYOPENMP
#pragma omp parallel for
#endif
for (i = 0; i < num_triplets; i++) {
for (i = 0; i < num_triplets; i++)
{
tpi_get_integration_weight_with_sigma(
iw + i * num_band_prod,
iw_zero + i * num_band_prod,
sigma,
cutoff,
frequency_points,
num_band0,
triplets[i],
const_adrs_shift,
frequencies,
num_band,
tp_type,
0);
iw + i * num_band_prod,
iw_zero + i * num_band_prod,
sigma,
cutoff,
frequency_points,
num_band0,
triplets[i],
const_adrs_shift,
frequencies,
num_band,
tp_type,
0);
}
}
long tpl_is_N(const long triplet[3], const long (*bz_grid_addresses)[3])
{
long i, j, sum_q, is_N;
is_N = 1;
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
sum_q = 0;
for (j = 0; j < 3; j++) { /* 1st, 2nd, 3rd triplet */
for (j = 0; j < 3; j++)
{ /* 1st, 2nd, 3rd triplet */
sum_q += bz_grid_addresses[triplet[j]][i];
}
if (sum_q) {
if (sum_q)
{
is_N = 0;
break;
}
@ -177,28 +180,33 @@ long tpl_is_N(const long triplet[3], const long (*bz_grid_addresses)[3])
}
void tpl_set_relative_grid_address(
long tp_relative_grid_address[2][24][4][3],
const long relative_grid_address[24][4][3],
const long tp_type)
long tp_relative_grid_address[2][24][4][3],
const long relative_grid_address[24][4][3],
const long tp_type)
{
long i, j, k, l;
long signs[2];
signs[0] = 1;
signs[1] = 1;
if ((tp_type == 2) || (tp_type == 3)) {
if ((tp_type == 2) || (tp_type == 3))
{
/* q1+q2+q3=G */
/* To set q2+1, q3-1 is needed to keep G */
signs[1] = -1;
}
/* tp_type == 4, q+k_i-k_f=G */
for (i = 0; i < 2; i++) {
for (j = 0; j < 24; j++) {
for (k = 0; k < 4; k++) {
for (l = 0; l < 3; l++) {
for (i = 0; i < 2; i++)
{
for (j = 0; j < 24; j++)
{
for (k = 0; k < 4; k++)
{
for (l = 0; l < 3; l++)
{
tp_relative_grid_address[i][j][k][l] =
relative_grid_address[j][k][l] * signs[i];
relative_grid_address[j][k][l] * signs[i];
}
}
}

256
c/triplet_grid.c
View File

@ -46,19 +46,19 @@ static long get_ir_triplets_at_q(long *map_triplets,
long *map_q,
const long grid_point,
const long D_diag[3],
const RotMats * rot_reciprocal,
const RotMats *rot_reciprocal,
const long swappable);
static long get_ir_triplets_at_q_perm_q1q2(long *map_triplets,
const long *map_q,
const long grid_point,
const long D_diag[3],
const RotMats * rot_reciprocal_q,
const RotMats *rot_reciprocal_q,
const long num_ir_q);
static long get_ir_triplets_at_q_noperm(long *map_triplets,
const long *map_q,
const long grid_point,
const long D_diag[3],
const RotMats * rot_reciprocal_q);
const RotMats *rot_reciprocal_q);
static long get_BZ_triplets_at_q(long (*triplets)[3],
const long grid_point,
const ConstBZGrid *bzgrid,
@ -76,7 +76,7 @@ static void get_BZ_triplets_at_q_type2(long (*triplets)[3],
static double get_squared_distance(const long G[3],
const double LQD_inv[3][3]);
static void get_LQD_inv(double LQD_inv[3][3], const ConstBZGrid *bzgrid);
static RotMats *get_reciprocal_point_group_with_q(const RotMats * rot_reciprocal,
static RotMats *get_reciprocal_point_group_with_q(const RotMats *rot_reciprocal,
const long D_diag[3],
const long grid_point);
static RotMats *get_reciprocal_point_group(const long (*rec_rotations_in)[3][3],
@ -100,7 +100,8 @@ long tpk_get_ir_triplets_at_q(long *map_triplets,
num_rot,
is_time_reversal,
0);
if (rotations == NULL) {
if (rotations == NULL)
{
return 0;
}
@ -131,7 +132,7 @@ static long get_ir_triplets_at_q(long *map_triplets,
long *map_q,
const long grid_point,
const long D_diag[3],
const RotMats * rot_reciprocal,
const RotMats *rot_reciprocal,
const long swappable)
{
long i, num_ir_q, num_ir_triplets;
@ -140,7 +141,8 @@ static long get_ir_triplets_at_q(long *map_triplets,
rot_reciprocal_q = NULL;
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
PS[i] = 0;
}
@ -155,20 +157,25 @@ static long get_ir_triplets_at_q(long *map_triplets,
D_diag,
PS);
num_ir_q = 0;
for (i = 0; i < D_diag[0] * D_diag[1] * D_diag[2]; i++) {
if (map_q[i] == i) {
for (i = 0; i < D_diag[0] * D_diag[1] * D_diag[2]; i++)
{
if (map_q[i] == i)
{
num_ir_q++;
}
}
if (swappable) {
if (swappable)
{
num_ir_triplets = get_ir_triplets_at_q_perm_q1q2(map_triplets,
map_q,
grid_point,
D_diag,
rot_reciprocal_q,
num_ir_q);
} else {
}
else
{
num_ir_triplets = get_ir_triplets_at_q_noperm(map_triplets,
map_q,
grid_point,
@ -186,7 +193,7 @@ static long get_ir_triplets_at_q_perm_q1q2(long *map_triplets,
const long *map_q,
const long grid_point,
const long D_diag[3],
const RotMats * rot_reciprocal_q,
const RotMats *rot_reciprocal_q,
const long num_ir_q)
{
long i, j, num_grid, num_ir_triplets, ir_gp, count;
@ -199,19 +206,23 @@ static long get_ir_triplets_at_q_perm_q1q2(long *map_triplets,
num_grid = D_diag[0] * D_diag[1] * D_diag[2];
if ((q_2 = (long*) malloc(sizeof(long) * num_ir_q)) == NULL) {
if ((q_2 = (long *)malloc(sizeof(long) * num_ir_q)) == NULL)
{
warning_print("Memory could not be allocated.");
goto ret;
}
if ((ir_gps_at_q = (long*) malloc(sizeof(long) * num_ir_q)) == NULL) {
if ((ir_gps_at_q = (long *)malloc(sizeof(long) * num_ir_q)) == NULL)
{
warning_print("Memory could not be allocated.");
goto ret;
}
count = 0;
for (i = 0; i < num_grid; i++) {
if (map_q[i] == i) {
for (i = 0; i < num_grid; i++)
{
if (map_q[i] == i)
{
ir_gps_at_q[count] = i;
count++;
}
@ -222,22 +233,28 @@ static long get_ir_triplets_at_q_perm_q1q2(long *map_triplets,
#ifdef PHPYOPENMP
#pragma omp parallel for private(j, adrs1, adrs2)
#endif
for (i = 0; i < num_ir_q; i++) {
for (i = 0; i < num_ir_q; i++)
{
grg_get_grid_address_from_index(adrs1, ir_gps_at_q[i], D_diag);
for (j = 0; j < 3; j++) { /* q'' */
adrs2[j] = - adrs0[j] - adrs1[j];
for (j = 0; j < 3; j++)
{ /* q'' */
adrs2[j] = -adrs0[j] - adrs1[j];
}
q_2[i] = grg_get_grid_index(adrs2, D_diag);
q_2[i] = grg_get_grid_index(adrs2, D_diag);
}
/* map_q[q_2[i]] is in ir_gps_at_q. */
/* If map_q[q_2[i]] < ir_gps_at_q[i], this should be already */
/* stored. So the counter is not incremented. */
for (i = 0; i < num_ir_q; i++) {
for (i = 0; i < num_ir_q; i++)
{
ir_gp = ir_gps_at_q[i];
if (map_q[q_2[i]] < ir_gp) {
if (map_q[q_2[i]] < ir_gp)
{
map_triplets[ir_gp] = map_q[q_2[i]];
} else {
}
else
{
map_triplets[ir_gp] = ir_gp;
num_ir_triplets++;
}
@ -246,16 +263,19 @@ static long get_ir_triplets_at_q_perm_q1q2(long *map_triplets,
#ifdef PHPYOPENMP
#pragma omp parallel for
#endif
for (i = 0; i < num_grid; i++) {
for (i = 0; i < num_grid; i++)
{
map_triplets[i] = map_triplets[map_q[i]];
}
ret:
if (q_2) {
if (q_2)
{
free(q_2);
q_2 = NULL;
}
if (ir_gps_at_q) {
if (ir_gps_at_q)
{
free(ir_gps_at_q);
ir_gps_at_q = NULL;
}
@ -266,18 +286,22 @@ static long get_ir_triplets_at_q_noperm(long *map_triplets,
const long *map_q,
const long grid_point,
const long D_diag[3],
const RotMats * rot_reciprocal_q)
const RotMats *rot_reciprocal_q)
{
long i, num_grid, num_ir_triplets;
num_ir_triplets = 0;
num_grid = D_diag[0] * D_diag[1] * D_diag[2];
for (i = 0; i < num_grid; i++) {
if (map_q[i] == i) {
for (i = 0; i < num_grid; i++)
{
if (map_q[i] == i)
{
map_triplets[i] = i;
num_ir_triplets++;
} else {
}
else
{
map_triplets[i] = map_triplets[map_q[i]];
}
}
@ -296,26 +320,31 @@ static long get_BZ_triplets_at_q(long (*triplets)[3],
ir_q1_gps = NULL;
num_ir = 0;
if ((ir_q1_gps = (long*) malloc(sizeof(long) * bzgrid->size))
== NULL) {
if ((ir_q1_gps = (long *)malloc(sizeof(long) * bzgrid->size)) == NULL)
{
warning_print("Memory could not be allocated.");
goto ret;
}
for (i = 0; i < bzgrid->size; i++) {
if (map_triplets[i] == i) {
for (i = 0; i < bzgrid->size; i++)
{
if (map_triplets[i] == i)
{
ir_q1_gps[num_ir] = i;
num_ir++;
}
}
if (bzgrid->type == 1) {
if (bzgrid->type == 1)
{
get_BZ_triplets_at_q_type1(triplets,
grid_point,
bzgrid,
ir_q1_gps,
num_ir);
} else {
}
else
{
get_BZ_triplets_at_q_type2(triplets,
grid_point,
bzgrid,
@ -340,7 +369,7 @@ static void get_BZ_triplets_at_q_type1(long (*triplets)[3],
long bzgp[3], G[3];
long bz_adrs0[3], bz_adrs1[3], bz_adrs2[3];
const long *gp_map;
const long (*bz_adrs)[3];
const long(*bz_adrs)[3];
double d2, min_d2, tolerance;
double LQD_inv[3][3];
@ -348,9 +377,10 @@ static void get_BZ_triplets_at_q_type1(long (*triplets)[3],
bz_adrs = bzgrid->addresses;
get_LQD_inv(LQD_inv, bzgrid);
/* This tolerance is used to be consistent to BZ reduction in bzgrid. */
tolerance = bzg_get_tolerance_for_BZ_reduction((BZGrid*)bzgrid);
tolerance = bzg_get_tolerance_for_BZ_reduction((BZGrid *)bzgrid);
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
bz_adrs0[i] = bz_adrs[grid_point][i];
}
num_gp = bzgrid->D_diag[0] * bzgrid->D_diag[1] * bzgrid->D_diag[2];
@ -359,61 +389,77 @@ static void get_BZ_triplets_at_q_type1(long (*triplets)[3],
#ifdef PHPYOPENMP
#pragma omp parallel for private(j, gp2, bzgp, G, bz_adrs1, bz_adrs2, d2, min_d2, bz0, bz1, bz2)
#endif
for (i = 0; i < num_ir; i++) {
for (j = 0; j < 3; j++) {
for (i = 0; i < num_ir; i++)
{
for (j = 0; j < 3; j++)
{
bz_adrs1[j] = bz_adrs[ir_q1_gps[i]][j];
bz_adrs2[j] = - bz_adrs0[j] - bz_adrs1[j];
bz_adrs2[j] = -bz_adrs0[j] - bz_adrs1[j];
}
gp2 = grg_get_grid_index(bz_adrs2, bzgrid->D_diag);
/* Negative value is the signal to initialize min_d2 later. */
min_d2 = -1;
for (bz0 = 0;
bz0 < gp_map[num_bzgp + grid_point + 1]
- gp_map[num_bzgp + grid_point] + 1;
bz0++) {
if (bz0 == 0) {
bz0 < gp_map[num_bzgp + grid_point + 1] - gp_map[num_bzgp + grid_point] + 1;
bz0++)
{
if (bz0 == 0)
{
bzgp[0] = grid_point;
} else {
}
else
{
bzgp[0] = num_gp + gp_map[num_bzgp + grid_point] + bz0 - 1;
}
for (bz1 = 0;
bz1 < gp_map[num_bzgp + ir_q1_gps[i] + 1]
- gp_map[num_bzgp + ir_q1_gps[i]] + 1;
bz1++) {
if (bz1 == 0) {
bz1 < gp_map[num_bzgp + ir_q1_gps[i] + 1] - gp_map[num_bzgp + ir_q1_gps[i]] + 1;
bz1++)
{
if (bz1 == 0)
{
bzgp[1] = ir_q1_gps[i];
} else {
}
else
{
bzgp[1] = num_gp + gp_map[num_bzgp + ir_q1_gps[i]] + bz1 - 1;
}
for (bz2 = 0;
bz2 < gp_map[num_bzgp + gp2 + 1] - gp_map[num_bzgp + gp2] + 1;
bz2++) {
if (bz2 == 0) {
bz2++)
{
if (bz2 == 0)
{
bzgp[2] = gp2;
} else {
}
else
{
bzgp[2] = num_gp + gp_map[num_bzgp + gp2] + bz2 - 1;
}
for (j = 0; j < 3; j++) {
for (j = 0; j < 3; j++)
{
G[j] = bz_adrs[bzgp[0]][j] + bz_adrs[bzgp[1]][j] + bz_adrs[bzgp[2]][j];
}
if (G[0] == 0 && G[1] == 0 && G[2] == 0) {
for (j = 0; j < 3; j++) {
if (G[0] == 0 && G[1] == 0 && G[2] == 0)
{
for (j = 0; j < 3; j++)
{
triplets[i][j] = bzgp[j];
}
goto found;
}
d2 = get_squared_distance(G, LQD_inv);
if (d2 < min_d2 - tolerance || min_d2 < 0) {
if (d2 < min_d2 - tolerance || min_d2 < 0)
{
min_d2 = d2;
for (j = 0; j < 3; j++) {
for (j = 0; j < 3; j++)
{
triplets[i][j] = bzgp[j];
}
}
}
}
}
found:
;
found:;
}
}
@ -427,7 +473,7 @@ static void get_BZ_triplets_at_q_type2(long (*triplets)[3],
long bzgp[3], G[3];
long bz_adrs0[3], bz_adrs1[3], bz_adrs2[3];
const long *gp_map;
const long (*bz_adrs)[3];
const long(*bz_adrs)[3];
double d2, min_d2, tolerance;
double LQD_inv[3][3];
@ -435,9 +481,10 @@ static void get_BZ_triplets_at_q_type2(long (*triplets)[3],
bz_adrs = bzgrid->addresses;
get_LQD_inv(LQD_inv, bzgrid);
/* This tolerance is used to be consistent to BZ reduction in bzgrid. */
tolerance = bzg_get_tolerance_for_BZ_reduction((BZGrid*)bzgrid);
tolerance = bzg_get_tolerance_for_BZ_reduction((BZGrid *)bzgrid);
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
bz_adrs0[i] = bz_adrs[grid_point][i];
}
gp0 = grg_get_grid_index(bz_adrs0, bzgrid->D_diag);
@ -445,39 +492,48 @@ static void get_BZ_triplets_at_q_type2(long (*triplets)[3],
#ifdef PHPYOPENMP
#pragma omp parallel for private(j, gp2, bzgp, G, bz_adrs1, bz_adrs2, d2, min_d2)
#endif
for (i = 0; i < num_ir; i++) {
for (j = 0; j < 3; j++) {
for (i = 0; i < num_ir; i++)
{
for (j = 0; j < 3; j++)
{
bz_adrs1[j] = bz_adrs[gp_map[ir_q1_gps[i]]][j];
bz_adrs2[j] = - bz_adrs0[j] - bz_adrs1[j];
bz_adrs2[j] = -bz_adrs0[j] - bz_adrs1[j];
}
gp2 = grg_get_grid_index(bz_adrs2, bzgrid->D_diag);
/* Negative value is the signal to initialize min_d2 later. */
min_d2 = -1;
for (bzgp[0] = gp_map[gp0]; bzgp[0] < gp_map[gp0 + 1]; bzgp[0]++) {
for (bzgp[0] = gp_map[gp0]; bzgp[0] < gp_map[gp0 + 1]; bzgp[0]++)
{
for (bzgp[1] = gp_map[ir_q1_gps[i]];
bzgp[1] < gp_map[ir_q1_gps[i] + 1]; bzgp[1]++) {
for (bzgp[2] = gp_map[gp2]; bzgp[2] < gp_map[gp2 + 1]; bzgp[2]++) {
for (j = 0; j < 3; j++) {
bzgp[1] < gp_map[ir_q1_gps[i] + 1]; bzgp[1]++)
{
for (bzgp[2] = gp_map[gp2]; bzgp[2] < gp_map[gp2 + 1]; bzgp[2]++)
{
for (j = 0; j < 3; j++)
{
G[j] = bz_adrs[bzgp[0]][j] + bz_adrs[bzgp[1]][j] + bz_adrs[bzgp[2]][j];
}
if (G[0] == 0 && G[1] == 0 && G[2] == 0) {
for (j = 0; j < 3; j++) {
if (G[0] == 0 && G[1] == 0 && G[2] == 0)
{
for (j = 0; j < 3; j++)
{
triplets[i][j] = bzgp[j];
}
goto found;
}
d2 = get_squared_distance(G, LQD_inv);
if (d2 < min_d2 - tolerance || min_d2 < 0) {
if (d2 < min_d2 - tolerance || min_d2 < 0)
{
min_d2 = d2;
for (j = 0; j < 3; j++) {
for (j = 0; j < 3; j++)
{
triplets[i][j] = bzgp[j];
}
}
}
}
}
found:
;
found:;
}
}
@ -488,7 +544,8 @@ static double get_squared_distance(const long G[3],
long i;
d2 = 0;
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
d = LQD_inv[i][0] * G[0] + LQD_inv[i][1] * G[1] + LQD_inv[i][2] * G[2];
d2 += d * d;
}
@ -501,25 +558,27 @@ static void get_LQD_inv(double LQD_inv[3][3], const ConstBZGrid *bzgrid)
long i, j, k;
/* LQD^-1 */
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
for (k = 0; k < 3; k++) {
LQD_inv[i][k]
= bzgrid->reclat[i][j] * bzgrid->Q[j][k] / bzgrid->D_diag[k];
for (i = 0; i < 3; i++)
{
for (j = 0; j < 3; j++)
{
for (k = 0; k < 3; k++)
{
LQD_inv[i][k] = bzgrid->reclat[i][j] * bzgrid->Q[j][k] / bzgrid->D_diag[k];
}
}
}
}
/* Return NULL if failed */
static RotMats *get_reciprocal_point_group_with_q(const RotMats * rot_reciprocal,
static RotMats *get_reciprocal_point_group_with_q(const RotMats *rot_reciprocal,
const long D_diag[3],
const long grid_point)
{
long i, num_rot, gp_rot;
long *ir_rot;
long adrs[3], adrs_rot[3];
RotMats * rot_reciprocal_q;
RotMats *rot_reciprocal_q;
ir_rot = NULL;
rot_reciprocal_q = NULL;
@ -527,26 +586,32 @@ static RotMats *get_reciprocal_point_group_with_q(const RotMats * rot_reciprocal
grg_get_grid_address_from_index(adrs, grid_point, D_diag);
if ((ir_rot = (long*)malloc(sizeof(long) * rot_reciprocal->size)) == NULL) {
if ((ir_rot = (long *)malloc(sizeof(long) * rot_reciprocal->size)) == NULL)
{
warning_print("Memory of ir_rot could not be allocated.");
return NULL;
}
for (i = 0; i < rot_reciprocal->size; i++) {
for (i = 0; i < rot_reciprocal->size; i++)
{
ir_rot[i] = -1;
}
for (i = 0; i < rot_reciprocal->size; i++) {
for (i = 0; i < rot_reciprocal->size; i++)
{
lagmat_multiply_matrix_vector_l3(adrs_rot, rot_reciprocal->mat[i], adrs);
gp_rot = grg_get_grid_index(adrs_rot, D_diag);
if (gp_rot == grid_point) {
if (gp_rot == grid_point)
{
ir_rot[num_rot] = i;
num_rot++;
}
}
if ((rot_reciprocal_q = bzg_alloc_RotMats(num_rot)) != NULL) {
for (i = 0; i < num_rot; i++) {
if ((rot_reciprocal_q = bzg_alloc_RotMats(num_rot)) != NULL)
{
for (i = 0; i < num_rot; i++)
{
lagmat_copy_matrix_l3(rot_reciprocal_q->mat[i],
rot_reciprocal->mat[ir_rot[i]]);
}
@ -558,7 +623,6 @@ static RotMats *get_reciprocal_point_group_with_q(const RotMats * rot_reciprocal
return rot_reciprocal_q;
}
static RotMats *get_reciprocal_point_group(const long (*rec_rotations_in)[3][3],
const long num_rot,
const long is_time_reversal,
@ -573,12 +637,14 @@ static RotMats *get_reciprocal_point_group(const long (*rec_rotations_in)[3][3],
num_rot,
is_time_reversal,
is_transpose);
if (num_rot_out == 0) {
if (num_rot_out == 0)
{
return NULL;
}
rec_rotations = bzg_alloc_RotMats(num_rot_out);
for (i = 0; i < num_rot_out; i++) {
for (i = 0; i < num_rot_out; i++)
{
lagmat_copy_matrix_l3(rec_rotations->mat[i], rec_rotations_out[i]);
}

196
c/triplet_iw.c
View File

@ -53,10 +53,10 @@ static long set_g(double g[3],
const double freq_vertices[3][24][4],
const long max_i);
static void get_triplet_tetrahedra_vertices(
long vertices[2][24][4],
const long tp_relative_grid_address[2][24][4][3],
const long triplet[3],
const ConstBZGrid *bzgrid);
long vertices[2][24][4],
const long tp_relative_grid_address[2][24][4][3],
const long triplet[3],
const ConstBZGrid *bzgrid);
static void
get_neighboring_grid_points_type1(long *neighboring_grid_points,
const long grid_point,
@ -70,21 +70,20 @@ get_neighboring_grid_points_type2(long *neighboring_grid_points,
const long num_relative_grid_address,
const ConstBZGrid *bzgrid);
void
tpi_get_integration_weight(double *iw,
char *iw_zero,
const double *frequency_points,
const long num_band0,
const long tp_relative_grid_address[2][24][4][3],
const long triplets[3],
const long num_triplets,
const ConstBZGrid *bzgrid,
const double *frequencies1,
const long num_band1,
const double *frequencies2,
const long num_band2,
const long tp_type,
const long openmp_per_bands)
void tpi_get_integration_weight(double *iw,
char *iw_zero,
const double *frequency_points,
const long num_band0,
const long tp_relative_grid_address[2][24][4][3],
const long triplets[3],
const long num_triplets,
const ConstBZGrid *bzgrid,
const double *frequencies1,
const long num_band1,
const double *frequencies2,
const long num_band2,
const long tp_type,
const long openmp_per_bands)
{
long max_i, j, b1, b2, b12, num_band_prod, adrs_shift;
long vertices[2][24][4];
@ -109,37 +108,44 @@ tpi_get_integration_weight(double *iw,
/* tp_type = 4: (g[0]) mainly for el-ph phonon decay, */
/* f0: ph, f1: el_i, f2: el_f */
if ((tp_type == 2) || (tp_type == 3)) {
if ((tp_type == 2) || (tp_type == 3))
{
max_i = 3;
}
if (tp_type == 4) {
if (tp_type == 4)
{
max_i = 1;
}
#ifdef PHPYOPENMP
#pragma omp parallel for private(j, b1, b2, adrs_shift, g, freq_vertices) if (openmp_per_bands)
#endif
for (b12 = 0; b12 < num_band1 * num_band2; b12++) {
for (b12 = 0; b12 < num_band1 * num_band2; b12++)
{
b1 = b12 / num_band2;
b2 = b12 % num_band2;
set_freq_vertices(freq_vertices, frequencies1, frequencies2,
vertices, num_band1, num_band2, b1, b2, tp_type);
for (j = 0; j < num_band0; j++) {
for (j = 0; j < num_band0; j++)
{
adrs_shift = j * num_band1 * num_band2 + b1 * num_band2 + b2;
iw_zero[adrs_shift] = set_g(g, frequency_points[j], freq_vertices, max_i);
if (tp_type == 2) {
if (tp_type == 2)
{
iw[adrs_shift] = g[2];
adrs_shift += num_band_prod;
iw[adrs_shift] = g[0] - g[1];
}
if (tp_type == 3) {
if (tp_type == 3)
{
iw[adrs_shift] = g[2];
adrs_shift += num_band_prod;
iw[adrs_shift] = g[0] - g[1];
adrs_shift += num_band_prod;
iw[adrs_shift] = g[0] + g[1] + g[2];
}
if (tp_type == 4) {
if (tp_type == 4)
{
iw[adrs_shift] = g[0];
}
}
@ -165,36 +171,44 @@ void tpi_get_integration_weight_with_sigma(double *iw,
#ifdef PHPYOPENMP
#pragma omp parallel for private(j, b1, b2, f0, f1, f2, g0, g1, g2, adrs_shift) if (openmp_per_bands)
#endif
for (b12 = 0; b12 < num_band * num_band; b12++) {
for (b12 = 0; b12 < num_band * num_band; b12++)
{
b1 = b12 / num_band;
b2 = b12 % num_band;
f1 = frequencies[triplet[1] * num_band + b1];
f2 = frequencies[triplet[2] * num_band + b2];
for (j = 0; j < num_band0; j++) {
for (j = 0; j < num_band0; j++)
{
f0 = frequency_points[j];
adrs_shift = j * num_band * num_band + b1 * num_band + b2;
if ((tp_type == 2) || (tp_type == 3)) {
if ((tp_type == 2) || (tp_type == 3))
{
if (cutoff > 0 &&
fabs(f0 + f1 - f2) > cutoff &&
fabs(f0 - f1 + f2) > cutoff &&
fabs(f0 - f1 - f2) > cutoff) {
fabs(f0 - f1 - f2) > cutoff)
{
iw_zero[adrs_shift] = 1;
g0 = 0;
g1 = 0;
g2 = 0;
} else {
}
else
{
iw_zero[adrs_shift] = 0;
g0 = phonoc_gaussian(f0 + f1 - f2, sigma);
g1 = phonoc_gaussian(f0 - f1 + f2, sigma);
g2 = phonoc_gaussian(f0 - f1 - f2, sigma);
}
if (tp_type == 2) {
if (tp_type == 2)
{
iw[adrs_shift] = g2;
adrs_shift += const_adrs_shift;
iw[adrs_shift] = g0 - g1;
}
if (tp_type == 3) {
if (tp_type == 3)
{
iw[adrs_shift] = g2;
adrs_shift += const_adrs_shift;
iw[adrs_shift] = g0 - g1;
@ -202,11 +216,15 @@ void tpi_get_integration_weight_with_sigma(double *iw,
iw[adrs_shift] = g0 + g1 + g2;
}
}
if (tp_type == 4) {
if (cutoff > 0 && fabs(f0 + f1 - f2) > cutoff) {
if (tp_type == 4)
{
if (cutoff > 0 && fabs(f0 + f1 - f2) > cutoff)
{
iw_zero[adrs_shift] = 1;
iw[adrs_shift] = 0;
} else {
}
else
{
iw_zero[adrs_shift] = 0;
iw[adrs_shift] = phonoc_gaussian(f0 + f1 - f2, sigma);
}
@ -215,22 +233,22 @@ void tpi_get_integration_weight_with_sigma(double *iw,
}
}
void
tpi_get_neighboring_grid_points(long *neighboring_grid_points,
const long grid_point,
const long (*relative_grid_address)[3],
const long num_relative_grid_address,
const ConstBZGrid *bzgrid)
void tpi_get_neighboring_grid_points(long *neighboring_grid_points,
const long grid_point,
const long (*relative_grid_address)[3],
const long num_relative_grid_address,
const ConstBZGrid *bzgrid)
{
if (bzgrid->type == 1) {
if (bzgrid->type == 1)
{
get_neighboring_grid_points_type1(neighboring_grid_points,
grid_point,
relative_grid_address,
num_relative_grid_address,
bzgrid);
} else {
}
else
{
get_neighboring_grid_points_type2(neighboring_grid_points,
grid_point,
relative_grid_address,
@ -252,17 +270,28 @@ static void set_freq_vertices(double freq_vertices[3][24][4],
long i, j;
double f1, f2;
for (i = 0; i < 24; i++) {
for (j = 0; j < 4; j++) {
for (i = 0; i < 24; i++)
{
for (j = 0; j < 4; j++)
{
f1 = frequencies1[vertices[0][i][j] * num_band1 + b1];
f2 = frequencies2[vertices[1][i][j] * num_band2 + b2];
if ((tp_type == 2) || (tp_type == 3)) {
if (f1 < 0) {f1 = 0;}
if (f2 < 0) {f2 = 0;}
if ((tp_type == 2) || (tp_type == 3))
{
if (f1 < 0)
{
f1 = 0;
}
if (f2 < 0)
{
f2 = 0;
}
freq_vertices[0][i][j] = -f1 + f2;
freq_vertices[1][i][j] = f1 - f2;
freq_vertices[2][i][j] = f1 + f2;
} else {
}
else
{
freq_vertices[0][i][j] = -f1 + f2;
}
}
@ -287,11 +316,15 @@ static long set_g(double g[3],
iw_zero = 1;
for (i = 0; i < max_i; i++) {
if (thm_in_tetrahedra(f0, freq_vertices[i])) {
for (i = 0; i < max_i; i++)
{
if (thm_in_tetrahedra(f0, freq_vertices[i]))
{
g[i] = thm_get_integration_weight(f0, freq_vertices[i], 'I');
iw_zero = 0;
} else {
}
else
{
g[i] = 0;
}
}
@ -300,14 +333,16 @@ static long set_g(double g[3],
}
static void get_triplet_tetrahedra_vertices(long vertices[2][24][4],
const long tp_relative_grid_address[2][24][4][3],
const long triplet[3],
const ConstBZGrid *bzgrid)
const long tp_relative_grid_address[2][24][4][3],
const long triplet[3],
const ConstBZGrid *bzgrid)
{
long i, j;
for (i = 0; i < 2; i++) {
for (j = 0; j < 24; j++) {
for (i = 0; i < 2; i++)
{
for (j = 0; j < 24; j++)
{
tpi_get_neighboring_grid_points(vertices[i][j],
triplet[i + 1],
tp_relative_grid_address[i][j],
@ -327,20 +362,25 @@ get_neighboring_grid_points_type1(long *neighboring_grid_points,
long bzmesh[3], bz_address[3];
long i, j, bz_gp, prod_bz_mesh;
for (i = 0; i < 3; i++) {
for (i = 0; i < 3; i++)
{
bzmesh[i] = bzgrid->D_diag[i] * 2;
}
prod_bz_mesh = bzmesh[0] * bzmesh[1] * bzmesh[2];
for (i = 0; i < num_relative_grid_address; i++) {
for (j = 0; j < 3; j++) {
bz_address[j] = bzgrid->addresses[grid_point][j]
+ relative_grid_address[i][j];
for (i = 0; i < num_relative_grid_address; i++)
{
for (j = 0; j < 3; j++)
{
bz_address[j] = bzgrid->addresses[grid_point][j] + relative_grid_address[i][j];
}
bz_gp = bzgrid->gp_map[grg_get_grid_index(bz_address, bzmesh)];
if (bz_gp == prod_bz_mesh) {
if (bz_gp == prod_bz_mesh)
{
neighboring_grid_points[i] =
grg_get_grid_index(bz_address, bzgrid->D_diag);
} else {
grg_get_grid_index(bz_address, bzgrid->D_diag);
}
else
{
neighboring_grid_points[i] = bz_gp;
}
}
@ -356,18 +396,20 @@ get_neighboring_grid_points_type2(long *neighboring_grid_points,
long bz_address[3];
long i, j, gp;
for (i = 0; i < num_relative_grid_address; i++) {
for (j = 0; j < 3; j++) {
bz_address[j] = bzgrid->addresses[grid_point][j]
+ relative_grid_address[i][j];
for (i = 0; i < num_relative_grid_address; i++)
{
for (j = 0; j < 3; j++)
{
bz_address[j] = bzgrid->addresses[grid_point][j] + relative_grid_address[i][j];
}
gp = grg_get_grid_index(bz_address, bzgrid->D_diag);
neighboring_grid_points[i] = bzgrid->gp_map[gp];
if (bzgrid->gp_map[gp + 1] - bzgrid->gp_map[gp] > 1) {
for (j = bzgrid->gp_map[gp]; j < bzgrid->gp_map[gp + 1]; j++) {
if (bz_address[0] == bzgrid->addresses[j][0]
&& bz_address[1] == bzgrid->addresses[j][1]
&& bz_address[2] == bzgrid->addresses[j][2]) {
if (bzgrid->gp_map[gp + 1] - bzgrid->gp_map[gp] > 1)
{
for (j = bzgrid->gp_map[gp]; j < bzgrid->gp_map[gp + 1]; j++)
{
if (bz_address[0] == bzgrid->addresses[j][0] && bz_address[1] == bzgrid->addresses[j][1] && bz_address[2] == bzgrid->addresses[j][2])
{
neighboring_grid_points[i] = j;
break;
}