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

Moved algebraic functions to lagrid.c

This commit is contained in:
Atsushi Togo 2021-03-05 12:07:54 +09:00
parent e5f5231d92
commit 9726060832
15 changed files with 494 additions and 502 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
CMakeLists.txt
View File

@ -50,6 +50,7 @@ set(SOURCES_PHONO3PY
${PROJECT_SOURCE_DIR}/c/imag_self_energy_with_g.c
${PROJECT_SOURCE_DIR}/c/interaction.c
${PROJECT_SOURCE_DIR}/c/isotope.c
${PROJECT_SOURCE_DIR}/c/lagrid.c
${PROJECT_SOURCE_DIR}/c/lapack_wrapper.c
${PROJECT_SOURCE_DIR}/c/phono3py.c
${PROJECT_SOURCE_DIR}/c/phonoc_utils.c

141
c/bzgrid.c
View File

@ -37,13 +37,7 @@
#include <stddef.h>
#include "bzgrid.h"
#include "grgrid.h"
#ifdef BZGWARNING
#include <stdio.h>
#define warning_print(...) fprintf(stderr,__VA_ARGS__)
#else
#define warning_print(...)
#endif
#include "lagrid.h"
#define BZG_NUM_BZ_SEARCH_SPACE 125
static long bz_search_space[BZG_NUM_BZ_SEARCH_SPACE][3] = {
@ -193,29 +187,24 @@ static long get_ir_reciprocal_mesh_distortion(long grid_address[][3],
const MatLONG *rot_reciprocal);
static long relocate_BZ_grid_address(long bz_grid_address[][3],
long bz_map[],
GRGCONST long grid_address[][3],
LAGCONST long grid_address[][3],
const long mesh[3],
GRGCONST double rec_lattice[3][3],
LAGCONST double rec_lattice[3][3],
const long is_shift[3]);
static long get_bz_grid_addresses(long bz_grid_address[][3],
long bz_map[][2],
GRGCONST long grid_address[][3],
LAGCONST long grid_address[][3],
const long mesh[3],
GRGCONST double rec_lattice[3][3],
LAGCONST double rec_lattice[3][3],
const long is_shift[3]);
static double get_tolerance_for_BZ_reduction(GRGCONST double rec_lattice[3][3],
static double get_tolerance_for_BZ_reduction(LAGCONST double rec_lattice[3][3],
const long mesh[3]);
static long get_num_ir(long ir_mapping_table[], const long mesh[3]);
static long check_mesh_symmetry(const long mesh[3],
const long is_shift[3],
const MatLONG *rot_reciprocal);
static void transpose_matrix_l3(long a[3][3], GRGCONST long b[3][3]);
static void multiply_matrix_l3(long m[3][3],
GRGCONST long a[3][3], GRGCONST long b[3][3]);
static long check_identity_matrix_l3(GRGCONST long a[3][3],
GRGCONST long b[3][3]);
static void multiply_matrix_vector_d3(double v[3],
GRGCONST double a[3][3],
LAGCONST double a[3][3],
const double b[3]);
static double norm_squared_d3(const double a[3]);
@ -268,9 +257,9 @@ long bzg_get_ir_reciprocal_mesh(long grid_address[][3],
long bzg_relocate_BZ_grid_address(long bz_grid_address[][3],
long bz_map[],
GRGCONST long grid_address[][3],
LAGCONST long grid_address[][3],
const long mesh[3],
GRGCONST double rec_lattice[3][3],
LAGCONST double rec_lattice[3][3],
const long is_shift[3])
{
return relocate_BZ_grid_address(bz_grid_address,
@ -283,9 +272,9 @@ long bzg_relocate_BZ_grid_address(long bz_grid_address[][3],
long bzg_get_bz_grid_addresses(long bz_grid_address[][3],
long bz_map[][2],
GRGCONST long grid_address[][3],
LAGCONST long grid_address[][3],
const long mesh[3],
GRGCONST double rec_lattice[3][3],
LAGCONST double rec_lattice[3][3],
const long is_shift[3])
{
return get_bz_grid_addresses(bz_grid_address,
@ -296,31 +285,6 @@ long bzg_get_bz_grid_addresses(long bz_grid_address[][3],
is_shift);
}
void bzg_copy_matrix_l3(long a[3][3], GRGCONST long b[3][3])
{
a[0][0] = b[0][0];
a[0][1] = b[0][1];
a[0][2] = b[0][2];
a[1][0] = b[1][0];
a[1][1] = b[1][1];
a[1][2] = b[1][2];
a[2][0] = b[2][0];
a[2][1] = b[2][1];
a[2][2] = b[2][2];
}
void bzg_multiply_matrix_vector_l3(long v[3],
GRGCONST long a[3][3],
const long b[3])
{
long i;
long 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];
}
MatLONG * bzg_alloc_MatLONG(const long size)
{
MatLONG *matlong;
@ -363,7 +327,7 @@ static MatLONG *get_point_group_reciprocal(const MatLONG * rotations,
long i, j, num_rot;
MatLONG *rot_reciprocal, *rot_return;
long *unique_rot;
GRGCONST long inversion[3][3] = {
LAGCONST long inversion[3][3] = {
{-1, 0, 0 },
{ 0,-1, 0 },
{ 0, 0,-1 }
@ -395,20 +359,20 @@ static MatLONG *get_point_group_reciprocal(const MatLONG * rotations,
}
for (i = 0; i < rotations->size; i++) {
transpose_matrix_l3(rot_reciprocal->mat[i], rotations->mat[i]);
lagmat_transpose_matrix_l3(rot_reciprocal->mat[i], rotations->mat[i]);
if (is_time_reversal) {
multiply_matrix_l3(rot_reciprocal->mat[rotations->size+i],
inversion,
rot_reciprocal->mat[i]);
lagmat_multiply_matrix_l3(rot_reciprocal->mat[rotations->size+i],
inversion,
rot_reciprocal->mat[i]);
}
}
num_rot = 0;
for (i = 0; i < rot_reciprocal->size; i++) {
for (j = 0; j < num_rot; j++) {
if (check_identity_matrix_l3(rot_reciprocal->mat[unique_rot[j]],
rot_reciprocal->mat[i])) {
if (lagmat_check_identity_matrix_l3(rot_reciprocal->mat[unique_rot[j]],
rot_reciprocal->mat[i])) {
goto escape;
}
}
@ -420,7 +384,7 @@ static MatLONG *get_point_group_reciprocal(const MatLONG * rotations,
if ((rot_return = bzg_alloc_MatLONG(num_rot)) != NULL) {
for (i = 0; i < num_rot; i++) {
bzg_copy_matrix_l3(rot_return->mat[i], rot_reciprocal->mat[unique_rot[i]]);
lagmat_copy_matrix_l3(rot_return->mat[i], rot_reciprocal->mat[unique_rot[i]]);
}
}
@ -481,9 +445,9 @@ static long get_ir_reciprocal_mesh_normal(long grid_address[][3],
is_shift);
ir_mapping_table[i] = i;
for (j = 0; j < rot_reciprocal->size; j++) {
bzg_multiply_matrix_vector_l3(address_double_rot,
rot_reciprocal->mat[j],
address_double);
lagmat_multiply_matrix_vector_l3(address_double_rot,
rot_reciprocal->mat[j],
address_double);
grid_point_rot = grg_get_double_grid_index(address_double_rot, mesh, is_shift);
if (grid_point_rot < ir_mapping_table[i]) {
#ifdef _OPENMP
@ -569,9 +533,9 @@ get_ir_reciprocal_mesh_distortion(long grid_address[][3],
static long relocate_BZ_grid_address(long bz_grid_address[][3],
long bz_map[],
GRGCONST long grid_address[][3],
LAGCONST long grid_address[][3],
const long mesh[3],
GRGCONST double rec_lattice[3][3],
LAGCONST double rec_lattice[3][3],
const long is_shift[3])
{
double tolerance, min_distance;
@ -641,9 +605,9 @@ static long relocate_BZ_grid_address(long bz_grid_address[][3],
static long get_bz_grid_addresses(long bz_grid_address[][3],
long bz_map[][2],
GRGCONST long grid_address[][3],
LAGCONST long grid_address[][3],
const long mesh[3],
GRGCONST double rec_lattice[3][3],
LAGCONST double rec_lattice[3][3],
const long is_shift[3])
{
double tolerance, min_distance;
@ -691,7 +655,7 @@ static long get_bz_grid_addresses(long bz_grid_address[][3],
return num_gp;
}
static double get_tolerance_for_BZ_reduction(GRGCONST double rec_lattice[3][3],
static double get_tolerance_for_BZ_reduction(LAGCONST double rec_lattice[3][3],
const long mesh[3])
{
long i, j;
@ -780,57 +744,8 @@ static long check_mesh_symmetry(const long mesh[3],
((eq[2] && mesh[2] == mesh[0] && is_shift[2] == is_shift[0]) || (!eq[2])));
}
static void transpose_matrix_l3(long a[3][3], GRGCONST long b[3][3])
{
long c[3][3];
c[0][0] = b[0][0];
c[0][1] = b[1][0];
c[0][2] = b[2][0];
c[1][0] = b[0][1];
c[1][1] = b[1][1];
c[1][2] = b[2][1];
c[2][0] = b[0][2];
c[2][1] = b[1][2];
c[2][2] = b[2][2];
bzg_copy_matrix_l3(a, c);
}
static void multiply_matrix_l3(long m[3][3],
GRGCONST long a[3][3],
GRGCONST long b[3][3])
{
long i, j; /* a_ij */
long c[3][3];
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];
}
}
bzg_copy_matrix_l3(m, c);
}
static long check_identity_matrix_l3(GRGCONST long a[3][3],
GRGCONST long b[3][3])
{
if ( a[0][0] - b[0][0] ||
a[0][1] - b[0][1] ||
a[0][2] - b[0][2] ||
a[1][0] - b[1][0] ||
a[1][1] - b[1][1] ||
a[1][2] - b[1][2] ||
a[2][0] - b[2][0] ||
a[2][1] - b[2][1] ||
a[2][2] - b[2][2]) {
return 0;
}
else {
return 1;
}
}
static void multiply_matrix_vector_d3(double v[3],
GRGCONST double a[3][3],
LAGCONST double a[3][3],
const double b[3])
{
long i;

14
c/bzgrid.h
View File

@ -35,7 +35,7 @@
#ifndef __bzgrid_H__
#define __bzgrid_H__
#include "grgrid.h"
#include "lagrid.h"
typedef struct {
long size;
@ -57,20 +57,16 @@ long bzg_get_ir_reciprocal_mesh(long grid_address[][3],
const MatLONG * rotations);
long bzg_relocate_BZ_grid_address(long bz_grid_address[][3],
long bz_map[],
GRGCONST long grid_address[][3],
LAGCONST long grid_address[][3],
const long mesh[3],
GRGCONST double rec_lattice[3][3],
LAGCONST double rec_lattice[3][3],
const long is_shift[3]);
long bzg_get_bz_grid_addresses(long bz_grid_address[][3],
long bz_map[][2],
GRGCONST long grid_address[][3],
LAGCONST long grid_address[][3],
const long mesh[3],
GRGCONST double rec_lattice[3][3],
LAGCONST double rec_lattice[3][3],
const long is_shift[3]);
void bzg_copy_matrix_l3(long a[3][3], GRGCONST long b[3][3]);
void bzg_multiply_matrix_vector_l3(long v[3],
GRGCONST long a[3][3],
const long b[3]);
MatLONG * bzg_alloc_MatLONG(const long size);
void bzg_free_MatLONG(MatLONG * matlong);

342
c/grgrid.c
View File

@ -36,6 +36,7 @@
#include <stddef.h>
#include <assert.h>
#include "grgrid.h"
#include "lagrid.h"
#include "snf3x3.h"
#define IDENTITY_TOL 1e-5
@ -60,62 +61,27 @@ 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,
GRGCONST long rotation[3][3],
LAGCONST long rotation[3][3],
const long D_diag[3],
const long PS[3]);
static void get_ir_grid_map(long ir_grid_indices[],
GRGCONST long (*rotations)[3][3],
LAGCONST long (*rotations)[3][3],
const long num_rot,
const long D_diag[3],
const long PS[3]);
static long mat_get_determinant_l3(GRGCONST long a[3][3]);
static double mat_get_determinant_d3(GRGCONST double a[3][3]);
static void mat_cast_matrix_3l_to_3d(double m[3][3], GRGCONST long a[3][3]);
static void mat_cast_matrix_3d_to_3l(long m[3][3], GRGCONST double a[3][3]);
static long mat_get_similar_matrix_ld3(double m[3][3],
GRGCONST long a[3][3],
GRGCONST double b[3][3],
const double precision);
static long mat_check_identity_matrix_l3(GRGCONST long a[3][3],
GRGCONST long b[3][3]);
static long mat_check_identity_matrix_ld3(GRGCONST long a[3][3],
GRGCONST double b[3][3],
const double symprec);
static long mat_inverse_matrix_d3(double m[3][3],
GRGCONST double a[3][3],
const double precision);
static void mat_transpose_matrix_l3(long a[3][3], GRGCONST long b[3][3]);
static void mat_multiply_matrix_vector_l3(long v[3],
GRGCONST long a[3][3],
const long b[3]);
static void mat_multiply_matrix_l3(long m[3][3],
GRGCONST long a[3][3],
GRGCONST long b[3][3]);
static void mat_multiply_matrix_ld3(double m[3][3],
GRGCONST long a[3][3],
GRGCONST double b[3][3]);
static void mat_multiply_matrix_d3(double m[3][3],
GRGCONST double a[3][3],
GRGCONST double b[3][3]);
static void mat_copy_matrix_l3(long a[3][3], GRGCONST long b[3][3]);
static void mat_copy_matrix_d3(double a[3][3], GRGCONST double b[3][3]);
static void mat_copy_vector_l3(long a[3], const long b[3]);
static long mat_modulo_l(const long a, const long b);
static long mat_Nint(const double a);
static double mat_Dabs(const double a);
long grg_get_snf3x3(long D_diag[3],
long P[3][3],
long Q[3][3],
GRGCONST long A[3][3])
LAGCONST long A[3][3])
{
long i, j, succeeded;
long D[3][3];
succeeded = 0;
if (mat_get_determinant_l3(A) == 0) {
if (lagmat_get_determinant_l3(A) == 0) {
goto err;
}
@ -143,10 +109,10 @@ err:
/* vectors. */
/* num_rot : Number of rotations */
long grg_transform_rotations(long (*transformed_rots)[3][3],
GRGCONST long (*rotations)[3][3],
LAGCONST long (*rotations)[3][3],
const long num_rot,
const long D_diag[3],
GRGCONST long Q[3][3])
LAGCONST long Q[3][3])
{
long i, j, k;
double r[3][3], Q_double[3][3];
@ -157,17 +123,17 @@ long grg_transform_rotations(long (*transformed_rots)[3][3],
/* 1. Compute (Q^-1)RQ */
/* 2. Compute D(Q^-1)RQ */
/* 3. Compute D(Q^-1)RQ(D^-1) */
mat_cast_matrix_3l_to_3d(Q_double, Q);
lagmat_cast_matrix_3l_to_3d(Q_double, Q);
for (i = 0; i < num_rot; i++) {
mat_get_similar_matrix_ld3(r, rotations[i], Q_double, 0);
lagmat_get_similar_matrix_ld3(r, rotations[i], Q_double, 0);
for (j = 0; j < 3; j++) {
for (k = 0; k < 3; k++) {
r[j][k] *= D_diag[j];
r[j][k] /= D_diag[k];
}
}
mat_cast_matrix_3d_to_3l(transformed_rots[i], r);
if (!mat_check_identity_matrix_ld3(transformed_rots[i], r, IDENTITY_TOL)) {
lagmat_cast_matrix_3d_to_3l(transformed_rots[i], r);
if (!lagmat_check_identity_matrix_ld3(transformed_rots[i], r, IDENTITY_TOL)) {
return 0;
}
}
@ -239,7 +205,7 @@ long grg_get_grid_index(const long address[3], const long D_diag[3])
{
long red_adrs[3];
mat_copy_vector_l3(red_adrs, address);
lagmat_copy_vector_l3(red_adrs, address);
reduce_grid_address(red_adrs, D_diag);
return get_grid_index_from_address(red_adrs, D_diag);
}
@ -261,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,
GRGCONST long rotation[3][3],
LAGCONST long rotation[3][3],
const long D_diag[3],
const long PS[3])
{
@ -272,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[],
GRGCONST long (*rotations)[3][3],
LAGCONST long (*rotations)[3][3],
const long num_rot,
const long D_diag[3],
const long PS[3])
@ -289,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],
GRGCONST long (*rotations)[3][3],
LAGCONST long (*rotations)[3][3],
const long num_rot,
const long is_time_reversal)
{
long i, j, num_rot_ret, inv_exist;
GRGCONST long inversion[3][3] = {
LAGCONST long inversion[3][3] = {
{-1, 0, 0 },
{ 0,-1, 0 },
{ 0, 0,-1 }
@ -303,14 +269,14 @@ long grg_get_reciprocal_point_group(long rec_rotations[48][3][3],
num_rot_ret = 0;
for (i = 0; i < num_rot; i++) {
for (j = 0; j < num_rot_ret; j++) {
if (mat_check_identity_matrix_l3(rotations[i], rec_rotations[j])) {
if (lagmat_check_identity_matrix_l3(rotations[i], rec_rotations[j])) {
goto escape;
}
}
if (num_rot_ret == 48) {
goto err;
}
mat_copy_matrix_l3(rec_rotations[num_rot_ret], rotations[i]);
lagmat_copy_matrix_l3(rec_rotations[num_rot_ret], rotations[i]);
num_rot_ret++;
escape:
;
@ -319,7 +285,7 @@ long grg_get_reciprocal_point_group(long rec_rotations[48][3][3],
inv_exist = 0;
if (is_time_reversal) {
for (i = 0; i < num_rot_ret; i++) {
if (mat_check_identity_matrix_l3(inversion, rec_rotations[i])) {
if (lagmat_check_identity_matrix_l3(inversion, rec_rotations[i])) {
inv_exist = 1;
break;
}
@ -331,15 +297,15 @@ long grg_get_reciprocal_point_group(long rec_rotations[48][3][3],
}
for (i = 0; i < num_rot_ret; i++) {
mat_multiply_matrix_l3(rec_rotations[num_rot_ret + i],
inversion, rec_rotations[i]);
lagmat_multiply_matrix_l3(rec_rotations[num_rot_ret + i],
inversion, rec_rotations[i]);
}
num_rot_ret *= 2;
}
}
for (i = 0; i < num_rot_ret; i++) {
mat_transpose_matrix_l3(rec_rotations[i], rec_rotations[i]);
lagmat_transpose_matrix_l3(rec_rotations[i], rec_rotations[i]);
}
return num_rot_ret;
@ -353,7 +319,7 @@ static void reduce_grid_address(long address[3], const long D_diag[3])
long i;
for (i = 0; i < 3; i++) {
address[i] = mat_modulo_l(address[i], D_diag[i]);
address[i] = lagmat_modulo_l(address[i], D_diag[i]);
}
}
@ -363,7 +329,7 @@ static void reduce_double_grid_address(long address_double[3],
long i;
for (i = 0; i < 3; i++) {
address_double[i] = mat_modulo_l(address_double[i], 2 * D_diag[i]);
address_double[i] = lagmat_modulo_l(address_double[i], 2 * D_diag[i]);
}
}
@ -409,7 +375,7 @@ static void get_all_grid_addresses(long grid_address[][3],
for (k = 0; k < D_diag[2]; k++) {
address[2] = k;
grid_index = get_grid_index_from_address(address, D_diag);
mat_copy_vector_l3(grid_address[grid_index], address);
lagmat_copy_vector_l3(grid_address[grid_index], address);
}
}
}
@ -462,7 +428,7 @@ static void get_double_grid_address(long address_double[3],
}
static long rotate_grid_index(const long grid_index,
GRGCONST long rotation[3][3],
LAGCONST long rotation[3][3],
const long D_diag[3],
const long PS[3])
{
@ -470,12 +436,12 @@ static long rotate_grid_index(const long grid_index,
get_grid_address_from_index(adrs, grid_index, D_diag);
get_double_grid_address(dadrs, adrs, PS);
mat_multiply_matrix_vector_l3(dadrs_rot, rotation, dadrs);
lagmat_multiply_matrix_vector_l3(dadrs_rot, rotation, dadrs);
return get_double_grid_index(dadrs_rot, D_diag, PS);
}
static void get_ir_grid_map(long ir_grid_indices[],
GRGCONST long (*rotations)[3][3],
LAGCONST long (*rotations)[3][3],
const long num_rot,
const long D_diag[3],
const long PS[3])
@ -505,255 +471,3 @@ static void get_ir_grid_map(long ir_grid_indices[],
}
}
static long mat_get_determinant_l3(GRGCONST 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]);
}
static double mat_get_determinant_d3(GRGCONST 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]);
}
static void mat_cast_matrix_3l_to_3d(double m[3][3], GRGCONST long a[3][3])
{
m[0][0] = a[0][0];
m[0][1] = a[0][1];
m[0][2] = a[0][2];
m[1][0] = a[1][0];
m[1][1] = a[1][1];
m[1][2] = a[1][2];
m[2][0] = a[2][0];
m[2][1] = a[2][1];
m[2][2] = a[2][2];
}
static void mat_cast_matrix_3d_to_3l(long m[3][3], GRGCONST double a[3][3])
{
m[0][0] = mat_Nint(a[0][0]);
m[0][1] = mat_Nint(a[0][1]);
m[0][2] = mat_Nint(a[0][2]);
m[1][0] = mat_Nint(a[1][0]);
m[1][1] = mat_Nint(a[1][1]);
m[1][2] = mat_Nint(a[1][2]);
m[2][0] = mat_Nint(a[2][0]);
m[2][1] = mat_Nint(a[2][1]);
m[2][2] = mat_Nint(a[2][2]);
}
static long mat_get_similar_matrix_ld3(double m[3][3],
GRGCONST long a[3][3],
GRGCONST double b[3][3],
const double precision)
{
double c[3][3];
if (!mat_inverse_matrix_d3(c, b, precision)) {
warning_print("No similar matrix due to 0 determinant.\n");
return 0;
}
mat_multiply_matrix_ld3(m, a, b);
mat_multiply_matrix_d3(m, c, m);
return 1;
}
static long mat_check_identity_matrix_l3(GRGCONST long a[3][3],
GRGCONST long b[3][3])
{
if (a[0][0] - b[0][0] ||
a[0][1] - b[0][1] ||
a[0][2] - b[0][2] ||
a[1][0] - b[1][0] ||
a[1][1] - b[1][1] ||
a[1][2] - b[1][2] ||
a[2][0] - b[2][0] ||
a[2][1] - b[2][1] ||
a[2][2] - b[2][2]) {
return 0;
}
else {
return 1;
}
}
static long mat_check_identity_matrix_ld3(GRGCONST long a[3][3],
GRGCONST double b[3][3],
const double symprec)
{
if (mat_Dabs(a[0][0] - b[0][0]) > symprec ||
mat_Dabs(a[0][1] - b[0][1]) > symprec ||
mat_Dabs(a[0][2] - b[0][2]) > symprec ||
mat_Dabs(a[1][0] - b[1][0]) > symprec ||
mat_Dabs(a[1][1] - b[1][1]) > symprec ||
mat_Dabs(a[1][2] - b[1][2]) > symprec ||
mat_Dabs(a[2][0] - b[2][0]) > symprec ||
mat_Dabs(a[2][1] - b[2][1]) > symprec ||
mat_Dabs(a[2][2] - b[2][2]) > symprec) {
return 0;
}
else {
return 1;
}
}
static long mat_inverse_matrix_d3(double m[3][3],
GRGCONST double a[3][3],
const double precision)
{
double det;
double c[3][3];
det = mat_get_determinant_d3(a);
if (mat_Dabs(det) < precision) {
warning_print("No inverse matrix (det=%f)\n", det);
return 0;
}
c[0][0] = (a[1][1] * a[2][2] - a[1][2] * a[2][1]) / det;
c[1][0] = (a[1][2] * a[2][0] - a[1][0] * a[2][2]) / det;
c[2][0] = (a[1][0] * a[2][1] - a[1][1] * a[2][0]) / det;
c[0][1] = (a[2][1] * a[0][2] - a[2][2] * a[0][1]) / det;
c[1][1] = (a[2][2] * a[0][0] - a[2][0] * a[0][2]) / det;
c[2][1] = (a[2][0] * a[0][1] - a[2][1] * a[0][0]) / det;
c[0][2] = (a[0][1] * a[1][2] - a[0][2] * a[1][1]) / det;
c[1][2] = (a[0][2] * a[1][0] - a[0][0] * a[1][2]) / det;
c[2][2] = (a[0][0] * a[1][1] - a[0][1] * a[1][0]) / det;
mat_copy_matrix_d3(m, c);
return 1;
}
static void mat_transpose_matrix_l3(long a[3][3], GRGCONST long b[3][3])
{
long c[3][3];
c[0][0] = b[0][0];
c[0][1] = b[1][0];
c[0][2] = b[2][0];
c[1][0] = b[0][1];
c[1][1] = b[1][1];
c[1][2] = b[2][1];
c[2][0] = b[0][2];
c[2][1] = b[1][2];
c[2][2] = b[2][2];
mat_copy_matrix_l3(a, c);
}
static void mat_multiply_matrix_vector_l3(long v[3],
GRGCONST long a[3][3],
const long b[3])
{
long i;
long 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 void mat_multiply_matrix_l3(long m[3][3],
GRGCONST long a[3][3],
GRGCONST long b[3][3])
{
long i, j; /* a_ij */
long c[3][3];
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];
}
}
mat_copy_matrix_l3(m, c);
}
static void mat_multiply_matrix_ld3(double m[3][3],
GRGCONST long a[3][3],
GRGCONST double b[3][3])
{
long i, j; /* a_ij */
double c[3][3];
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];
}
}
mat_copy_matrix_d3(m, c);
}
static void mat_multiply_matrix_d3(double m[3][3],
GRGCONST double a[3][3],
GRGCONST double b[3][3])
{
long i, j; /* a_ij */
double c[3][3];
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];
}
}
mat_copy_matrix_d3(m, c);
}
static void mat_copy_matrix_l3(long a[3][3], GRGCONST long b[3][3])
{
a[0][0] = b[0][0];
a[0][1] = b[0][1];
a[0][2] = b[0][2];
a[1][0] = b[1][0];
a[1][1] = b[1][1];
a[1][2] = b[1][2];
a[2][0] = b[2][0];
a[2][1] = b[2][1];
a[2][2] = b[2][2];
}
static void mat_copy_matrix_d3(double a[3][3], GRGCONST double b[3][3])
{
a[0][0] = b[0][0];
a[0][1] = b[0][1];
a[0][2] = b[0][2];
a[1][0] = b[1][0];
a[1][1] = b[1][1];
a[1][2] = b[1][2];
a[2][0] = b[2][0];
a[2][1] = b[2][1];
a[2][2] = b[2][2];
}
static void mat_copy_vector_l3(long a[3], const long b[3])
{
a[0] = b[0];
a[1] = b[1];
a[2] = b[2];
}
static long mat_modulo_l(const long a, const long b)
{
long c;
c = a % b;
if (c < 0) {
c += b;
}
return c;
}
static long mat_Nint(const double a)
{
if (a < 0.0)
return (long) (a - 0.5);
else
return (long) (a + 0.5);
}
static double mat_Dabs(const double a)
{
if (a < 0.0)
return -a;
else
return a;
}

23
c/grgrid.h
View File

@ -36,26 +36,17 @@
#define __grgrid_H__
#include <stddef.h>
#ifndef GRGCONST
#define GRGCONST
#endif
#ifdef MATWARNING
#define warning_print(...) fprintf(stderr,__VA_ARGS__)
#else
#define warning_print(...)
#endif
#include "lagrid.h"
long grg_get_snf3x3(long D_diag[3],
long P[3][3],
long Q[3][3],
GRGCONST long A[3][3]);
LAGCONST long A[3][3]);
long grg_transform_rotations(long (*transformed_rots)[3][3],
GRGCONST long (*rotations)[3][3],
LAGCONST long (*rotations)[3][3],
const long num_rot,
const long D_diag[3],
GRGCONST long Q[3][3]);
LAGCONST 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],
@ -74,16 +65,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,
GRGCONST long rotations[3][3],
LAGCONST long rotations[3][3],
const long D_diag[3],
const long PS[3]);
void grg_get_ir_grid_map(long ir_grid_indices[],
GRGCONST long (*rotations)[3][3],
LAGCONST 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],
GRGCONST long (*rotations)[3][3],
LAGCONST long (*rotations)[3][3],
const long num_rot,
const long is_time_reversal);

287
c/lagrid.c Normal file
View File

@ -0,0 +1,287 @@
/* Copyright (C) 2021 Atsushi Togo */
/* All rights reserved. */
/* This file is part of kspclib. */
/* Redistribution and use in source and binary forms, with or without */
/* modification, are permitted provided that the following conditions */
/* are met: */
/* * Redistributions of source code must retain the above copyright */
/* notice, this list of conditions and the following disclaimer. */
/* * Redistributions in binary form must reproduce the above copyright */
/* notice, this list of conditions and the following disclaimer in */
/* the documentation and/or other materials provided with the */
/* distribution. */
/* * Neither the name of the kspclib project nor the names of its */
/* contributors may be used to endorse or promote products derived */
/* from this software without specific prior written permission. */
/* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS */
/* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT */
/* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS */
/* FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE */
/* COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, */
/* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, */
/* BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; */
/* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER */
/* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT */
/* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN */
/* ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE */
/* POSSIBILITY OF SUCH DAMAGE. */
#include "lagrid.h"
long lagmat_get_determinant_l3(LAGCONST 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])
{
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])
{
m[0][0] = a[0][0];
m[0][1] = a[0][1];
m[0][2] = a[0][2];
m[1][0] = a[1][0];
m[1][1] = a[1][1];
m[1][2] = a[1][2];
m[2][0] = a[2][0];
m[2][1] = a[2][1];
m[2][2] = a[2][2];
}
void lagmat_cast_matrix_3d_to_3l(long m[3][3], LAGCONST double a[3][3])
{
m[0][0] = lagmat_Nint(a[0][0]);
m[0][1] = lagmat_Nint(a[0][1]);
m[0][2] = lagmat_Nint(a[0][2]);
m[1][0] = lagmat_Nint(a[1][0]);
m[1][1] = lagmat_Nint(a[1][1]);
m[1][2] = lagmat_Nint(a[1][2]);
m[2][0] = lagmat_Nint(a[2][0]);
m[2][1] = lagmat_Nint(a[2][1]);
m[2][2] = lagmat_Nint(a[2][2]);
}
long lagmat_get_similar_matrix_ld3(double m[3][3],
LAGCONST long a[3][3],
LAGCONST double b[3][3],
const double precision)
{
double c[3][3];
if (!lagmat_inverse_matrix_d3(c, b, precision)) {
warning_print("No similar matrix due to 0 determinant.\n");
return 0;
}
lagmat_multiply_matrix_ld3(m, a, b);
lagmat_multiply_matrix_d3(m, c, m);
return 1;
}
long lagmat_check_identity_matrix_l3(LAGCONST long a[3][3],
LAGCONST long b[3][3])
{
if (a[0][0] - b[0][0] ||
a[0][1] - b[0][1] ||
a[0][2] - b[0][2] ||
a[1][0] - b[1][0] ||
a[1][1] - b[1][1] ||
a[1][2] - b[1][2] ||
a[2][0] - b[2][0] ||
a[2][1] - b[2][1] ||
a[2][2] - b[2][2]) {
return 0;
}
else {
return 1;
}
}
long lagmat_check_identity_matrix_ld3(LAGCONST long a[3][3],
LAGCONST double b[3][3],
const double symprec)
{
if (lagmat_Dabs(a[0][0] - b[0][0]) > symprec ||
lagmat_Dabs(a[0][1] - b[0][1]) > symprec ||
lagmat_Dabs(a[0][2] - b[0][2]) > symprec ||
lagmat_Dabs(a[1][0] - b[1][0]) > symprec ||
lagmat_Dabs(a[1][1] - b[1][1]) > symprec ||
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) {
return 0;
}
else {
return 1;
}
}
long lagmat_inverse_matrix_d3(double m[3][3],
LAGCONST double a[3][3],
const double precision)
{
double det;
double c[3][3];
det = lagmat_get_determinant_d3(a);
if (lagmat_Dabs(det) < precision) {
warning_print("No inverse matrix (det=%f)\n", det);
return 0;
}
c[0][0] = (a[1][1] * a[2][2] - a[1][2] * a[2][1]) / det;
c[1][0] = (a[1][2] * a[2][0] - a[1][0] * a[2][2]) / det;
c[2][0] = (a[1][0] * a[2][1] - a[1][1] * a[2][0]) / det;
c[0][1] = (a[2][1] * a[0][2] - a[2][2] * a[0][1]) / det;
c[1][1] = (a[2][2] * a[0][0] - a[2][0] * a[0][2]) / det;
c[2][1] = (a[2][0] * a[0][1] - a[2][1] * a[0][0]) / det;
c[0][2] = (a[0][1] * a[1][2] - a[0][2] * a[1][1]) / det;
c[1][2] = (a[0][2] * a[1][0] - a[0][0] * a[1][2]) / det;
c[2][2] = (a[0][0] * a[1][1] - a[0][1] * a[1][0]) / det;
lagmat_copy_matrix_d3(m, c);
return 1;
}
void lagmat_transpose_matrix_l3(long a[3][3], LAGCONST long b[3][3])
{
long c[3][3];
c[0][0] = b[0][0];
c[0][1] = b[1][0];
c[0][2] = b[2][0];
c[1][0] = b[0][1];
c[1][1] = b[1][1];
c[1][2] = b[2][1];
c[2][0] = b[0][2];
c[2][1] = b[1][2];
c[2][2] = b[2][2];
lagmat_copy_matrix_l3(a, c);
}
void lagmat_multiply_matrix_vector_l3(long v[3],
LAGCONST long a[3][3],
const long b[3])
{
long i;
long 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];
}
}
void lagmat_multiply_matrix_l3(long m[3][3],
LAGCONST long a[3][3],
LAGCONST long b[3][3])
{
long i, j; /* a_ij */
long c[3][3];
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];
}
}
lagmat_copy_matrix_l3(m, c);
}
void lagmat_multiply_matrix_ld3(double m[3][3],
LAGCONST long a[3][3],
LAGCONST double b[3][3])
{
long i, j; /* a_ij */
double c[3][3];
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];
}
}
lagmat_copy_matrix_d3(m, c);
}
void lagmat_multiply_matrix_d3(double m[3][3],
LAGCONST double a[3][3],
LAGCONST double b[3][3])
{
long i, j; /* a_ij */
double c[3][3];
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];
}
}
lagmat_copy_matrix_d3(m, c);
}
void lagmat_copy_matrix_l3(long a[3][3], LAGCONST long b[3][3])
{
a[0][0] = b[0][0];
a[0][1] = b[0][1];
a[0][2] = b[0][2];
a[1][0] = b[1][0];
a[1][1] = b[1][1];
a[1][2] = b[1][2];
a[2][0] = b[2][0];
a[2][1] = b[2][1];
a[2][2] = b[2][2];
}
void lagmat_copy_matrix_d3(double a[3][3], LAGCONST double b[3][3])
{
a[0][0] = b[0][0];
a[0][1] = b[0][1];
a[0][2] = b[0][2];
a[1][0] = b[1][0];
a[1][1] = b[1][1];
a[1][2] = b[1][2];
a[2][0] = b[2][0];
a[2][1] = b[2][1];
a[2][2] = b[2][2];
}
void lagmat_copy_vector_l3(long a[3], const long b[3])
{
a[0] = b[0];
a[1] = b[1];
a[2] = b[2];
}
long lagmat_modulo_l(const long a, const long b)
{
long c;
c = a % b;
if (c < 0) {
c += b;
}
return c;
}
long lagmat_Nint(const double a)
{
if (a < 0.0)
return (long) (a - 0.5);
else
return (long) (a + 0.5);
}
double lagmat_Dabs(const double a)
{
if (a < 0.0)
return -a;
else
return a;
}

84
c/lagrid.h Normal file
View File

@ -0,0 +1,84 @@
/* Copyright (C) 2021 Atsushi Togo */
/* All rights reserved. */
/* This file is part of kspclib. */
/* Redistribution and use in source and binary forms, with or without */
/* modification, are permitted provided that the following conditions */
/* are met: */
/* * Redistributions of source code must retain the above copyright */
/* notice, this list of conditions and the following disclaimer. */
/* * Redistributions in binary form must reproduce the above copyright */
/* notice, this list of conditions and the following disclaimer in */
/* the documentation and/or other materials provided with the */
/* distribution. */
/* * Neither the name of the kspclib project nor the names of its */
/* contributors may be used to endorse or promote products derived */
/* from this software without specific prior written permission. */
/* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS */
/* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT */
/* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS */
/* FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE */
/* COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, */
/* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, */
/* BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; */
/* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER */
/* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT */
/* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN */
/* ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE */
/* POSSIBILITY OF SUCH DAMAGE. */
#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_similar_matrix_ld3(double m[3][3],
LAGCONST long a[3][3],
LAGCONST 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],
const double symprec);
long lagmat_inverse_matrix_d3(double m[3][3],
LAGCONST double a[3][3],
const double precision);
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 b[3]);
void lagmat_multiply_matrix_l3(long m[3][3],
LAGCONST long a[3][3],
LAGCONST long b[3][3]);
void lagmat_multiply_matrix_ld3(double m[3][3],
LAGCONST long a[3][3],
LAGCONST 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]);
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);
double lagmat_Dabs(const double a);
#endif

3
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"
@ -623,7 +624,7 @@ long ph3py_get_ir_reciprocal_mesh(long grid_address[][3],
rotations = bzg_alloc_MatLONG(num_rot);
for (i = 0; i < num_rot; i++) {
bzg_copy_matrix_l3(rotations->mat[i], rotations_in[i]);
lagmat_copy_matrix_l3(rotations->mat[i], rotations_in[i]);
}
num_ir = bzg_get_ir_reciprocal_mesh(grid_address,
ir_mapping_table,

20
c/triplet.c
View File

@ -34,7 +34,7 @@
/* ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE */
/* POSSIBILITY OF SUCH DAMAGE. */
#include "grgrid.h"
#include "lagrid.h"
#include "triplet.h"
#include "triplet_iw.h"
#include "triplet_kpoint.h"
@ -46,12 +46,12 @@ static long get_triplets_reciprocal_mesh_at_q(long *map_triplets,
const long mesh[3],
const long is_time_reversal,
const long num_rot,
GRGCONST long (*rotations)[3][3],
LAGCONST long (*rotations)[3][3],
const long swappable);
long tpl_get_BZ_triplets_at_q(long (*triplets)[3],
const long grid_point,
GRGCONST long (*bz_grid_address)[3],
LAGCONST long (*bz_grid_address)[3],
const long *bz_map,
const long *map_triplets,
const long num_map_triplets,
@ -73,7 +73,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,
GRGCONST long (*rotations)[3][3],
LAGCONST long (*rotations)[3][3],
const long swappable)
{
return get_triplets_reciprocal_mesh_at_q(map_triplets,
@ -91,11 +91,11 @@ void tpl_get_integration_weight(double *iw,
char *iw_zero,
const double *frequency_points,
const long num_band0,
GRGCONST long relative_grid_address[24][4][3],
LAGCONST long relative_grid_address[24][4][3],
const long mesh[3],
GRGCONST long (*triplets)[3],
LAGCONST long (*triplets)[3],
const long num_triplets,
GRGCONST long (*bz_grid_address)[3],
LAGCONST long (*bz_grid_address)[3],
const long *bz_map,
const double *frequencies1,
const long num_band1,
@ -141,7 +141,7 @@ void tpl_get_integration_weight_with_sigma(double *iw,
const double sigma_cutoff,
const double *frequency_points,
const long num_band0,
GRGCONST long (*triplets)[3],
LAGCONST long (*triplets)[3],
const long num_triplets,
const double *frequencies,
const long num_band,
@ -193,7 +193,7 @@ long tpl_is_N(const long triplet[3], const long *grid_address)
void tpl_set_relative_grid_address(
long tp_relative_grid_address[2][24][4][3],
GRGCONST long relative_grid_address[24][4][3],
LAGCONST long relative_grid_address[24][4][3],
const long tp_type)
{
long i, j, k, l;
@ -227,7 +227,7 @@ static long get_triplets_reciprocal_mesh_at_q(long *map_triplets,
const long mesh[3],
const long is_time_reversal,
const long num_rot,
GRGCONST long (*rotations)[3][3],
LAGCONST long (*rotations)[3][3],
const long swappable)
{
long num_ir;

16
c/triplet.h
View File

@ -38,7 +38,7 @@
#define __triplet_H__
#include <stddef.h>
#include "grgrid.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 +53,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,
GRGCONST long (*rotations)[3][3],
LAGCONST long (*rotations)[3][3],
const long swappable);
/* Irreducible grid-point-triplets in BZ are stored. */
/* triplets are recovered from grid_point and triplet_weights. */
@ -63,7 +63,7 @@ long tpl_get_triplets_reciprocal_mesh_at_q(long *map_triplets,
/* Number of ir-triplets is returned. */
long tpl_get_BZ_triplets_at_q(long (*triplets)[3],
const long grid_point,
GRGCONST long (*bz_grid_address)[3],
LAGCONST long (*bz_grid_address)[3],
const long *bz_map,
const long *map_triplets,
const long num_map_triplets,
@ -72,11 +72,11 @@ void tpl_get_integration_weight(double *iw,
char *iw_zero,
const double *frequency_points,
const long num_band0,
GRGCONST long relative_grid_address[24][4][3],
LAGCONST long relative_grid_address[24][4][3],
const long mesh[3],
GRGCONST long (*triplets)[3],
LAGCONST long (*triplets)[3],
const long num_triplets,
GRGCONST long (*bz_grid_address)[3],
LAGCONST long (*bz_grid_address)[3],
const long *bz_map,
const double *frequencies1,
const long num_band1,
@ -91,7 +91,7 @@ void tpl_get_integration_weight_with_sigma(double *iw,
const double sigma_cutoff,
const double *frequency_points,
const long num_band0,
GRGCONST long (*triplets)[3],
LAGCONST long (*triplets)[3],
const long num_triplets,
const double *frequencies,
const long num_band,
@ -100,7 +100,7 @@ void tpl_get_integration_weight_with_sigma(double *iw,
long tpl_is_N(const long triplet[3], const long *grid_address);
void tpl_set_relative_grid_address(
long tp_relative_grid_address[2][24][4][3],
GRGCONST long relative_grid_address[24][4][3],
LAGCONST long relative_grid_address[24][4][3],
const long tp_type);
#endif

29
c/triplet_iw.c
View File

@ -34,6 +34,7 @@
#include <math.h>
#include "grgrid.h"
#include "lagrid.h"
#include "phonoc_utils.h"
#include "triplet.h"
#include "triplet_iw.h"
@ -42,7 +43,7 @@
static void set_freq_vertices(double freq_vertices[3][24][4],
const double *frequencies1,
const double *frequencies2,
GRGCONST long vertices[2][24][4],
LAGCONST long vertices[2][24][4],
const long num_band1,
const long num_band2,
const long b1,
@ -50,15 +51,15 @@ 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,
GRGCONST double freq_vertices[3][24][4],
LAGCONST double freq_vertices[3][24][4],
const long max_i);
static long in_tetrahedra(const double f0, GRGCONST double freq_vertices[24][4]);
static long in_tetrahedra(const double f0, LAGCONST double freq_vertices[24][4]);
static void get_triplet_tetrahedra_vertices(
long vertices[2][24][4],
GRGCONST long tp_relative_grid_address[2][24][4][3],
LAGCONST long tp_relative_grid_address[2][24][4][3],
const long mesh[3],
const long triplet[3],
GRGCONST long (*bz_grid_address)[3],
LAGCONST long (*bz_grid_address)[3],
const long *bz_map);
void
@ -66,11 +67,11 @@ tpi_get_integration_weight(double *iw,
char *iw_zero,
const double *frequency_points,
const long num_band0,
GRGCONST long tp_relative_grid_address[2][24][4][3],
LAGCONST long tp_relative_grid_address[2][24][4][3],
const long mesh[3],
const long triplets[3],
const long num_triplets,
GRGCONST long (*bz_grid_address)[3],
LAGCONST long (*bz_grid_address)[3],
const long *bz_map,
const double *frequencies1,
const long num_band1,
@ -211,10 +212,10 @@ void tpi_get_integration_weight_with_sigma(double *iw,
void
tpi_get_neighboring_grid_points(long neighboring_grid_points[],
const long grid_point,
GRGCONST long relative_grid_address[][3],
LAGCONST long relative_grid_address[][3],
const long num_relative_grid_address,
const long mesh[3],
GRGCONST long bz_grid_address[][3],
LAGCONST long bz_grid_address[][3],
const long bz_map[])
{
long bzmesh[3], address_double[3], bz_address_double[3], PS[3];
@ -245,7 +246,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,
GRGCONST long vertices[2][24][4],
LAGCONST long vertices[2][24][4],
const long num_band1,
const long num_band2,
const long b1,
@ -283,7 +284,7 @@ static void set_freq_vertices(double freq_vertices[3][24][4],
/* calculation. */
static long set_g(double g[3],
const double f0,
GRGCONST double freq_vertices[3][24][4],
LAGCONST double freq_vertices[3][24][4],
const long max_i)
{
long i, iw_zero;
@ -302,7 +303,7 @@ static long set_g(double g[3],
return iw_zero;
}
static long in_tetrahedra(const double f0, GRGCONST double freq_vertices[24][4])
static long in_tetrahedra(const double f0, LAGCONST double freq_vertices[24][4])
{
long i, j;
double fmin, fmax;
@ -330,10 +331,10 @@ static long in_tetrahedra(const double f0, GRGCONST double freq_vertices[24][4])
static void get_triplet_tetrahedra_vertices(
long vertices[2][24][4],
GRGCONST long tp_relative_grid_address[2][24][4][3],
LAGCONST long tp_relative_grid_address[2][24][4][3],
const long mesh[3],
const long triplet[3],
GRGCONST long (*bz_grid_address)[3],
LAGCONST long (*bz_grid_address)[3],
const long *bz_map)
{
long i, j;

10
c/triplet_iw.h
View File

@ -35,18 +35,18 @@
#ifndef __triplet_iw_H__
#define __triplet_iw_H__
#include "grgrid.h"
#include "lagrid.h"
void
tpi_get_integration_weight(double *iw,
char *iw_zero,
const double *frequency_points,
const long num_band0,
GRGCONST long tp_relative_grid_address[2][24][4][3],
LAGCONST long tp_relative_grid_address[2][24][4][3],
const long mesh[3],
const long triplets[3],
const long num_triplets,
GRGCONST long (*bz_grid_address)[3],
LAGCONST long (*bz_grid_address)[3],
const long *bz_map,
const double *frequencies1,
const long num_band1,
@ -69,10 +69,10 @@ void tpi_get_integration_weight_with_sigma(double *iw,
void
tpi_get_neighboring_grid_points(long neighboring_grid_points[],
const long grid_point,
GRGCONST long relative_grid_address[][3],
LAGCONST long relative_grid_address[][3],
const long num_relative_grid_address,
const long mesh[3],
GRGCONST long bz_grid_address[][3],
LAGCONST long bz_grid_address[][3],
const long bz_map[]);
#endif

17
c/triplet_kpoint.c
View File

@ -38,6 +38,7 @@
#include <stdlib.h>
#include "bzgrid.h"
#include "grgrid.h"
#include "lagrid.h"
#include "triplet.h"
#include "triplet_kpoint.h"
@ -179,7 +180,7 @@ static long get_ir_triplets_at_q(long *map_triplets,
const long swappable);
static long get_BZ_triplets_at_q(long (*triplets)[3],
const long grid_point,
GRGCONST long (*bz_grid_address)[3],
LAGCONST long (*bz_grid_address)[3],
const long *bz_map,
const long *map_triplets,
const long num_map_triplets,
@ -200,7 +201,7 @@ long tpk_get_ir_triplets_at_q(long *map_triplets,
const long grid_point,
const long mesh[3],
const long is_time_reversal,
GRGCONST long (*rotations)[3][3],
LAGCONST long (*rotations)[3][3],
const long num_rot,
const long swappable)
{
@ -209,7 +210,7 @@ long tpk_get_ir_triplets_at_q(long *map_triplets,
rot_real = bzg_alloc_MatLONG(num_rot);
for (i = 0; i < num_rot; i++) {
bzg_copy_matrix_l3(rot_real->mat[i], rotations[i]);
lagmat_copy_matrix_l3(rot_real->mat[i], rotations[i]);
}
rot_reciprocal = bzg_get_point_group_reciprocal(rot_real, is_time_reversal);
bzg_free_MatLONG(rot_real);
@ -228,7 +229,7 @@ long tpk_get_ir_triplets_at_q(long *map_triplets,
long tpk_get_BZ_triplets_at_q(long (*triplets)[3],
const long grid_point,
GRGCONST long (*bz_grid_address)[3],
LAGCONST long (*bz_grid_address)[3],
const long *bz_map,
const long *map_triplets,
const long num_map_triplets,
@ -351,7 +352,7 @@ ret:
static long get_BZ_triplets_at_q(long (*triplets)[3],
const long grid_point,
GRGCONST long (*bz_grid_address)[3],
LAGCONST long (*bz_grid_address)[3],
const long *bz_map,
const long *map_triplets,
const long num_map_triplets,
@ -502,7 +503,7 @@ static MatLONG *get_point_group_reciprocal_with_q(const MatLONG * rot_reciprocal
ir_rot[i] = -1;
}
for (i = 0; i < rot_reciprocal->size; i++) {
bzg_multiply_matrix_vector_l3(adrs_rot, rot_reciprocal->mat[i], adrs);
lagmat_multiply_matrix_vector_l3(adrs_rot, rot_reciprocal->mat[i], adrs);
gp_rot = grg_get_grid_index(adrs_rot, mesh);
if (gp_rot == grid_point) {
@ -513,8 +514,8 @@ static MatLONG *get_point_group_reciprocal_with_q(const MatLONG * rot_reciprocal
if ((rot_reciprocal_q = bzg_alloc_MatLONG(num_rot)) != NULL) {
for (i = 0; i < num_rot; i++) {
bzg_copy_matrix_l3(rot_reciprocal_q->mat[i],
rot_reciprocal->mat[ir_rot[i]]);
lagmat_copy_matrix_l3(rot_reciprocal_q->mat[i],
rot_reciprocal->mat[ir_rot[i]]);
}
}

6
c/triplet_kpoint.h
View File

@ -37,7 +37,7 @@
#ifndef __triplet_kpoint_H__
#define __triplet_kpoint_H__
#include "grgrid.h"
#include "lagrid.h"
long tpk_get_ir_triplets_at_q(long *map_triplets,
long *map_q,
@ -45,12 +45,12 @@ long tpk_get_ir_triplets_at_q(long *map_triplets,
const long grid_point,
const long mesh[3],
const long is_time_reversal,
GRGCONST long (*rotations)[3][3],
LAGCONST long (*rotations)[3][3],
const long num_rot,
const long swappable);
long tpk_get_BZ_triplets_at_q(long (*triplets)[3],
const long grid_point,
GRGCONST long (*bz_grid_address)[3],
LAGCONST long (*bz_grid_address)[3],
const long *bz_map,
const long *map_triplets,
const long num_map_triplets,

3
setup.py
View File

@ -209,11 +209,12 @@ print("extra_link_args", extra_link_args)
sources_phono3py = ['c/_phono3py.c',
'c/bzgrid.c',
'c/collision_matrix.c',
'c/grgrid.c',
'c/fc3.c',
'c/grgrid.c',
'c/imag_self_energy_with_g.c',
'c/interaction.c',
'c/isotope.c',
'c/lagrid.c',
'c/lapack_wrapper.c',
'c/phono3py.c',
'c/phonoc_utils.c',