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

Replace rgrid by grgrid

This commit is contained in:
Atsushi Togo 2021-03-02 16:50:15 +09:00
parent 3ab893614c
commit 4deb9af767
12 changed files with 1572 additions and 279 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
CMakeLists.txt
View File

@ -45,6 +45,7 @@ include_directories("${PROJECT_SOURCE_DIR}/c")
set(SOURCES_PHONO3PY
${PROJECT_SOURCE_DIR}/c/collision_matrix.c
${PROJECT_SOURCE_DIR}/c/fc3.c
${PROJECT_SOURCE_DIR}/c/grgrid.c
${PROJECT_SOURCE_DIR}/c/imag_self_energy_with_g.c
${PROJECT_SOURCE_DIR}/c/interaction.c
${PROJECT_SOURCE_DIR}/c/isotope.c
@ -56,7 +57,7 @@ set(SOURCES_PHONO3PY
${PROJECT_SOURCE_DIR}/c/real_self_energy.c
${PROJECT_SOURCE_DIR}/c/real_to_reciprocal.c
${PROJECT_SOURCE_DIR}/c/reciprocal_to_normal.c
${PROJECT_SOURCE_DIR}/c/rgrid.c
${PROJECT_SOURCE_DIR}/c/snf3x3.c
${PROJECT_SOURCE_DIR}/c/tetrahedron_method.c
${PROJECT_SOURCE_DIR}/c/triplet.c
${PROJECT_SOURCE_DIR}/c/triplet_iw.c

759
c/grgrid.c Normal file
View File

@ -0,0 +1,759 @@
/* Copyright (C) 2020 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 <stdio.h>
#include <stddef.h>
#include <assert.h>
#include "grgrid.h"
#include "snf3x3.h"
#define IDENTITY_TOL 1e-5
static void reduce_grid_address(long address[3], const long D_diag[3]);
static void reduce_double_grid_address(long address_double[3],
const long D_diag[3]);
static long get_double_grid_index(const long address_double[3],
const long D_diag[3],
const long PS[3]);
static long get_grid_index_from_address(const long address[3],
const long D_diag[3]);
static void get_all_grid_addresses(long grid_address[][3],
const long D_diag[3]);
static void get_grid_address_from_index(long address[3],
const long grid_index,
const long D_diag[3]);
static void get_grid_address(long address[3],
const long address_double[3],
const long PS[3]);
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,
MATCONST long rotation[3][3],
const long D_diag[3],
const long PS[3]);
static void get_ir_grid_map(long ir_grid_indices[],
MATCONST long (*rotations)[3][3],
const long num_rot,
const long D_diag[3],
const long PS[3]);
static long mat_get_determinant_l3(MATCONST long a[3][3]);
static double mat_get_determinant_d3(MATCONST double a[3][3]);
static void mat_cast_matrix_3l_to_3d(double m[3][3], MATCONST long a[3][3]);
static void mat_cast_matrix_3d_to_3l(long m[3][3], MATCONST double a[3][3]);
static long mat_get_similar_matrix_ld3(double m[3][3],
MATCONST long a[3][3],
MATCONST double b[3][3],
const double precision);
static long mat_check_identity_matrix_l3(MATCONST long a[3][3],
MATCONST long b[3][3]);
static long mat_check_identity_matrix_ld3(MATCONST long a[3][3],
MATCONST double b[3][3],
const double symprec);
static long mat_inverse_matrix_d3(double m[3][3],
MATCONST double a[3][3],
const double precision);
static void mat_transpose_matrix_l3(long a[3][3], MATCONST long b[3][3]);
static void mat_multiply_matrix_vector_l3(long v[3],
MATCONST long a[3][3],
const long b[3]);
static void mat_multiply_matrix_l3(long m[3][3],
MATCONST long a[3][3],
MATCONST long b[3][3]);
static void mat_multiply_matrix_ld3(double m[3][3],
MATCONST long a[3][3],
MATCONST double b[3][3]);
static void mat_multiply_matrix_d3(double m[3][3],
MATCONST double a[3][3],
MATCONST double b[3][3]);
static void mat_copy_matrix_l3(long a[3][3], MATCONST long b[3][3]);
static void mat_copy_matrix_d3(double a[3][3], MATCONST 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],
MATCONST long A[3][3])
{
long i, j, succeeded;
long D[3][3];
succeeded = 0;
if (mat_get_determinant_l3(A) == 0) {
goto err;
}
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
D[i][j] = A[i][j];
}
}
succeeded = snf3x3(D, P, Q);
for (i = 0; i < 3; i++) {
D_diag[i] = D[i][i];
}
err:
return succeeded;
}
/*----------------------------------------*/
/* Transform rotations by D(Q^-1)RQ(D^-1) */
/*----------------------------------------*/
/* transformed_rots : D(Q^-1)RQ(D^-1) */
/* rotations : [num_rot][3][3] */
/* Defined as q' = Rq where q is in the reciprocal primitive basis */
/* vectors. */
/* num_rot : Number of rotations */
long grg_transform_rotations(long (*transformed_rots)[3][3],
MATCONST long (*rotations)[3][3],
const long num_rot,
const long D_diag[3],
MATCONST long Q[3][3])
{
long i, j, k;
double r[3][3], Q_double[3][3];
/* Compute D(Q^-1)RQ(D^-1) by three steps */
/* It is assumed that |det(Q)|=1 and Q^-1 has relatively small round-off */
/* error, and we want to divide by D carefully. */
/* 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);
for (i = 0; i < num_rot; i++) {
mat_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)) {
return 0;
}
}
return 1;
}
/* -------------------------------*/
/* Get all address in single grid */
/* -------------------------------*/
/* 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])
{
get_all_grid_addresses(grid_address, D_diag);
}
/* -------------------------------------------------------*/
/* Get address in double grid from address in single grid */
/* -------------------------------------------------------*/
/* address_double : Double grid address. */
/* address : Single grid address. */
/* D_diag : Diagnal elements of D. */
/* PS : Shifts transformed by P. s_i is 0 or 1. */
void grg_get_double_grid_address(long address_double[3],
const long address[3],
const long D_diag[3],
const long PS[3])
{
get_double_grid_address(address_double, address, PS);
reduce_double_grid_address(address_double, D_diag);
}
/* -------------------------------------------------------*/
/* Get address in single grid from address in double grid */
/* -------------------------------------------------------*/
/* address : Single grid address. */
/* address_double : Double grid address. */
/* D_diag : Diagnal elements of D. */
/* PS : Shifts transformed by P. s_i is 0 or 1. */
void grg_get_grid_address(long address[3],
const long address_double[3],
const long D_diag[3],
const long PS[3])
{
get_grid_address(address, address_double, PS);
reduce_grid_address(address, D_diag);
}
/* -------------------------------------------------*/
/* Get grid point index from address in double grid */
/* -------------------------------------------------*/
/* address_double : Double grid address. */
/* D_diag : Diagnal elements of D. */
/* PS : Shifts transformed by P. s_i is 0 or 1. */
long grg_get_double_grid_index(const long address_double[3],
const long D_diag[3],
const long PS[3])
{
return get_double_grid_index(address_double, D_diag, PS);
}
/* -------------------------------------------------*/
/* Get grid point index from address in single grid */
/* -------------------------------------------------*/
/* address : Single grid address. */
/* D_diag : Diagnal elements of D. */
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);
reduce_grid_address(red_adrs, D_diag);
return get_grid_index_from_address(red_adrs, D_diag);
}
/* ---------------------------------------*/
/* Get grid address from grid point index */
/* ---------------------------------------*/
/* address : Single grid address. */
/* D_diag : Diagnal elements of D. */
void grg_get_grid_address_from_index(long address[3],
const long grid_index,
const long D_diag[3])
{
get_grid_address_from_index(address, grid_index, D_diag);
}
/* ---------------------------*/
/* Rotate grid point by index */
/* ---------------------------*/
long grg_rotate_grid_index(const long grid_index,
MATCONST long rotation[3][3],
const long D_diag[3],
const long PS[3])
{
return rotate_grid_index(grid_index, rotation, D_diag, PS);
}
/* -----------------------------*/
/* Find irreducible grid points */
/* -----------------------------*/
void grg_get_ir_grid_map(long ir_grid_indices[],
MATCONST long (*rotations)[3][3],
const long num_rot,
const long D_diag[3],
const long PS[3])
{
get_ir_grid_map(ir_grid_indices,
rotations,
num_rot,
D_diag,
PS);
}
/* Extract unique rotation matrices and transpose them. */
/* When is_time_reversal == 1, inverse of the extracted matrices are */
/* included. */
/* Return 0 if failed */
long grg_get_reciprocal_point_group(long rec_rotations[48][3][3],
MATCONST long (*rotations)[3][3],
const long num_rot,
const long is_time_reversal)
{
long i, j, num_rot_ret, inv_exist;
MATCONST long inversion[3][3] = {
{-1, 0, 0 },
{ 0,-1, 0 },
{ 0, 0,-1 }
};
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])) {
goto escape;
}
}
if (num_rot_ret == 48) {
goto err;
}
mat_copy_matrix_l3(rec_rotations[num_rot_ret], rotations[i]);
num_rot_ret++;
escape:
;
}
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])) {
inv_exist = 1;
break;
}
}
if (!inv_exist) {
if (num_rot_ret > 24) {
goto err;
}
for (i = 0; i < num_rot_ret; i++) {
mat_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]);
}
return num_rot_ret;
err:
return 0;
}
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]);
}
}
static void reduce_double_grid_address(long address_double[3],
const long D_diag[3])
{
long i;
for (i = 0; i < 3; i++) {
address_double[i] = mat_modulo_l(address_double[i], 2 * D_diag[i]);
}
}
static long get_double_grid_index(const long address_double[3],
const long D_diag[3],
const long PS[3])
{
long address[3];
grg_get_grid_address(address,
address_double,
D_diag,
PS);
return get_grid_index_from_address(address, D_diag);
}
/* Here address elements have to be zero or positive. */
/* Therefore reduction to interval [0, D_diag[i]) has to be */
/* done outside of this function. */
/* See kgrid.h about GRID_ORDER_XYZ information. */
static long get_grid_index_from_address(const long address[3],
const long D_diag[3])
{
#ifndef GRID_ORDER_XYZ
return (address[2] * D_diag[0] * D_diag[1]
+ address[1] * D_diag[0] + address[0]);
#else
return (address[0] * D_diag[1] * D_diag[2]
+ address[1] * D_diag[2] + address[2]);
#endif
}
static void get_all_grid_addresses(long grid_address[][3],
const long D_diag[3])
{
long i, j, k, grid_index;
long address[3];
for (i = 0; i < D_diag[0]; i++) {
address[0] = i;
for (j = 0; j < D_diag[1]; j++) {
address[1] = j;
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);
}
}
}
}
/* See grg_get_grid_address_from_index */
static void get_grid_address_from_index(long address[3],
const long grid_index,
const long D_diag[3])
{
long nn;
#ifndef GRID_ORDER_XYZ
nn = D_diag[0] * D_diag[1];
address[0] = grid_index % D_diag[0];
address[2] = grid_index / nn;
address[1] = (grid_index - address[2] * nn) / D_diag[0];
#else
nn = D_diag[1] * D_diag[2];
address[2] = grid_index % D_diag[2];
address[0] = grid_index / nn;
address[1] = (grid_index - address[0] * nn) / D_diag[2];
#endif
}
/* Usually address has to be reduced to [0, D_diag[i]) */
/* by calling reduce_grid_address after this operation. */
static void get_grid_address(long address[3],
const long address_double[3],
const long PS[3])
{
long i;
for (i = 0; i < 3; i++) {
address[i] = (address_double[i] - PS[i]) / 2;
}
}
/* Usually address_double has to be reduced to [0, 2*D_diag[i]) */
/* by calling reduce_double_grid_address after this operation. */
static void get_double_grid_address(long address_double[3],
const long address[3],
const long PS[3])
{
long i;
for (i = 0; i < 3; i++) {
address_double[i] = address[i] * 2 + PS[i];
}
}
static long rotate_grid_index(const long grid_index,
MATCONST long rotation[3][3],
const long D_diag[3],
const long PS[3])
{
long adrs[3], dadrs[3], dadrs_rot[3];
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);
return get_double_grid_index(dadrs_rot, D_diag, PS);
}
static void get_ir_grid_map(long ir_grid_indices[],
MATCONST long (*rotations)[3][3],
const long num_rot,
const long D_diag[3],
const long PS[3])
{
long gp, num_gp, r_gp;
long i;
num_gp = D_diag[0] * D_diag[1] * D_diag[2];
for (gp = 0; gp < num_gp; gp++) {
ir_grid_indices[gp] = num_gp;
}
/* Do not simply multithreaded this for-loop. */
/* This algorithm contains race condition in different gp's. */
for (gp = 0; gp < num_gp; gp++) {
for (i = 0; i < num_rot; i++) {
r_gp = rotate_grid_index(gp, rotations[i], D_diag, PS);
if (r_gp < gp) {
ir_grid_indices[gp] = ir_grid_indices[r_gp];
break;
}
}
if (ir_grid_indices[gp] == num_gp) {
ir_grid_indices[gp] = gp;
}
}
}
static long mat_get_determinant_l3(MATCONST 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(MATCONST 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], MATCONST 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], MATCONST 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],
MATCONST long a[3][3],
MATCONST 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(MATCONST long a[3][3],
MATCONST 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(MATCONST long a[3][3],
MATCONST 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],
MATCONST 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], MATCONST 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],
MATCONST 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],
MATCONST long a[3][3],
MATCONST 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],
MATCONST long a[3][3],
MATCONST 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],
MATCONST double a[3][3],
MATCONST 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], MATCONST 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], MATCONST 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;
}

90
c/grgrid.h Normal file
View File

@ -0,0 +1,90 @@
/* Copyright (C) 2020 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 __grgrid_H__
#define __grgrid_H__
#include <stddef.h>
#ifndef MATCONST
#define MATCONST
#endif
#ifdef MATWARNING
#define warning_print(...) fprintf(stderr,__VA_ARGS__)
#else
#define warning_print(...)
#endif
long grg_get_snf3x3(long D_diag[3],
long P[3][3],
long Q[3][3],
MATCONST long A[3][3]);
long grg_transform_rotations(long (*transformed_rots)[3][3],
MATCONST long (*rotations)[3][3],
const long num_rot,
const long D_diag[3],
MATCONST 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],
const long PS[3]);
void grg_get_grid_address(long address[3],
const long address_double[3],
const long D_diag[3],
const long PS[3]);
long grg_get_grid_index(const long address[3],
const long D_diag[3]);
long grg_get_double_grid_index(const long address_double[3],
const long D_diag[3],
const long PS[3]);
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,
MATCONST long rotations[3][3],
const long D_diag[3],
const long PS[3]);
void grg_get_ir_grid_map(long ir_grid_indices[],
MATCONST 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],
MATCONST long (*rotations)[3][3],
const long num_rot,
const long is_time_reversal);
#endif

16
c/kpoint.c
View File

@ -36,7 +36,7 @@
#include <stdlib.h>
#include <stddef.h>
#include "kpoint.h"
#include "rgrid.h"
#include "grgrid.h"
#ifdef KPTWARNING
#include <stdio.h>
@ -542,11 +542,11 @@ static long get_ir_reciprocal_mesh_normal(long grid_address[][3],
long j;
long address_double[3], address_double_rot[3];
rgd_get_all_grid_addresses(grid_address, mesh);
grg_get_all_grid_addresses(grid_address, mesh);
#pragma omp parallel for private(j, grid_point_rot, address_double, address_double_rot)
for (i = 0; i < mesh[0] * mesh[1] * (long)(mesh[2]); i++) {
rgd_get_double_grid_address(address_double,
grg_get_double_grid_address(address_double,
grid_address[i],
mesh,
is_shift);
@ -555,7 +555,7 @@ static long get_ir_reciprocal_mesh_normal(long grid_address[][3],
multiply_matrix_vector_l3(address_double_rot,
rot_reciprocal->mat[j],
address_double);
grid_point_rot = rgd_get_double_grid_index(address_double_rot, mesh);
grid_point_rot = grg_get_double_grid_index(address_double_rot, mesh, is_shift);
if (grid_point_rot < ir_mapping_table[i]) {
#ifdef _OPENMP
ir_mapping_table[i] = grid_point_rot;
@ -585,7 +585,7 @@ get_ir_reciprocal_mesh_distortion(long grid_address[][3],
/* divisor, long_address_double, and long_address_double_rot have */
/* long integer type to treat dense mesh. */
rgd_get_all_grid_addresses(grid_address, mesh);
grg_get_all_grid_addresses(grid_address, mesh);
for (j = 0; j < 3; j++) {
divisor[j] = mesh[(j + 1) % 3] * mesh[(j + 2) % 3];
@ -593,7 +593,7 @@ get_ir_reciprocal_mesh_distortion(long grid_address[][3],
#pragma omp parallel for private(j, k, grid_point_rot, address_double, address_double_rot, long_address_double, long_address_double_rot)
for (i = 0; i < mesh[0] * mesh[1] * (long)(mesh[2]); i++) {
rgd_get_double_grid_address(address_double,
grg_get_double_grid_address(address_double,
grid_address[i],
mesh,
is_shift);
@ -623,7 +623,7 @@ get_ir_reciprocal_mesh_distortion(long grid_address[][3],
}
if (indivisible) {continue;}
grid_point_rot =
rgd_get_double_grid_index(address_double_rot, mesh);
grg_get_double_grid_index(address_double_rot, mesh, is_shift);
if (grid_point_rot < ir_mapping_table[i]) {
#ifdef _OPENMP
ir_mapping_table[i] = grid_point_rot;
@ -698,7 +698,7 @@ static long relocate_BZ_grid_address(long bz_grid_address[][3],
grid_address[i][k] + bz_search_space[j][k] * mesh[k];
bz_address_double[k] = bz_grid_address[gp][k] * 2 + is_shift[k];
}
bzgp = rgd_get_double_grid_index(bz_address_double, bzmesh);
bzgp = grg_get_double_grid_index(bz_address_double, bzmesh, is_shift);
bz_map[bzgp] = gp;
if (j != min_index) {
boundary_num_gp++;

4
c/phono3py.c
View File

@ -44,7 +44,7 @@
#include "kpoint.h"
#include "pp_collision.h"
#include "real_self_energy.h"
#include "rgrid.h"
#include "grgrid.h"
#include "tetrahedron_method.h"
#include "triplet.h"
#include "triplet_iw.h"
@ -605,7 +605,7 @@ void ph3py_get_integration_weight_with_sigma(double *iw,
long ph3py_get_grid_index_from_address(const long address[3],
const long mesh[3])
{
return rgd_get_single_grid_index(address, mesh);
return grg_get_grid_index(address, mesh);
}

174
c/rgrid.c
View File

@ -1,174 +0,0 @@
/* Copyright (C) 2015 Atsushi Togo */
/* All rights reserved. */
/* This file was originally part of spglib and 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 <stdio.h>
#include <stddef.h>
#include <assert.h>
#include "rgrid.h"
static void get_all_grid_addresses(long grid_address[][3], const long mesh[3]);
static long get_double_grid_index(const long address_double[3],
const long mesh[3]);
static long get_single_grid_index(const long address[3],
const long mesh[3]);
static void reduce_grid_address(long address[3], const long mesh[3]);
static void reduce_double_grid_address(long address[3], const long mesh[3]);
static long mat_modulo_l(const long a, const long b);
void rgd_get_all_grid_addresses(long grid_address[][3], const long mesh[3])
{
get_all_grid_addresses(grid_address, mesh);
}
long rgd_get_double_grid_index(const long address_double[3],
const long mesh[3])
{
return get_double_grid_index(address_double, mesh);
}
long rgd_get_single_grid_index(const long address[3],
const long mesh[3])
{
long i;
long address_mod[3];
for (i = 0; i < 3; i++) {
address_mod[i] = mat_modulo_l(address[i], mesh[i]);
}
return get_single_grid_index(address_mod, mesh);
}
void rgd_get_double_grid_address(long address_double[3],
const long address[3],
const long mesh[3],
const long is_shift[3])
{
long i;
for (i = 0; i < 3; i++) {
address_double[i] = address[i] * 2 + (is_shift[i] != 0);
}
reduce_double_grid_address(address_double, mesh);
}
static void get_all_grid_addresses(long grid_address[][3], const long mesh[3])
{
long i, j, k;
long grid_index;
long address[3];
for (i = 0; i < mesh[0]; i++) {
address[0] = i;
for (j = 0; j < mesh[1]; j++) {
address[1] = j;
for (k = 0; k < mesh[2]; k++) {
address[2] = k;
grid_index = get_single_grid_index(address, mesh);
assert(mesh[0] * mesh[1] * mesh[2] > grid_index);
grid_address[grid_index][0] = address[0];
grid_address[grid_index][1] = address[1];
grid_address[grid_index][2] = address[2];
reduce_grid_address(grid_address[grid_index], mesh);
}
}
}
}
static long get_double_grid_index(const long address_double[3],
const long mesh[3])
{
long i;
long address[3];
for (i = 0; i < 3; i++) {
if (address_double[i] % 2 == 0) {
address[i] = address_double[i] / 2;
} else {
address[i] = (address_double[i] - 1) / 2;
}
}
return rgd_get_single_grid_index(address, mesh);
}
static long get_single_grid_index(const long address[3],
const long mesh[3])
{
#ifndef GRID_ORDER_XYZ
return (address[2] * mesh[0] * (long)(mesh[1])
+ address[1] * mesh[0] + address[0]);
#else
return (address[0] * mesh[1] * (long)(mesh[2])
+ address[1] * mesh[2] + address[2]);
#endif
}
static void reduce_grid_address(long address[3], const long mesh[3])
{
long i;
for (i = 0; i < 3; i++) {
#ifndef GRID_BOUNDARY_AS_NEGATIVE
address[i] -= mesh[i] * (address[i] > mesh[i] / 2);
#else
address[i] -= mesh[i] * (address[i] > (mesh[i] - 1) / 2);
#endif
}
}
static void reduce_double_grid_address(long address[3], const long mesh[3])
{
long i;
for (i = 0; i < 3; i++) {
#ifndef GRID_BOUNDARY_AS_NEGATIVE
address[i] -= 2 * mesh[i] * (address[i] > mesh[i]);
#else
address[i] -= 2 * mesh[i] * (address[i] > mesh[i] - 1);
#endif
}
}
static long mat_modulo_l(const long a, const long b)
{
long c;
c = a % b;
if (c < 0) {
c += b;
}
return c;
}

83
c/rgrid.h
View File

@ -1,83 +0,0 @@
/* Copyright (C) 2015 Atsushi Togo */
/* All rights reserved. */
/* This file was originally part of spglib and 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 __rgrid_H__
#define __rgrid_H__
/* #define GRID_ORDER_XYZ */
/* This changes behaviour of index order of address. */
/* Without GRID_ORDER_XYZ, left most element of address runs first. */
/* grid_address (e.g. 4x4x4 mesh, unless GRID_ORDER_XYZ is defined) */
/* [[ 0 0 0] */
/* [ 1 0 0] */
/* [ 2 0 0] */
/* [-1 0 0] */
/* [ 0 1 0] */
/* [ 1 1 0] */
/* [ 2 1 0] */
/* [-1 1 0] */
/* .... ] */
/* */
/* With GRID_ORDER_XYZ, right most element of address runs first. */
/* grid_address (e.g. 4x4x4 mesh, if GRID_ORDER_XYZ is defined) */
/* [[ 0 0 0] */
/* [ 0 0 1] */
/* [ 0 0 2] */
/* [ 0 0 -1] */
/* [ 0 1 0] */
/* [ 0 1 1] */
/* [ 0 1 2] */
/* [ 0 1 -1] */
/* .... ] */
/* #define GRID_BOUNDARY_AS_NEGATIVE */
/* This changes the behaviour of address elements on the surface of */
/* parallelepiped. */
/* For odd mesh number, this affects nothing, e.g., [-2, -1, 0, 1, 2]. */
/* regardless of with and without GRID_BOUNDARY_AS_NEGATIVE. */
/* For even mesh number, this affects as follows: */
/* without GRID_BOUNDARY_AS_NEGATIVE, e.g., [-2, -1, 0, 1, 2, 3]. */
/* with GRID_BOUNDARY_AS_NEGATIVE, e.g., [-3, -2, -1, 0, 1, 2]. */
void rgd_get_all_grid_addresses(long grid_address[][3], const long mesh[3]);
long rgd_get_double_grid_index(const long address_double[3],
const long mesh[3]);
long rgd_get_single_grid_index(const long address[3],
const long mesh[3]);
void rgd_get_double_grid_address(long address_double[3],
const long address[3],
const long mesh[3],
const long is_shift[3]);
#endif

679
c/snf3x3.c Normal file
View File

@ -0,0 +1,679 @@
#include <assert.h>
/* #include <stdio.h> */
#ifndef SNF3x3CONST
#define SNF3x3CONST
#endif
static void initialize_PQ(long P[3][3], long Q[3][3]);
static int first(long A[3][3], long P[3][3], long Q[3][3]);
static void first_one_loop(long A[3][3], long P[3][3], long Q[3][3]);
static void first_column(long A[3][3], long P[3][3]);
static void zero_first_column(long L[3][3], const int j,
SNF3x3CONST long A[3][3]);
static int search_first_pivot(SNF3x3CONST long A[3][3]);
static void first_finalize(long L[3][3], SNF3x3CONST long A[3][3]);
static int second(long A[3][3], long P[3][3], long Q[3][3]);
static void second_one_loop(long A[3][3], long P[3][3], long Q[3][3]);
static void second_column(long A[3][3], long P[3][3]);
static void zero_second_column(long L[3][3], SNF3x3CONST long A[3][3]);
static void second_finalize(long L[3][3], SNF3x3CONST long A[3][3]);
static void finalize(long A[3][3], long P[3][3], long Q[3][3]);
static void finalize_sort(long A[3][3], long P[3][3], long Q[3][3]);
static void finalize_disturb(long A[3][3], long Q[3][3], const int i, const int j);
static void disturb_rows(long L[3][3], const int i, const int j);
static void swap_diag_elems(long A[3][3], long P[3][3], long Q[3][3], const int i, const int j);
static void make_diagA_positive(long A[3][3], long P[3][3]);
static void flip_PQ(long P[3][3], long Q[3][3]);
static void swap_rows(long L[3][3], const int i, const int j);
static void set_zero(long L[3][3],
const int i, const int j, const long a, const long b,
const long r, const long s, const long t);
static void extended_gcd(long retvals[3], const long a, const long b);
static void extended_gcd_step(long vals[6]);
static void flip_sign_row(long L[3][3], const int i);
static void transpose(long m[3][3]);
static void matmul(long m[3][3],
SNF3x3CONST long a[3][3],
SNF3x3CONST long b[3][3]);
static long det(SNF3x3CONST long m[3][3]);
/* static void test_set_A(long A[3][3]);
* static void test_show_A(SNF3x3CONST long A[3][3]);
* static void test_extended_gcd(void);
* static void test_transpose(void);
* static void test_swap_rows(void);
* static void test_zero_first_column(void);
* static void test_first_column(void);
* static void test_first_one_loop(void);
* static void test_first(void);
* static void test_second_column(void);
* static void test_second_one_loop(void);
* static void test_second(void); */
int snf3x3(long A[3][3], long P[3][3], long Q[3][3])
{
int i;
initialize_PQ(P, Q);
for (i = 0; i < 100; i++) {
if (first(A, P, Q)) {
if (second(A, P, Q)) {
finalize(A, P, Q);
transpose(Q);
goto succeeded;
}
}
}
return 0;
succeeded:
return 1;
}
static void initialize_PQ(long P[3][3], long Q[3][3])
{
int i, j;
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
if (i == j) {
P[i][j] = 1;
Q[i][j] = 1;
} else {
P[i][j] = 0;
Q[i][j] = 0;
}
}
}
}
static int first(long A[3][3], long P[3][3], long Q[3][3])
{
long L[3][3];
first_one_loop(A, P, Q);
/* rows and columns are all zero except for the pivot */
if ((A[1][0] == 0) && (A[2][0] == 0)) {
return 1;
}
/* columns of the pivot are assumed zero because of first_one_loop. */
/* rows of the pivot are non-zero, but divisible by the pivot. */
/* first_finalize makes the rows be zero. */
if ((A[1][0] % A[0][0] == 0) && (A[2][0] % A[0][0] == 0)) {
first_finalize(L, A);
matmul(A, L, A);
matmul(P, L, P);
return 1;
}
return 0;
}
static void first_one_loop(long A[3][3], long P[3][3], long Q[3][3])
{
first_column(A, P);
transpose(A);
first_column(A, Q);
transpose(A);
}
static void first_column(long A[3][3], long P[3][3])
{
int i;
long L[3][3];
i = search_first_pivot(A);
if (i > 0) {
swap_rows(L, 0, i);
matmul(A, L, A);
matmul(P, L, P);
}
if (i < 0) {
goto err;
}
if (A[1][0] != 0) {
zero_first_column(L, 1, A);
matmul(A, L, A);
matmul(P, L, P);
}
if (A[2][0] != 0) {
zero_first_column(L, 2, A);
matmul(A, L, A);
matmul(P, L, P);
}
err:
;
}
static void zero_first_column(long L[3][3], const int j,
SNF3x3CONST long A[3][3])
{
long vals[3];
extended_gcd(vals, A[0][0], A[j][0]);
set_zero(L, 0, j, A[0][0], A[j][0], vals[0], vals[1], vals[2]);
}
static int search_first_pivot(SNF3x3CONST long A[3][3])
{
int i;
for (i = 0; i < 3; i++) {
if (A[i][0] != 0) {
return i;
}
}
return -1;
}
static void first_finalize(long L[3][3], SNF3x3CONST long A[3][3])
{
L[0][0] = 1;
L[0][1] = 0;
L[0][2] = 0;
L[1][0] = -A[1][0] / A[0][0];
L[1][1] = 1;
L[1][2] = 0;
L[2][0] = -A[2][0] / A[0][0];
L[2][1] = 0;
L[2][2] = 1;
}
static int second(long A[3][3], long P[3][3], long Q[3][3])
{
long L[3][3];
second_one_loop(A, P, Q);
if (A[2][1] == 0) {
return 1;
}
if (A[2][1] % A[1][1] == 0) {
second_finalize(L, A);
matmul(A, L, A);
matmul(P, L, P);
return 1;
}
return 0;
}
static void second_one_loop(long A[3][3], long P[3][3], long Q[3][3])
{
second_column(A, P);
transpose(A);
second_column(A, Q);
transpose(A);
}
static void second_column(long A[3][3], long P[3][3])
{
long L[3][3];
if ((A[1][1] == 0) && (A[2][1] != 0)) {
swap_rows(L, 1, 2);
matmul(A, L, A);
matmul(P, L, P);
}
if (A[2][1] != 0) {
zero_second_column(L, A);
matmul(A, L, A);
matmul(P, L, P);
}
}
static void zero_second_column(long L[3][3], SNF3x3CONST long A[3][3])
{
long vals[3];
extended_gcd(vals, A[1][1], A[2][1]);
set_zero(L, 1, 2, A[1][1], A[2][1], vals[0], vals[1], vals[2]);
}
static void second_finalize(long L[3][3], SNF3x3CONST long A[3][3])
{
L[0][0] = 1;
L[0][1] = 0;
L[0][2] = 0;
L[1][0] = 0;
L[1][1] = 1;
L[1][2] = 0;
L[2][0] = 0;
L[2][1] = -A[2][1] / A[1][1];
L[2][2] = 1;
}
static void finalize(long A[3][3], long P[3][3], long Q[3][3])
{
make_diagA_positive(A, P);
finalize_sort(A, P, Q);
finalize_disturb(A, Q, 0, 1);
first(A, P, Q);
finalize_sort(A, P, Q);
finalize_disturb(A, Q, 1, 2);
second(A, P, Q);
flip_PQ(P, Q);
}
static void finalize_sort(long A[3][3], long P[3][3], long Q[3][3])
{
if (A[0][0] > A[1][1]) {
swap_diag_elems(A, P, Q, 0, 1);
}
if (A[1][1] > A[2][2]) {
swap_diag_elems(A, P, Q, 1, 2);
}
if (A[0][0] > A[1][1]) {
swap_diag_elems(A, P, Q, 0, 1);
}
}
static void finalize_disturb(long A[3][3], long Q[3][3], const int i, const int j)
{
long L[3][3];
if (A[j][j] % A[i][i] != 0) {
transpose(A);
disturb_rows(L, i, j);
matmul(A, L, A);
matmul(Q, L, Q);
transpose(A);
}
}
static void disturb_rows(long L[3][3], const int i, const int j)
{
L[0][0] = 1;
L[0][1] = 0;
L[0][2] = 0;
L[1][0] = 0;
L[1][1] = 1;
L[1][2] = 0;
L[2][0] = 0;
L[2][1] = 0;
L[2][2] = 1;
L[i][i] = 1;
L[i][j] = 1;
L[j][i] = 0;
L[j][j] = 1;
}
static void swap_diag_elems(long A[3][3], long P[3][3], long Q[3][3], const int i, const int j)
{
long L[3][3];
swap_rows(L, i, j);
matmul(A, L, A);
matmul(P, L, P);
transpose(A);
swap_rows(L, i, j);
matmul(A, L, A);
matmul(Q, L, Q);
transpose(A);
}
static void make_diagA_positive(long A[3][3], long P[3][3])
{
int i;
long L[3][3];
for (i = 0; i < 3; i++) {
if (A[i][i] < 0) {
flip_sign_row(L, i);
matmul(A, L, A);
matmul(P, L, P);
}
}
}
static void flip_PQ(long P[3][3], long Q[3][3])
{
int i, j;
if (det(P) < 0) {
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
P[i][j] *= -1;
Q[i][j] *= -1;
}
}
}
}
static void swap_rows(long L[3][3], const int r1, const int r2)
{
L[0][0] = 1;
L[0][1] = 0;
L[0][2] = 0;
L[1][0] = 0;
L[1][1] = 1;
L[1][2] = 0;
L[2][0] = 0;
L[2][1] = 0;
L[2][2] = 1;
L[r1][r1] = 0;
L[r2][r2] = 0;
L[r1][r2] = 1;
L[r2][r1] = 1;
}
static void set_zero(long L[3][3],
const int i, const int j, const long a, const long b,
const long r, const long s, const long t)
{
L[0][0] = 1;
L[0][1] = 0;
L[0][2] = 0;
L[1][0] = 0;
L[1][1] = 1;
L[1][2] = 0;
L[2][0] = 0;
L[2][1] = 0;
L[2][2] = 1;
L[i][i] = s;
L[i][j] = t;
L[j][i] = -b / r;
L[j][j] = a / r;
}
/**
* Extended Euclidean algorithm
* See https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm
*/
static void extended_gcd(long retvals[3], const long a, const long b)
{
int i;
long vals[6];
vals[0] = a; /* r0 */
vals[1] = b; /* r1 */
vals[2] = 1; /* s0 */
vals[3] = 0; /* s1 */
vals[4] = 0; /* t0 */
vals[5] = 1; /* t1 */
for (i = 0; i < 1000; i++) {
extended_gcd_step(vals);
if (vals[1] == 0) {
break;
}
}
retvals[0] = vals[0];
retvals[1] = vals[2];
retvals[2] = vals[4];
assert(vals[0] == a * vals[2] + b * vals[4]);
}
static void extended_gcd_step(long vals[6])
{
long q, r2, s2, t2;
q = vals[0] / vals[1];
r2 = vals[0] % vals[1];
if (r2 < 0) {
if (vals[1] > 0) {
r2 += vals[1];
q -= 1;
}
if (vals[1] < 0) {
r2 -= vals[1];
q += 1;
}
}
s2 = vals[2] - q * vals[3];
t2 = vals[4] - q * vals[5];
vals[0] = vals[1];
vals[1] = r2;
vals[2] = vals[3];
vals[3] = s2;
vals[4] = vals[5];
vals[5] = t2;
}
static void flip_sign_row(long L[3][3], const int i)
{
L[0][0] = 1;
L[0][1] = 0;
L[0][2] = 0;
L[1][0] = 0;
L[1][1] = 1;
L[1][2] = 0;
L[2][0] = 0;
L[2][1] = 0;
L[2][2] = 1;
L[i][i] = -1;
}
/**
* Matrix operation utils
*/
static void transpose(long m[3][3])
{
long tmp;
int i, j;
for (i = 0; i < 3; i++) {
for (j = i; j < 3; j++) {
tmp = m[i][j];
m[i][j] = m[j][i];
m[j][i] = tmp;
}
}
}
static void matmul(long m[3][3],
SNF3x3CONST long a[3][3],
SNF3x3CONST long b[3][3])
{
int i, j;
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];
}
}
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
m[i][j] = c[i][j];
}
}
}
static long det(SNF3x3CONST long m[3][3])
{
return m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1])
+ m[0][1] * (m[1][2] * m[2][0] - m[1][0] * m[2][2])
+ m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0]);
}
/* int main()
* {
* test_extended_gcd();
* test_transpose();
* test_swap_rows();
* test_zero_first_column();
* test_first_column();
* test_first_one_loop();
* test_first();
* test_second_column();
* test_second_one_loop();
* test_second();
* }
*
* static void test_set_A(long A[3][3])
* {
* int i, j;
*
* for (i = 0; i < 3; i++) {
* for (j = 0; j < 3; j++) {
* A[i][j] = i * 3 + j;
* }
* }
* A[0][0] = 1; /\* to avoid det(A) = 0 *\/
* }
*
* static void test_show_A(SNF3x3CONST long A[3][3])
* {
* int i, j;
*
* for (i = 0; i < 3; i++) {
* for (j = 0; j < 3; j++) {
* printf("%d ", A[i][j]);
* }
* printf("\n");
* }
* }
*
* static void test_extended_gcd(void)
* {
* int vals[3];
*
* printf("Test extended_gcd\n");
* extended_gcd(vals, 4 , -3);
* printf("%d %d %d\n", vals[0], vals[1], vals[2]);
* }
*
* static void test_transpose(void)
* {
* int i, j;
* long A[3][3];
*
* printf("Test transpose\n");
*
* test_set_A(A);
* test_show_A(A);
* transpose(A);
* test_show_A(A);
* }
*
* static void test_swap_rows(void)
* {
* int i, j;
* long A[3][3], L[3][3];
*
* printf("Test swap_rows 1 <-> 2\n");
*
* test_set_A(A);
* test_show_A(A);
* swap_rows(L, 0, 1);
* matmul(A, L, A);
* test_show_A(A);
* }
*
* static void test_zero_first_column(void)
* {
* int i, j;
* long A[3][3], L[3][3];
*
* printf("Test zero_first_column\n");
*
* test_set_A(A);
* test_show_A(A);
* zero_first_column(L, 1, A);
* matmul(A, L, A);
* test_show_A(A);
* zero_first_column(L, 2, A);
* matmul(A, L, A);
* test_show_A(A);
* }
*
* static void test_first_column(void)
* {
* int i, j;
* long A[3][3], P[3][3];
*
* printf("Test first_column\n");
*
* test_set_A(A);
* test_show_A(A);
* first_column(A, P);
* test_show_A(A);
* transpose(A);
* first_column(A, P);
* transpose(A);
* test_show_A(A);
* }
*
* static void test_first_one_loop(void)
* {
* int i, j;
* long A[3][3], P[3][3], Q[3][3];
*
* printf("Test first_one_loop\n");
*
* test_set_A(A);
* test_show_A(A);
* first_one_loop(A, P, Q);
* test_show_A(A);
* }
*
* static void test_first(void)
* {
* int i, j;
* long A[3][3], P[3][3], Q[3][3];
*
* printf("Test first\n");
*
* test_set_A(A);
* test_show_A(A);
* first(A, P, Q);
* test_show_A(A);
* }
*
* static void test_second_column(void)
* {
* int i, j;
* long A[3][3], P[3][3], Q[3][3];
*
* printf("Test second_column\n");
*
* test_set_A(A);
* test_show_A(A);
* first(A, P, Q);
* test_show_A(A);
* second_column(A, P);
* test_show_A(A);
* }
*
* static void test_second_one_loop(void)
* {
* int i, j;
* long A[3][3], P[3][3], Q[3][3];
*
* printf("Test second_one_loop\n");
*
* test_set_A(A);
* test_show_A(A);
* first(A, P, Q);
* test_show_A(A);
* second_one_loop(A, P, Q);
* test_show_A(A);
* }
*
* static void test_second(void)
* {
* int i, j;
* long A[3][3], P[3][3], Q[3][3];
*
* printf("Test second\n");
*
* test_set_A(A);
* test_show_A(A);
* first(A, P, Q);
* test_show_A(A);
* second(A, P, Q);
* test_show_A(A);
* } */

17
c/snf3x3.h Normal file
View File

@ -0,0 +1,17 @@
#ifndef __snf3x3_H__
#define __snf3x3_H__
#ifdef __cplusplus
extern "C" {
#endif
#define SNF3X3_MAJOR_VERSION 0
#define SNF3X3_MINOR_VERSION 1
#define SNF3X3_MICRO_VERSION 0
int snf3x3(long A[3][3], long P[3][3], long Q[3][3]);
#ifdef __cplusplus
}
#endif
#endif

9
c/triplet_iw.c
View File

@ -33,7 +33,7 @@
/* POSSIBILITY OF SUCH DAMAGE. */
#include <math.h>
#include "rgrid.h"
#include "grgrid.h"
#include "phonoc_utils.h"
#include "triplet.h"
#include "triplet_iw.h"
@ -217,13 +217,14 @@ tpi_get_neighboring_grid_points(long neighboring_grid_points[],
TPLCONST long bz_grid_address[][3],
const long bz_map[])
{
long bzmesh[3], address_double[3], bz_address_double[3];
long bzmesh[3], address_double[3], bz_address_double[3], PS[3];
long i, j, bz_gp, prod_bz_mesh;
prod_bz_mesh = 1;
for (i = 0; i < 3; i++) {
bzmesh[i] = mesh[i] * 2;
prod_bz_mesh *= bzmesh[i];
PS[i] = 0;
}
for (i = 0; i < num_relative_grid_address; i++) {
for (j = 0; j < 3; j++) {
@ -231,10 +232,10 @@ tpi_get_neighboring_grid_points(long neighboring_grid_points[],
relative_grid_address[i][j]) * 2;
bz_address_double[j] = address_double[j];
}
bz_gp = bz_map[rgd_get_double_grid_index(bz_address_double, bzmesh)];
bz_gp = bz_map[grg_get_double_grid_index(bz_address_double, bzmesh, PS)];
if (bz_gp == prod_bz_mesh) {
neighboring_grid_points[i] =
rgd_get_double_grid_index(address_double, mesh);
grg_get_double_grid_index(address_double, mesh, PS);
} else {
neighboring_grid_points[i] = bz_gp;
}

14
c/triplet_kpoint.c
View File

@ -37,7 +37,7 @@
#include <stddef.h>
#include <stdlib.h>
#include "kpoint.h"
#include "rgrid.h"
#include "grgrid.h"
#include "triplet.h"
#include "triplet_kpoint.h"
@ -308,7 +308,7 @@ static long get_ir_triplets_at_q(long *map_triplets,
for (j = 0; j < 3; j++) { /* q'' */
address_double2[j] = - address_double0[j] - address_double1[j];
}
third_q[i] = rgd_get_double_grid_index(address_double2, mesh);
third_q[i] = grg_get_double_grid_index(address_double2, mesh, is_shift);
}
num_ir_triplets = 0;
@ -355,12 +355,13 @@ static long get_BZ_triplets_at_q(long (*triplets)[3],
{
long i, num_ir;
long j, k;
long bz_address[3][3], bz_address_double[3], bzmesh[3];
long bz_address[3][3], bz_address_double[3], bzmesh[3], PS[3];
long *ir_grid_points;
ir_grid_points = NULL;
for (i = 0; i < 3; i++) {
PS[i] = 0;
bzmesh[i] = mesh[i] * 2;
}
@ -394,7 +395,7 @@ static long get_BZ_triplets_at_q(long (*triplets)[3],
bz_address_double[k] = bz_address[j][k] * 2;
}
triplets[i][j] =
bz_map[rgd_get_double_grid_index(bz_address_double, bzmesh)];
bz_map[grg_get_double_grid_index(bz_address_double, bzmesh, PS)];
}
}
@ -413,12 +414,13 @@ static long get_third_q_of_triplets_at_q(long bz_address[3][3],
long i, j, smallest_g, smallest_index, sum_g, delta_g[3];
long prod_bzmesh;
long bzgp[KPT_NUM_BZ_SEARCH_SPACE];
long bz_address_double[3];
long bz_address_double[3], PS[3];
prod_bzmesh = (long)bzmesh[0] * bzmesh[1] * bzmesh[2];
modulo_l3(bz_address[q_index], mesh);
for (i = 0; i < 3; i++) {
PS[i] = 0;
delta_g[i] = 0;
for (j = 0; j < 3; j++) {
delta_g[i] += bz_address[j][i];
@ -431,7 +433,7 @@ static long get_third_q_of_triplets_at_q(long bz_address[3][3],
bz_address_double[j] = (bz_address[q_index][j] +
bz_search_space[i][j] * mesh[j]) * 2;
}
bzgp[i] = bz_map[rgd_get_double_grid_index(bz_address_double, bzmesh)];
bzgp[i] = bz_map[grg_get_double_grid_index(bz_address_double, bzmesh, PS)];
}
for (i = 0; i < KPT_NUM_BZ_SEARCH_SPACE; i++) {

3
setup.py
View File

@ -208,6 +208,7 @@ include_dirs += include_dirs_lapacke
print("extra_link_args", extra_link_args)
sources_phono3py = ['c/_phono3py.c',
'c/collision_matrix.c',
'c/grgrid.c',
'c/fc3.c',
'c/imag_self_energy_with_g.c',
'c/interaction.c',
@ -220,7 +221,7 @@ sources_phono3py = ['c/_phono3py.c',
'c/real_self_energy.c',
'c/real_to_reciprocal.c',
'c/reciprocal_to_normal.c',
'c/rgrid.c',
'c/snf3x3.c',
'c/tetrahedron_method.c',
'c/triplet.c',
'c/triplet_iw.c',