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

Changed to 64bit grid point addressing by size_t

This commit is contained in:
Atsushi Togo 2019-02-04 11:54:28 +09:00
parent 11851b9662
commit 3f4e578f33
18 changed files with 1029 additions and 1222 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

50
c/_phonopy.c
View File

@ -34,6 +34,7 @@
#include <Python.h>
#include <stdio.h>
#include <stddef.h>
#include <math.h>
#include <float.h>
#include <numpy/arrayobject.h>
@ -1011,7 +1012,6 @@ static PyObject * py_get_thermal_properties(PyObject *self, PyObject *args)
int num_temp;
int i, j, k;
long sum_weights;
double f;
double *tp;
@ -1140,16 +1140,16 @@ static PyObject *py_thm_neighboring_grid_points(PyObject *self, PyObject *args)
PyArrayObject* py_mesh;
PyArrayObject* py_bz_grid_address;
PyArrayObject* py_bz_map;
int grid_point;
long grid_point;
int* relative_grid_points;
int (*relative_grid_address)[3];
int num_relative_grid_address;
int *mesh;
int (*bz_grid_address)[3];
int *bz_map;
size_t *bz_map_size_t;
size_t *relative_grid_points_size_t;
if (!PyArg_ParseTuple(args, "OiOOOO",
if (!PyArg_ParseTuple(args, "OlOOOO",
&py_relative_grid_points,
&grid_point,
&py_relative_grid_address,
@ -1159,20 +1159,20 @@ static PyObject *py_thm_neighboring_grid_points(PyObject *self, PyObject *args)
return NULL;
}
relative_grid_points = (int*)PyArray_DATA(py_relative_grid_points);
relative_grid_address = (int(*)[3])PyArray_DATA(py_relative_grid_address);
num_relative_grid_address = PyArray_DIMS(py_relative_grid_address)[0];
mesh = (int*)PyArray_DATA(py_mesh);
bz_grid_address = (int(*)[3])PyArray_DATA(py_bz_grid_address);
bz_map = (int*)PyArray_DATA(py_bz_map);
bz_map_size_t = (size_t*)PyArray_DATA(py_bz_map);
relative_grid_points_size_t = (size_t*)PyArray_DATA(py_relative_grid_points);
thm_get_neighboring_grid_points(relative_grid_points,
grid_point,
relative_grid_address,
num_relative_grid_address,
mesh,
bz_grid_address,
bz_map);
thm_get_dense_neighboring_grid_points(relative_grid_points_size_t,
grid_point,
relative_grid_address,
num_relative_grid_address,
mesh,
bz_grid_address,
bz_map_size_t);
Py_RETURN_NONE;
}
@ -1364,20 +1364,20 @@ static PyObject * py_tetrahedron_method_dos(PyObject *self, PyObject *args)
double* frequencies;
double* coef;
int (*grid_address)[3];
int num_gp;
int num_ir_gp;
size_t num_gp, num_ir_gp;
int num_coef;
int num_band;
int* grid_mapping_table;
size_t *grid_mapping_table;
int (*relative_grid_address)[4][3];
int is_shift[3] = {0, 0, 0};
int i, j, k, l, m, q, r, count;
size_t i, j, k, l, m, q, r, count;
size_t ir_gps[24][4];
int g_addr[3];
int ir_gps[24][4];
double tetrahedra[24][4];
int address_double[3];
int *gp2ir, *ir_grid_points, *weights;
size_t *gp2ir, *ir_grid_points;
int *weights;
double iw;
gp2ir = NULL;
@ -1402,17 +1402,17 @@ static PyObject * py_tetrahedron_method_dos(PyObject *self, PyObject *args)
freq_points = (double*)PyArray_DATA(py_freq_points);
num_freq_points = (int)PyArray_DIMS(py_freq_points)[0];
frequencies = (double*)PyArray_DATA(py_frequencies);
num_ir_gp = (int)PyArray_DIMS(py_frequencies)[0];
num_ir_gp = (size_t)PyArray_DIMS(py_frequencies)[0];
num_band = (int)PyArray_DIMS(py_frequencies)[1];
coef = (double*)PyArray_DATA(py_coef);
num_coef = (int)PyArray_DIMS(py_coef)[1];
grid_address = (int(*)[3])PyArray_DATA(py_grid_address);
num_gp = (int)PyArray_DIMS(py_grid_address)[0];
grid_mapping_table = (int*)PyArray_DATA(py_grid_mapping_table);
num_gp = (size_t)PyArray_DIMS(py_grid_address)[0];
grid_mapping_table = (size_t*)PyArray_DATA(py_grid_mapping_table);
relative_grid_address = (int(*)[4][3])PyArray_DATA(py_relative_grid_address);
gp2ir = (int*)malloc(sizeof(int) * num_gp);
ir_grid_points = (int*)malloc(sizeof(int) * num_ir_gp);
gp2ir = (size_t*)malloc(sizeof(size_t) * num_gp);
ir_grid_points = (size_t*)malloc(sizeof(size_t) * num_ir_gp);
weights = (int*)malloc(sizeof(int) * num_ir_gp);
count = 0;

565
c/_spglib.c
View File

@ -211,10 +211,10 @@ static PyObject * py_get_dataset(PyObject *self, PyObject *args)
{
int hall_number;
double symprec, angle_tolerance;
SpglibDataset *dataset;
PyArrayObject* lattice;
PyArrayObject* position;
PyArrayObject* atom_type;
PyArrayObject* py_lattice;
PyArrayObject* py_positions;
PyArrayObject* py_atom_types;
PyObject *array, *vec, *mat, *rot, *trans, *wyckoffs, *equiv_atoms;
PyObject *site_symmetry_symbols, *mapping_to_primitive;
PyObject *std_lattice, *std_types, *std_positions, *std_mapping_to_primitive;
@ -225,21 +225,22 @@ static PyObject * py_get_dataset(PyObject *self, PyObject *args)
double (*pos)[3];
int num_atom, len_list;
int* typat;
SpglibDataset *dataset;
if (!PyArg_ParseTuple(args, "OOOidd",
&lattice,
&position,
&atom_type,
&py_lattice,
&py_positions,
&py_atom_types,
&hall_number,
&symprec,
&angle_tolerance)) {
return NULL;
}
lat = (double(*)[3])PyArray_DATA(lattice);
pos = (double(*)[3])PyArray_DATA(position);
num_atom = PyArray_DIMS(position)[0];
typat = (int*)PyArray_DATA(atom_type);
lat = (double(*)[3])PyArray_DATA(py_lattice);
pos = (double(*)[3])PyArray_DATA(py_positions);
num_atom = PyArray_DIMS(py_positions)[0];
typat = (int*)PyArray_DATA(py_atom_types);
if ((dataset = spgat_get_dataset_with_hall_number(lat,
pos,
@ -400,26 +401,26 @@ static PyObject * py_get_dataset(PyObject *self, PyObject *args)
static PyObject * py_get_symmetry_from_database(PyObject *self, PyObject *args)
{
int hall_number;
PyArrayObject* rotation;
PyArrayObject* translation;
PyArrayObject* py_rotations;
PyArrayObject* py_translations;
int (*rot)[3][3];
double (*trans)[3];
int num_sym;
if (!PyArg_ParseTuple(args, "OOi",
&rotation,
&translation,
&py_rotations,
&py_translations,
&hall_number)) {
return NULL;
}
if (PyArray_DIMS(rotation)[0] < 192 || PyArray_DIMS(translation)[0] < 192) {
if (PyArray_DIMS(py_rotations)[0] < 192 || PyArray_DIMS(py_translations)[0] < 192) {
Py_RETURN_NONE;
}
rot = (int(*)[3][3])PyArray_DATA(rotation);
trans = (double(*)[3])PyArray_DATA(translation);
rot = (int(*)[3][3])PyArray_DATA(py_rotations);
trans = (double(*)[3])PyArray_DATA(py_translations);
num_sym = spg_get_symmetry_from_database(rot, trans, hall_number);
@ -472,7 +473,7 @@ static PyObject * py_get_spacegroup_type(PyObject *self, PyObject *args)
static PyObject * py_get_pointgroup(PyObject *self, PyObject *args)
{
PyArrayObject* rotations;
PyArrayObject* py_rotations;
int i, j;
int trans_mat[3][3];
@ -482,12 +483,12 @@ static PyObject * py_get_pointgroup(PyObject *self, PyObject *args)
int num_rot;
int ptg_num;
if (! PyArg_ParseTuple(args, "O", &rotations)) {
if (! PyArg_ParseTuple(args, "O", &py_rotations)) {
return NULL;
}
rot = (int(*)[3][3])PyArray_DATA(rotations);
num_rot = PyArray_DIMS(rotations)[0];
rot = (int(*)[3][3])PyArray_DATA(py_rotations);
num_rot = PyArray_DIMS(py_rotations)[0];
ptg_num = spg_get_pointgroup(symbol, trans_mat, rot, num_rot);
/* Transformation matrix */
@ -512,9 +513,9 @@ static PyObject * py_standardize_cell(PyObject *self, PyObject *args)
{
int num_atom, to_primitive, no_idealize;
double symprec, angle_tolerance;
PyArrayObject* lattice;
PyArrayObject* position;
PyArrayObject* atom_type;
PyArrayObject* py_lattice;
PyArrayObject* py_positions;
PyArrayObject* py_atom_types;
double (*lat)[3];
double (*pos)[3];
@ -522,9 +523,9 @@ static PyObject * py_standardize_cell(PyObject *self, PyObject *args)
int num_atom_std;
if (!PyArg_ParseTuple(args, "OOOiiidd",
&lattice,
&position,
&atom_type,
&py_lattice,
&py_positions,
&py_atom_types,
&num_atom,
&to_primitive,
&no_idealize,
@ -533,9 +534,9 @@ static PyObject * py_standardize_cell(PyObject *self, PyObject *args)
return NULL;
}
lat = (double(*)[3])PyArray_DATA(lattice);
pos = (double(*)[3])PyArray_DATA(position);
typat = (int*)PyArray_DATA(atom_type);
lat = (double(*)[3])PyArray_DATA(py_lattice);
pos = (double(*)[3])PyArray_DATA(py_positions);
typat = (int*)PyArray_DATA(py_atom_types);
num_atom_std = spgat_standardize_cell(lat,
pos,
@ -553,9 +554,9 @@ static PyObject * py_refine_cell(PyObject *self, PyObject *args)
{
int num_atom;
double symprec, angle_tolerance;
PyArrayObject* lattice;
PyArrayObject* position;
PyArrayObject* atom_type;
PyArrayObject* py_lattice;
PyArrayObject* py_positions;
PyArrayObject* py_atom_types;
double (*lat)[3];
double (*pos)[3];
@ -563,18 +564,18 @@ static PyObject * py_refine_cell(PyObject *self, PyObject *args)
int num_atom_std;
if (!PyArg_ParseTuple(args, "OOOidd",
&lattice,
&position,
&atom_type,
&py_lattice,
&py_positions,
&py_atom_types,
&num_atom,
&symprec,
&angle_tolerance)) {
return NULL;
}
lat = (double(*)[3])PyArray_DATA(lattice);
pos = (double(*)[3])PyArray_DATA(position);
typat = (int*)PyArray_DATA(atom_type);
lat = (double(*)[3])PyArray_DATA(py_lattice);
pos = (double(*)[3])PyArray_DATA(py_positions);
typat = (int*)PyArray_DATA(py_atom_types);
num_atom_std = spgat_refine_cell(lat,
pos,
@ -590,9 +591,9 @@ static PyObject * py_refine_cell(PyObject *self, PyObject *args)
static PyObject * py_find_primitive(PyObject *self, PyObject *args)
{
double symprec, angle_tolerance;
PyArrayObject* lattice;
PyArrayObject* position;
PyArrayObject* atom_type;
PyArrayObject* py_lattice;
PyArrayObject* py_positions;
PyArrayObject* py_atom_types;
double (*lat)[3];
double (*pos)[3];
@ -601,18 +602,18 @@ static PyObject * py_find_primitive(PyObject *self, PyObject *args)
int num_atom_prim;
if (!PyArg_ParseTuple(args, "OOOdd",
&lattice,
&position,
&atom_type,
&py_lattice,
&py_positions,
&py_atom_types,
&symprec,
&angle_tolerance)) {
return NULL;
}
lat = (double(*)[3])PyArray_DATA(lattice);
pos = (double(*)[3])PyArray_DATA(position);
num_atom = PyArray_DIMS(position)[0];
types = (int*)PyArray_DATA(atom_type);
lat = (double(*)[3])PyArray_DATA(py_lattice);
pos = (double(*)[3])PyArray_DATA(py_positions);
num_atom = PyArray_DIMS(py_positions)[0];
types = (int*)PyArray_DATA(py_atom_types);
num_atom_prim = spgat_find_primitive(lat,
pos,
@ -627,11 +628,11 @@ static PyObject * py_find_primitive(PyObject *self, PyObject *args)
static PyObject * py_get_symmetry(PyObject *self, PyObject *args)
{
double symprec, angle_tolerance;
PyArrayObject* lattice;
PyArrayObject* position;
PyArrayObject* rotation;
PyArrayObject* translation;
PyArrayObject* atom_type;
PyArrayObject* py_lattice;
PyArrayObject* py_positions;
PyArrayObject* py_rotations;
PyArrayObject* py_translations;
PyArrayObject* py_atom_types;
double (*lat)[3];
double (*pos)[3];
@ -643,23 +644,23 @@ static PyObject * py_get_symmetry(PyObject *self, PyObject *args)
int num_sym;
if (!PyArg_ParseTuple(args, "OOOOOdd",
&rotation,
&translation,
&lattice,
&position,
&atom_type,
&py_rotations,
&py_translations,
&py_lattice,
&py_positions,
&py_atom_types,
&symprec,
&angle_tolerance)) {
return NULL;
}
lat = (double(*)[3])PyArray_DATA(lattice);
pos = (double(*)[3])PyArray_DATA(position);
types = (int*)PyArray_DATA(atom_type);
num_atom = PyArray_DIMS(position)[0];
rot = (int(*)[3][3])PyArray_DATA(rotation);
trans = (double(*)[3])PyArray_DATA(translation);
num_sym_from_array_size = PyArray_DIMS(rotation)[0];
lat = (double(*)[3])PyArray_DATA(py_lattice);
pos = (double(*)[3])PyArray_DATA(py_positions);
types = (int*)PyArray_DATA(py_atom_types);
num_atom = PyArray_DIMS(py_positions)[0];
rot = (int(*)[3][3])PyArray_DATA(py_rotations);
trans = (double(*)[3])PyArray_DATA(py_translations);
num_sym_from_array_size = PyArray_DIMS(py_rotations)[0];
/* num_sym has to be larger than num_sym_from_array_size. */
num_sym = spgat_get_symmetry(rot,
@ -678,13 +679,13 @@ static PyObject * py_get_symmetry_with_collinear_spin(PyObject *self,
PyObject *args)
{
double symprec, angle_tolerance;
PyArrayObject* lattice;
PyArrayObject* position;
PyArrayObject* rotation;
PyArrayObject* translation;
PyArrayObject* atom_type;
PyArrayObject* magmom;
PyArrayObject* equiv_atoms_py;
PyArrayObject* py_lattice;
PyArrayObject* py_positions;
PyArrayObject* py_rotations;
PyArrayObject* py_translations;
PyArrayObject* py_atom_types;
PyArrayObject* py_magmoms;
PyArrayObject* py_equiv_atoms;
double (*lat)[3];
double (*pos)[3];
@ -698,27 +699,27 @@ static PyObject * py_get_symmetry_with_collinear_spin(PyObject *self,
int num_sym;
if (!PyArg_ParseTuple(args, "OOOOOOOdd",
&rotation,
&translation,
&equiv_atoms_py,
&lattice,
&position,
&atom_type,
&magmom,
&py_rotations,
&py_translations,
&py_equiv_atoms,
&py_lattice,
&py_positions,
&py_atom_types,
&py_magmoms,
&symprec,
&angle_tolerance)) {
return NULL;
}
lat = (double(*)[3])PyArray_DATA(lattice);
pos = (double(*)[3])PyArray_DATA(position);
spins = (double*)PyArray_DATA(magmom);
types = (int*)PyArray_DATA(atom_type);
num_atom = PyArray_DIMS(position)[0];
rot = (int(*)[3][3])PyArray_DATA(rotation);
trans = (double(*)[3])PyArray_DATA(translation);
equiv_atoms = (int*)PyArray_DATA(equiv_atoms_py);
num_sym_from_array_size = PyArray_DIMS(rotation)[0];
lat = (double(*)[3])PyArray_DATA(py_lattice);
pos = (double(*)[3])PyArray_DATA(py_positions);
spins = (double*)PyArray_DATA(py_magmoms);
types = (int*)PyArray_DATA(py_atom_types);
num_atom = PyArray_DIMS(py_positions)[0];
rot = (int(*)[3][3])PyArray_DATA(py_rotations);
trans = (double(*)[3])PyArray_DATA(py_translations);
equiv_atoms = (int*)PyArray_DATA(py_equiv_atoms);
num_sym_from_array_size = PyArray_DIMS(py_rotations)[0];
/* num_sym has to be larger than num_sym_from_array_size. */
num_sym = spgat_get_symmetry_with_collinear_spin(rot,
@ -739,8 +740,8 @@ static PyObject *
py_get_hall_number_from_symmetry(PyObject *self, PyObject *args)
{
double symprec;
PyArrayObject* rotation;
PyArrayObject* translation;
PyArrayObject* py_rotations;
PyArrayObject* py_translations;
int (*rot)[3][3];
double (*trans)[3];
@ -748,15 +749,15 @@ py_get_hall_number_from_symmetry(PyObject *self, PyObject *args)
int hall_number;
if (!PyArg_ParseTuple(args, "OOd",
&rotation,
&translation,
&py_rotations,
&py_translations,
&symprec)) {
return NULL;
}
rot = (int(*)[3][3])PyArray_DATA(rotation);
trans = (double(*)[3])PyArray_DATA(translation);
num_sym = PyArray_DIMS(rotation)[0];
rot = (int(*)[3][3])PyArray_DATA(py_rotations);
trans = (double(*)[3])PyArray_DATA(py_translations);
num_sym = PyArray_DIMS(py_rotations)[0];
hall_number = spg_get_hall_number_from_symmetry(rot,
trans,
@ -768,150 +769,188 @@ py_get_hall_number_from_symmetry(PyObject *self, PyObject *args)
static PyObject * py_get_grid_point_from_address(PyObject *self, PyObject *args)
{
PyArrayObject* grid_address_py;
PyArrayObject* mesh_py;
PyArrayObject* py_grid_address;
PyArrayObject* py_mesh;
int* grid_address;
int* mesh;
int gp;
size_t gp;
if (!PyArg_ParseTuple(args, "OO",
&grid_address_py,
&mesh_py)) {
&py_grid_address,
&py_mesh)) {
return NULL;
}
grid_address = (int*)PyArray_DATA(grid_address_py);
mesh = (int*)PyArray_DATA(mesh_py);
grid_address = (int*)PyArray_DATA(py_grid_address);
mesh = (int*)PyArray_DATA(py_mesh);
gp = spg_get_grid_point_from_address(grid_address, mesh);
gp = spg_get_dense_grid_point_from_address(grid_address, mesh);
return PyLong_FromLong((long) gp);
return PyLong_FromSize_t(gp);
}
static PyObject * py_get_ir_reciprocal_mesh(PyObject *self, PyObject *args)
{
double symprec;
PyArrayObject* grid_address_py;
PyArrayObject* map;
PyArrayObject* mesh;
PyArrayObject* is_shift;
PyArrayObject* py_grid_address;
PyArrayObject* py_grid_mapping_table;
PyArrayObject* py_mesh;
PyArrayObject* py_is_shift;
int is_time_reversal;
PyArrayObject* lattice;
PyArrayObject* position;
PyArrayObject* atom_type;
PyArrayObject* py_lattice;
PyArrayObject* py_positions;
PyArrayObject* py_atom_types;
double (*lat)[3];
double (*pos)[3];
int* types;
int* mesh_int;
int* is_shift_int;
int* mesh;
int* is_shift;
int num_atom;
int (*grid_address)[3];
int *map_int;
int num_ir;
int *grid_mapping_table_int;
size_t *grid_mapping_table_size_t;
int num_ir_int;
size_t num_ir_size_t;
if (!PyArg_ParseTuple(args, "OOOOiOOOd",
&grid_address_py,
&map,
&mesh,
&is_shift,
&py_grid_address,
&py_grid_mapping_table,
&py_mesh,
&py_is_shift,
&is_time_reversal,
&lattice,
&position,
&atom_type,
&py_lattice,
&py_positions,
&py_atom_types,
&symprec)) {
return NULL;
}
lat = (double(*)[3])PyArray_DATA(lattice);
pos = (double(*)[3])PyArray_DATA(position);
types = (int*)PyArray_DATA(atom_type);
mesh_int = (int*)PyArray_DATA(mesh);
is_shift_int = (int*)PyArray_DATA(is_shift);
num_atom = PyArray_DIMS(position)[0];
grid_address = (int(*)[3])PyArray_DATA(grid_address_py);
map_int = (int*)PyArray_DATA(map);
lat = (double(*)[3])PyArray_DATA(py_lattice);
pos = (double(*)[3])PyArray_DATA(py_positions);
types = (int*)PyArray_DATA(py_atom_types);
mesh = (int*)PyArray_DATA(py_mesh);
is_shift = (int*)PyArray_DATA(py_is_shift);
num_atom = PyArray_DIMS(py_positions)[0];
grid_address = (int(*)[3])PyArray_DATA(py_grid_address);
if (PyArray_TYPE(py_grid_mapping_table) == NPY_UINTP) {
grid_mapping_table_size_t = (size_t*)PyArray_DATA(py_grid_mapping_table);
num_ir_size_t = spg_get_dense_ir_reciprocal_mesh(grid_address,
grid_mapping_table_size_t,
mesh,
is_shift,
is_time_reversal,
lat,
pos,
types,
num_atom,
symprec);
return PyLong_FromSize_t(num_ir_size_t);
}
if (PyArray_TYPE(py_grid_mapping_table) == NPY_INT) {
grid_mapping_table_int = (int*)PyArray_DATA(py_grid_mapping_table);
/* num_sym has to be larger than num_sym_from_array_size. */
num_ir_int = spg_get_ir_reciprocal_mesh(grid_address,
grid_mapping_table_int,
mesh,
is_shift,
is_time_reversal,
lat,
pos,
types,
num_atom,
symprec);
return PyLong_FromLong((long) num_ir_int);
}
/* num_sym has to be larger than num_sym_from_array_size. */
num_ir = spg_get_ir_reciprocal_mesh(grid_address,
map_int,
mesh_int,
is_shift_int,
is_time_reversal,
lat,
pos,
types,
num_atom,
symprec);
return PyLong_FromLong((long) num_ir);
Py_RETURN_NONE;
}
static PyObject *
py_get_stabilized_reciprocal_mesh(PyObject *self, PyObject *args)
{
PyArrayObject* grid_address_py;
PyArrayObject* map;
PyArrayObject* mesh;
PyArrayObject* is_shift;
PyArrayObject* py_grid_address;
PyArrayObject* py_grid_mapping_table;
PyArrayObject* py_mesh;
PyArrayObject* py_is_shift;
int is_time_reversal;
PyArrayObject* rotations;
PyArrayObject* qpoints;
PyArrayObject* py_rotations;
PyArrayObject* py_qpoints;
int (*grid_address)[3];
int *map_int;
int* mesh_int;
int* is_shift_int;
int* mesh;
int* is_shift;
int (*rot)[3][3];
int num_rot;
double (*q)[3];
int num_q;
int num_ir;
int *grid_mapping_table_int;
size_t *grid_mapping_table_size_t;
int num_ir_int;
size_t num_ir_size_t;
if (!PyArg_ParseTuple(args, "OOOOiOO",
&grid_address_py,
&map,
&mesh,
&is_shift,
&py_grid_address,
&py_grid_mapping_table,
&py_mesh,
&py_is_shift,
&is_time_reversal,
&rotations,
&qpoints)) {
&py_rotations,
&py_qpoints)) {
return NULL;
}
grid_address = (int(*)[3])PyArray_DATA(grid_address_py);
map_int = (int*)PyArray_DATA(map);
mesh_int = (int*)PyArray_DATA(mesh);
is_shift_int = (int*)PyArray_DATA(is_shift);
rot = (int(*)[3][3])PyArray_DATA(rotations);
num_rot = PyArray_DIMS(rotations)[0];
q = (double(*)[3])PyArray_DATA(qpoints);
num_q = PyArray_DIMS(qpoints)[0];
grid_address = (int(*)[3])PyArray_DATA(py_grid_address);
mesh = (int*)PyArray_DATA(py_mesh);
is_shift = (int*)PyArray_DATA(py_is_shift);
rot = (int(*)[3][3])PyArray_DATA(py_rotations);
num_rot = PyArray_DIMS(py_rotations)[0];
q = (double(*)[3])PyArray_DATA(py_qpoints);
num_q = PyArray_DIMS(py_qpoints)[0];
num_ir = spg_get_stabilized_reciprocal_mesh(grid_address,
map_int,
mesh_int,
is_shift_int,
is_time_reversal,
num_rot,
rot,
num_q,
q);
if (PyArray_TYPE(py_grid_mapping_table) == NPY_UINTP) {
grid_mapping_table_size_t = (size_t*)PyArray_DATA(py_grid_mapping_table);
num_ir_size_t =
spg_get_dense_stabilized_reciprocal_mesh(grid_address,
grid_mapping_table_size_t,
mesh,
is_shift,
is_time_reversal,
num_rot,
rot,
num_q,
q);
return PyLong_FromSize_t(num_ir_size_t);
}
if (PyArray_TYPE(py_grid_mapping_table) == NPY_INT) {
grid_mapping_table_int = (int*)PyArray_DATA(py_grid_mapping_table);
num_ir_int = spg_get_stabilized_reciprocal_mesh(grid_address,
grid_mapping_table_int,
mesh,
is_shift,
is_time_reversal,
num_rot,
rot,
num_q,
q);
return PyLong_FromLong((long) num_ir_int);
}
return PyLong_FromLong((long) num_ir);
Py_RETURN_NONE;
}
static PyObject *
py_get_grid_points_by_rotations(PyObject *self, PyObject *args)
{
PyArrayObject* rot_grid_points_py;
PyArrayObject* address_orig_py;
PyArrayObject* rot_reciprocal_py;
PyArrayObject* mesh_py;
PyArrayObject* is_shift_py;
PyArrayObject* py_rot_grid_points;
PyArrayObject* py_address_orig;
PyArrayObject* py_rot_reciprocal;
PyArrayObject* py_mesh;
PyArrayObject* py_is_shift;
int *rot_grid_points;
size_t *rot_grid_points;
int *address_orig;
int (*rot_reciprocal)[3][3];
int num_rot;
@ -919,133 +958,133 @@ py_get_grid_points_by_rotations(PyObject *self, PyObject *args)
int* is_shift;
if (!PyArg_ParseTuple(args, "OOOOO",
&rot_grid_points_py,
&address_orig_py,
&rot_reciprocal_py,
&mesh_py,
&is_shift_py)) {
&py_rot_grid_points,
&py_address_orig,
&py_rot_reciprocal,
&py_mesh,
&py_is_shift)) {
return NULL;
}
rot_grid_points = (int*)PyArray_DATA(rot_grid_points_py);
address_orig = (int*)PyArray_DATA(address_orig_py);
rot_reciprocal = (int(*)[3][3])PyArray_DATA(rot_reciprocal_py);
num_rot = PyArray_DIMS(rot_reciprocal_py)[0];
mesh = (int*)PyArray_DATA(mesh_py);
is_shift = (int*)PyArray_DATA(is_shift_py);
rot_grid_points = (size_t*)PyArray_DATA(py_rot_grid_points);
address_orig = (int*)PyArray_DATA(py_address_orig);
rot_reciprocal = (int(*)[3][3])PyArray_DATA(py_rot_reciprocal);
num_rot = PyArray_DIMS(py_rot_reciprocal)[0];
mesh = (int*)PyArray_DATA(py_mesh);
is_shift = (int*)PyArray_DATA(py_is_shift);
spg_get_grid_points_by_rotations(rot_grid_points,
address_orig,
num_rot,
rot_reciprocal,
mesh,
is_shift);
spg_get_dense_grid_points_by_rotations(rot_grid_points,
address_orig,
num_rot,
rot_reciprocal,
mesh,
is_shift);
Py_RETURN_NONE;
}
static PyObject *
py_get_BZ_grid_points_by_rotations(PyObject *self, PyObject *args)
{
PyArrayObject* rot_grid_points_py;
PyArrayObject* address_orig_py;
PyArrayObject* rot_reciprocal_py;
PyArrayObject* mesh_py;
PyArrayObject* is_shift_py;
PyArrayObject* bz_map_py;
PyArrayObject* py_rot_grid_points;
PyArrayObject* py_address_orig;
PyArrayObject* py_rot_reciprocal;
PyArrayObject* py_mesh;
PyArrayObject* py_is_shift;
PyArrayObject* py_bz_map;
int *rot_grid_points;
size_t *rot_grid_points;
int *address_orig;
int (*rot_reciprocal)[3][3];
int num_rot;
int* mesh;
int* is_shift;
int* bz_map;
size_t* bz_map;
if (!PyArg_ParseTuple(args, "OOOOOO",
&rot_grid_points_py,
&address_orig_py,
&rot_reciprocal_py,
&mesh_py,
&is_shift_py,
&bz_map_py)) {
&py_rot_grid_points,
&py_address_orig,
&py_rot_reciprocal,
&py_mesh,
&py_is_shift,
&py_bz_map)) {
return NULL;
}
rot_grid_points = (int*)PyArray_DATA(rot_grid_points_py);
address_orig = (int*)PyArray_DATA(address_orig_py);
rot_reciprocal = (int(*)[3][3])PyArray_DATA(rot_reciprocal_py);
num_rot = PyArray_DIMS(rot_reciprocal_py)[0];
mesh = (int*)PyArray_DATA(mesh_py);
is_shift = (int*)PyArray_DATA(is_shift_py);
bz_map = (int*)PyArray_DATA(bz_map_py);
rot_grid_points = (size_t*)PyArray_DATA(py_rot_grid_points);
address_orig = (int*)PyArray_DATA(py_address_orig);
rot_reciprocal = (int(*)[3][3])PyArray_DATA(py_rot_reciprocal);
num_rot = PyArray_DIMS(py_rot_reciprocal)[0];
mesh = (int*)PyArray_DATA(py_mesh);
is_shift = (int*)PyArray_DATA(py_is_shift);
bz_map = (size_t*)PyArray_DATA(py_bz_map);
spg_get_BZ_grid_points_by_rotations(rot_grid_points,
address_orig,
num_rot,
rot_reciprocal,
mesh,
is_shift,
bz_map);
spg_get_dense_BZ_grid_points_by_rotations(rot_grid_points,
address_orig,
num_rot,
rot_reciprocal,
mesh,
is_shift,
bz_map);
Py_RETURN_NONE;
}
static PyObject * py_relocate_BZ_grid_address(PyObject *self, PyObject *args)
{
PyArrayObject* bz_grid_address_py;
PyArrayObject* bz_map_py;
PyArrayObject* grid_address_py;
PyArrayObject* mesh_py;
PyArrayObject* is_shift_py;
PyArrayObject* reciprocal_lattice_py;
PyArrayObject* py_bz_grid_address;
PyArrayObject* py_bz_map;
PyArrayObject* py_grid_address;
PyArrayObject* py_mesh;
PyArrayObject* py_is_shift;
PyArrayObject* py_reciprocal_lattice;
int (*bz_grid_address)[3];
int *bz_map;
size_t *bz_map;
int (*grid_address)[3];
int* mesh;
int* is_shift;
double (*reciprocal_lattice)[3];
int num_ir_gp;
size_t num_ir_gp;
if (!PyArg_ParseTuple(args, "OOOOOO",
&bz_grid_address_py,
&bz_map_py,
&grid_address_py,
&mesh_py,
&reciprocal_lattice_py,
&is_shift_py)) {
&py_bz_grid_address,
&py_bz_map,
&py_grid_address,
&py_mesh,
&py_reciprocal_lattice,
&py_is_shift)) {
return NULL;
}
bz_grid_address = (int(*)[3])PyArray_DATA(bz_grid_address_py);
bz_map = (int*)PyArray_DATA(bz_map_py);
grid_address = (int(*)[3])PyArray_DATA(grid_address_py);
mesh = (int*)PyArray_DATA(mesh_py);
is_shift = (int*)PyArray_DATA(is_shift_py);
reciprocal_lattice = (double(*)[3])PyArray_DATA(reciprocal_lattice_py);
bz_grid_address = (int(*)[3])PyArray_DATA(py_bz_grid_address);
bz_map = (size_t*)PyArray_DATA(py_bz_map);
grid_address = (int(*)[3])PyArray_DATA(py_grid_address);
mesh = (int*)PyArray_DATA(py_mesh);
is_shift = (int*)PyArray_DATA(py_is_shift);
reciprocal_lattice = (double(*)[3])PyArray_DATA(py_reciprocal_lattice);
num_ir_gp = spg_relocate_BZ_grid_address(bz_grid_address,
bz_map,
grid_address,
mesh,
reciprocal_lattice,
is_shift);
num_ir_gp = spg_relocate_dense_BZ_grid_address(bz_grid_address,
bz_map,
grid_address,
mesh,
reciprocal_lattice,
is_shift);
return PyLong_FromLong((long) num_ir_gp);
return PyLong_FromSize_t(num_ir_gp);
}
static PyObject * py_delaunay_reduce(PyObject *self, PyObject *args)
{
PyArrayObject* lattice_py;
PyArrayObject* py_lattice;
double symprec;
double (*lattice)[3];
int result;
if (!PyArg_ParseTuple(args, "Od", &lattice_py, &symprec)) {
if (!PyArg_ParseTuple(args, "Od", &py_lattice, &symprec)) {
return NULL;
}
lattice = (double(*)[3])PyArray_DATA(lattice_py);
lattice = (double(*)[3])PyArray_DATA(py_lattice);
result = spg_delaunay_reduce(lattice, symprec);
@ -1054,17 +1093,17 @@ static PyObject * py_delaunay_reduce(PyObject *self, PyObject *args)
static PyObject * py_niggli_reduce(PyObject *self, PyObject *args)
{
PyArrayObject* lattice_py;
PyArrayObject* py_lattice;
double eps;
double (*lattice)[3];
int result;
if (!PyArg_ParseTuple(args, "Od", &lattice_py, &eps)) {
if (!PyArg_ParseTuple(args, "Od", &py_lattice, &eps)) {
return NULL;
}
lattice = (double(*)[3])PyArray_DATA(lattice_py);
lattice = (double(*)[3])PyArray_DATA(py_lattice);
result = spg_niggli_reduce(lattice, eps);

154
c/kspclib/kgrid.c
View File

@ -1,154 +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 phonopy 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 "kgrid.h"
static void get_all_grid_addresses(int grid_address[][3], const int mesh[3]);
static int get_grid_point_double_mesh(const int address_double[3],
const int mesh[3]);
static int get_grid_point_single_mesh(const int address[3], const int mesh[3]);
static void modulo_i3(int v[3], const int m[3]);
static void reduce_grid_address(int address[3], const int mesh[3]);
static void reduce_grid_address_double(int address[3], const int mesh[3]);
void kgd_get_all_grid_addresses(int grid_address[][3], const int mesh[3])
{
get_all_grid_addresses(grid_address, mesh);
}
int kgd_get_grid_point_double_mesh(const int address_double[3],
const int mesh[3])
{
return get_grid_point_double_mesh(address_double, mesh);
}
void kgd_get_grid_address_double_mesh(int address_double[3],
const int address[3],
const int mesh[3],
const int is_shift[3])
{
int i;
for (i = 0; i < 3; i++) {
address_double[i] = address[i] * 2 + (is_shift[i] != 0);
}
reduce_grid_address_double(address_double, mesh);
}
static void get_all_grid_addresses(int grid_address[][3], const int mesh[3])
{
int i, j, k, grid_point;
int 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_point = get_grid_point_single_mesh(address, mesh);
grid_address[grid_point][0] = address[0];
grid_address[grid_point][1] = address[1];
grid_address[grid_point][2] = address[2];
reduce_grid_address(grid_address[grid_point], mesh);
}
}
}
}
static int get_grid_point_double_mesh(const int address_double[3],
const int mesh[3])
{
int i, 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;
}
}
modulo_i3(address, mesh);
return get_grid_point_single_mesh(address, mesh);
}
static int get_grid_point_single_mesh(const int address[3],
const int mesh[3])
{
#ifndef GRID_ORDER_XYZ
return address[2] * mesh[0] * mesh[1] + address[1] * mesh[0] + address[0];
#else
return address[0] * mesh[1] * mesh[2] + address[1] * mesh[2] + address[2];
#endif
}
static void modulo_i3(int v[3], const int m[3])
{
int i;
for (i = 0; i < 3; i++) {
v[i] = v[i] % m[i];
if (v[i] < 0) {
v[i] += m[i];
}
}
}
static void reduce_grid_address(int address[3], const int mesh[3])
{
int 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_grid_address_double(int address[3], const int mesh[3])
{
int 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
}
}

656
c/kspclib/tetrahedron_method.c
View File

File diff suppressed because it is too large Load Diff

81
c/kspclib_h/kgrid.h
View File

@ -1,81 +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 phonopy 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 __kgrid_H__
#define __kgrid_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 kgd_get_all_grid_addresses(int grid_address[][3], const int mesh[3]);
int kgd_get_grid_point_double_mesh(const int address_double[3],
const int mesh[3]);
void kgd_get_grid_address_double_mesh(int address_double[3],
const int address[3],
const int mesh[3],
const int is_shift[3]);
#endif

36
c/kspclib_h/tetrahedron_method.h
View File

@ -39,24 +39,34 @@
#define THMCONST
#endif
#include <stddef.h>
void thm_get_relative_grid_address(int relative_grid_address[24][4][3],
THMCONST double rec_lattice[3][3]);
THMCONST double rec_lattice[3][3]);
void thm_get_all_relative_grid_address(int relative_grid_address[4][24][4][3]);
double thm_get_integration_weight(const double omega,
THMCONST double tetrahedra_omegas[24][4],
const char function);
THMCONST double tetrahedra_omegas[24][4],
const char function);
void
thm_get_integration_weight_at_omegas(double *integration_weights,
const int num_omegas,
const double *omegas,
THMCONST double tetrahedra_omegas[24][4],
const char function);
const int num_omegas,
const double *omegas,
THMCONST double tetrahedra_omegas[24][4],
const char function);
void thm_get_neighboring_grid_points(int neighboring_grid_points[],
const int grid_point,
THMCONST int relative_grid_address[][3],
const int num_relative_grid_address,
const int mesh[3],
THMCONST int bz_grid_address[][3],
const int bz_map[]);
const int grid_point,
THMCONST int relative_grid_address[][3],
const int num_relative_grid_address,
const int mesh[3],
THMCONST int bz_grid_address[][3],
const int bz_map[]);
void
thm_get_dense_neighboring_grid_points(size_t neighboring_grid_points[],
const size_t grid_point,
THMCONST int relative_grid_address[][3],
const int num_relative_grid_address,
const int mesh[3],
THMCONST int bz_grid_address[][3],
const size_t bz_map[]);
#endif

6
c/spglib/kgrid.c
View File

@ -33,6 +33,7 @@
/* POSSIBILITY OF SUCH DAMAGE. */
#include <stddef.h>
#include <assert.h>
#include "kgrid.h"
static void get_all_grid_addresses(int grid_address[][3], const int mesh[3]);
@ -77,7 +78,7 @@ void kgd_get_grid_address_double_mesh(int address_double[3],
static void get_all_grid_addresses(int grid_address[][3], const int mesh[3])
{
int i, j, k;
long grid_point;
size_t grid_point;
int address[3];
for (i = 0; i < mesh[0]; i++) {
@ -87,6 +88,9 @@ static void get_all_grid_addresses(int grid_address[][3], const int mesh[3])
for (k = 0; k < mesh[2]; k++) {
address[2] = k;
grid_point = get_grid_point_single_mesh(address, mesh);
assert(mesh[0] * mesh[1] * mesh[2] > grid_point);
grid_address[grid_point][0] = address[0];
grid_address[grid_point][1] = address[1];
grid_address[grid_point][2] = address[2];

117
c/spglib/kpoint.c
View File

@ -345,43 +345,13 @@ size_t kpt_get_dense_stabilized_reciprocal_mesh(int grid_address[][3],
return num_ir;
}
void kpt_get_grid_points_by_rotations(int rot_grid_points[],
const int address_orig[3],
const MatINT * rot_reciprocal,
const int mesh[3],
const int is_shift[3])
{
int i;
size_t *dense_rot_grid_points;
if ((dense_rot_grid_points =
(size_t*)malloc(sizeof(size_t) * rot_reciprocal->size)) == NULL) {
warning_print("spglib: Memory of unique_rot could not be allocated.");
goto err;
}
kpt_get_dense_grid_points_by_rotations(dense_rot_grid_points,
address_orig,
rot_reciprocal,
mesh,
is_shift);
for (i = 0; i < rot_reciprocal->size; i++) {
rot_grid_points[i] = dense_rot_grid_points[i];
}
free(dense_rot_grid_points);
dense_rot_grid_points = NULL;
err:
;
}
void kpt_get_dense_grid_points_by_rotations(size_t rot_grid_points[],
const int address_orig[3],
const MatINT * rot_reciprocal,
const int mesh[3],
const int is_shift[3])
void
kpt_get_dense_grid_points_by_rotations(size_t rot_grid_points[],
const int address_orig[3],
SPGCONST int (*rot_reciprocal)[3][3],
const int num_rot,
const int mesh[3],
const int is_shift[3])
{
int i;
int address_double_orig[3], address_double[3];
@ -389,71 +359,22 @@ void kpt_get_dense_grid_points_by_rotations(size_t rot_grid_points[],
for (i = 0; i < 3; i++) {
address_double_orig[i] = address_orig[i] * 2 + is_shift[i];
}
for (i = 0; i < rot_reciprocal->size; i++) {
for (i = 0; i < num_rot; i++) {
mat_multiply_matrix_vector_i3(address_double,
rot_reciprocal->mat[i],
rot_reciprocal[i],
address_double_orig);
rot_grid_points[i] = kgd_get_dense_grid_point_double_mesh(address_double, mesh);
}
}
void kpt_get_BZ_grid_points_by_rotations(int rot_grid_points[],
const int address_orig[3],
const MatINT * rot_reciprocal,
const int mesh[3],
const int is_shift[3],
const int bz_map[])
{
int i, num_bz_map;
size_t *dense_rot_grid_points, *dense_bz_map;
if ((dense_rot_grid_points =
(size_t*)malloc(sizeof(size_t) * rot_reciprocal->size)) == NULL) {
warning_print("spglib: Memory of unique_rot could not be allocated.");
goto err;
}
num_bz_map = mesh[0] * mesh[1] * mesh[2] * 8;
if ((dense_bz_map =
(size_t*)malloc(sizeof(size_t) * num_bz_map)) == NULL) {
warning_print("spglib: Memory of unique_rot could not be allocated.");
free(dense_rot_grid_points);
dense_rot_grid_points = NULL;
goto err;
}
for (i = 0; i < num_bz_map; i++) {
dense_bz_map[i] = bz_map[i];
}
kpt_get_dense_BZ_grid_points_by_rotations(dense_rot_grid_points,
address_orig,
rot_reciprocal,
mesh,
is_shift,
dense_bz_map);
free(dense_bz_map);
dense_bz_map = NULL;
for (i = 0; i < rot_reciprocal->size; i++) {
rot_grid_points[i] = dense_rot_grid_points[i];
}
free(dense_rot_grid_points);
dense_rot_grid_points = NULL;
err:
;
}
void kpt_get_dense_BZ_grid_points_by_rotations(size_t rot_grid_points[],
const int address_orig[3],
const MatINT * rot_reciprocal,
const int mesh[3],
const int is_shift[3],
const size_t bz_map[])
void
kpt_get_dense_BZ_grid_points_by_rotations(size_t rot_grid_points[],
const int address_orig[3],
SPGCONST int (*rot_reciprocal)[3][3],
const int num_rot,
const int mesh[3],
const int is_shift[3],
const size_t bz_map[])
{
int i;
int address_double_orig[3], address_double[3], bzmesh[3];
@ -462,9 +383,9 @@ void kpt_get_dense_BZ_grid_points_by_rotations(size_t rot_grid_points[],
bzmesh[i] = mesh[i] * 2;
address_double_orig[i] = address_orig[i] * 2 + is_shift[i];
}
for (i = 0; i < rot_reciprocal->size; i++) {
for (i = 0; i < num_rot; i++) {
mat_multiply_matrix_vector_i3(address_double,
rot_reciprocal->mat[i],
rot_reciprocal[i],
address_double_orig);
rot_grid_points[i] =
bz_map[kgd_get_dense_grid_point_double_mesh(address_double, bzmesh)];

4
c/spglib/refinement.c
View File

@ -301,6 +301,10 @@ void ref_free_exact_structure(ExactStructure *exstr)
free(exstr->std_mapping_to_primitive);
exstr->std_mapping_to_primitive = NULL;
}
if (exstr->site_symmetry_symbols != NULL) {
free(exstr->site_symmetry_symbols);
exstr->site_symmetry_symbols = NULL;
}
free(exstr);
}
}

3
c/spglib/spacegroup.c
View File

@ -439,6 +439,9 @@ int spa_search_spacegroup_with_symmetry(const Symmetry *symmetry,
symmetry,
symprec,
-1.0);
cel_free_cell(primitive);
primitive = NULL;
if (spacegroup != NULL) {
hall_number = spacegroup->hall_number;
free(spacegroup);

166
c/spglib/spglib.c
View File

@ -357,6 +357,11 @@ void spg_free_dataset(SpglibDataset *dataset)
dataset->n_std_atoms = 0;
}
if (dataset->site_symmetry_symbols != NULL) {
free(dataset->site_symmetry_symbols);
dataset->site_symmetry_symbols = NULL;
}
dataset->spacegroup_number = 0;
dataset->hall_number = 0;
strcpy(dataset->international_symbol, "");
@ -471,23 +476,37 @@ int spg_get_hall_number_from_symmetry(SPGCONST int rotation[][3][3],
symmetry = NULL;
prim_symmetry = NULL;
hall_number = 0;
if ((symmetry = sym_alloc_symmetry(num_operations)) == NULL) {
return 0;
}
symmetry = sym_alloc_symmetry(num_operations);
for (i = 0; i < num_operations; i++) {
mat_copy_matrix_i3(symmetry->rot[i], rotation[i]);
mat_copy_vector_d3(symmetry->trans[i], translation[i]);
}
prim_symmetry = prm_get_primitive_symmetry(symmetry, symprec);
sym_free_symmetry(symmetry);
symmetry = NULL;
if (prim_symmetry == NULL) {
return 0;
}
hall_number = spa_search_spacegroup_with_symmetry(prim_symmetry, symprec);
if (hall_number) {
spglib_error_code = SPGLIB_SUCCESS;
return hall_number;
} else {
spglib_error_code = SPGERR_SPACEGROUP_SEARCH_FAILED;
return 0;
}
sym_free_symmetry(prim_symmetry);
prim_symmetry = NULL;
return hall_number;
}
/* Return 0 if failed */
@ -968,143 +987,36 @@ size_t spg_get_dense_stabilized_reciprocal_mesh(int grid_address[][3],
qpoints);
}
int spg_get_grid_points_by_rotations(int rot_grid_points[],
const int address_orig[3],
const int num_rot,
SPGCONST int rot_reciprocal[][3][3],
const int mesh[3],
const int is_shift[3])
void spg_get_dense_grid_points_by_rotations(size_t rot_grid_points[],
const int address_orig[3],
const int num_rot,
SPGCONST int rot_reciprocal[][3][3],
const int mesh[3],
const int is_shift[3])
{
int i;
MatINT *rot;
rot = NULL;
if ((rot = mat_alloc_MatINT(num_rot)) == NULL) {
return 0;
}
for (i = 0; i < num_rot; i++) {
mat_copy_matrix_i3(rot->mat[i], rot_reciprocal[i]);
}
kpt_get_grid_points_by_rotations(rot_grid_points,
address_orig,
rot,
mesh,
is_shift);
mat_free_MatINT(rot);
rot = NULL;
return 1;
}
int spg_get_dense_grid_points_by_rotations(size_t rot_grid_points[],
const int address_orig[3],
const int num_rot,
SPGCONST int rot_reciprocal[][3][3],
const int mesh[3],
const int is_shift[3])
{
int i;
MatINT *rot;
rot = NULL;
if ((rot = mat_alloc_MatINT(num_rot)) == NULL) {
return 0;
}
for (i = 0; i < num_rot; i++) {
mat_copy_matrix_i3(rot->mat[i], rot_reciprocal[i]);
}
kpt_get_dense_grid_points_by_rotations(rot_grid_points,
address_orig,
rot,
rot_reciprocal,
num_rot,
mesh,
is_shift);
mat_free_MatINT(rot);
rot = NULL;
return 1;
}
int spg_get_BZ_grid_points_by_rotations(int rot_grid_points[],
const int address_orig[3],
const int num_rot,
SPGCONST int rot_reciprocal[][3][3],
const int mesh[3],
const int is_shift[3],
const int bz_map[])
void spg_get_dense_BZ_grid_points_by_rotations(size_t rot_grid_points[],
const int address_orig[3],
const int num_rot,
SPGCONST int rot_reciprocal[][3][3],
const int mesh[3],
const int is_shift[3],
const size_t bz_map[])
{
int i;
MatINT *rot;
rot = NULL;
if ((rot = mat_alloc_MatINT(num_rot)) == NULL) {
return 0;
}
for (i = 0; i < num_rot; i++) {
mat_copy_matrix_i3(rot->mat[i], rot_reciprocal[i]);
}
kpt_get_BZ_grid_points_by_rotations(rot_grid_points,
address_orig,
rot,
mesh,
is_shift,
bz_map);
mat_free_MatINT(rot);
rot = NULL;
return 1;
}
int spg_get_dense_BZ_grid_points_by_rotations(size_t rot_grid_points[],
const int address_orig[3],
const int num_rot,
SPGCONST int rot_reciprocal[][3][3],
const int mesh[3],
const int is_shift[3],
const size_t bz_map[])
{
int i;
MatINT *rot;
rot = NULL;
if ((rot = mat_alloc_MatINT(num_rot)) == NULL) {
return 0;
}
for (i = 0; i < num_rot; i++) {
mat_copy_matrix_i3(rot->mat[i], rot_reciprocal[i]);
}
kpt_get_dense_BZ_grid_points_by_rotations(rot_grid_points,
address_orig,
rot,
rot_reciprocal,
num_rot,
mesh,
is_shift,
bz_map);
mat_free_MatINT(rot);
rot = NULL;
return 1;
}
int spg_relocate_BZ_grid_address(int bz_grid_address[][3],
int bz_map[],
SPGCONST int grid_address[][3],
const int mesh[3],
SPGCONST double rec_lattice[3][3],
const int is_shift[3])
{
return kpt_relocate_BZ_grid_address(bz_grid_address,
bz_map,
grid_address,
mesh,
rec_lattice,
is_shift);
}
size_t spg_relocate_dense_BZ_grid_address(int bz_grid_address[][3],

37
c/spglib_h/kpoint.h
View File

@ -64,28 +64,21 @@ size_t kpt_get_dense_stabilized_reciprocal_mesh(int grid_address[][3],
const MatINT * rotations,
const size_t num_q,
SPGCONST double qpoints[][3]);
void kpt_get_grid_points_by_rotations(int rot_grid_points[],
const int address_orig[3],
const MatINT * rot_reciprocal,
const int mesh[3],
const int is_shift[3]);
void kpt_get_dense_grid_points_by_rotations(size_t rot_grid_points[],
const int address_orig[3],
const MatINT * rot_reciprocal,
const int mesh[3],
const int is_shift[3]);
void kpt_get_BZ_grid_points_by_rotations(int rot_grid_points[],
const int address_orig[3],
const MatINT * rot_reciprocal,
const int mesh[3],
const int is_shift[3],
const int bz_map[]);
void kpt_get_dense_BZ_grid_points_by_rotations(size_t rot_grid_points[],
const int address_orig[3],
const MatINT * rot_reciprocal,
const int mesh[3],
const int is_shift[3],
const size_t bz_map[]);
void
kpt_get_dense_grid_points_by_rotations(size_t rot_grid_points[],
const int address_orig[3],
SPGCONST int (*rot_reciprocal)[3][3],
const int num_rot,
const int mesh[3],
const int is_shift[3]);
void
kpt_get_dense_BZ_grid_points_by_rotations(size_t rot_grid_points[],
const int address_orig[3],
SPGCONST int (*rot_reciprocal)[3][3],
const int num_rot,
const int mesh[3],
const int is_shift[3],
const size_t bz_map[]);
int kpt_relocate_BZ_grid_address(int bz_grid_address[][3],
int bz_map[],
SPGCONST int grid_address[][3],

66
c/spglib_h/spglib.h
View File

@ -451,33 +451,33 @@ extern "C" {
/* Rotation operations in reciprocal space ``rot_reciprocal`` are applied */
/* to a grid address ``address_orig`` and resulting grid points are stored in */
/* ``rot_grid_points``. Return 0 if failed. */
int spg_get_grid_points_by_rotations(int rot_grid_points[],
const int address_orig[3],
const int num_rot,
SPGCONST int rot_reciprocal[][3][3],
const int mesh[3],
const int is_shift[3]);
int spg_get_dense_grid_points_by_rotations(size_t rot_grid_points[],
const int address_orig[3],
const int num_rot,
SPGCONST int rot_reciprocal[][3][3],
const int mesh[3],
const int is_shift[3]);
void spg_get_grid_points_by_rotations(int rot_grid_points[],
const int address_orig[3],
const int num_rot,
SPGCONST int rot_reciprocal[][3][3],
const int mesh[3],
const int is_shift[3]);
void spg_get_dense_grid_points_by_rotations(size_t rot_grid_points[],
const int address_orig[3],
const int num_rot,
SPGCONST int rot_reciprocal[][3][3],
const int mesh[3],
const int is_shift[3]);
int spg_get_BZ_grid_points_by_rotations(int rot_grid_points[],
const int address_orig[3],
const int num_rot,
SPGCONST int rot_reciprocal[][3][3],
const int mesh[3],
const int is_shift[3],
const int bz_map[]);
int spg_get_dense_BZ_grid_points_by_rotations(size_t rot_grid_points[],
const int address_orig[3],
const int num_rot,
SPGCONST int rot_reciprocal[][3][3],
const int mesh[3],
const int is_shift[3],
const size_t bz_map[]);
void spg_get_BZ_grid_points_by_rotations(int rot_grid_points[],
const int address_orig[3],
const int num_rot,
SPGCONST int rot_reciprocal[][3][3],
const int mesh[3],
const int is_shift[3],
const int bz_map[]);
void spg_get_dense_BZ_grid_points_by_rotations(size_t rot_grid_points[],
const int address_orig[3],
const int num_rot,
SPGCONST int rot_reciprocal[][3][3],
const int mesh[3],
const int is_shift[3],
const size_t bz_map[]);
/* Grid addresses are relocated inside Brillouin zone. */
/* Number of ir-grid-points inside Brillouin zone is returned. */
@ -501,12 +501,6 @@ extern "C" {
/* bz_map is used to recover grid point index expanded to include BZ */
/* surface from grid address. The grid point indices are mapped to */
/* (mesh[0] * 2) x (mesh[1] * 2) x (mesh[2] * 2) space (bz_map). */
int spg_relocate_BZ_grid_address(int bz_grid_address[][3],
int bz_map[],
SPGCONST int grid_address[][3],
const int mesh[3],
SPGCONST double rec_lattice[3][3],
const int is_shift[3]);
size_t spg_relocate_dense_BZ_grid_address(int bz_grid_address[][3],
size_t bz_map[],
SPGCONST int grid_address[][3],
@ -514,14 +508,6 @@ extern "C" {
SPGCONST double rec_lattice[3][3],
const int is_shift[3]);
void spg_get_neighboring_grid_points(int relative_grid_points[],
const int grid_point,
SPGCONST int relative_grid_address[][3],
const int num_relative_grid_address,
const int mesh[3],
SPGCONST int bz_grid_address[][3],
const int bz_map[]);
/*--------*/
/* Niggli */
/*--------*/

4
phonopy/api_phonopy.py
View File

@ -1707,8 +1707,8 @@ class Phonopy(object):
"Use Phonopy.plot_projected_dos (lowercase on DOS).",
DeprecationWarning)
return self.plot_partial_dos(pdos_indices=pdos_indices,
legend=legend)
return self.plot_projected_dos(pdos_indices=pdos_indices,
legend=legend)
def plot_projected_dos(self, pdos_indices=None, legend=None):
"""Plot projected DOS

12
phonopy/structure/grid_points.py
View File

@ -60,7 +60,8 @@ def get_qpoints(mesh_numbers,
def extract_ir_grid_points(grid_mapping_table):
ir_grid_points = np.array(np.unique(grid_mapping_table), dtype='intc')
ir_grid_points = np.array(np.unique(grid_mapping_table),
dtype=grid_mapping_table.dtype)
weights = np.zeros_like(grid_mapping_table)
for i, gp in enumerate(grid_mapping_table):
weights[gp] += 1
@ -269,16 +270,19 @@ class GridPoints(object):
self._mesh,
rotations,
is_shift=self._is_shift,
is_time_reversal=is_time_reversal)
is_time_reversal=is_time_reversal,
is_dense=True)
shift = np.array(self._is_shift, dtype='intc') * 0.5
if self._fit_in_BZ:
self._grid_address = relocate_BZ_grid_address(
grid_address, _ = relocate_BZ_grid_address(
grid_address,
self._mesh,
self._rec_lat,
is_shift=self._is_shift)[0][:np.prod(self._mesh)]
is_shift=self._is_shift,
is_dense=True)
self._grid_address = grid_address[:np.prod(self._mesh)]
else:
self._grid_address = grid_address

271
phonopy/structure/spglib.py
View File

@ -495,27 +495,40 @@ def get_ir_reciprocal_mesh(mesh,
cell,
is_shift=None,
is_time_reversal=True,
symprec=1e-5):
symprec=1e-5,
is_dense=False):
"""Return k-points mesh and k-point map to the irreducible k-points.
The symmetry is serched from the input cell.
Args:
mesh:
int array (3,): Uniform sampling mesh numbers
cell, symprec:
See the docstring of get_symmetry.
is_shift:
int array (3,): [0, 0, 0] gives Gamma center mesh and value 1 gives
half mesh shift.
is_time_reversal:
bool: Time reversal symmetry is included or not.
Parameters
----------
mesh : array_like
Uniform sampling mesh numbers.
dtype='intc', shape=(3,)
cell : spglib cell tuple
Crystal structure.
is_shift : array_like, optional
[0, 0, 0] gives Gamma center mesh and value 1 gives half mesh shift.
Default is None which equals to [0, 0, 0].
dtype='intc', shape=(3,)
is_time_reversal : bool, optional
Whether time reversal symmetry is included or not. Default is True.
symprec : float, optional
Symmetry tolerance in distance. Default is 1e-5.
is_dense : bool, optional
grid_mapping_table is returned with dtype='uintp' if True. Otherwise
its dtype='intc'. Default is False.
Returns
-------
grid_mapping_table : ndarray
Grid point mapping table to ir-gird-points.
dtype='intc' or 'uintp', shape=(prod(mesh),)
grid_address : ndarray
Address of all grid points.
dtype='intc', shspe=(prod(mesh), 3)
Returns:
mapping_table:
int array (N,): Grid point mapping table to ir-gird-points
grid_address:
int array (N, 3): Address of all grid points
"""
_set_no_error()
@ -523,13 +536,17 @@ def get_ir_reciprocal_mesh(mesh,
if lattice is None:
return None
mapping = np.zeros(np.prod(mesh), dtype='intc')
if is_dense:
dtype = 'uintp'
else:
dtype = 'intc'
grid_mapping_table = np.zeros(np.prod(mesh), dtype=dtype)
grid_address = np.zeros((np.prod(mesh), 3), dtype='intc')
if is_shift is None:
is_shift = [0, 0, 0]
if spg.ir_reciprocal_mesh(
grid_address,
mapping,
grid_mapping_table,
np.array(mesh, dtype='intc'),
np.array(is_shift, dtype='intc'),
is_time_reversal * 1,
@ -537,7 +554,7 @@ def get_ir_reciprocal_mesh(mesh,
positions,
numbers,
symprec) > 0:
return mapping, grid_address
return grid_mapping_table, grid_address
else:
return None
@ -546,32 +563,51 @@ def get_stabilized_reciprocal_mesh(mesh,
rotations,
is_shift=None,
is_time_reversal=True,
qpoints=None):
"""Return k-point map to the irreducible k-points and k-point grid points .
qpoints=None,
is_dense=False):
"""Return k-point map to the irreducible k-points and k-point grid points.
The symmetry is searched from the input rotation matrices in real space.
Args:
mesh:
int array (3,): Uniform sampling mesh numbers
is_shift:
int array (3,): [0, 0, 0] gives Gamma center mesh and value 1 gives
half mesh shift.
is_time_reversal:
bool: Time reversal symmetry is included or not.
qpoints:
float array (N ,3) or (3,):
Stabilizer(s) in the fractional coordinates.
Parameters
----------
mesh : array_like
Uniform sampling mesh numbers.
dtype='intc', shape=(3,)
rotations : array_like
Rotation matrices with respect to real space basis vectors.
dtype='intc', shape=(rotations, 3)
is_shift : array_like
[0, 0, 0] gives Gamma center mesh and value 1 gives half mesh shift.
dtype='intc', shape=(3,)
is_time_reversal : bool
Time reversal symmetry is included or not.
qpoints : array_like
q-points used as stabilizer(s) given in reciprocal space with respect
to reciprocal basis vectors.
dtype='double', shape=(qpoints ,3) or (3,)
is_dense : bool, optional
grid_mapping_table is returned with dtype='uintp' if True. Otherwise
its dtype='intc'. Default is False.
Returns
-------
grid_mapping_table : ndarray
Grid point mapping table to ir-gird-points.
dtype='intc', shape=(prod(mesh),)
grid_address : ndarray
Address of all grid points. Each address is given by three unsigned
integers.
dtype='intc', shape=(prod(mesh), 3)
Returns:
mapping_table:
int array (N,): Grid point mapping table to ir-gird-points
grid_address:
int array (N, 3): Address of all grid points
"""
_set_no_error()
mapping_table = np.zeros(np.prod(mesh), dtype='intc')
if is_dense:
dtype = 'uintp'
else:
dtype = 'intc'
mapping_table = np.zeros(np.prod(mesh), dtype=dtype)
grid_address = np.zeros((np.prod(mesh), 3), dtype='intc')
if is_shift is None:
is_shift = [0, 0, 0]
@ -598,87 +634,176 @@ def get_stabilized_reciprocal_mesh(mesh,
def get_grid_points_by_rotations(address_orig,
reciprocal_rotations,
mesh,
is_shift=np.zeros(3, dtype='intc')):
"""Rotation operations in reciprocal space ``reciprocal_rotations`` are applied
to a grid point ``grid_point`` and resulting grid points are returned.
is_shift=None,
is_dense=False):
"""Returns grid points obtained after rotating input grid address
Parameters
----------
address_orig : array_like
Grid point address to be rotated.
dtype='intc', shape=(3,)
reciprocal_rotations : array_like
Rotation matrices {R} with respect to reciprocal basis vectors.
Defined by q'=Rq.
dtype='intc', shape=(rotations, 3, 3)
mesh : array_like
dtype='intc', shape=(3,)
is_shift : array_like, optional
With (1) or without (0) half grid shifts with respect to grid intervals
sampled along reciprocal basis vectors. Default is None, which
gives [0, 0, 0].
is_dense : bool, optional
rot_grid_points is returned with dtype='uintp' if True. Otherwise
its dtype='intc'. Default is False.
Returns
-------
rot_grid_points : ndarray
Grid points obtained after rotating input grid address
dtype='intc' or 'uintp', shape=(rotations,)
"""
_set_no_error()
rot_grid_points = np.zeros(len(reciprocal_rotations), dtype='intc')
if is_shift is None:
_is_shift = np.zeros(3, dtype='intc')
else:
_is_shift = np.array(is_shift, dtype='intc')
rot_grid_points = np.zeros(len(reciprocal_rotations), dtype='uintp')
spg.grid_points_by_rotations(
rot_grid_points,
np.array(address_orig, dtype='intc'),
np.array(reciprocal_rotations, dtype='intc', order='C'),
np.array(mesh, dtype='intc'),
np.array(is_shift, dtype='intc'))
_is_shift)
return rot_grid_points
if is_dense:
return rot_grid_points
else:
return np.array(rot_grid_points, dtype='intc')
def get_BZ_grid_points_by_rotations(address_orig,
reciprocal_rotations,
mesh,
bz_map,
is_shift=np.zeros(3, dtype='intc')):
"""Rotation operations in reciprocal space ``reciprocal_rotations`` are applied
to a grid point ``grid_point`` and resulting grid points are returned.
is_shift=None,
is_dense=False):
"""Returns grid points obtained after rotating input grid address
Parameters
----------
address_orig : array_like
Grid point address to be rotated.
dtype='intc', shape=(3,)
reciprocal_rotations : array_like
Rotation matrices {R} with respect to reciprocal basis vectors.
Defined by q'=Rq.
dtype='intc', shape=(rotations, 3, 3)
mesh : array_like
dtype='intc', shape=(3,)
is_shift : array_like, optional
With (1) or without (0) half grid shifts with respect to grid intervals
sampled along reciprocal basis vectors. Default is None, which
gives [0, 0, 0].
is_dense : bool, optional
rot_grid_points is returned with dtype='uintp' if True. Otherwise
its dtype='intc'. Default is False.
Returns
-------
rot_grid_points : ndarray
Grid points obtained after rotating input grid address
dtype='intc' or 'uintp', shape=(rotations,)
"""
_set_no_error()
rot_grid_points = np.zeros(len(reciprocal_rotations), dtype='intc')
if is_shift is None:
_is_shift = np.zeros(3, dtype='intc')
else:
_is_shift = np.array(is_shift, dtype='intc')
if bz_map.dtype == 'uintp' and bz_map.flags.c_contiguous:
_bz_map = bz_map
else:
_bz_map = np.array(bz_map, dtype='uintp')
rot_grid_points = np.zeros(len(reciprocal_rotations), dtype='uintp')
spg.BZ_grid_points_by_rotations(
rot_grid_points,
np.array(address_orig, dtype='intc'),
np.array(reciprocal_rotations, dtype='intc', order='C'),
np.array(mesh, dtype='intc'),
np.array(is_shift, dtype='intc'),
bz_map)
_is_shift,
_bz_map)
return rot_grid_points
if is_dense:
return rot_grid_points
else:
return np.array(rot_grid_points, dtype='intc')
def relocate_BZ_grid_address(grid_address,
mesh,
reciprocal_lattice, # column vectors
is_shift=np.zeros(3, dtype='intc')):
"""Grid addresses are relocated inside Brillouin zone.
is_shift=None,
is_dense=False):
"""Grid addresses are relocated to be inside first Brillouin zone.
Number of ir-grid-points inside Brillouin zone is returned.
It is assumed that the following arrays have the shapes of
bz_grid_address[prod(mesh + 1)][3]
bz_map[prod(mesh * 2)]
where grid_address[prod(mesh)][3].
Each element of grid_address is mapped to each element of
bz_grid_address with keeping element order. bz_grid_address has
larger memory space to represent BZ surface even if some points
on a surface are translationally equivalent to the other points
on the other surface. Those equivalent points are added successively
as grid point numbers to bz_grid_address. Those added grid points
are stored after the address of end point of grid_address, i.e.
bz_grid_address : (num_grid_points_in_FBZ, 3)
bz_map (prod(mesh * 2), )
|-----------------array size of bz_grid_address---------------------|
|--grid addresses similar to grid_address--|--newly added ones--|xxx|
Note that the shape of grid_address is (prod(mesh), 3) and the
addresses in grid_address are arranged to be in parallelepiped
made of reciprocal basis vectors. The addresses in bz_grid_address
are inside the first Brillouin zone or on its surface. Each
address in grid_address is mapped to one of those in
bz_grid_address by a reciprocal lattice vector (including zero
vector) with keeping element order. For those inside first
Brillouin zone, the mapping is one-to-one. For those on the first
Brillouin zone surface, more than one addresses in bz_grid_address
that are equivalent by the reciprocal lattice translations are
mapped to one address in grid_address. In this case, those grid
points except for one of them are appended to the tail of this array,
for which bz_grid_address has the following data storing:
where xxx means the memory space that may not be used. Number of grid
points stored in bz_grid_address is returned.
|------------------array size of bz_grid_address-------------------------|
|--those equivalent to grid_address--|--those on surface except for one--|
|-----array size of grid_address-----|
Number of grid points stored in bz_grid_address is returned.
bz_map is used to recover grid point index expanded to include BZ
surface from grid address. The grid point indices are mapped to
(mesh[0] * 2) x (mesh[1] * 2) x (mesh[2] * 2) space (bz_map).
"""
_set_no_error()
bz_grid_address = np.zeros(
((mesh[0] + 1) * (mesh[1] + 1) * (mesh[2] + 1), 3), dtype='intc')
bz_map = np.zeros(
(2 * mesh[0]) * (2 * mesh[1]) * (2 * mesh[2]), dtype='intc')
if is_shift is None:
_is_shift = np.zeros(3, dtype='intc')
else:
_is_shift = np.array(is_shift, dtype='intc')
bz_grid_address = np.zeros((np.prod(np.add(mesh, 1)), 3), dtype='intc')
bz_map = np.zeros(np.prod(np.multiply(mesh, 2)), dtype='uintp')
num_bz_ir = spg.BZ_grid_address(
bz_grid_address,
bz_map,
grid_address,
np.array(mesh, dtype='intc'),
np.array(reciprocal_lattice, dtype='double', order='C'),
np.array(is_shift, dtype='intc'))
_is_shift)
return bz_grid_address[:num_bz_ir], bz_map
if is_dense:
return bz_grid_address[:num_bz_ir], bz_map
else:
return bz_grid_address[:num_bz_ir], np.array(bz_map, dtype='intc')
def delaunay_reduce(lattice, eps=1e-5):

18
phonopy/structure/tetrahedron_method.py
View File

@ -49,20 +49,23 @@ parallelepiped_vertices = np.array([[0, 0, 0],
[0, 1, 1],
[1, 1, 1]], dtype='intc', order='C')
def get_neighboring_grid_points(grid_point,
relative_grid_address,
mesh,
bz_grid_address,
bz_map):
relative_grid_points = np.zeros(len(relative_grid_address), dtype='intc')
relative_grid_points = np.zeros(len(relative_grid_address),
dtype=bz_map.dtype)
phonoc.neighboring_grid_points(relative_grid_points,
grid_point,
relative_grid_address,
mesh,
bz_grid_address,
bz_map)
grid_point,
relative_grid_address,
mesh,
bz_grid_address,
bz_map)
return relative_grid_points
def get_tetrahedra_relative_grid_address(microzone_lattice):
"""Returns relative (differences of) grid addresses from the central
@ -81,12 +84,14 @@ def get_tetrahedra_relative_grid_address(microzone_lattice):
return relative_grid_address
def get_all_tetrahedra_relative_grid_address():
relative_grid_address = np.zeros((4, 24, 4, 3), dtype='intc')
phonoc.all_tetrahedra_relative_grid_address(relative_grid_address)
return relative_grid_address
def get_tetrahedra_integration_weight(omegas,
tetrahedra_omegas,
function='I'):
@ -119,6 +124,7 @@ def get_tetrahedra_integration_weight(omegas,
function)
return integration_weights
class TetrahedronMethod(object):
def __init__(self,
primitive_vectors=None, # column vectors

5
setup.py
View File

@ -45,11 +45,12 @@ if (config_var is not None and
######################
# _phonopy extension #
######################
include_dirs_phonopy = ['c/harmonic_h', 'c/kspclib_h'] + include_dirs_numpy
include_dirs_phonopy = (['c/harmonic_h', 'c/kspclib_h', 'c/spglib_h']
+ include_dirs_numpy)
sources_phonopy = ['c/_phonopy.c',
'c/harmonic/dynmat.c',
'c/harmonic/derivative_dynmat.c',
'c/kspclib/kgrid.c',
'c/spglib/kgrid.c',
'c/kspclib/tetrahedron_method.c']
if with_openmp: