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

Removed LAGCONST from phono3c

This commit is contained in:
Atsushi Togo 2021-03-15 13:02:18 +09:00
parent eeae56ca25
commit 5115478996
37 changed files with 315 additions and 280 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
c/_phono3py.c
View File

@ -322,6 +322,7 @@ static PyObject * py_get_interaction(PyObject *self, PyObject *args)
PyArrayObject *py_triplets;
PyArrayObject *py_bz_grid_addresses;
PyArrayObject *py_D_diag;
PyArrayObject *py_Q;
PyArrayObject *py_shortest_vectors;
PyArrayObject *py_multiplicities;
PyArrayObject *py_fc3;
@ -340,6 +341,7 @@ static PyObject * py_get_interaction(PyObject *self, PyObject *args)
char* g_zero;
long (*bz_grid_addresses)[3];
long *D_diag;
long (*Q)[3];
double *fc3;
double *svecs;
long *multi;
@ -351,7 +353,7 @@ static PyObject * py_get_interaction(PyObject *self, PyObject *args)
long i;
long is_compact_fc3;
if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOOld",
if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOOOld",
&py_fc3_normal_squared,
&py_g_zero,
&py_frequencies,
@ -359,6 +361,7 @@ static PyObject * py_get_interaction(PyObject *self, PyObject *args)
&py_triplets,
&py_bz_grid_addresses,
&py_D_diag,
&py_Q,
&py_fc3,
&py_shortest_vectors,
&py_multiplicities,
@ -382,6 +385,7 @@ static PyObject * py_get_interaction(PyObject *self, PyObject *args)
g_zero = (char*)PyArray_DATA(py_g_zero);
bz_grid_addresses = (long(*)[3])PyArray_DATA(py_bz_grid_addresses);
D_diag = (long*)PyArray_DATA(py_D_diag);
Q = (long(*)[3])PyArray_DATA(py_Q);
fc3 = (double*)PyArray_DATA(py_fc3);
if (PyArray_DIMS(py_fc3)[0] == PyArray_DIMS(py_fc3)[1]) {
is_compact_fc3 = 0;
@ -406,6 +410,7 @@ static PyObject * py_get_interaction(PyObject *self, PyObject *args)
num_triplets,
bz_grid_addresses,
D_diag,
Q,
fc3,
is_compact_fc3,
svecs,

54
c/bzgrid.c
View File

@ -175,18 +175,15 @@ static long get_ir_grid_map(long ir_mapping_table[],
const long PS[3],
const RotMats *rot_reciprocal);
static long get_bz_grid_addresses_type1(BZGrid *bzgrid,
LAGCONST long grid_address[][3]);
const long (*grid_address)[3]);
static long get_bz_grid_addresses_type2(BZGrid *bzgrid,
LAGCONST long grid_address[][3]);
const long (*grid_address)[3]);
static void multiply_matrix_vector_d3(double v[3],
LAGCONST double a[3][3],
const double a[3][3],
const double b[3]);
static void multiply_matrix_vector_ld3(double v[3],
LAGCONST long a[3][3],
const double b[3]);
static double norm_squared_d3(const double a[3]);
static long inverse_unimodular_matrix_l3(long m[3][3],
LAGCONST long a[3][3]);
const long a[3][3]);
long bzg_get_ir_grid_map(long ir_mapping_table[],
const long D_diag[3],
@ -212,7 +209,7 @@ long bzg_get_ir_reciprocal_mesh(long *ir_mapping_table,
const long mesh[3],
const long is_shift[3],
const long is_time_reversal,
LAGCONST long (*rotations_in)[3][3],
const long (*rotations_in)[3][3],
const long num_rot)
{
long i, num_ir;
@ -239,7 +236,7 @@ long bzg_get_ir_reciprocal_mesh(long *ir_mapping_table,
}
long bzg_get_bz_grid_addresses(BZGrid *bzgrid,
LAGCONST long grid_address[][3])
const long (*grid_address)[3])
{
if (bzgrid->type == 1) {
return get_bz_grid_addresses_type1(bzgrid, grid_address);
@ -318,6 +315,19 @@ void bzg_free_RotMats(RotMats * rotmats)
free(rotmats);
}
void bzg_multiply_matrix_vector_ld3(double v[3],
const long a[3][3],
const double b[3])
{
long i;
double c[3];
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++) {
v[i] = c[i];
}
}
/* Return NULL if failed */
static RotMats *get_point_group_reciprocal(const RotMats * rotations,
@ -326,7 +336,7 @@ static RotMats *get_point_group_reciprocal(const RotMats * rotations,
long i, j, num_rot;
RotMats *rot_reciprocal, *rot_return;
long *unique_rot;
LAGCONST long inversion[3][3] = {
const long inversion[3][3] = {
{-1, 0, 0 },
{ 0,-1, 0 },
{ 0, 0,-1 }
@ -420,7 +430,7 @@ static long get_ir_grid_map(long ir_mapping_table[],
}
static long get_bz_grid_addresses_type1(BZGrid *bzgrid,
LAGCONST long grid_address[][3])
const long (*grid_address)[3])
{
double tolerance, min_distance;
double q_vector[3], distance[BZG_NUM_BZ_SEARCH_SPACE];
@ -496,7 +506,7 @@ static long get_bz_grid_addresses_type1(BZGrid *bzgrid,
}
static long get_bz_grid_addresses_type2(BZGrid *bzgrid,
LAGCONST long grid_address[][3])
const long (*grid_address)[3])
{
double tolerance, min_distance;
double q_vec[3], q_red[3], distance[BZG_NUM_BZ_SEARCH_SPACE];
@ -520,7 +530,7 @@ static long get_bz_grid_addresses_type2(BZGrid *bzgrid,
q_red[j] = grid_address[i][j] + bzgrid->PS[j] / 2.0;
q_red[j] /= bzgrid->D_diag[j];
}
multiply_matrix_vector_ld3(q_red, bzgrid->Q, q_red);
bzg_multiply_matrix_vector_ld3(q_red, bzgrid->Q, q_red);
for (j = 0; j < 3; j++) {
nint[j] = lagmat_Nint(q_red[j]);
q_red[j] -= nint[j];
@ -562,7 +572,7 @@ static long get_bz_grid_addresses_type2(BZGrid *bzgrid,
}
static void multiply_matrix_vector_d3(double v[3],
LAGCONST double a[3][3],
const double a[3][3],
const double b[3])
{
long i;
@ -575,27 +585,13 @@ static void multiply_matrix_vector_d3(double v[3],
}
}
static void multiply_matrix_vector_ld3(double v[3],
LAGCONST long a[3][3],
const double b[3])
{
long i;
double c[3];
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++) {
v[i] = c[i];
}
}
static double norm_squared_d3(const double a[3])
{
return a[0] * a[0] + a[1] * a[1] + a[2] * a[2];
}
static long inverse_unimodular_matrix_l3(long m[3][3],
LAGCONST long a[3][3])
const long a[3][3])
{
long det;
long c[3][3];

9
c/bzgrid.h
View File

@ -35,8 +35,6 @@
#ifndef __bzgrid_H__
#define __bzgrid_H__
#include "lagrid.h"
typedef struct {
long size;
long (*mat)[3][3];
@ -104,12 +102,15 @@ long bzg_get_ir_reciprocal_mesh(long *ir_mapping_table,
const long mesh[3],
const long is_shift[3],
const long is_time_reversal,
LAGCONST long (*rotations_in)[3][3],
const long (*rotations_in)[3][3],
const long num_rot);
long bzg_get_bz_grid_addresses(BZGrid *bzgrid,
LAGCONST long grid_address[][3]);
const long grid_address[][3]);
double bzg_get_tolerance_for_BZ_reduction(const BZGrid *bzgrid);
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]);
#endif

2
c/collision_matrix.c
View File

@ -313,7 +313,7 @@ static void get_inv_sinh(double *inv_sinh,
for (i = 0; i < num_band; i++) {
f = frequencies[gp2 * num_band + i];
if (f > cutoff_frequency) {
inv_sinh[i] = inv_sinh_occupation(f, temperature);
inv_sinh[i] = phonoc_inv_sinh_occupation(f, temperature);
} else {
inv_sinh[i] = 0;
}

21
c/fc3.c
View File

@ -33,21 +33,20 @@
/* POSSIBILITY OF SUCH DAMAGE. */
#include <stdlib.h>
#include "phonoc_const.h"
#include "fc3.h"
static void rotate_delta_fc2s(double (*rot_delta_fc2s)[3][3],
const long i_atom,
const long j_atom,
PHPYCONST double (*delta_fc2s)[3][3],
PHPYCONST double (*site_sym_cart)[3][3],
const double (*delta_fc2s)[3][3],
const double (*site_sym_cart)[3][3],
const long *rot_map_sym,
const long num_atom,
const long num_site_sym,
const long num_disp);
static void tensor2_rotation(double rot_tensor[3][3],
PHPYCONST double tensor[3][3],
PHPYCONST double r[3][3]);
const double tensor[3][3],
const double r[3][3]);
static void tensor3_rotation(double *rot_tensor,
const double *tensor,
const double *rot_cartesian);
@ -110,9 +109,9 @@ void fc3_distribute_fc3(double *fc3,
}
void fc3_rotate_delta_fc2(double (*fc3)[3][3][3],
PHPYCONST double (*delta_fc2s)[3][3],
const double (*delta_fc2s)[3][3],
const double *inv_U,
PHPYCONST double (*site_sym_cart)[3][3],
const double (*site_sym_cart)[3][3],
const long *rot_map_syms,
const long num_atom,
const long num_site_sym,
@ -225,8 +224,8 @@ void fc3_transpose_compact_fc3(double * fc3,
static void rotate_delta_fc2s(double (*rot_delta_fc2s)[3][3],
const long i_atom,
const long j_atom,
PHPYCONST double (*delta_fc2s)[3][3],
PHPYCONST double (*site_sym_cart)[3][3],
const double (*delta_fc2s)[3][3],
const double (*site_sym_cart)[3][3],
const long *rot_map_sym,
const long num_atom,
const long num_site_sym,
@ -246,8 +245,8 @@ static void rotate_delta_fc2s(double (*rot_delta_fc2s)[3][3],
}
static void tensor2_rotation(double rot_tensor[3][3],
PHPYCONST double tensor[3][3],
PHPYCONST double r[3][3])
const double tensor[3][3],
const double r[3][3])
{
long i, j, k, l;

4
c/fc3.h
View File

@ -42,9 +42,9 @@ void fc3_distribute_fc3(double *fc3,
const long num_atom,
const double *rot_cart);
void fc3_rotate_delta_fc2(double (*fc3)[3][3][3],
PHPYCONST double (*delta_fc2s)[3][3],
const double (*delta_fc2s)[3][3],
const double *inv_U,
PHPYCONST double (*site_sym_cart)[3][3],
const double (*site_sym_cart)[3][3],
const long *rot_map_syms,
const long num_atom,
const long num_site_sym,

24
c/grgrid.c
View File

@ -61,11 +61,11 @@ static void get_double_grid_address(long address_double[3],
const long address[3],
const long PS[3]);
static long rotate_grid_index(const long grid_index,
LAGCONST long rotation[3][3],
const long rotation[3][3],
const long D_diag[3],
const long PS[3]);
static void get_ir_grid_map(long ir_grid_indices[],
LAGCONST long (*rotations)[3][3],
const long (*rotations)[3][3],
const long num_rot,
const long D_diag[3],
const long PS[3]);
@ -74,7 +74,7 @@ static void get_ir_grid_map(long ir_grid_indices[],
long grg_get_snf3x3(long D_diag[3],
long P[3][3],
long Q[3][3],
LAGCONST long A[3][3])
const long A[3][3])
{
long i, j, succeeded;
long D[3][3];
@ -109,10 +109,10 @@ err:
/* vectors. */
/* num_rot : Number of rotations */
long grg_transform_rotations(long (*transformed_rots)[3][3],
LAGCONST long (*rotations)[3][3],
const long (*rotations)[3][3],
const long num_rot,
const long D_diag[3],
LAGCONST long Q[3][3])
const long Q[3][3])
{
long i, j, k;
double r[3][3], Q_double[3][3];
@ -146,7 +146,7 @@ long grg_transform_rotations(long (*transformed_rots)[3][3],
/* -------------------------------*/
/* address : Single grid address. */
/* D_diag : Diagnal elements of D. */
void grg_get_all_grid_addresses(long grid_address[][3], const long D_diag[3])
void grg_get_all_grid_addresses(long (*grid_address)[3], const long D_diag[3])
{
get_all_grid_addresses(grid_address, D_diag);
}
@ -227,7 +227,7 @@ void grg_get_grid_address_from_index(long address[3],
/* Rotate grid point by index */
/* ---------------------------*/
long grg_rotate_grid_index(const long grid_index,
LAGCONST long rotation[3][3],
const long rotation[3][3],
const long D_diag[3],
const long PS[3])
{
@ -238,7 +238,7 @@ long grg_rotate_grid_index(const long grid_index,
/* Find irreducible grid points */
/* -----------------------------*/
void grg_get_ir_grid_map(long ir_grid_indices[],
LAGCONST long (*rotations)[3][3],
const long (*rotations)[3][3],
const long num_rot,
const long D_diag[3],
const long PS[3])
@ -255,12 +255,12 @@ void grg_get_ir_grid_map(long ir_grid_indices[],
/* included. */
/* Return 0 if failed */
long grg_get_reciprocal_point_group(long rec_rotations[48][3][3],
LAGCONST long (*rotations)[3][3],
const long (*rotations)[3][3],
const long num_rot,
const long is_time_reversal)
{
long i, j, num_rot_ret, inv_exist;
LAGCONST long inversion[3][3] = {
const long inversion[3][3] = {
{-1, 0, 0 },
{ 0,-1, 0 },
{ 0, 0,-1 }
@ -428,7 +428,7 @@ static void get_double_grid_address(long address_double[3],
}
static long rotate_grid_index(const long grid_index,
LAGCONST long rotation[3][3],
const long rotation[3][3],
const long D_diag[3],
const long PS[3])
{
@ -441,7 +441,7 @@ static long rotate_grid_index(const long grid_index,
}
static void get_ir_grid_map(long ir_grid_indices[],
LAGCONST long (*rotations)[3][3],
const long (*rotations)[3][3],
const long num_rot,
const long D_diag[3],
const long PS[3])

15
c/grgrid.h
View File

@ -36,18 +36,17 @@
#define __grgrid_H__
#include <stddef.h>
#include "lagrid.h"
long grg_get_snf3x3(long D_diag[3],
long P[3][3],
long Q[3][3],
LAGCONST long A[3][3]);
const long A[3][3]);
long grg_transform_rotations(long (*transformed_rots)[3][3],
LAGCONST long (*rotations)[3][3],
const long (*rotations)[3][3],
const long num_rot,
const long D_diag[3],
LAGCONST long Q[3][3]);
void grg_get_all_grid_addresses(long grid_address[][3], const long D_diag[3]);
const long Q[3][3]);
void grg_get_all_grid_addresses(long (*grid_address)[3], const long D_diag[3]);
void grg_get_double_grid_address(long address_double[3],
const long address[3],
const long D_diag[3],
@ -65,16 +64,16 @@ void grg_get_grid_address_from_index(long address[3],
const long grid_index,
const long D_diag[3]);
long grg_rotate_grid_index(const long grid_index,
LAGCONST long rotations[3][3],
const long rotations[3][3],
const long D_diag[3],
const long PS[3]);
void grg_get_ir_grid_map(long ir_grid_indices[],
LAGCONST long (*rotations)[3][3],
const long (*rotations)[3][3],
const long num_rot,
const long D_diag[3],
const long PS[3]);
long grg_get_reciprocal_point_group(long rec_rotations[48][3][3],
LAGCONST long (*rotations)[3][3],
const long (*rotations)[3][3],
const long num_rot,
const long is_time_reversal);

8
c/imag_self_energy_with_g.c
View File

@ -35,9 +35,9 @@
#include <stdio.h>
#include <stdlib.h>
#include <stddef.h>
#include "lagrid.h"
#include "phonoc_array.h"
#include "phonoc_utils.h"
#include "phonoc_const.h"
#include "imag_self_energy_with_g.h"
#include "triplet.h"
@ -263,7 +263,7 @@ void ise_imag_self_energy_at_triplet(double *imag_self_energy,
const long triplet_weight,
const double *g1,
const double *g2_3,
PHPYCONST long (*g_pos)[4],
const long (*g_pos)[4],
const long num_g_pos,
const double *temperatures,
const long num_temps,
@ -518,12 +518,12 @@ static void set_occupations(double *n1,
f1 = frequencies[triplet[1] * num_band + j];
f2 = frequencies[triplet[2] * num_band + j];
if (f1 > cutoff_frequency) {
n1[j] = bose_einstein(f1, temperature);
n1[j] = phonoc_bose_einstein(f1, temperature);
} else {
n1[j] = -1;
}
if (f2 > cutoff_frequency) {
n2[j] = bose_einstein(f2, temperature);
n2[j] = phonoc_bose_einstein(f2, temperature);
} else {
n2[j] = -1;
}

2
c/imag_self_energy_with_g.h
View File

@ -71,7 +71,7 @@ void ise_imag_self_energy_at_triplet(double *imag_self_energy,
const long triplet_weight,
const double *g1,
const double *g2_3,
PHPYCONST long (*g_pos)[4],
const long (*g_pos)[4],
const long num_g_pos,
const double *temperatures,
const long num_temps,

41
c/interaction.c
View File

@ -35,10 +35,9 @@
#include <stdio.h>
#include <stdlib.h>
#include "bzgrid.h"
#include "phonoc_array.h"
#include "phonoc_const.h"
#include "interaction.h"
#include "imag_self_energy_with_g.h"
#include "phonoc_array.h"
#include "real_to_reciprocal.h"
#include "reciprocal_to_normal.h"
#include "lapack_wrapper.h"
@ -50,7 +49,7 @@ static const long index_exchange[6][3] = {{0, 1, 2},
{0, 2, 1},
{1, 0, 2}};
static void real_to_normal(double *fc3_normal_squared,
PHPYCONST long (*g_pos)[4],
const long (*g_pos)[4],
const long num_g_pos,
const double *freqs0,
const double *freqs1,
@ -60,7 +59,7 @@ static void real_to_normal(double *fc3_normal_squared,
const lapack_complex_double *eigvecs2,
const double *fc3,
const long is_compact_fc3,
const double q[9], /* q0, q1, q2 */
const double q_vecs[9], /* q0, q1, q2 */
const double *shortest_vectors,
const long svecs_dims[3],
const long *multiplicity,
@ -75,13 +74,13 @@ static void real_to_normal(double *fc3_normal_squared,
const long num_triplets,
const long openmp_at_bands);
static void real_to_normal_sym_q(double *fc3_normal_squared,
PHPYCONST long (*g_pos)[4],
const long (*g_pos)[4],
const long num_g_pos,
PHPYCONST double *freqs[3],
PHPYCONST lapack_complex_double *eigvecs[3],
double * const freqs[3],
lapack_complex_double * const eigvecs[3],
const double *fc3,
const long is_compact_fc3,
const double q[9], /* q0, q1, q2 */
const double q_vecs[9], /* q0, q1, q2 */
const double *shortest_vectors,
const long svecs_dims[3],
const long *multiplicity,
@ -174,7 +173,7 @@ void itr_get_interaction(Darray *fc3_normal_squared,
void itr_get_interaction_at_triplet(double *fc3_normal_squared,
const long num_band0,
const long num_band,
PHPYCONST long (*g_pos)[4],
const long (*g_pos)[4],
const long num_g_pos,
const double *frequencies,
const lapack_complex_double *eigenvectors,
@ -198,11 +197,11 @@ void itr_get_interaction_at_triplet(double *fc3_normal_squared,
long j, k;
double *freqs[3];
lapack_complex_double *eigvecs[3];
double q[9];
double q_vecs[9];
for (j = 0; j < 3; j++) {
for (k = 0; k < 3; k++) {
q[j * 3 + k] = ((double)bzgrid->addresses[triplet[j]][k]) / bzgrid->D_diag[k];
q_vecs[j * 3 + k] = ((double)bzgrid->addresses[triplet[j]][k]) / bzgrid->D_diag[k];
}
}
@ -225,7 +224,7 @@ void itr_get_interaction_at_triplet(double *fc3_normal_squared,
eigvecs,
fc3,
is_compact_fc3,
q, /* q0, q1, q2 */
q_vecs, /* q0, q1, q2 */
shortest_vectors,
svecs_dims,
multiplicity,
@ -257,7 +256,7 @@ void itr_get_interaction_at_triplet(double *fc3_normal_squared,
eigenvectors + triplet[2] * num_band * num_band,
fc3,
is_compact_fc3,
q, /* q0, q1, q2 */
q_vecs, /* q0, q1, q2 */
shortest_vectors,
svecs_dims,
multiplicity,
@ -275,7 +274,7 @@ void itr_get_interaction_at_triplet(double *fc3_normal_squared,
}
static void real_to_normal(double *fc3_normal_squared,
PHPYCONST long (*g_pos)[4],
const long (*g_pos)[4],
const long num_g_pos,
const double *freqs0,
const double *freqs1,
@ -285,7 +284,7 @@ static void real_to_normal(double *fc3_normal_squared,
const lapack_complex_double *eigvecs2,
const double *fc3,
const long is_compact_fc3,
const double q[9], /* q0, q1, q2 */
const double q_vecs[9], /* q0, q1, q2 */
const double *shortest_vectors,
const long svecs_dims[3],
const long *multiplicity,
@ -309,7 +308,7 @@ static void real_to_normal(double *fc3_normal_squared,
(lapack_complex_double*)malloc(sizeof(lapack_complex_double) *
num_patom * num_patom * num_patom * 27);
r2r_real_to_reciprocal(fc3_reciprocal,
q,
q_vecs,
fc3,
is_compact_fc3,
shortest_vectors,
@ -347,13 +346,13 @@ static void real_to_normal(double *fc3_normal_squared,
}
static void real_to_normal_sym_q(double *fc3_normal_squared,
PHPYCONST long (*g_pos)[4],
const long (*g_pos)[4],
const long num_g_pos,
PHPYCONST double *freqs[3],
PHPYCONST lapack_complex_double *eigvecs[3],
double * const freqs[3],
lapack_complex_double * const eigvecs[3],
const double *fc3,
const long is_compact_fc3,
const double q[9], /* q0, q1, q2 */
const double q_vecs[9], /* q0, q1, q2 */
const double *shortest_vectors,
const long svecs_dims[3],
const long *multiplicity,
@ -383,7 +382,7 @@ static void real_to_normal_sym_q(double *fc3_normal_squared,
for (i = 0; i < 6; i++) {
for (j = 0; j < 3; j ++) {
for (k = 0; k < 3; k ++) {
q_ex[j * 3 + k] = q[index_exchange[i][j] * 3 + k];
q_ex[j * 3 + k] = q_vecs[index_exchange[i][j] * 3 + k];
}
}
real_to_normal(fc3_normal_squared_ex,

5
c/interaction.h
View File

@ -35,8 +35,9 @@
#ifndef __interaction_H__
#define __interaction_H__
#include "phonoc_array.h"
#include "bzgrid.h"
#include "lapack_wrapper.h"
#include "phonoc_array.h"
void itr_get_interaction(Darray *fc3_normal_squared,
const char *g_zero,
@ -59,7 +60,7 @@ void itr_get_interaction(Darray *fc3_normal_squared,
void itr_get_interaction_at_triplet(double *fc3_normal_squared,
const long num_band0,
const long num_band,
PHPYCONST long (*g_pos)[4],
const long (*g_pos)[4],
const long num_g_pos,
const double *frequencies,
const lapack_complex_double *eigenvectors,

2
c/isotope.c
View File

@ -87,7 +87,7 @@ iso_get_isotope_scattering_strength(double *gamma,
if (f < cutoff_frequency) {
continue;
}
dist = gaussian(f - f0[i], sigma);
dist = phonoc_gaussian(f - f0[i], sigma);
for (l = 0; l < num_band / 3; l++) { /* elements */
a = 0;
b = 0;

42
c/lagrid.c
View File

@ -34,21 +34,21 @@
#include "lagrid.h"
long lagmat_get_determinant_l3(LAGCONST long a[3][3])
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]);
}
double lagmat_get_determinant_d3(LAGCONST double a[3][3])
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]);
}
void lagmat_cast_matrix_3l_to_3d(double m[3][3], LAGCONST long a[3][3])
void lagmat_cast_matrix_3l_to_3d(double m[3][3], const long a[3][3])
{
m[0][0] = a[0][0];
m[0][1] = a[0][1];
@ -61,7 +61,7 @@ void lagmat_cast_matrix_3l_to_3d(double m[3][3], LAGCONST long a[3][3])
m[2][2] = a[2][2];
}
void lagmat_cast_matrix_3d_to_3l(long m[3][3], LAGCONST double a[3][3])
void lagmat_cast_matrix_3d_to_3l(long m[3][3], const double a[3][3])
{
m[0][0] = lagmat_Nint(a[0][0]);
m[0][1] = lagmat_Nint(a[0][1]);
@ -75,8 +75,8 @@ void lagmat_cast_matrix_3d_to_3l(long m[3][3], LAGCONST double a[3][3])
}
long lagmat_get_similar_matrix_ld3(double m[3][3],
LAGCONST long a[3][3],
LAGCONST double b[3][3],
const long a[3][3],
const double b[3][3],
const double precision)
{
double c[3][3];
@ -89,8 +89,8 @@ long lagmat_get_similar_matrix_ld3(double m[3][3],
return 1;
}
long lagmat_check_identity_matrix_l3(LAGCONST long a[3][3],
LAGCONST long b[3][3])
long lagmat_check_identity_matrix_l3(const long a[3][3],
const long b[3][3])
{
if (a[0][0] - b[0][0] ||
a[0][1] - b[0][1] ||
@ -108,8 +108,8 @@ long lagmat_check_identity_matrix_l3(LAGCONST long a[3][3],
}
}
long lagmat_check_identity_matrix_ld3(LAGCONST long a[3][3],
LAGCONST double b[3][3],
long lagmat_check_identity_matrix_ld3(const long a[3][3],
const double b[3][3],
const double symprec)
{
if (lagmat_Dabs(a[0][0] - b[0][0]) > symprec ||
@ -129,7 +129,7 @@ long lagmat_check_identity_matrix_ld3(LAGCONST long a[3][3],
}
long lagmat_inverse_matrix_d3(double m[3][3],
LAGCONST double a[3][3],
const double a[3][3],
const double precision)
{
double det;
@ -153,7 +153,7 @@ long lagmat_inverse_matrix_d3(double m[3][3],
return 1;
}
void lagmat_transpose_matrix_l3(long a[3][3], LAGCONST long b[3][3])
void lagmat_transpose_matrix_l3(long a[3][3], const long b[3][3])
{
long c[3][3];
c[0][0] = b[0][0];
@ -169,7 +169,7 @@ void lagmat_transpose_matrix_l3(long a[3][3], LAGCONST long b[3][3])
}
void lagmat_multiply_matrix_vector_l3(long v[3],
LAGCONST long a[3][3],
const long a[3][3],
const long b[3])
{
long i;
@ -183,8 +183,8 @@ void lagmat_multiply_matrix_vector_l3(long v[3],
}
void lagmat_multiply_matrix_l3(long m[3][3],
LAGCONST long a[3][3],
LAGCONST long b[3][3])
const long a[3][3],
const long b[3][3])
{
long i, j; /* a_ij */
long c[3][3];
@ -198,8 +198,8 @@ void lagmat_multiply_matrix_l3(long m[3][3],
}
void lagmat_multiply_matrix_ld3(double m[3][3],
LAGCONST long a[3][3],
LAGCONST double b[3][3])
const long a[3][3],
const double b[3][3])
{
long i, j; /* a_ij */
double c[3][3];
@ -213,8 +213,8 @@ void lagmat_multiply_matrix_ld3(double m[3][3],
}
void lagmat_multiply_matrix_d3(double m[3][3],
LAGCONST double a[3][3],
LAGCONST double b[3][3])
const double a[3][3],
const double b[3][3])
{
long i, j; /* a_ij */
double c[3][3];
@ -227,7 +227,7 @@ void lagmat_multiply_matrix_d3(double m[3][3],
lagmat_copy_matrix_d3(m, c);
}
void lagmat_copy_matrix_l3(long a[3][3], LAGCONST long b[3][3])
void lagmat_copy_matrix_l3(long a[3][3], const long b[3][3])
{
a[0][0] = b[0][0];
a[0][1] = b[0][1];
@ -240,7 +240,7 @@ void lagmat_copy_matrix_l3(long a[3][3], LAGCONST long b[3][3])
a[2][2] = b[2][2];
}
void lagmat_copy_matrix_d3(double a[3][3], LAGCONST double b[3][3])
void lagmat_copy_matrix_d3(double a[3][3], const double b[3][3])
{
a[0][0] = b[0][0];
a[0][1] = b[0][1];

46
c/lagrid.h
View File

@ -35,47 +35,43 @@
#ifndef __lagrid_H__
#define __lagrid_H__
#ifndef LAGCONST
#define LAGCONST
#endif
#ifdef LAGWARNING
#define warning_print(...) fprintf(stderr,__VA_ARGS__)
#else
#define warning_print(...)
#endif
long lagmat_get_determinant_l3(LAGCONST long a[3][3]);
double lagmat_get_determinant_d3(LAGCONST double a[3][3]);
void lagmat_cast_matrix_3l_to_3d(double m[3][3], LAGCONST long a[3][3]);
void lagmat_cast_matrix_3d_to_3l(long m[3][3], LAGCONST double a[3][3]);
long lagmat_get_determinant_l3(const long a[3][3]);
double lagmat_get_determinant_d3(const double a[3][3]);
void lagmat_cast_matrix_3l_to_3d(double m[3][3], const long a[3][3]);
void lagmat_cast_matrix_3d_to_3l(long m[3][3], const double a[3][3]);
long lagmat_get_similar_matrix_ld3(double m[3][3],
LAGCONST long a[3][3],
LAGCONST double b[3][3],
const long a[3][3],
const double b[3][3],
const double precision);
long lagmat_check_identity_matrix_l3(LAGCONST long a[3][3],
LAGCONST long b[3][3]);
long lagmat_check_identity_matrix_ld3(LAGCONST long a[3][3],
LAGCONST double b[3][3],
long lagmat_check_identity_matrix_l3(const long a[3][3],
const long b[3][3]);
long lagmat_check_identity_matrix_ld3(const long a[3][3],
const double b[3][3],
const double symprec);
long lagmat_inverse_matrix_d3(double m[3][3],
LAGCONST double a[3][3],
const double a[3][3],
const double precision);
void lagmat_transpose_matrix_l3(long a[3][3], LAGCONST long b[3][3]);
void lagmat_transpose_matrix_l3(long a[3][3], const long b[3][3]);
void lagmat_multiply_matrix_vector_l3(long v[3],
LAGCONST long a[3][3],
const long a[3][3],
const long b[3]);
void lagmat_multiply_matrix_l3(long m[3][3],
LAGCONST long a[3][3],
LAGCONST long b[3][3]);
const long a[3][3],
const long b[3][3]);
void lagmat_multiply_matrix_ld3(double m[3][3],
LAGCONST long a[3][3],
LAGCONST double b[3][3]);
const long a[3][3],
const double b[3][3]);
void lagmat_multiply_matrix_d3(double m[3][3],
LAGCONST double a[3][3],
LAGCONST double b[3][3]);
void lagmat_copy_matrix_l3(long a[3][3], LAGCONST long b[3][3]);
void lagmat_copy_matrix_d3(double a[3][3], LAGCONST double b[3][3]);
const double a[3][3],
const double b[3][3]);
void lagmat_copy_matrix_l3(long a[3][3], const long b[3][3]);
void lagmat_copy_matrix_d3(double a[3][3], const double b[3][3]);
void lagmat_copy_vector_l3(long a[3], const long b[3]);
long lagmat_modulo_l(const long a, const long b);
long lagmat_Nint(const double a);

39
c/phono3py.c
View File

@ -42,6 +42,7 @@
#include "interaction.h"
#include "imag_self_energy_with_g.h"
#include "isotope.h"
#include "lagrid.h"
#include "pp_collision.h"
#include "real_self_energy.h"
#include "grgrid.h"
@ -59,7 +60,8 @@ long ph3py_get_interaction(Darray *fc3_normal_squared,
const long (*triplets)[3],
const long num_triplets,
const long (*bz_grid_addresses)[3],
const long *D_diag,
const long D_diag[3],
const long Q[3][3],
const double *fc3,
const long is_compact_fc3,
const double *shortest_vectors,
@ -73,7 +75,7 @@ long ph3py_get_interaction(Darray *fc3_normal_squared,
const double cutoff_frequency)
{
ConstBZGrid *bzgrid;
long i;
long i, j;
if ((bzgrid = (ConstBZGrid*) malloc(sizeof(ConstBZGrid))) == NULL) {
warning_print("Memory could not be allocated.");
@ -84,6 +86,9 @@ long ph3py_get_interaction(Darray *fc3_normal_squared,
for (i = 0; i < 3; i++) {
bzgrid->D_diag[i] = D_diag[i];
bzgrid->PS[i] = 0;
for (j = 0; j < 3; j++) {
bzgrid->Q[i][j] = Q[i][j];
}
}
itr_get_interaction(fc3_normal_squared,
@ -112,7 +117,7 @@ long ph3py_get_interaction(Darray *fc3_normal_squared,
long ph3py_get_pp_collision(double *imag_self_energy,
PHPYCONST long relative_grid_address[24][4][3], /* thm */
const long relative_grid_address[24][4][3], /* thm */
const double *frequencies,
const lapack_complex_double *eigenvectors,
const long (*triplets)[3],
@ -491,9 +496,9 @@ void ph3py_distribute_fc3(double *fc3,
void ph3py_rotate_delta_fc2(double (*fc3)[3][3][3],
PHPYCONST double (*delta_fc2s)[3][3],
const double (*delta_fc2s)[3][3],
const double *inv_U,
PHPYCONST double (*site_sym_cart)[3][3],
const double (*site_sym_cart)[3][3],
const long *rot_map_syms,
const long num_atom,
const long num_site_sym,
@ -559,7 +564,7 @@ long ph3py_get_triplets_reciprocal_mesh_at_q(long *map_triplets,
const long D_diag[3],
const long is_time_reversal,
const long num_rot,
PHPYCONST long (*rotations)[3][3],
const long (*rotations)[3][3],
const long swappable)
{
return tpl_get_triplets_reciprocal_mesh_at_q(map_triplets,
@ -618,11 +623,11 @@ long ph3py_get_integration_weight(double *iw,
char *iw_zero,
const double *frequency_points,
const long num_band0,
PHPYCONST long relative_grid_address[24][4][3],
const long relative_grid_address[24][4][3],
const long D_diag[3],
PHPYCONST long (*triplets)[3],
const long (*triplets)[3],
const long num_triplets,
PHPYCONST long (*bz_grid_addresses)[3],
const long (*bz_grid_addresses)[3],
const long *bz_map,
const long bz_grid_type,
const double *frequencies1,
@ -676,7 +681,7 @@ void ph3py_get_integration_weight_with_sigma(double *iw,
const double sigma_cutoff,
const double *frequency_points,
const long num_band0,
PHPYCONST long (*triplets)[3],
const long (*triplets)[3],
const long num_triplets,
const double *frequencies,
const long num_band,
@ -714,7 +719,7 @@ long ph3py_get_ir_reciprocal_mesh(long *ir_mapping_table,
const long D_diag[3],
const long is_shift[3],
const long is_time_reversal,
PHPYCONST long (*rotations)[3][3],
const long (*rotations)[3][3],
const long num_rot)
{
long num_ir;
@ -731,11 +736,11 @@ long ph3py_get_ir_reciprocal_mesh(long *ir_mapping_table,
long ph3py_get_bz_grid_address(long (*bz_grid_addresses)[3],
long *bz_map,
long *bzg2grg,
PHPYCONST long grid_address[][3],
const long (*grid_address)[3],
const long D_diag[3],
const long Q[3][3],
const long PS[3],
PHPYCONST double rec_lattice[3][3],
const double rec_lattice[3][3],
const long type)
{
BZGrid *bzgrid;
@ -870,9 +875,9 @@ void ph3py_expand_collision_matrix(double *collision_matrix,
long ph3py_get_neighboring_gird_points(long *relative_grid_points,
const long *grid_points,
PHPYCONST long (*relative_grid_address)[3],
const long (*relative_grid_address)[3],
const long D_diag[3],
PHPYCONST long (*bz_grid_addresses)[3],
const long (*bz_grid_addresses)[3],
const long *bz_map,
const long bz_grid_type,
const long num_grid_points,
@ -916,10 +921,10 @@ long ph3py_set_integration_weights(double *iw,
const long num_band0,
const long num_band,
const long num_gp,
PHPYCONST long (*relative_grid_address)[4][3],
const long (*relative_grid_address)[4][3],
const long D_diag[3],
const long *grid_points,
PHPYCONST long (*bz_grid_addresses)[3],
const long (*bz_grid_addresses)[3],
const long *bz_map,
const long bz_grid_type,
const double *frequencies)

33
c/phono3py.h
View File

@ -53,7 +53,8 @@ long ph3py_get_interaction(Darray *fc3_normal_squared,
const long (*triplets)[3],
const long num_triplets,
const long (*bz_grid_addresses)[3],
const long *mesh,
const long D_diag[3],
const long Q[3][3],
const double *fc3,
const long is_compact_fc3,
const double *shortest_vectors,
@ -66,7 +67,7 @@ long ph3py_get_interaction(Darray *fc3_normal_squared,
const long symmetrize_fc3_q,
const double cutoff_frequency);
long ph3py_get_pp_collision(double *imag_self_energy,
PHPYCONST long relative_grid_address[24][4][3], /* thm */
const long relative_grid_address[24][4][3], /* thm */
const double *frequencies,
const lapack_complex_double *eigenvectors,
const long (*triplets)[3],
@ -219,9 +220,9 @@ void ph3py_distribute_fc3(double *fc3,
const long num_atom,
const double *rot_cart);
void ph3py_rotate_delta_fc2(double (*fc3)[3][3][3],
PHPYCONST double (*delta_fc2s)[3][3],
const double (*delta_fc2s)[3][3],
const double *inv_U,
PHPYCONST double (*site_sym_cart)[3][3],
const double (*site_sym_cart)[3][3],
const long *rot_map_syms,
const long num_atom,
const long num_site_sym,
@ -248,7 +249,7 @@ long ph3py_get_triplets_reciprocal_mesh_at_q(long *map_triplets,
const long mesh[3],
const long is_time_reversal,
const long num_rot,
PHPYCONST long (*rotations)[3][3],
const long (*rotations)[3][3],
const long swappable);
long ph3py_get_BZ_triplets_at_q(long (*triplets)[3],
const long grid_point,
@ -263,11 +264,11 @@ long ph3py_get_integration_weight(double *iw,
char *iw_zero,
const double *frequency_points,
const long num_band0,
PHPYCONST long relative_grid_address[24][4][3],
const long relative_grid_address[24][4][3],
const long mesh[3],
PHPYCONST long (*triplets)[3],
const long (*triplets)[3],
const long num_triplets,
PHPYCONST long (*bz_grid_addresses)[3],
const long (*bz_grid_addresses)[3],
const long *bz_map,
const long bz_grid_type,
const double *frequencies1,
@ -283,7 +284,7 @@ void ph3py_get_integration_weight_with_sigma(double *iw,
const double sigma_cutoff,
const double *frequency_points,
const long num_band0,
PHPYCONST long (*triplets)[3],
const long (*triplets)[3],
const long num_triplets,
const double *frequencies,
const long num_band,
@ -296,16 +297,16 @@ long ph3py_get_ir_reciprocal_mesh(long *ir_mapping_table,
const long D_diag[3],
const long is_shift[3],
const long is_time_reversal,
PHPYCONST long (*rotations)[3][3],
const long (*rotations)[3][3],
const long num_rot);
long ph3py_get_bz_grid_address(long (*bz_grid_addresses)[3],
long *bz_map,
long *bzg2grg,
PHPYCONST long grid_address[][3],
const long (*grid_address)[3],
const long D_diag[3],
const long Q[3][3],
const long PS[3],
PHPYCONST double rec_lattice[3][3],
const double rec_lattice[3][3],
const long type);
void ph3py_symmetrize_collision_matrix(double *collision_matrix,
@ -323,9 +324,9 @@ void ph3py_expand_collision_matrix(double *collision_matrix,
const long num_band);
long ph3py_get_neighboring_gird_points(long *relative_grid_points,
const long *grid_points,
PHPYCONST long (*relative_grid_address)[3],
const long (*relative_grid_address)[3],
const long mesh[3],
PHPYCONST long (*bz_grid_addresses)[3],
const long (*bz_grid_addresses)[3],
const long *bz_map,
const long bz_grid_type,
const long num_grid_points,
@ -335,10 +336,10 @@ long ph3py_set_integration_weights(double *iw,
const long num_band0,
const long num_band,
const long num_gp,
PHPYCONST long (*relative_grid_address)[4][3],
const long (*relative_grid_address)[4][3],
const long mesh[3],
const long *grid_points,
PHPYCONST long (*bz_grid_addresses)[3],
const long (*bz_grid_addresses)[3],
const long *bz_map,
const long bz_grid_type,
const double *frequencies);

7
c/phonoc_const.h
View File

@ -38,11 +38,4 @@
#define M_2PI 6.283185307179586
#define KB 8.6173382568083159E-05
/* PHPYCONST is used instead of 'const' so to avoid gcc warning. */
/* However there should be better way than this way.... */
#ifndef PHPYCONST
#define PHPYCONST
#endif
#endif

6
c/phonoc_utils.c
View File

@ -40,17 +40,17 @@
#define THZTOEVPARKB 47.992398658977166
#define INVSQRT2PI 0.3989422804014327
double bose_einstein(const double x, const double t)
double phonoc_bose_einstein(const double x, const double t)
{
return 1.0 / (exp(THZTOEVPARKB * x / t) - 1);
}
double gaussian(const double x, const double sigma)
double phonoc_gaussian(const double x, const double sigma)
{
return INVSQRT2PI / sigma * exp(-x * x / 2 / sigma / sigma);
}
double inv_sinh_occupation(const double x, const double t)
double phonoc_inv_sinh_occupation(const double x, const double t)
{
return 1.0 / sinh(x * THZTOEVPARKB / 2 / t);
}

6
c/phonoc_utils.h
View File

@ -37,9 +37,9 @@
#include "lapack_wrapper.h"
double bose_einstein(const double x, const double t);
double gaussian(const double x, const double sigma);
double inv_sinh_occupation(const double x, const double t);
double phonoc_bose_einstein(const double x, const double t);
double phonoc_gaussian(const double x, const double sigma);
double phonoc_inv_sinh_occupation(const double x, const double t);
lapack_complex_double
phonoc_complex_prod(const lapack_complex_double a,
const lapack_complex_double b);

9
c/pp_collision.c
View File

@ -34,12 +34,11 @@
#include <stdio.h>
#include <stdlib.h>
#include "phonoc_array.h"
#include "phonoc_const.h"
#include "phonoc_utils.h"
#include "imag_self_energy_with_g.h"
#include "pp_collision.h"
#include "interaction.h"
#include "phonoc_array.h"
#include "phonoc_utils.h"
#include "pp_collision.h"
#include "triplet.h"
#include "triplet_iw.h"
#include "lapack_wrapper.h"
@ -78,7 +77,7 @@ static void finalize_ise(double *imag_self_energy,
const long is_NU);
void ppc_get_pp_collision(double *imag_self_energy,
PHPYCONST long relative_grid_address[24][4][3], /* thm */
const long relative_grid_address[24][4][3], /* thm */
const double *frequencies,
const lapack_complex_double *eigenvectors,
const long (*triplets)[3],

5
c/pp_collision.h
View File

@ -35,13 +35,12 @@
#ifndef __pp_collision_H__
#define __pp_collision_H__
#include "phonoc_array.h"
#include "phonoc_const.h"
#include "bzgrid.h"
#include "lapack_wrapper.h"
#include "phonoc_array.h"
void ppc_get_pp_collision(double *imag_self_energy,
PHPYCONST long relative_grid_address[24][4][3], /* thm */
const long relative_grid_address[24][4][3], /* thm */
const double *frequencies,
const lapack_complex_double *eigenvectors,
const long (*triplets)[3],

4
c/real_self_energy.c
View File

@ -210,10 +210,10 @@ static double sum_real_self_energy_at_band(const long num_band,
shift = 0;
for (i = 0; i < num_band; i++) {
if (freqs1[i] > cutoff_frequency) {
n1 = bose_einstein(freqs1[i], temperature);
n1 = phonoc_bose_einstein(freqs1[i], temperature);
for (j = 0; j < num_band; j++) {
if (freqs2[j] > cutoff_frequency) {
n2 = bose_einstein(freqs2[j], temperature);
n2 = phonoc_bose_einstein(freqs2[j], temperature);
f1 = fpoint + freqs1[i] + freqs2[j];
f2 = fpoint - freqs1[i] - freqs2[j];
f3 = fpoint - freqs1[i] + freqs2[j];

36
c/real_to_reciprocal.c
View File

@ -44,7 +44,7 @@
static void
real_to_reciprocal_single_thread(lapack_complex_double *fc3_reciprocal,
const double q[9],
const double q_vecs[9],
const double *fc3,
const long is_compact_fc3,
const double *shortest_vectors,
@ -54,7 +54,7 @@ real_to_reciprocal_single_thread(lapack_complex_double *fc3_reciprocal,
const long *s2p_map);
static void
real_to_reciprocal_openmp(lapack_complex_double *fc3_reciprocal,
const double q[9],
const double q_vecs[9],
const double *fc3,
const long is_compact_fc3,
const double *shortest_vectors,
@ -63,7 +63,7 @@ real_to_reciprocal_openmp(lapack_complex_double *fc3_reciprocal,
const long *p2s_map,
const long *s2p_map);
static void real_to_reciprocal_elements(lapack_complex_double *fc3_rec_elem,
const double q[9],
const double q_vecs[9],
const double *fc3,
const long is_compact_fc3,
const double *shortest_vectors,
@ -80,7 +80,7 @@ static lapack_complex_double get_phase_factor(const double q[],
const long multi);
static lapack_complex_double
get_pre_phase_factor(const long i,
const double q[9],
const double q_vecs[9],
const double *shortest_vectors,
const long svecs_dims[3],
const long *multiplicity,
@ -88,7 +88,7 @@ get_pre_phase_factor(const long i,
/* fc3_reciprocal[num_patom, num_patom, num_patom, 3, 3, 3] */
void r2r_real_to_reciprocal(lapack_complex_double *fc3_reciprocal,
const double q[9],
const double q_vecs[9],
const double *fc3,
const long is_compact_fc3,
const double *shortest_vectors,
@ -100,7 +100,7 @@ void r2r_real_to_reciprocal(lapack_complex_double *fc3_reciprocal,
{
if (openmp_at_bands) {
real_to_reciprocal_openmp(fc3_reciprocal,
q,
q_vecs,
fc3,
is_compact_fc3,
shortest_vectors,
@ -110,7 +110,7 @@ void r2r_real_to_reciprocal(lapack_complex_double *fc3_reciprocal,
s2p_map);
} else {
real_to_reciprocal_single_thread(fc3_reciprocal,
q,
q_vecs,
fc3,
is_compact_fc3,
shortest_vectors,
@ -124,7 +124,7 @@ void r2r_real_to_reciprocal(lapack_complex_double *fc3_reciprocal,
static void
real_to_reciprocal_single_thread(lapack_complex_double *fc3_reciprocal,
const double q[9],
const double q_vecs[9],
const double *fc3,
const long is_compact_fc3,
const double *shortest_vectors,
@ -146,7 +146,7 @@ real_to_reciprocal_single_thread(lapack_complex_double *fc3_reciprocal,
i * 27 * num_patom * num_patom +
j * 27 * num_patom +
k * 27,
q,
q_vecs,
fc3,
is_compact_fc3,
shortest_vectors,
@ -159,7 +159,7 @@ real_to_reciprocal_single_thread(lapack_complex_double *fc3_reciprocal,
}
}
pre_phase_factor = get_pre_phase_factor(
i, q, shortest_vectors, svecs_dims, multiplicity, p2s_map);
i, q_vecs, shortest_vectors, svecs_dims, multiplicity, p2s_map);
adrs_shift = i * num_patom * num_patom * 27;
for (j = 0; j < num_patom * num_patom * 27; j++) {
fc3_reciprocal[adrs_shift + j] =
@ -170,7 +170,7 @@ real_to_reciprocal_single_thread(lapack_complex_double *fc3_reciprocal,
static void
real_to_reciprocal_openmp(lapack_complex_double *fc3_reciprocal,
const double q[9],
const double q_vecs[9],
const double *fc3,
const long is_compact_fc3,
const double *shortest_vectors,
@ -194,7 +194,7 @@ real_to_reciprocal_openmp(lapack_complex_double *fc3_reciprocal,
i * 27 * num_patom * num_patom +
j * 27 * num_patom +
k * 27,
q,
q_vecs,
fc3,
is_compact_fc3,
shortest_vectors,
@ -206,7 +206,7 @@ real_to_reciprocal_openmp(lapack_complex_double *fc3_reciprocal,
}
pre_phase_factor = get_pre_phase_factor(
i, q, shortest_vectors, svecs_dims, multiplicity, p2s_map);
i, q_vecs, shortest_vectors, svecs_dims, multiplicity, p2s_map);
adrs_shift = i * num_patom * num_patom * 27;
#pragma omp parallel for
for (j = 0; j < num_patom * num_patom * 27; j++) {
@ -217,7 +217,7 @@ real_to_reciprocal_openmp(lapack_complex_double *fc3_reciprocal,
}
static void real_to_reciprocal_elements(lapack_complex_double *fc3_rec_elem,
const double q[9],
const double q_vecs[9],
const double *fc3,
const long is_compact_fc3,
const double *shortest_vectors,
@ -253,7 +253,7 @@ static void real_to_reciprocal_elements(lapack_complex_double *fc3_rec_elem,
}
adrs_vec1 = j * svecs_dims[1] + pi0;
phase_factor1 = get_phase_factor(q,
phase_factor1 = get_phase_factor(q_vecs,
1,
shortest_vectors +
adrs_vec1 * svecs_dims[2] * 3,
@ -263,7 +263,7 @@ static void real_to_reciprocal_elements(lapack_complex_double *fc3_rec_elem,
continue;
}
adrs_vec2 = k * svecs_dims[1] + pi0;
phase_factor2 = get_phase_factor(q,
phase_factor2 = get_phase_factor(q_vecs,
2,
shortest_vectors +
adrs_vec2 * svecs_dims[2] * 3,
@ -287,7 +287,7 @@ static void real_to_reciprocal_elements(lapack_complex_double *fc3_rec_elem,
static lapack_complex_double
get_pre_phase_factor(const long i,
const double q[9],
const double q_vecs[9],
const double *shortest_vectors,
const long svecs_dims[3],
const long *multiplicity,
@ -301,7 +301,7 @@ get_pre_phase_factor(const long i,
for (j = 0; j < 3; j++) {
pre_phase += shortest_vectors[
p2s_map[i] * svecs_dims[1] * svecs_dims[2] * 3 + j] *
(q[j] + q[3 + j] + q[6 + j]);
(q_vecs[j] + q_vecs[3 + j] + q_vecs[6 + j]);
}
assert(multiplicity[p2s_map[i] * svecs_dims[1]] == 1);

2
c/real_to_reciprocal.h
View File

@ -39,7 +39,7 @@
#include "lapack_wrapper.h"
void r2r_real_to_reciprocal(lapack_complex_double *fc3_reciprocal,
const double q[9],
const double q_vecs[9],
const double *fc3,
const long is_compact_fc3,
const double *shortest_vectors,

4
c/reciprocal_to_normal.c
View File

@ -35,8 +35,6 @@
#include <stdlib.h>
#include <math.h>
#include "phonoc_utils.h"
#include "phonoc_const.h"
#include "phonoc_array.h"
#include "reciprocal_to_normal.h"
#include "lapack_wrapper.h"
@ -73,7 +71,7 @@ static double get_fc3_sum
void reciprocal_to_normal_squared
(double *fc3_normal_squared,
PHPYCONST long (*g_pos)[4],
const long (*g_pos)[4],
const long num_g_pos,
const lapack_complex_double *fc3_reciprocal,
const double *freqs0,

4
c/reciprocal_to_normal.h
View File

@ -35,13 +35,11 @@
#ifndef __reciprocal_to_normal_H__
#define __reciprocal_to_normal_H__
#include "phonoc_array.h"
#include "phonoc_const.h"
#include "lapack_wrapper.h"
void reciprocal_to_normal_squared
(double *fc3_normal_squared,
PHPYCONST long (*g_pos)[4],
const long (*g_pos)[4],
const long num_g_pos,
const lapack_complex_double *fc3_reciprocal,
const double *freqs0,

24
c/tetrahedron_method.c
View File

@ -178,7 +178,7 @@ static long db_relative_grid_address[4][24][4][3] = {
static double
get_integration_weight(const double omega,
THMCONST double tetrahedra_omegas[24][4],
const double tetrahedra_omegas[24][4],
double (*gn)(const long,
const double,
const double[4]),
@ -186,11 +186,11 @@ get_integration_weight(const double omega,
const long,
const double,
const double[4]));
static long get_main_diagonal(THMCONST double rec_lattice[3][3]);
static long get_main_diagonal(const double rec_lattice[3][3]);
static long sort_omegas(double v[4]);
static double norm_squared_d3(const double a[3]);
static void multiply_matrix_vector_di3(double v[3],
THMCONST double a[3][3],
static void multiply_matrix_vector_dl3(double v[3],
const double a[3][3],
const long b[3]);
static double _f(const long n,
const long m,
@ -281,7 +281,7 @@ static double _I_4(void);
void thm_get_relative_grid_address(long relative_grid_address[24][4][3],
THMCONST double rec_lattice[3][3])
const double rec_lattice[3][3])
{
long i, j, k, main_diag_index;
@ -314,7 +314,7 @@ void thm_get_all_relative_grid_address(long relative_grid_address[4][24][4][3])
}
double thm_get_integration_weight(const double omega,
THMCONST double tetrahedra_omegas[24][4],
const double tetrahedra_omegas[24][4],
const char function)
{
if (function == 'I') {
@ -330,7 +330,7 @@ double thm_get_integration_weight(const double omega,
static double
get_integration_weight(const double omega,
THMCONST double tetrahedra_omegas[24][4],
const double tetrahedra_omegas[24][4],
double (*gn)(const long,
const double,
const double[4]),
@ -442,17 +442,17 @@ static long sort_omegas(double v[4])
return i;
}
static long get_main_diagonal(THMCONST double rec_lattice[3][3])
static long get_main_diagonal(const double rec_lattice[3][3])
{
long i, shortest;
double length, min_length;
double main_diag[3];
shortest = 0;
multiply_matrix_vector_di3(main_diag, rec_lattice, main_diagonals[0]);
multiply_matrix_vector_dl3(main_diag, rec_lattice, main_diagonals[0]);
min_length = norm_squared_d3(main_diag);
for (i = 1; i < 4; i++) {
multiply_matrix_vector_di3(main_diag, rec_lattice, main_diagonals[i]);
multiply_matrix_vector_dl3(main_diag, rec_lattice, main_diagonals[i]);
length = norm_squared_d3(main_diag);
if (min_length > length) {
min_length = length;
@ -467,8 +467,8 @@ static double norm_squared_d3(const double a[3])
return a[0] * a[0] + a[1] * a[1] + a[2] * a[2];
}
static void multiply_matrix_vector_di3(double v[3],
THMCONST double a[3][3],
static void multiply_matrix_vector_dl3(double v[3],
const double a[3][3],
const long b[3])
{
long i;

8
c/tetrahedron_method.h
View File

@ -35,17 +35,13 @@
#ifndef __tetrahedron_method_H__
#define __tetrahedron_method_H__
#ifndef THMCONST
#define THMCONST
#endif
#include <stddef.h>
void thm_get_relative_grid_address(long relative_grid_address[24][4][3],
THMCONST double rec_lattice[3][3]);
const double rec_lattice[3][3]);
void thm_get_all_relative_grid_address(long relative_grid_address[4][24][4][3]);
double thm_get_integration_weight(const double omega,
THMCONST double tetrahedra_omegas[24][4],
const double tetrahedra_omegas[24][4],
const char function);
#endif

11
c/triplet.c
View File

@ -35,7 +35,6 @@
/* POSSIBILITY OF SUCH DAMAGE. */
#include "bzgrid.h"
#include "lagrid.h"
#include "triplet.h"
#include "triplet_iw.h"
#include "triplet_grid.h"
@ -57,7 +56,7 @@ long tpl_get_triplets_reciprocal_mesh_at_q(long *map_triplets,
const long mesh[3],
const long is_time_reversal,
const long num_rot,
LAGCONST long (*rotations)[3][3],
const long (*rotations)[3][3],
const long swappable)
{
long num_ir;
@ -77,8 +76,8 @@ void tpl_get_integration_weight(double *iw,
char *iw_zero,
const double *frequency_points,
const long num_band0,
LAGCONST long relative_grid_address[24][4][3],
LAGCONST long (*triplets)[3],
const long relative_grid_address[24][4][3],
const long (*triplets)[3],
const long num_triplets,
const ConstBZGrid *bzgrid,
const double *frequencies1,
@ -123,7 +122,7 @@ void tpl_get_integration_weight_with_sigma(double *iw,
const double sigma_cutoff,
const double *frequency_points,
const long num_band0,
LAGCONST long (*triplets)[3],
const long (*triplets)[3],
const long num_triplets,
const double *frequencies,
const long num_band,
@ -175,7 +174,7 @@ 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],
LAGCONST long relative_grid_address[24][4][3],
const long relative_grid_address[24][4][3],
const long tp_type)
{
long i, j, k, l;

11
c/triplet.h
View File

@ -39,7 +39,6 @@
#include <stddef.h>
#include "bzgrid.h"
#include "lagrid.h"
/* Irreducible triplets of k-points are searched under conservation of */
/* :math:``\mathbf{k}_1 + \mathbf{k}_2 + \mathbf{k}_3 = \mathbf{G}``. */
@ -53,7 +52,7 @@ long tpl_get_triplets_reciprocal_mesh_at_q(long *map_triplets,
const long mesh[3],
const long is_time_reversal,
const long num_rot,
LAGCONST long (*rotations)[3][3],
const long (*rotations)[3][3],
const long swappable);
/* Irreducible grid-point-triplets in BZ are stored. */
/* triplets are recovered from grid_point and triplet_weights. */
@ -69,8 +68,8 @@ void tpl_get_integration_weight(double *iw,
char *iw_zero,
const double *frequency_points,
const long num_band0,
LAGCONST long relative_grid_address[24][4][3],
LAGCONST long (*triplets)[3],
const long relative_grid_address[24][4][3],
const long (*triplets)[3],
const long num_triplets,
const ConstBZGrid *bzgrid,
const double *frequencies1,
@ -86,7 +85,7 @@ void tpl_get_integration_weight_with_sigma(double *iw,
const double sigma_cutoff,
const double *frequency_points,
const long num_band0,
LAGCONST long (*triplets)[3],
const long (*triplets)[3],
const long num_triplets,
const double *frequencies,
const long num_band,
@ -95,7 +94,7 @@ void tpl_get_integration_weight_with_sigma(double *iw,
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],
LAGCONST long relative_grid_address[24][4][3],
const long relative_grid_address[24][4][3],
const long tp_type);
#endif

2
c/triplet_grid.c
View File

@ -207,7 +207,7 @@ long tpk_get_ir_triplets_at_q(long *map_triplets,
const long grid_point,
const long D_diag[3],
const long is_time_reversal,
LAGCONST long (*rotations)[3][3],
const long (*rotations)[3][3],
const long num_rot,
const long swappable)
{

2
c/triplet_grid.h
View File

@ -44,7 +44,7 @@ long tpk_get_ir_triplets_at_q(long *map_triplets,
const long grid_point,
const long D_diag[3],
const long is_time_reversal,
LAGCONST long (*rotations)[3][3],
const long (*rotations)[3][3],
const long num_rot,
const long swappable);
long tpk_get_BZ_triplets_at_q(long (*triplets)[3],

37
c/triplet_iw.c
View File

@ -34,7 +34,6 @@
#include <math.h>
#include "grgrid.h"
#include "lagrid.h"
#include "phonoc_utils.h"
#include "triplet.h"
#include "triplet_iw.h"
@ -43,7 +42,7 @@
static void set_freq_vertices(double freq_vertices[3][24][4],
const double *frequencies1,
const double *frequencies2,
LAGCONST long vertices[2][24][4],
const long vertices[2][24][4],
const long num_band1,
const long num_band2,
const long b1,
@ -51,24 +50,24 @@ static void set_freq_vertices(double freq_vertices[3][24][4],
const long tp_type);
static long set_g(double g[3],
const double f0,
LAGCONST double freq_vertices[3][24][4],
const double freq_vertices[3][24][4],
const long max_i);
static long in_tetrahedra(const double f0, LAGCONST double freq_vertices[24][4]);
static long in_tetrahedra(const double f0, const double freq_vertices[24][4]);
static void get_triplet_tetrahedra_vertices(
long vertices[2][24][4],
LAGCONST long tp_relative_grid_address[2][24][4][3],
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,
LAGCONST long relative_grid_address[][3],
const long (*relative_grid_address)[3],
const long num_relative_grid_address,
const ConstBZGrid *bzgrid);
static void
get_neighboring_grid_points_type2(long neighboring_grid_points[],
const long grid_point,
LAGCONST long relative_grid_address[][3],
const long (*relative_grid_address)[3],
const long num_relative_grid_address,
const ConstBZGrid *bzgrid);
@ -77,7 +76,7 @@ tpi_get_integration_weight(double *iw,
char *iw_zero,
const double *frequency_points,
const long num_band0,
LAGCONST long tp_relative_grid_address[2][24][4][3],
const long tp_relative_grid_address[2][24][4][3],
const long triplets[3],
const long num_triplets,
const ConstBZGrid *bzgrid,
@ -183,9 +182,9 @@ void tpi_get_integration_weight_with_sigma(double *iw,
g2 = 0;
} else {
iw_zero[adrs_shift] = 0;
g0 = gaussian(f0 + f1 - f2, sigma);
g1 = gaussian(f0 - f1 + f2, sigma);
g2 = gaussian(f0 - f1 - f2, sigma);
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) {
iw[adrs_shift] = g2;
@ -206,7 +205,7 @@ void tpi_get_integration_weight_with_sigma(double *iw,
iw[adrs_shift] = 0;
} else {
iw_zero[adrs_shift] = 0;
iw[adrs_shift] = gaussian(f0 + f1 - f2, sigma);
iw[adrs_shift] = phonoc_gaussian(f0 + f1 - f2, sigma);
}
}
}
@ -218,7 +217,7 @@ void tpi_get_integration_weight_with_sigma(double *iw,
void
tpi_get_neighboring_grid_points(long neighboring_grid_points[],
const long grid_point,
LAGCONST long relative_grid_address[][3],
const long (*relative_grid_address)[3],
const long num_relative_grid_address,
const ConstBZGrid *bzgrid)
{
@ -240,7 +239,7 @@ tpi_get_neighboring_grid_points(long neighboring_grid_points[],
static void set_freq_vertices(double freq_vertices[3][24][4],
const double *frequencies1,
const double *frequencies2,
LAGCONST long vertices[2][24][4],
const long vertices[2][24][4],
const long num_band1,
const long num_band2,
const long b1,
@ -278,7 +277,7 @@ static void set_freq_vertices(double freq_vertices[3][24][4],
/* calculation. */
static long set_g(double g[3],
const double f0,
LAGCONST double freq_vertices[3][24][4],
const double freq_vertices[3][24][4],
const long max_i)
{
long i, iw_zero;
@ -297,7 +296,7 @@ static long set_g(double g[3],
return iw_zero;
}
static long in_tetrahedra(const double f0, LAGCONST double freq_vertices[24][4])
static long in_tetrahedra(const double f0, const double freq_vertices[24][4])
{
long i, j;
double fmin, fmax;
@ -324,7 +323,7 @@ static long in_tetrahedra(const double f0, LAGCONST double freq_vertices[24][4])
}
static void get_triplet_tetrahedra_vertices(long vertices[2][24][4],
LAGCONST long tp_relative_grid_address[2][24][4][3],
const long tp_relative_grid_address[2][24][4][3],
const long triplet[3],
const ConstBZGrid *bzgrid)
{
@ -344,7 +343,7 @@ static void get_triplet_tetrahedra_vertices(long vertices[2][24][4],
static void
get_neighboring_grid_points_type1(long neighboring_grid_points[],
const long grid_point,
LAGCONST long relative_grid_address[][3],
const long (*relative_grid_address)[3],
const long num_relative_grid_address,
const ConstBZGrid *bzgrid)
{
@ -373,7 +372,7 @@ get_neighboring_grid_points_type1(long neighboring_grid_points[],
static void
get_neighboring_grid_points_type2(long neighboring_grid_points[],
const long grid_point,
LAGCONST long relative_grid_address[][3],
const long (*relative_grid_address)[3],
const long num_relative_grid_address,
const ConstBZGrid *bzgrid)
{

6
c/triplet_iw.h
View File

@ -35,14 +35,12 @@
#ifndef __triplet_iw_H__
#define __triplet_iw_H__
#include "lagrid.h"
void
tpi_get_integration_weight(double *iw,
char *iw_zero,
const double *frequency_points,
const long num_band0,
LAGCONST long tp_relative_grid_address[2][24][4][3],
const long tp_relative_grid_address[2][24][4][3],
const long triplets[3],
const long num_triplets,
const ConstBZGrid *bzgrid,
@ -67,7 +65,7 @@ void tpi_get_integration_weight_with_sigma(double *iw,
void
tpi_get_neighboring_grid_points(long neighboring_grid_points[],
const long grid_point,
LAGCONST long relative_grid_address[][3],
const long (*relative_grid_address)[3],
const long num_relative_grid_address,
const ConstBZGrid *bzgrid);

3
phono3py/phonon3/interaction.py
View File

@ -459,7 +459,8 @@ class Interaction(object):
self._eigenvectors,
self._triplets_at_q,
self._bz_grid.addresses,
self._mesh,
self._bz_grid.D_diag,
self._bz_grid.Q,
self._fc3,
self._smallest_vectors,
self._multiplicity,

54
phono3py/phonon3/triplets.py
View File

@ -119,10 +119,64 @@ class BZGrid(object):
self._addresses = None
self._gp_map = None
self._Q = np.eye(3, dtype='int_')
self._P = np.eye(3, dtype='int_')
@property
def mesh_numbers(self):
"""Mesh numbers of conventional regular grid"""
return self._mesh_numbers
@property
def D_diag(self):
"""Diagonal elements of diagonal matrix after SNF: D=PAQ
This corresponds to the mesh numbers in transformed reciprocal
basis vectors. q-points with respect to the original recirpocal
basis vectors are given by
q = np.dot(Q, addresses[gp] / D_diag.astype('double'))
for the Gamma cetnred grid. With shifted, where only half grid shifts
that satisfy the symmetry are considered,
q = np.dot(Q, (addresses[gp] + np.dot(P, s)) / D_diag.astype('double'))
where s is the shift vectors that are 0 or 1/2. But it is more
convenient to use the integer shift vectors S by 0 or 1, which gives
q = (np.dot(Q, (2 * addresses[gp] + np.dot(P, S))
/ D_diag.astype('double'))) / 2
and this is the definition of PS in this class.
"""
return self._mesh_numbers
@property
def P(self):
"""Left unimodular matrix after SNF: D=PAQ"""
return self._P
@property
def Q(self):
"""Right unimodular matrix after SNF: D=PAQ"""
return self._Q
@property
def QDinv(self):
"""QD^-1"""
return np.array(self.Q * (1 / self.D_diag.astype('double')),
dtype='double', order='C')
@property
def PS(self):
"""Integer shift vectors of GRGrid"""
if self._is_shift is None:
return np.zeros(3, dtype='int_')
else:
return np.array(np.dot(self.P, self._is_shift), dtype='int_')
@property
def addresses(self):
return self._addresses