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

Spglib update to version 1.9.2

This commit is contained in:
Atsushi Togo 2016-04-22 10:41:52 +09:00
parent 893ec5470d
commit 2c903fe1c5
22 changed files with 865 additions and 1085 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
c/_spglib.c
View File

@ -60,6 +60,7 @@ static PyObject * get_grid_points_by_rotations(PyObject *self, PyObject *args);
static PyObject * get_BZ_grid_points_by_rotations(PyObject *self, PyObject *args); static PyObject * get_BZ_grid_points_by_rotations(PyObject *self, PyObject *args);
static PyObject * relocate_BZ_grid_address(PyObject *self, PyObject *args); static PyObject * relocate_BZ_grid_address(PyObject *self, PyObject *args);
static PyObject * get_symmetry_from_database(PyObject *self, PyObject *args); static PyObject * get_symmetry_from_database(PyObject *self, PyObject *args);
static PyObject * py_niggli_reduce(PyObject *self, PyObject *args);
struct module_state { struct module_state {
PyObject *error; PyObject *error;
@ -107,6 +108,7 @@ static PyMethodDef _spglib_methods[] = {
"Rotated grid points in BZ are returned"}, "Rotated grid points in BZ are returned"},
{"BZ_grid_address", relocate_BZ_grid_address, METH_VARARGS, {"BZ_grid_address", relocate_BZ_grid_address, METH_VARARGS,
"Relocate grid addresses inside Brillouin zone"}, "Relocate grid addresses inside Brillouin zone"},
{"niggli_reduce", py_niggli_reduce, METH_VARARGS, "Niggli reduction"},
{NULL, NULL, 0, NULL} {NULL, NULL, 0, NULL}
}; };
@ -363,7 +365,7 @@ static PyObject * get_spacegroup_type(PyObject *self, PyObject *args)
PyObject *array; PyObject *array;
SpglibSpacegroupType symbols; SpglibSpacegroupType symbols;
if (!PyArg_ParseTuple(args, "i",&hall_number)) { if (!PyArg_ParseTuple(args, "i", &hall_number)) {
return NULL; return NULL;
} }
@ -819,3 +821,18 @@ static PyObject * relocate_BZ_grid_address(PyObject *self, PyObject *args)
return PyLong_FromLong((long) num_ir_gp); return PyLong_FromLong((long) num_ir_gp);
} }
static PyObject * py_niggli_reduce(PyObject *self, PyObject *args)
{
PyArrayObject* lattice_py;
double eps;
if (!PyArg_ParseTuple(args, "Od", &lattice_py, &eps)) {
return NULL;
}
double (*lattice)[3] = (double(*)[3])PyArray_DATA(lattice_py);
int result = spg_niggli_reduce(lattice, eps);
return PyLong_FromLong((long) result);
}

96
c/spglib/cell.c
View File

@ -58,7 +58,8 @@ static VecDBL *
translate_atoms_in_trimmed_lattice(SPGCONST Cell * cell, translate_atoms_in_trimmed_lattice(SPGCONST Cell * cell,
SPGCONST double prim_lat[3][3]); SPGCONST double prim_lat[3][3]);
static int * get_overlap_table(const VecDBL * position, static int * get_overlap_table(const VecDBL * position,
SPGCONST Cell * cell, const int cell_size,
const int * cell_types,
SPGCONST Cell * trimmed_cell, SPGCONST Cell * trimmed_cell,
const double symprec); const double symprec);
@ -67,7 +68,13 @@ Cell * cel_alloc_cell(const int size)
{ {
Cell *cell; Cell *cell;
int i, j; int i, j;
cell = NULL;
if (size < 1) {
return NULL;
}
cell = NULL; cell = NULL;
if ((cell = (Cell*) malloc(sizeof(Cell))) == NULL) { if ((cell = (Cell*) malloc(sizeof(Cell))) == NULL) {
@ -81,23 +88,21 @@ Cell * cel_alloc_cell(const int size)
} }
} }
cell->size = size; cell->size = size;
if (size > 0) { if ((cell->types = (int *) malloc(sizeof(int) * size)) == NULL) {
if ((cell->types = (int *) malloc(sizeof(int) * size)) == NULL) { warning_print("spglib: Memory could not be allocated.");
warning_print("spglib: Memory could not be allocated."); free(cell);
free(cell); cell = NULL;
cell = NULL; return NULL;
return NULL; }
} if ((cell->position =
if ((cell->position = (double (*)[3]) malloc(sizeof(double[3]) * size)) == NULL) {
(double (*)[3]) malloc(sizeof(double[3]) * size)) == NULL) { warning_print("spglib: Memory could not be allocated.");
warning_print("spglib: Memory could not be allocated."); free(cell->types);
free(cell->types); cell->types = NULL;
cell->types = NULL; free(cell);
free(cell); cell = NULL;
cell = NULL; return NULL;
return NULL;
}
} }
return cell; return cell;
@ -105,14 +110,17 @@ Cell * cel_alloc_cell(const int size)
void cel_free_cell(Cell * cell) void cel_free_cell(Cell * cell)
{ {
if (cell->size > 0) { if (cell != NULL) {
free(cell->position); if (cell->position != NULL) {
cell->position = NULL; free(cell->position);
free(cell->types); cell->position = NULL;
cell->types = NULL; }
if (cell->types != NULL) {
free(cell->types);
cell->types = NULL;
}
free(cell);
} }
free (cell);
cell = NULL;
} }
void cel_set_cell(Cell * cell, void cel_set_cell(Cell * cell,
@ -135,7 +143,7 @@ Cell * cel_copy_cell(SPGCONST Cell * cell)
Cell * cell_new; Cell * cell_new;
cell_new = NULL; cell_new = NULL;
if ((cell_new = cel_alloc_cell(cell->size)) == NULL) { if ((cell_new = cel_alloc_cell(cell->size)) == NULL) {
return NULL; return NULL;
} }
@ -162,7 +170,7 @@ int cel_is_overlap(const double a[3],
} }
mat_multiply_matrix_vector_d3(v_diff, lattice, v_diff); mat_multiply_matrix_vector_d3(v_diff, lattice, v_diff);
if ( mat_norm_squared_d3(v_diff) < symprec * symprec) { if (mat_norm_squared_d3(v_diff) < symprec * symprec) {
return 1; return 1;
} else { } else {
return 0; return 0;
@ -196,8 +204,8 @@ static Cell * trim_cell(int * mapping_table,
overlap_table = NULL; overlap_table = NULL;
trimmed_cell = NULL; trimmed_cell = NULL;
ratio = mat_Nint(mat_get_determinant_d3(cell->lattice) / ratio = abs(mat_Nint(mat_get_determinant_d3(cell->lattice) /
mat_get_determinant_d3(trimmed_lattice)); mat_get_determinant_d3(trimmed_lattice)));
/* Check if cell->size is dividable by ratio */ /* Check if cell->size is dividable by ratio */
if ((cell->size / ratio) * ratio != cell->size) { if ((cell->size / ratio) * ratio != cell->size) {
@ -212,17 +220,21 @@ static Cell * trim_cell(int * mapping_table,
trimmed_lattice)) trimmed_lattice))
== NULL) { == NULL) {
cel_free_cell(trimmed_cell); cel_free_cell(trimmed_cell);
trimmed_cell = NULL;
goto err; goto err;
} }
mat_copy_matrix_d3(trimmed_cell->lattice, trimmed_lattice); mat_copy_matrix_d3(trimmed_cell->lattice, trimmed_lattice);
if ((overlap_table = get_overlap_table(position, if ((overlap_table = get_overlap_table(position,
cell, cell->size,
cell->types,
trimmed_cell, trimmed_cell,
symprec)) == NULL) { symprec)) == NULL) {
mat_free_VecDBL(position); mat_free_VecDBL(position);
position = NULL;
cel_free_cell(trimmed_cell); cel_free_cell(trimmed_cell);
trimmed_cell = NULL;
goto err; goto err;
} }
@ -243,6 +255,7 @@ static Cell * trim_cell(int * mapping_table,
overlap_table); overlap_table);
mat_free_VecDBL(position); mat_free_VecDBL(position);
position = NULL;
free(overlap_table); free(overlap_table);
return trimmed_cell; return trimmed_cell;
@ -281,7 +294,7 @@ static void set_positions(Cell * trimmed_cell,
trimmed_cell->position[j][l] += position->vec[i][l]; trimmed_cell->position[j][l] += position->vec[i][l];
} }
} }
} }
multi = position->size / trimmed_cell->size; multi = position->size / trimmed_cell->size;
@ -327,7 +340,8 @@ translate_atoms_in_trimmed_lattice(SPGCONST Cell * cell,
/* Return NULL if failed */ /* Return NULL if failed */
static int * get_overlap_table(const VecDBL * position, static int * get_overlap_table(const VecDBL * position,
SPGCONST Cell * cell, const int cell_size,
const int * cell_types,
SPGCONST Cell * trimmed_cell, SPGCONST Cell * trimmed_cell,
const double symprec) const double symprec)
{ {
@ -337,18 +351,18 @@ static int * get_overlap_table(const VecDBL * position,
trim_tolerance = symprec; trim_tolerance = symprec;
ratio = cell->size / trimmed_cell->size; ratio = cell_size / trimmed_cell->size;
if ((overlap_table = (int*)malloc(sizeof(int) * cell->size)) == NULL) { if ((overlap_table = (int*)malloc(sizeof(int) * cell_size)) == NULL) {
return NULL; return NULL;
} }
for (attempt = 0; attempt < 100; attempt++) { for (attempt = 0; attempt < 100; attempt++) {
for (i = 0; i < cell->size; i++) { for (i = 0; i < cell_size; i++) {
overlap_table[i] = -1; overlap_table[i] = -1;
num_overlap = 0; num_overlap = 0;
for (j = 0; j < cell->size; j++) { for (j = 0; j < cell_size; j++) {
if (cell->types[i] == cell->types[j]) { if (cell_types[i] == cell_types[j]) {
if (cel_is_overlap(position->vec[i], if (cel_is_overlap(position->vec[i],
position->vec[j], position->vec[j],
trimmed_cell->lattice, trimmed_cell->lattice,
@ -379,12 +393,12 @@ static int * get_overlap_table(const VecDBL * position,
} }
} }
for (i = 0; i < cell->size; i++) { for (i = 0; i < cell_size; i++) {
if (overlap_table[i] != i) { if (overlap_table[i] != i) {
continue; continue;
} }
count = 0; count = 0;
for (j = 0; j < cell->size; j++) { for (j = 0; j < cell_size; j++) {
if (i == overlap_table[j]) { if (i == overlap_table[j]) {
count++; count++;
} }

3
c/spglib/kpoint.c
View File

@ -252,7 +252,9 @@ int kpt_get_stabilized_reciprocal_mesh(int grid_address[][3],
rot_reciprocal_q); rot_reciprocal_q);
mat_free_MatINT(rot_reciprocal_q); mat_free_MatINT(rot_reciprocal_q);
rot_reciprocal_q = NULL;
mat_free_MatINT(rot_reciprocal); mat_free_MatINT(rot_reciprocal);
rot_reciprocal = NULL;
return num_ir; return num_ir;
} }
@ -362,6 +364,7 @@ static MatINT *get_point_group_reciprocal(const MatINT * rotations,
if ((unique_rot = (int*)malloc(sizeof(int) * rot_reciprocal->size)) == NULL) { if ((unique_rot = (int*)malloc(sizeof(int) * rot_reciprocal->size)) == NULL) {
warning_print("spglib: Memory of unique_rot could not be allocated."); warning_print("spglib: Memory of unique_rot could not be allocated.");
mat_free_MatINT(rot_reciprocal); mat_free_MatINT(rot_reciprocal);
rot_reciprocal = NULL;
return NULL; return NULL;
} }

396
c/spglib/lattice.c
View File

@ -1,396 +0,0 @@
/* Copyright (C) 2010 Atsushi Togo */
/* All rights reserved. */
/* This file is part of spglib. */
/* 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 <stdio.h>
#include <stdlib.h>
#include "lattice.h"
#include "mathfunc.h"
#include "debug.h"
static int get_Delaunay_reduction(double red_lattice[3][3],
SPGCONST double lattice[3][3],
SPGCONST double symprec);
static int get_Delaunay_reduction_basis(double basis[4][3],
const double symprec);
static void get_Delaunay_shortest_vectors(double basis[4][3],
const double symprec);
static void get_exteneded_basis(double basis[4][3],
SPGCONST double lattice[3][3]);
static int get_Delaunay_reduction_2D(double red_lattice[3][3],
SPGCONST double lattice[3][3],
const int unique_axis,
const double symprec);
static int get_Delaunay_reduction_basis_2D(double basis[3][3],
const double symprec);
static void get_Delaunay_shortest_vectors_2D(double basis[3][3],
const double unique_vec[3],
const double symprec);
static void get_exteneded_basis_2D(double basis[3][3],
SPGCONST double lattice[3][2]);
/* Return 0 if failed */
int lat_smallest_lattice_vector(double min_lattice[3][3],
SPGCONST double lattice[3][3],
const double symprec)
{
debug_print("lat_smallest_lattice_vector (tolerance = %f):\n", symprec);
return get_Delaunay_reduction(min_lattice, lattice, symprec);
}
int lat_smallest_lattice_vector_2D(double min_lattice[3][3],
SPGCONST double lattice[3][3],
const int unique_axis,
const double symprec)
{
debug_print("lat_smallest_lattice_vector_2D:\n");
return get_Delaunay_reduction_2D(min_lattice, lattice, unique_axis, symprec);
}
/* Delaunay reduction */
/* Reference can be found in International table A. */
/* Return 0 if failed */
static int get_Delaunay_reduction(double red_lattice[3][3],
SPGCONST double lattice[3][3],
const double symprec)
{
int i, j;
double volume, sum;
double basis[4][3];
get_exteneded_basis(basis, lattice);
sum = 0;
for (i = 0; i < 4; i++) {
for (j = 0; j < 3; j++) {
sum += basis[i][j] * basis[i][j];
}
}
while (1) {
if (get_Delaunay_reduction_basis(basis, symprec)) {
break;
}
}
get_Delaunay_shortest_vectors(basis, symprec);
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
red_lattice[i][j] = basis[j][i];
}
}
volume = mat_get_determinant_d3(red_lattice);
if (mat_Dabs(volume) < symprec) {
warning_print("spglib: Minimum lattice has no volume (line %d, %s).\n", __LINE__, __FILE__);
goto err;
}
if (volume < 0) {
/* Flip axes */
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
red_lattice[i][j] = -red_lattice[i][j];
}
}
}
return 1;
err:
return 0;
}
static void get_Delaunay_shortest_vectors(double basis[4][3],
const double symprec)
{
int i, j;
double tmpmat[3][3], b[7][3], tmpvec[3];
/* Search in the set {b1, b2, b3, b4, b1+b2, b2+b3, b3+b1} */
for (i = 0; i < 4; i++) {
for (j = 0; j < 3; j++) {
b[i][j] = basis[i][j];
}
}
for (i = 0; i < 3; i++) {
b[4][i] = basis[0][i] + basis[1][i];
}
for (i = 0; i < 3; i++) {
b[5][i] = basis[1][i] + basis[2][i];
}
for (i = 0; i < 3; i++) {
b[6][i] = basis[2][i] + basis[0][i];
}
/* Bubble sort */
for (i = 0; i < 6; i++) {
for (j = 0; j < 6; j++) {
if (mat_norm_squared_d3(b[j]) > mat_norm_squared_d3(b[j+1])) {
mat_copy_vector_d3(tmpvec, b[j]);
mat_copy_vector_d3(b[j], b[j+1]);
mat_copy_vector_d3(b[j+1], tmpvec);
}
}
}
for (i = 2; i < 7; i++) {
for (j = 0; j < 3; j++) {
tmpmat[j][0] = b[0][j];
tmpmat[j][1] = b[1][j];
tmpmat[j][2] = b[i][j];
}
if (mat_Dabs(mat_get_determinant_d3(tmpmat)) > symprec) {
for (j = 0; j < 3; j++) {
basis[0][j] = b[0][j];
basis[1][j] = b[1][j];
basis[2][j] = b[i][j];
}
break;
}
}
}
static int get_Delaunay_reduction_basis(double basis[4][3],
const double symprec)
{
int i, j, k, l;
double dot_product;
for (i = 0; i < 4; i++) {
for (j = i+1; j < 4; j++) {
dot_product = 0.0;
for (k = 0; k < 3; k++) {
dot_product += basis[i][k] * basis[j][k];
}
if (dot_product > symprec) {
for (k = 0; k < 4; k++) {
if (! (k == i || k == j)) {
for (l = 0; l < 3; l++) {
basis[k][l] += basis[i][l];
}
}
}
for (k = 0; k < 3; k++) {
basis[i][k] = -basis[i][k];
}
return 0;
}
}
}
return 1;
}
static void get_exteneded_basis(double basis[4][3],
SPGCONST double lattice[3][3])
{
int i, j;
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
basis[i][j] = lattice[j][i];
}
}
for (i = 0; i < 3; i++) {
basis[3][i] = -lattice[i][0] -lattice[i][1] -lattice[i][2];
}
}
static int get_Delaunay_reduction_2D(double red_lattice[3][3],
SPGCONST double lattice[3][3],
const int unique_axis,
const double symprec)
{
int i, j, k;
double volume;
double basis[3][3], lattice_2D[3][2], unique_vec[3];
k = 0;
for (i = 0; i < 3; i++) {
unique_vec[i] = lattice[i][unique_axis];
}
for (i = 0; i < 3; i++) {
if (i != unique_axis) {
for (j = 0; j < 3; j++) {
lattice_2D[j][k] = lattice[j][i];
}
k++;
}
}
get_exteneded_basis_2D(basis, lattice_2D);
while (1) {
if (get_Delaunay_reduction_basis_2D(basis, symprec)) {
break;
}
}
get_Delaunay_shortest_vectors_2D(basis, unique_vec, symprec);
k = 0;
for (i = 0; i < 3; i++) {
if (i == unique_axis) {
for (j = 0; j < 3; j++) {
red_lattice[j][i] = lattice[j][i];
}
} else {
for (j = 0; j < 3; j++) {
red_lattice[j][i] = basis[k][j];
}
k++;
}
}
volume = mat_get_determinant_d3(red_lattice);
if (mat_Dabs(volume) < symprec) {
warning_print("spglib: Minimum lattice has no volume (line %d, %s).\n", __LINE__, __FILE__);
goto err;
}
if (volume < 0) {
for (i = 0; i < 3; i++) {
red_lattice[i][unique_axis] = -red_lattice[i][unique_axis];
}
}
return 1;
err:
return 0;
}
static int get_Delaunay_reduction_basis_2D(double basis[3][3],
const double symprec)
{
int i, j, k, l;
double dot_product;
for (i = 0; i < 3; i++) {
for (j = i + 1; j < 3; j++) {
dot_product = 0.0;
for (k = 0; k < 3; k++) {
dot_product += basis[i][k] * basis[j][k];
}
if (dot_product > symprec) {
for (k = 0; k < 3; k++) {
if (! (k == i || k == j)) {
for (l = 0; l < 3; l++) {
basis[k][l] += 2 * basis[i][l];
}
break;
}
}
for (k = 0; k < 3; k++) {
basis[i][k] = -basis[i][k];
}
return 0;
}
}
}
return 1;
}
static void get_Delaunay_shortest_vectors_2D(double basis[3][3],
const double unique_vec[3],
const double symprec)
{
int i, j;
double b[4][3], tmpmat[3][3];
double tmpvec[3];
/* Search in the set {b1, b2, b3, b1+b2} */
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
b[i][j] = basis[i][j];
}
}
for (i = 0; i < 3; i++) {
b[3][i] = basis[0][i] + basis[1][i];
}
/* Bubble sort */
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
if (mat_norm_squared_d3(b[j]) > mat_norm_squared_d3(b[j + 1])) {
mat_copy_vector_d3(tmpvec, b[j]);
mat_copy_vector_d3(b[j], b[j + 1]);
mat_copy_vector_d3(b[j + 1], tmpvec);
}
}
}
for (i = 0; i < 3; i++) {
tmpmat[i][0] = b[0][i];
tmpmat[i][1] = unique_vec[i];
}
for (i = 1; i < 4; i++) {
for (j = 0; j < 3; j++) {
tmpmat[j][2] = b[i][j];
}
if (mat_Dabs(mat_get_determinant_d3(tmpmat)) > symprec) {
for (j = 0; j < 3; j++) {
basis[0][j] = b[0][j];
basis[1][j] = b[i][j];
}
break;
}
}
}
static void get_exteneded_basis_2D(double basis[3][3],
SPGCONST double lattice[3][2])
{
int i, j;
for (i = 0; i < 2; i++) {
for (j = 0; j < 3; j++) {
basis[i][j] = lattice[j][i];
}
}
for (i = 0; i < 3; i++) {
basis[2][i] = -lattice[i][0] -lattice[i][1];
}
}

2
c/spglib/mathfunc.c
View File

@ -455,7 +455,6 @@ void mat_free_MatINT(MatINT * matint)
matint->mat = NULL; matint->mat = NULL;
} }
free(matint); free(matint);
matint = NULL;
} }
VecDBL * mat_alloc_VecDBL(const int size) VecDBL * mat_alloc_VecDBL(const int size)
@ -490,7 +489,6 @@ void mat_free_VecDBL(VecDBL * vecdbl)
vecdbl->vec = NULL; vecdbl->vec = NULL;
} }
free(vecdbl); free(vecdbl);
vecdbl = NULL;
} }

583
c/spglib/niggli.c
View File

@ -38,41 +38,59 @@
#include <string.h> #include <string.h>
#include "niggli.h" #include "niggli.h"
static double A, B, C, eta, xi, zeta, eps; #define NIGGLI_MAX_NUM_LOOP 100
static int l, m, n;
static double *tmat = NULL;
static double *lattice = NULL;
static int initialize(const double *lattice_, const double eps_); typedef struct {
static void finalize(double *lattice_); double A;
static void reset(void); double B;
static void step0(void); double C;
static int step1(void); double eta;
static int step2(void); double xi;
static int step3(void); double zeta;
static int step4(void); double eps;
static int step5(void); int l;
static int step6(void); int m;
static int step7(void); int n;
static int step8(void); double *tmat;
static void set_parameters(void); double *lattice;
static void set_angle_types(void); } NiggliParams;
static NiggliParams * initialize(const double *lattice_, const double eps_);
static void finalize(double *lattice_, NiggliParams *p);
static int reset(NiggliParams *p);
static int step1(NiggliParams *p);
static int step2(NiggliParams *p);
static int step3(NiggliParams *p);
static int step4(NiggliParams *p);
static int step5(NiggliParams *p);
static int step6(NiggliParams *p);
static int step7(NiggliParams *p);
static int step8(NiggliParams *p);
static int set_parameters(NiggliParams *p);
static void set_angle_types(NiggliParams *p);
static double * get_transpose(const double *M); static double * get_transpose(const double *M);
static double * get_metric(const double *M); static double * get_metric(const double *M);
static double * multiply_matrices(const double *A, const double *B); static double * multiply_matrices(const double *A, const double *B);
#ifdef NIGGLI_DEBUG #ifdef NIGGLI_DEBUG
#define debug_print(...) printf(__VA_ARGS__) #define debug_print(...) printf(__VA_ARGS__)
static void debug_show(void); static void debug_show(const int j, const NiggliParams *p);
static void debug_show(void) static void debug_show(const int j, const NiggliParams *p)
{ {
int i; /* int i; */
printf("%f %f %f %f %f %f\n", A, B, C, xi, eta, zeta);
printf("%d %d %d\n", l, m, n); if (j < 0) {
printf("Finish: ");
for (i = 0; i < 3; i++) { } else {
printf("%f %f %f\n", lattice[i * 3], lattice[i * 3 + 1], lattice[i * 3 + 2]); printf("Step %d: ", j);
} }
printf("%f %f %f %f %f %f\n", p->A, p->B, p->C, p->xi, p->eta, p->zeta);
/* printf("%d %d %d\n", p->l, p->m, p->n); */
/* for (i = 0; i < 3; i++) { */
/* printf("%f %f %f\n", */
/* p->lattice[i * 3], p->lattice[i * 3 + 1], p->lattice[i * 3 + 2]); */
/* } */
} }
#else #else
#define debug_print(...) #define debug_print(...)
@ -86,280 +104,303 @@ static void debug_show(void)
#define warning_print(...) #define warning_print(...)
#endif #endif
/*--------------------------------------------*/
/* Version: niggli-[major].[minor].[micro] */
/*--------------------------------------------*/
int niggli_get_major_version(void)
{
return NIGGLI_MAJOR_VERSION;
}
int niggli_get_minor_version(void)
{
return NIGGLI_MINOR_VERSION;
}
int niggli_get_micro_version(void)
{
return NIGGLI_MICRO_VERSION;
}
/* return 0 if failed */ /* return 0 if failed */
int niggli_reduce(double *lattice_, const double eps_) int niggli_reduce(double *lattice_, const double eps_)
{ {
int i; int i, j, succeeded;
NiggliParams *p;
int (*steps[8])(NiggliParams *p) = {step1, step2, step3, step4,
step5, step6, step7, step8};
if (! initialize(lattice_, eps_)) { p = NULL;
succeeded = 0;
if ((p = initialize(lattice_, eps_)) == NULL) {
return 0; return 0;
} }
step0(); /* Step 0 */
if (! set_parameters(p)) {
goto ret;
}
for (i = 0; i < NIGGLI_MAX_NUM_LOOP; i++) {
for (j = 0; j < 8; j++) {
if ((*steps[j])(p)) {
debug_show(j + 1, p);
if (! reset(p)) {goto ret;}
if (j == 1 || j == 4 || j == 5 || j == 6 || j == 7) {break;}
}
}
if (j == 8) {
succeeded = 1;
break;
}
}
debug_show(-1, p);
ret:
finalize(lattice_, p);
return succeeded;
}
static NiggliParams * initialize(const double *lattice_, const double eps_)
{
NiggliParams * p;
p = NULL;
if ((p = (NiggliParams*)malloc(sizeof(NiggliParams))) == NULL) {
warning_print("niggli: Memory could not be allocated.");
return 0;
}
p->A = 0;
p->B = 0;
p->C = 0;
p->eta = 0;
p->xi = 0;
p->zeta = 0;
p->eps = 0;
p->l = 0;
p->m = 0;
p->n = 0;
p->tmat = NULL;
p->lattice = NULL;
if ((p->tmat = (double*)malloc(sizeof(double) * 9)) == NULL) {
warning_print("niggli: Memory could not be allocated.");
free(p);
p = NULL;
return NULL;
}
p->eps = eps_;
if ((p->lattice = (double*)malloc(sizeof(double) * 9)) == NULL) {
warning_print("niggli: Memory could not be allocated.");
free(p->tmat);
p->tmat = NULL;
free(p);
p = NULL;
return NULL;
}
for (i = 0; i < 10; i++) { memcpy(p->lattice, lattice_, sizeof(double) * 9);
if (step1()) {
debug_print("step1\n");
debug_show();
debug_print("\n");
}
if (step2()) { return p;
debug_print("step2\n");
debug_show();
debug_print("\n");
continue;
}
if (step3()) {
debug_print("step3\n");
debug_show();
debug_print("\n");
}
if (step4()) {
debug_print("step4\n");
debug_show();
debug_print("\n");
}
if (step5()) {
debug_print("step5\n");
debug_show();
debug_print("\n");
continue;
}
if (step6()) {
debug_print("step6\n");
debug_show();
debug_print("\n");
continue;
}
if (step7()) {
debug_print("step7\n");
debug_show();
debug_print("\n");
continue;
}
if (step8()) {
debug_print("step7\n");
debug_show();
debug_print("\n");
continue;
}
break;
}
finalize(lattice_);
return 1;
} }
static int initialize(const double *lattice_, const double eps_) static void finalize(double *lattice_, NiggliParams *p)
{ {
if ((tmat = (double*)malloc(sizeof(double) * 9)) == NULL) { free(p->tmat);
warning_print("niggli: Memory could not be allocated."); p->tmat = NULL;
return 0; memcpy(lattice_, p->lattice, sizeof(double) * 9);
} free(p->lattice);
p->lattice = NULL;
eps = eps_; free(p);
if ((lattice = (double*)malloc(sizeof(double) * 9)) == NULL) { p = NULL;
warning_print("niggli: Memory could not be allocated.");
free(tmat);
tmat = NULL;
return 0;
}
memcpy(lattice, lattice_, sizeof(double) * 9);
return 1;
} }
static void finalize(double *lattice_) static int reset(NiggliParams *p)
{
free(tmat);
tmat = NULL;
memcpy(lattice_, lattice, sizeof(double) * 9);
free(lattice);
lattice = NULL;
}
static void reset(void)
{ {
double *lat_tmp; double *lat_tmp;
lat_tmp = NULL; lat_tmp = NULL;
lat_tmp = multiply_matrices(lattice, tmat); if ((lat_tmp = multiply_matrices(p->lattice, p->tmat)) == NULL) {return 0;}
memcpy(p->lattice, lat_tmp, sizeof(double) * 9);
memcpy(lattice, lat_tmp, sizeof(double) * 9);
step0();
free(lat_tmp); free(lat_tmp);
lat_tmp = NULL; lat_tmp = NULL;
return set_parameters(p);
} }
static void step0(void) static int set_parameters(NiggliParams *p)
{
set_parameters();
set_angle_types();
}
static int step1(void)
{
if (A > B + eps ||
(! (fabs(A -B) > eps) && fabs(xi) > fabs(eta) + eps)) {
tmat[0] = 0, tmat[1] = -1, tmat[2] = 0;
tmat[3] = -1, tmat[4] = 0, tmat[5] = 0;
tmat[6] = 0, tmat[7] = 0, tmat[8] = -1;
reset();
return 1;
}
else {return 0;}
}
static int step2(void)
{
if (B > C + eps ||
(! (fabs(B - C) > eps) && fabs(eta) > fabs(zeta) + eps)) {
tmat[0] = -1, tmat[1] = 0, tmat[2] = 0;
tmat[3] = 0, tmat[4] = 0, tmat[5] = -1;
tmat[6] = 0, tmat[7] = -1, tmat[8] = 0;
reset();
return 1;
}
else {return 0;}
}
static int step3(void)
{
int i, j, k;
if (l * m * n == 1) {
if (l == -1) {i = -1;} else {i = 1;}
if (m == -1) {j = -1;} else {j = 1;}
if (n == -1) {k = -1;} else {k = 1;}
tmat[0] = i, tmat[1] = 0, tmat[2] = 0;
tmat[3] = 0, tmat[4] = j, tmat[5] = 0;
tmat[6] = 0, tmat[7] = 0, tmat[8] = k;
reset();
return 1;
}
else {return 0;}
}
static int step4(void)
{
int i, j, k;
if (l * m * n == 0 || l * m * n == -1) {
if (l == -1) {i = -1;} else {i = 1;}
if (m == -1) {j = -1;} else {j = 1;}
if (n == -1) {k = -1;} else {k = 1;}
if (i * j * k == -1) {
if (l == 0) {i = -1;}
if (m == 0) {j = -1;}
if (n == 0) {k = -1;}
}
tmat[0] = i, tmat[1] = 0, tmat[2] = 0;
tmat[3] = 0, tmat[4] = j, tmat[5] = 0;
tmat[6] = 0, tmat[7] = 0, tmat[8] = k;
reset();
return 1;
}
else {return 0;}
}
static int step5(void)
{
if (fabs(xi) > B + eps ||
(! (fabs(B - xi) > eps) && 2 * eta < zeta - eps) ||
(! (fabs(B + xi) > eps) && zeta < -eps)) {
tmat[0] = 1, tmat[1] = 0, tmat[2] = 0;
tmat[3] = 0, tmat[4] = 1, tmat[5] = 0;
tmat[6] = 0, tmat[7] = 0, tmat[8] = 1;
if (xi > 0) {tmat[5] = -1;}
if (xi < 0) {tmat[5] = 1;}
reset();
return 1;
}
else {return 0;}
}
static int step6(void)
{
if (fabs(eta) > A + eps ||
(! (fabs(A - eta) > eps) && 2 * xi < zeta - eps) ||
(! (fabs(A + eta) > eps) && zeta < -eps)) {
tmat[0] = 1, tmat[1] = 0, tmat[2] = 0;
tmat[3] = 0, tmat[4] = 1, tmat[5] = 0;
tmat[6] = 0, tmat[7] = 0, tmat[8] = 1;
if (eta > 0) {tmat[2] = -1;}
if (eta < 0) {tmat[2] = 1;}
reset();
return 1;
}
else {return 0;}
}
static int step7(void)
{
if (fabs(zeta) > A + eps ||
(! (fabs(A - zeta) > eps) && 2 * xi < eta - eps) ||
(! (fabs(A + zeta) > eps) && eta < -eps)) {
tmat[0] = 1, tmat[1] = 0, tmat[2] = 0;
tmat[3] = 0, tmat[4] = 1, tmat[5] = 0;
tmat[6] = 0, tmat[7] = 0, tmat[8] = 1;
if (zeta > 0) {tmat[1] = -1;}
if (zeta < 0) {tmat[1] = 1;}
reset();
return 1;
}
else {return 0;}
}
static int step8(void)
{
if (xi + eta + zeta + A + B < -eps ||
(! (fabs(xi + eta + zeta + A + B) > eps) && 2 * (A + eta) + zeta > eps)) {
tmat[0] = 1, tmat[1] = 0, tmat[2] = 1;
tmat[3] = 0, tmat[4] = 1, tmat[5] = 1;
tmat[6] = 0, tmat[7] = 0, tmat[8] = 1;
reset();
return 1;
}
else {return 0;}
}
static void set_angle_types(void)
{
l = 0, m = 0, n = 0;
if (xi < -eps) {l = -1;}
if (xi > eps) {l = 1;}
if (eta < -eps) {m = -1;}
if (eta > eps) {m = 1;}
if (zeta < -eps) {n = -1;}
if (zeta > eps) {n = 1;}
}
static void set_parameters(void)
{ {
double *G; double *G;
G = NULL; G = NULL;
G = get_metric(lattice); if ((G = get_metric(p->lattice)) == NULL) {return 0;}
A = G[0]; p->A = G[0];
B = G[4]; p->B = G[4];
C = G[8]; p->C = G[8];
xi = G[5] * 2; p->xi = G[5] * 2;
eta = G[2] * 2; p->eta = G[2] * 2;
zeta = G[1] * 2; p->zeta = G[1] * 2;
free(G); free(G);
G = NULL; G = NULL;
set_angle_types(p);
return 1;
}
static void set_angle_types(NiggliParams *p)
{
p->l = 0;
p->m = 0;
p->n = 0;
if (p->xi < -p->eps) {p->l = -1;}
if (p->xi > p->eps) {p->l = 1;}
if (p->eta < -p->eps) {p->m = -1;}
if (p->eta > p->eps) {p->m = 1;}
if (p->zeta < -p->eps) {p->n = -1;}
if (p->zeta > p->eps) {p->n = 1;}
}
static int step1(NiggliParams *p)
{
if (p->A > p->B + p->eps ||
(! (fabs(p->A - p->B) > p->eps) &&
fabs(p->xi) > fabs(p->eta) + p->eps)) {
p->tmat[0] = 0, p->tmat[1] = -1, p->tmat[2] = 0;
p->tmat[3] = -1, p->tmat[4] = 0, p->tmat[5] = 0;
p->tmat[6] = 0, p->tmat[7] = 0, p->tmat[8] = -1;
return 1;
}
else {return 0;}
}
static int step2(NiggliParams *p)
{
if (p->B > p->C + p->eps ||
(! (fabs(p->B - p->C) > p->eps)
&& fabs(p->eta) > fabs(p->zeta) + p->eps)) {
p->tmat[0] = -1, p->tmat[1] = 0, p->tmat[2] = 0;
p->tmat[3] = 0, p->tmat[4] = 0, p->tmat[5] = -1;
p->tmat[6] = 0, p->tmat[7] = -1, p->tmat[8] = 0;
return 1;
}
else {return 0;}
}
static int step3(NiggliParams *p)
{
int i, j, k;
if (p->l * p->m * p->n == 1) {
if (p->l == -1) {i = -1;} else {i = 1;}
if (p->m == -1) {j = -1;} else {j = 1;}
if (p->n == -1) {k = -1;} else {k = 1;}
p->tmat[0] = i, p->tmat[1] = 0, p->tmat[2] = 0;
p->tmat[3] = 0, p->tmat[4] = j, p->tmat[5] = 0;
p->tmat[6] = 0, p->tmat[7] = 0, p->tmat[8] = k;
return 1;
}
else {return 0;}
}
static int step4(NiggliParams *p)
{
int i, j, k, r;
if (p->l == -1 && p->m == -1 && p->n == -1) {
return 0;
}
if (p->l * p->m * p->n == 0 || p->l * p->m * p->n == -1) {
i = 1;
j = 1;
k = 1;
r = -1; /* 0: i, 1: j, 2: k */
if (p->l == 1) {i = -1;}
if (p->l == 0) {r = 0;}
if (p->m == 1) {j = -1;}
if (p->m == 0) {r = 1;}
if (p->n == 1) {k = -1;}
if (p->n == 0) {r = 2;}
if (i * j * k == -1) {
if (r == 0) {i = -1;}
if (r == 1) {j = -1;}
if (r == 2) {k = -1;}
}
p->tmat[0] = i, p->tmat[1] = 0, p->tmat[2] = 0;
p->tmat[3] = 0, p->tmat[4] = j, p->tmat[5] = 0;
p->tmat[6] = 0, p->tmat[7] = 0, p->tmat[8] = k;
return 1;
}
else {return 0;}
}
static int step5(NiggliParams *p)
{
if (fabs(p->xi) > p->B + p->eps ||
(! (fabs(p->B - p->xi) > p->eps) && 2 * p->eta < p->zeta - p->eps) ||
(! (fabs(p->B + p->xi) > p->eps) && p->zeta < -p->eps)) {
p->tmat[0] = 1, p->tmat[1] = 0, p->tmat[2] = 0;
p->tmat[3] = 0, p->tmat[4] = 1, p->tmat[5] = 0;
p->tmat[6] = 0, p->tmat[7] = 0, p->tmat[8] = 1;
if (p->xi > 0) {p->tmat[5] = -1;}
if (p->xi < 0) {p->tmat[5] = 1;}
return 1;
}
else {return 0;}
}
static int step6(NiggliParams *p)
{
if (fabs(p->eta) > p->A + p->eps ||
(! (fabs(p->A - p->eta) > p->eps) && 2 * p->xi < p->zeta - p->eps) ||
(! (fabs(p->A + p->eta) > p->eps) && p->zeta < -p->eps)) {
p->tmat[0] = 1, p->tmat[1] = 0, p->tmat[2] = 0;
p->tmat[3] = 0, p->tmat[4] = 1, p->tmat[5] = 0;
p->tmat[6] = 0, p->tmat[7] = 0, p->tmat[8] = 1;
if (p->eta > 0) {p->tmat[2] = -1;}
if (p->eta < 0) {p->tmat[2] = 1;}
return 1;
}
else {return 0;}
}
static int step7(NiggliParams *p)
{
if (fabs(p->zeta) > p->A + p->eps ||
(! (fabs(p->A - p->zeta) > p->eps) && 2 * p->xi < p->eta - p->eps) ||
(! (fabs(p->A + p->zeta) > p->eps) && p->eta < -p->eps)) {
p->tmat[0] = 1, p->tmat[1] = 0, p->tmat[2] = 0;
p->tmat[3] = 0, p->tmat[4] = 1, p->tmat[5] = 0;
p->tmat[6] = 0, p->tmat[7] = 0, p->tmat[8] = 1;
if (p->zeta > 0) {p->tmat[1] = -1;}
if (p->zeta < 0) {p->tmat[1] = 1;}
return 1;
}
else {return 0;}
}
static int step8(NiggliParams *p)
{
if (p->xi + p->eta + p->zeta + p->A + p->B < -p->eps ||
(! (fabs(p->xi + p->eta + p->zeta + p->A + p->B) > p->eps) &&
2 * (p->A + p->eta) + p->zeta > p->eps)) {
p->tmat[0] = 1, p->tmat[1] = 0, p->tmat[2] = 1;
p->tmat[3] = 0, p->tmat[4] = 1, p->tmat[5] = 1;
p->tmat[6] = 0, p->tmat[7] = 0, p->tmat[8] = 1;
return 1;
}
else {return 0;}
} }
static double * get_transpose(const double *M) static double * get_transpose(const double *M)
@ -379,6 +420,7 @@ static double * get_transpose(const double *M)
M_T[i * 3 + j] = M[j * 3 + i]; M_T[i * 3 + j] = M[j * 3 + i];
} }
} }
return M_T; return M_T;
} }
@ -389,9 +431,9 @@ static double * get_metric(const double *M)
G = NULL; G = NULL;
M_T = NULL; M_T = NULL;
M_T = get_transpose(M); if ((M_T = get_transpose(M)) == NULL) {return NULL;}
if ((G = multiply_matrices(M_T, M)) == NULL) {return NULL;}
G = multiply_matrices(M_T, M);
free(M_T); free(M_T);
M_T = NULL; M_T = NULL;
@ -418,5 +460,6 @@ static double * multiply_matrices(const double *L, const double *R)
} }
} }
} }
return M; return M;
} }

31
c/spglib/pointgroup.c
View File

@ -35,7 +35,6 @@
#include <string.h> #include <string.h>
#include <stdio.h> #include <stdio.h>
#include <stdlib.h> #include <stdlib.h>
#include "lattice.h"
#include "pointgroup.h" #include "pointgroup.h"
#include "symmetry.h" #include "symmetry.h"
#include "mathfunc.h" #include "mathfunc.h"
@ -376,8 +375,6 @@ static int rot_axes[][3] = {
static int get_pointgroup_number_by_rotations(SPGCONST int rotations[][3][3], static int get_pointgroup_number_by_rotations(SPGCONST int rotations[][3][3],
const int num_rotations); const int num_rotations);
static PointSymmetry get_pointsymmetry(SPGCONST int rotations[][3][3],
const int num_rotations);
static int get_pointgroup_number(SPGCONST PointSymmetry * pointsym); static int get_pointgroup_number(SPGCONST PointSymmetry * pointsym);
static int get_pointgroup_class_table(int table[10], static int get_pointgroup_class_table(int table[10],
SPGCONST PointSymmetry * pointsym); SPGCONST PointSymmetry * pointsym);
@ -442,7 +439,7 @@ Pointgroup ptg_get_transformation_matrix(int transform_mat[3][3],
if (pg_num > 0) { if (pg_num > 0) {
pointgroup = ptg_get_pointgroup(pg_num); pointgroup = ptg_get_pointgroup(pg_num);
pointsym = get_pointsymmetry(rotations, num_rotations); pointsym = ptg_get_pointsymmetry(rotations, num_rotations);
get_axes(axes, pointgroup.laue, &pointsym); get_axes(axes, pointgroup.laue, &pointsym);
set_transformation_matrix(transform_mat, axes); set_transformation_matrix(transform_mat, axes);
} else { } else {
@ -472,17 +469,8 @@ Pointgroup ptg_get_pointgroup(const int pointgroup_number)
return pointgroup; return pointgroup;
} }
static int get_pointgroup_number_by_rotations(SPGCONST int rotations[][3][3], PointSymmetry ptg_get_pointsymmetry(SPGCONST int rotations[][3][3],
const int num_rotations) const int num_rotations)
{
PointSymmetry pointsym;
pointsym = get_pointsymmetry(rotations, num_rotations);
return get_pointgroup_number(&pointsym);
}
static PointSymmetry get_pointsymmetry(SPGCONST int rotations[][3][3],
const int num_rotations)
{ {
int i, j; int i, j;
PointSymmetry pointsym; PointSymmetry pointsym;
@ -503,6 +491,15 @@ static PointSymmetry get_pointsymmetry(SPGCONST int rotations[][3][3],
return pointsym; return pointsym;
} }
static int get_pointgroup_number_by_rotations(SPGCONST int rotations[][3][3],
const int num_rotations)
{
PointSymmetry pointsym;
pointsym = ptg_get_pointsymmetry(rotations, num_rotations);
return get_pointgroup_number(&pointsym);
}
static int get_pointgroup_number(SPGCONST PointSymmetry * pointsym) static int get_pointgroup_number(SPGCONST PointSymmetry * pointsym)
{ {
int i, j, pg_num, counter; int i, j, pg_num, counter;
@ -739,7 +736,7 @@ static int laue2m(int axes[3],
static int lauemmm(int axes[3], static int lauemmm(int axes[3],
SPGCONST PointSymmetry * pointsym) SPGCONST PointSymmetry * pointsym)
{ {
int i, count, axis, tmpval; int i, count, axis;
int prop_rot[3][3]; int prop_rot[3][3];
@ -1024,7 +1021,7 @@ static int laue3m(int axes[3],
static int lauem3m(int axes[3], static int lauem3m(int axes[3],
SPGCONST PointSymmetry * pointsym) SPGCONST PointSymmetry * pointsym)
{ {
int i, count, axis, tmpval; int i, count, axis;
int prop_rot[3][3]; int prop_rot[3][3];
for (i = 0; i < 3; i++) { axes[i] = -1; } for (i = 0; i < 3; i++) { axes[i] = -1; }

48
c/spglib/primitive.c
View File

@ -35,7 +35,7 @@
#include <stdio.h> #include <stdio.h>
#include <stdlib.h> #include <stdlib.h>
#include "cell.h" #include "cell.h"
#include "lattice.h" #include "delaunay.h"
#include "mathfunc.h" #include "mathfunc.h"
#include "primitive.h" #include "primitive.h"
#include "symmetry.h" #include "symmetry.h"
@ -94,6 +94,10 @@ Primitive * prm_alloc_primitive(const int size)
} }
} }
for (i = 0; i < size; i++) {
primitive->mapping_table[i] = -1;
}
return primitive; return primitive;
} }
@ -107,9 +111,9 @@ void prm_free_primitive(Primitive * primitive)
if (primitive->cell != NULL) { if (primitive->cell != NULL) {
cel_free_cell(primitive->cell); cel_free_cell(primitive->cell);
primitive->cell = NULL;
} }
free(primitive); free(primitive);
primitive = NULL;
} }
} }
@ -139,6 +143,7 @@ static Primitive * get_primitive(SPGCONST Cell * cell, const double symprec)
tolerance = symprec; tolerance = symprec;
for (attempt = 0; attempt < 100; attempt++) { for (attempt = 0; attempt < 100; attempt++) {
debug_print("get_primitive (attempt = %d):\n", attempt);
if ((pure_trans = sym_get_pure_translation(cell, tolerance)) == NULL) { if ((pure_trans = sym_get_pure_translation(cell, tolerance)) == NULL) {
goto cont; goto cont;
} }
@ -161,6 +166,7 @@ static Primitive * get_primitive(SPGCONST Cell * cell, const double symprec)
} }
mat_free_VecDBL(pure_trans); mat_free_VecDBL(pure_trans);
pure_trans = NULL;
cont: cont:
tolerance *= REDUCE_RATE; tolerance *= REDUCE_RATE;
@ -176,6 +182,7 @@ static Primitive * get_primitive(SPGCONST Cell * cell, const double symprec)
mat_inverse_matrix_d3(inv_lat, cell->lattice, 0); mat_inverse_matrix_d3(inv_lat, cell->lattice, 0);
mat_multiply_matrix_d3(primitive->t_mat, primitive->cell->lattice, inv_lat); mat_multiply_matrix_d3(primitive->t_mat, primitive->cell->lattice, inv_lat);
mat_free_VecDBL(pure_trans); mat_free_VecDBL(pure_trans);
pure_trans = NULL;
return primitive; return primitive;
} }
@ -188,20 +195,18 @@ static Cell * get_cell_with_smallest_lattice(SPGCONST Cell * cell,
Cell * smallest_cell; Cell * smallest_cell;
debug_print("get_cell_with_smallest_lattice:\n"); debug_print("get_cell_with_smallest_lattice:\n");
smallest_cell = NULL; smallest_cell = NULL;
if (!lat_smallest_lattice_vector(min_lat, if (!del_delaunay_reduce(min_lat, cell->lattice, symprec)) {
cell->lattice, return NULL;
symprec)) {
goto err;
} }
mat_inverse_matrix_d3(inv_lat, min_lat, 0); mat_inverse_matrix_d3(inv_lat, min_lat, 0);
mat_multiply_matrix_d3(trans_mat, inv_lat, cell->lattice); mat_multiply_matrix_d3(trans_mat, inv_lat, cell->lattice);
if ((smallest_cell = cel_alloc_cell(cell->size)) == NULL) { if ((smallest_cell = cel_alloc_cell(cell->size)) == NULL) {
goto err; return NULL;
} }
mat_copy_matrix_d3(smallest_cell->lattice, min_lat); mat_copy_matrix_d3(smallest_cell->lattice, min_lat);
@ -210,14 +215,11 @@ static Cell * get_cell_with_smallest_lattice(SPGCONST Cell * cell,
mat_multiply_matrix_vector_d3(smallest_cell->position[i], mat_multiply_matrix_vector_d3(smallest_cell->position[i],
trans_mat, cell->position[i]); trans_mat, cell->position[i]);
for (j = 0; j < 3; j++) { for (j = 0; j < 3; j++) {
cell->position[i][j] = mat_Dmod1(cell->position[i][j]); smallest_cell->position[i][j] = mat_Dmod1(smallest_cell->position[i][j]);
} }
} }
return smallest_cell; return smallest_cell;
err:
return NULL;
} }
/* Return NULL if failed */ /* Return NULL if failed */
@ -230,7 +232,7 @@ static Cell * get_primitive_cell(int * mapping_table,
double prim_lat[3][3], smallest_lat[3][3]; double prim_lat[3][3], smallest_lat[3][3];
Cell * primitive_cell; Cell * primitive_cell;
debug_print("get_primitive:\n"); debug_print("get_primitive_cell:\n");
primitive_cell = NULL; primitive_cell = NULL;
@ -245,7 +247,7 @@ static Cell * get_primitive_cell(int * mapping_table,
goto not_found; goto not_found;
} }
if (! lat_smallest_lattice_vector(smallest_lat, prim_lat, symprec)) { if (! del_delaunay_reduce(smallest_lat, prim_lat, symprec)) {
goto not_found; goto not_found;
} }
@ -290,12 +292,13 @@ static int get_primitive_lattice_vectors_iterative(double prim_lattice[3][3],
for (i = 0; i < pure_trans->size; i++) { for (i = 0; i < pure_trans->size; i++) {
mat_copy_vector_d3(pure_trans_reduced->vec[i], pure_trans->vec[i]); mat_copy_vector_d3(pure_trans_reduced->vec[i], pure_trans->vec[i]);
} }
for (attempt = 0; attempt < 100; attempt++) { for (attempt = 0; attempt < 100; attempt++) {
multi = pure_trans_reduced->size; multi = pure_trans_reduced->size;
if ((vectors = get_translation_candidates(pure_trans_reduced)) == NULL) { if ((vectors = get_translation_candidates(pure_trans_reduced)) == NULL) {
mat_free_VecDBL(pure_trans_reduced); mat_free_VecDBL(pure_trans_reduced);
pure_trans_reduced = NULL;
goto fail; goto fail;
} }
@ -306,7 +309,9 @@ static int get_primitive_lattice_vectors_iterative(double prim_lattice[3][3],
tolerance)) { tolerance)) {
mat_free_VecDBL(vectors); mat_free_VecDBL(vectors);
vectors = NULL;
mat_free_VecDBL(pure_trans_reduced); mat_free_VecDBL(pure_trans_reduced);
pure_trans_reduced = NULL;
goto found; goto found;
@ -314,7 +319,9 @@ static int get_primitive_lattice_vectors_iterative(double prim_lattice[3][3],
if ((tmp_vec = mat_alloc_VecDBL(multi)) == NULL) { if ((tmp_vec = mat_alloc_VecDBL(multi)) == NULL) {
mat_free_VecDBL(vectors); mat_free_VecDBL(vectors);
vectors = NULL;
mat_free_VecDBL(pure_trans_reduced); mat_free_VecDBL(pure_trans_reduced);
pure_trans_reduced = NULL;
goto fail; goto fail;
} }
@ -322,19 +329,22 @@ static int get_primitive_lattice_vectors_iterative(double prim_lattice[3][3],
mat_copy_vector_d3(tmp_vec->vec[i], pure_trans_reduced->vec[i]); mat_copy_vector_d3(tmp_vec->vec[i], pure_trans_reduced->vec[i]);
} }
mat_free_VecDBL(pure_trans_reduced); mat_free_VecDBL(pure_trans_reduced);
pure_trans_reduced = NULL;
pure_trans_reduced = sym_reduce_pure_translation(cell, pure_trans_reduced = sym_reduce_pure_translation(cell,
tmp_vec, tmp_vec,
tolerance); tolerance);
mat_free_VecDBL(tmp_vec); mat_free_VecDBL(tmp_vec);
tmp_vec = NULL;
mat_free_VecDBL(vectors); mat_free_VecDBL(vectors);
vectors = NULL;
if (pure_trans_reduced == NULL) { if (pure_trans_reduced == NULL) {
goto fail; goto fail;
} }
warning_print("Tolerance is reduced to %f (%d), size = %d\n", warning_print("Tolerance is reduced to %f (%d), num_pure_trans = %d\n",
tolerance, attempt, pure_trans_reduced->size); tolerance, attempt, pure_trans_reduced->size);
tolerance *= REDUCE_RATE; tolerance *= REDUCE_RATE;
@ -342,6 +352,7 @@ static int get_primitive_lattice_vectors_iterative(double prim_lattice[3][3],
} }
mat_free_VecDBL(pure_trans_reduced); mat_free_VecDBL(pure_trans_reduced);
pure_trans_reduced = NULL;
fail: fail:
return 0; return 0;
@ -417,7 +428,7 @@ static int get_primitive_lattice_vectors(double prim_lattice[3][3],
} }
mat_multiply_matrix_d3(prim_lattice, cell->lattice, relative_lattice); mat_multiply_matrix_d3(prim_lattice, cell->lattice, relative_lattice);
return 1; return 1;
} }
static VecDBL * get_translation_candidates(const VecDBL * pure_trans) static VecDBL * get_translation_candidates(const VecDBL * pure_trans)
@ -432,7 +443,7 @@ static VecDBL * get_translation_candidates(const VecDBL * pure_trans)
return NULL; return NULL;
} }
/* store pure translations in original cell */ /* store pure translations in original cell */
/* as trial primitive lattice vectors */ /* as trial primitive lattice vectors */
for (i = 0; i < multi - 1; i++) { for (i = 0; i < multi - 1; i++) {
mat_copy_vector_d3(vectors->vec[i], pure_trans->vec[i + 1]); mat_copy_vector_d3(vectors->vec[i], pure_trans->vec[i + 1]);
@ -452,4 +463,3 @@ static VecDBL * get_translation_candidates(const VecDBL * pure_trans)
return vectors; return vectors;
} }

84
c/spglib/refinement.c
View File

@ -116,12 +116,13 @@ static int search_equivalent_atom(const int atom_index,
SPGCONST Cell * cell, SPGCONST Cell * cell,
const Symmetry *symmetry, const Symmetry *symmetry,
const double symprec); const double symprec);
static Symmetry * reduce_symmetry_in_frame(const int frame[3], static Symmetry *
SPGCONST Symmetry *prim_sym, recover_symmetry_in_original_cell(const int frame[3],
SPGCONST int t_mat[3][3], SPGCONST Symmetry *prim_sym,
SPGCONST double lattice[3][3], SPGCONST int t_mat[3][3],
const int multiplicity, SPGCONST double lattice[3][3],
const double symprec); const int multiplicity,
const double symprec);
static VecDBL * get_lattice_translations(const int frame[3], static VecDBL * get_lattice_translations(const int frame[3],
SPGCONST double inv_tmat[3][3]); SPGCONST double inv_tmat[3][3]);
static VecDBL * static VecDBL *
@ -212,7 +213,7 @@ Cell * get_Wyckoff_positions(int * wyckoffs,
wyckoffs_bravais = NULL; wyckoffs_bravais = NULL;
return NULL; return NULL;
} }
if ((bravais = get_bravais_exact_positions_and_lattice if ((bravais = get_bravais_exact_positions_and_lattice
(wyckoffs_bravais, (wyckoffs_bravais,
equiv_atoms_bravais, equiv_atoms_bravais,
@ -225,7 +226,7 @@ Cell * get_Wyckoff_positions(int * wyckoffs,
for (i = 0; i < cell->size; i++) { for (i = 0; i < cell->size; i++) {
wyckoffs[i] = wyckoffs_bravais[mapping_table[i]]; wyckoffs[i] = wyckoffs_bravais[mapping_table[i]];
} }
spgdb_get_operation_index(operation_index, spacegroup->hall_number); spgdb_get_operation_index(operation_index, spacegroup->hall_number);
num_prim_sym = operation_index[0] / (bravais->size / primitive->size); num_prim_sym = operation_index[0] / (bravais->size / primitive->size);
@ -247,12 +248,12 @@ Cell * get_Wyckoff_positions(int * wyckoffs,
} }
} }
ret: ret:
free(equiv_atoms_bravais); free(equiv_atoms_bravais);
equiv_atoms_bravais = NULL; equiv_atoms_bravais = NULL;
free(wyckoffs_bravais); free(wyckoffs_bravais);
wyckoffs_bravais = NULL; wyckoffs_bravais = NULL;
return bravais; return bravais;
} }
@ -323,6 +324,7 @@ get_bravais_exact_positions_and_lattice(int * wyckoffs,
spacegroup->hall_number, spacegroup->hall_number,
symprec)) == NULL) { symprec)) == NULL) {
sym_free_symmetry(conv_sym); sym_free_symmetry(conv_sym);
conv_sym = NULL;
goto err; goto err;
} }
@ -338,13 +340,16 @@ get_bravais_exact_positions_and_lattice(int * wyckoffs,
equiv_atoms_prim); equiv_atoms_prim);
mat_free_VecDBL(exact_positions); mat_free_VecDBL(exact_positions);
exact_positions = NULL;
sym_free_symmetry(conv_sym); sym_free_symmetry(conv_sym);
conv_sym = NULL;
err: err:
free(wyckoffs_prim); free(wyckoffs_prim);
wyckoffs_prim = NULL; wyckoffs_prim = NULL;
free(equiv_atoms_prim); free(equiv_atoms_prim);
equiv_atoms_prim = NULL; equiv_atoms_prim = NULL;
cel_free_cell(conv_prim); cel_free_cell(conv_prim);
conv_prim = NULL;
return bravais; return bravais;
} }
@ -379,7 +384,7 @@ static Cell * expand_positions(int * wyckoffs,
conv_prim->position[j]); conv_prim->position[j]);
for (k = 0; k < 3; k++) { for (k = 0; k < 3; k++) {
bravais->position[num_atom][k] += conv_sym->trans[i][k]; bravais->position[num_atom][k] += conv_sym->trans[i][k];
bravais->position[num_atom][k] = bravais->position[num_atom][k] =
mat_Dmod1(bravais->position[num_atom][k]); mat_Dmod1(bravais->position[num_atom][k]);
} }
wyckoffs[num_atom] = wyckoffs_prim[j]; wyckoffs[num_atom] = wyckoffs_prim[j];
@ -398,7 +403,7 @@ static Cell * expand_positions(int * wyckoffs,
static int get_number_of_pure_translation(SPGCONST Symmetry * conv_sym) static int get_number_of_pure_translation(SPGCONST Symmetry * conv_sym)
{ {
int i, num_pure_trans = 0; int i, num_pure_trans = 0;
for (i = 0; i < conv_sym->size; i++) { for (i = 0; i < conv_sym->size; i++) {
if (mat_check_identity_matrix_i3(identity, conv_sym->rot[i])) { if (mat_check_identity_matrix_i3(identity, conv_sym->rot[i])) {
num_pure_trans++; num_pure_trans++;
@ -423,7 +428,7 @@ static Cell * get_conventional_primitive(SPGCONST Spacegroup * spacegroup,
mat_inverse_matrix_d3(inv_brv, spacegroup->bravais_lattice, 0); mat_inverse_matrix_d3(inv_brv, spacegroup->bravais_lattice, 0);
mat_multiply_matrix_d3(trans_mat, inv_brv, primitive->lattice); mat_multiply_matrix_d3(trans_mat, inv_brv, primitive->lattice);
for (i = 0; i < primitive->size; i++) { for (i = 0; i < primitive->size; i++) {
conv_prim->types[i] = primitive->types[i]; conv_prim->types[i] = primitive->types[i];
mat_multiply_matrix_vector_d3(conv_prim->position[i], mat_multiply_matrix_vector_d3(conv_prim->position[i],
@ -446,7 +451,7 @@ static void get_conventional_lattice(double lattice[3][3],
Pointgroup pointgroup; Pointgroup pointgroup;
pointgroup = ptg_get_pointgroup(spacegroup->pointgroup_number); pointgroup = ptg_get_pointgroup(spacegroup->pointgroup_number);
for (i = 0; i < 3; i++) { for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) { for (j = 0; j < 3; j++) {
lattice[i][j] = 0; lattice[i][j] = 0;
@ -539,7 +544,7 @@ static void set_monocli(double lattice[3][3],
debug_print("beta %f\n", beta); debug_print("beta %f\n", beta);
debug_print_matrix_d3(lattice); debug_print_matrix_d3(lattice);
} }
static void set_ortho(double lattice[3][3], static void set_ortho(double lattice[3][3],
SPGCONST double metric[3][3]) SPGCONST double metric[3][3])
{ {
@ -663,24 +668,27 @@ get_refined_symmetry_operations(SPGCONST Cell * cell,
if ((prim_sym = get_primitive_db_symmetry(t_mat, conv_sym)) == NULL) { if ((prim_sym = get_primitive_db_symmetry(t_mat, conv_sym)) == NULL) {
sym_free_symmetry(conv_sym); sym_free_symmetry(conv_sym);
conv_sym = NULL;
return NULL; return NULL;
} }
sym_free_symmetry(conv_sym); sym_free_symmetry(conv_sym);
conv_sym = NULL;
/* Input cell symmetry from primitive symmetry */ /* Input cell symmetry from primitive symmetry */
mat_multiply_matrix_d3(t_mat, inv_prim_lat, cell->lattice); mat_multiply_matrix_d3(t_mat, inv_prim_lat, cell->lattice);
mat_cast_matrix_3d_to_3i(t_mat_int, t_mat); mat_cast_matrix_3d_to_3i(t_mat_int, t_mat);
get_surrounding_frame(frame, t_mat_int); get_surrounding_frame(frame, t_mat_int);
symmetry = reduce_symmetry_in_frame(frame, symmetry = recover_symmetry_in_original_cell(frame,
prim_sym, prim_sym,
t_mat_int, t_mat_int,
cell->lattice, cell->lattice,
cell->size / primitive->size, cell->size / primitive->size,
symprec); symprec);
sym_free_symmetry(prim_sym); sym_free_symmetry(prim_sym);
prim_sym = NULL;
return symmetry; return symmetry;
} }
@ -802,13 +810,14 @@ static Symmetry * get_primitive_db_symmetry(SPGCONST double t_mat[3][3],
r_prim = NULL; r_prim = NULL;
t_prim = NULL; t_prim = NULL;
prim_sym = NULL; prim_sym = NULL;
if ((r_prim = mat_alloc_MatINT(conv_sym->size)) == NULL) { if ((r_prim = mat_alloc_MatINT(conv_sym->size)) == NULL) {
return NULL; return NULL;
} }
if ((t_prim = mat_alloc_VecDBL(conv_sym->size)) == NULL) { if ((t_prim = mat_alloc_VecDBL(conv_sym->size)) == NULL) {
mat_free_MatINT(r_prim); mat_free_MatINT(r_prim);
r_prim = NULL;
return NULL; return NULL;
} }
@ -850,7 +859,9 @@ static Symmetry * get_primitive_db_symmetry(SPGCONST double t_mat[3][3],
ret: ret:
mat_free_MatINT(r_prim); mat_free_MatINT(r_prim);
r_prim = NULL;
mat_free_VecDBL(t_prim); mat_free_VecDBL(t_prim);
t_prim = NULL;
return prim_sym; return prim_sym;
} }
@ -908,12 +919,13 @@ static void get_corners(int corners[3][8],
} }
} }
static Symmetry * reduce_symmetry_in_frame(const int frame[3], static Symmetry *
SPGCONST Symmetry *prim_sym, recover_symmetry_in_original_cell(const int frame[3],
SPGCONST int t_mat[3][3], SPGCONST Symmetry *prim_sym,
SPGCONST double lattice[3][3], SPGCONST int t_mat[3][3],
const int multiplicity, SPGCONST double lattice[3][3],
const double symprec) const int multiplicity,
const double symprec)
{ {
Symmetry *symmetry, *t_sym; Symmetry *symmetry, *t_sym;
double inv_tmat[3][3], tmp_mat[3][3]; double inv_tmat[3][3], tmp_mat[3][3];
@ -935,6 +947,7 @@ static Symmetry * reduce_symmetry_in_frame(const int frame[3],
lattice_trans, lattice_trans,
symprec)) == NULL) { symprec)) == NULL) {
mat_free_VecDBL(lattice_trans); mat_free_VecDBL(lattice_trans);
lattice_trans = NULL;
return NULL; return NULL;
} }
@ -944,7 +957,9 @@ static Symmetry * reduce_symmetry_in_frame(const int frame[3],
prim_sym, prim_sym,
symprec)) == NULL) { symprec)) == NULL) {
mat_free_VecDBL(pure_trans); mat_free_VecDBL(pure_trans);
pure_trans = NULL;
mat_free_VecDBL(lattice_trans); mat_free_VecDBL(lattice_trans);
lattice_trans = NULL;
return NULL; return NULL;
} }
@ -953,8 +968,11 @@ static Symmetry * reduce_symmetry_in_frame(const int frame[3],
} }
mat_free_VecDBL(lattice_trans); mat_free_VecDBL(lattice_trans);
lattice_trans = NULL;
mat_free_VecDBL(pure_trans); mat_free_VecDBL(pure_trans);
pure_trans = NULL;
sym_free_symmetry(t_sym); sym_free_symmetry(t_sym);
t_sym = NULL;
return symmetry; return symmetry;
} }
@ -987,7 +1005,7 @@ static VecDBL * get_lattice_translations(const int frame[3],
lattice_trans->vec[num_trans]); lattice_trans->vec[num_trans]);
for (l = 0; l < 3; l++) { for (l = 0; l < 3; l++) {
lattice_trans->vec[num_trans][l] = lattice_trans->vec[num_trans][l] =
mat_Dmod1(lattice_trans->vec[num_trans][l]); mat_Dmod1(lattice_trans->vec[num_trans][l]);
} }
num_trans++; num_trans++;
} }
@ -1004,7 +1022,7 @@ remove_overlapping_lattice_points(SPGCONST double lattice[3][3],
{ {
int i, j, is_found, num_pure_trans; int i, j, is_found, num_pure_trans;
VecDBL *pure_trans, *t; VecDBL *pure_trans, *t;
pure_trans = NULL; pure_trans = NULL;
t = NULL; t = NULL;
@ -1030,6 +1048,7 @@ remove_overlapping_lattice_points(SPGCONST double lattice[3][3],
if ((pure_trans = mat_alloc_VecDBL(num_pure_trans)) == NULL) { if ((pure_trans = mat_alloc_VecDBL(num_pure_trans)) == NULL) {
mat_free_VecDBL(t); mat_free_VecDBL(t);
t = NULL;
return NULL; return NULL;
} }
@ -1037,6 +1056,7 @@ remove_overlapping_lattice_points(SPGCONST double lattice[3][3],
mat_copy_vector_d3(pure_trans->vec[i], t->vec[i]); mat_copy_vector_d3(pure_trans->vec[i], t->vec[i]);
} }
mat_free_VecDBL(t); mat_free_VecDBL(t);
t = NULL;
return pure_trans; return pure_trans;
} }
@ -1048,7 +1068,7 @@ get_symmetry_in_original_cell(SPGCONST int t_mat[3][3],
SPGCONST double lattice[3][3], SPGCONST double lattice[3][3],
SPGCONST Symmetry *prim_sym, SPGCONST Symmetry *prim_sym,
const double symprec) const double symprec)
{ {
int i, size_sym_orig; int i, size_sym_orig;
double tmp_rot_d[3][3], tmp_lat_d[3][3], tmp_lat_i[3][3], tmp_mat[3][3]; double tmp_rot_d[3][3], tmp_lat_d[3][3], tmp_lat_i[3][3], tmp_mat[3][3];
int tmp_rot_i[3][3]; int tmp_rot_i[3][3];
@ -1091,6 +1111,7 @@ get_symmetry_in_original_cell(SPGCONST int t_mat[3][3],
if ((t_red_sym = sym_alloc_symmetry(size_sym_orig)) == NULL) { if ((t_red_sym = sym_alloc_symmetry(size_sym_orig)) == NULL) {
sym_free_symmetry(t_sym); sym_free_symmetry(t_sym);
t_sym = NULL;
return NULL; return NULL;
} }
@ -1100,6 +1121,7 @@ get_symmetry_in_original_cell(SPGCONST int t_mat[3][3],
} }
sym_free_symmetry(t_sym); sym_free_symmetry(t_sym);
t_sym = NULL;
t_sym = t_red_sym; t_sym = t_red_sym;
t_red_sym = NULL; t_red_sym = NULL;

1
c/spglib/site_symmetry.c
View File

@ -288,5 +288,6 @@ static int get_Wyckoff_notation(double position[3],
end: end:
mat_free_VecDBL(pos_rot); mat_free_VecDBL(pos_rot);
pos_rot = NULL;
return wyckoff_letter; return wyckoff_letter;
} }

71
c/spglib/spacegroup.c
View File

@ -36,8 +36,8 @@
#include <stdlib.h> #include <stdlib.h>
#include <string.h> #include <string.h>
#include "cell.h" #include "cell.h"
#include "delaunay.h"
#include "hall_symbol.h" #include "hall_symbol.h"
#include "lattice.h"
#include "mathfunc.h" #include "mathfunc.h"
#include "niggli.h" #include "niggli.h"
#include "pointgroup.h" #include "pointgroup.h"
@ -316,11 +316,12 @@ Primitive * spa_get_spacegroup(Spacegroup * spacegroup,
primitive->tolerance); primitive->tolerance);
if (spacegroup->number > 0) { if (spacegroup->number > 0) {
break; break;
} else {
prm_free_primitive(primitive);
primitive = NULL;
} }
prm_free_primitive(primitive); cont:
cont:
warning_print("spglib: Attempt %d tolerance = %f failed.", warning_print("spglib: Attempt %d tolerance = %f failed.",
attempt, tolerance); attempt, tolerance);
warning_print(" (line %d, %s).\n", __LINE__, __FILE__); warning_print(" (line %d, %s).\n", __LINE__, __FILE__);
@ -343,13 +344,13 @@ Spacegroup spa_get_spacegroup_with_hall_number(SPGCONST Primitive * primitive,
int num_candidates; int num_candidates;
int candidate[1]; int candidate[1];
Spacegroup spacegroup; Spacegroup spacegroup;
spacegroup.number = 0; spacegroup.number = 0;
if (hall_number < 1 || hall_number > 530) { if (hall_number < 1 || hall_number > 530) {
goto err; goto err;
} }
num_candidates = 1; num_candidates = 1;
candidate[0] = hall_number; candidate[0] = hall_number;
spacegroup = search_spacegroup(primitive->cell, spacegroup = search_spacegroup(primitive->cell,
@ -434,6 +435,7 @@ static Spacegroup search_spacegroup(SPGCONST Cell * primitive,
double origin_shift[3]; double origin_shift[3];
Spacegroup spacegroup; Spacegroup spacegroup;
Symmetry *symmetry; Symmetry *symmetry;
PointSymmetry pointsym;
debug_print("search_spacegroup (tolerance = %f):\n", symprec); debug_print("search_spacegroup (tolerance = %f):\n", symprec);
@ -445,6 +447,15 @@ static Spacegroup search_spacegroup(SPGCONST Cell * primitive,
goto ret; goto ret;
} }
pointsym = ptg_get_pointsymmetry(symmetry->rot, symmetry->size);
if (pointsym.size < symmetry->size) {
warning_print("spglib: Point symmetry of primitive cell is broken. ");
warning_print("(line %d, %s).\n", __LINE__, __FILE__);
sym_free_symmetry(symmetry);
symmetry = NULL;
goto ret;
}
hall_number = iterative_search_hall_number(origin_shift, hall_number = iterative_search_hall_number(origin_shift,
conv_lattice, conv_lattice,
candidates, candidates,
@ -453,6 +464,7 @@ static Spacegroup search_spacegroup(SPGCONST Cell * primitive,
symmetry, symmetry,
symprec); symprec);
sym_free_symmetry(symmetry); sym_free_symmetry(symmetry);
symmetry = NULL;
spacegroup = get_spacegroup(hall_number, origin_shift, conv_lattice); spacegroup = get_spacegroup(hall_number, origin_shift, conv_lattice);
ret: ret:
@ -466,7 +478,7 @@ static Spacegroup get_spacegroup(const int hall_number,
{ {
Spacegroup spacegroup; Spacegroup spacegroup;
SpacegroupType spacegroup_type; SpacegroupType spacegroup_type;
spacegroup.number = 0; spacegroup.number = 0;
spacegroup_type = spgdb_get_spacegroup_type(hall_number); spacegroup_type = spgdb_get_spacegroup_type(hall_number);
@ -540,6 +552,7 @@ static int iterative_search_hall_number(double origin_shift[3],
sym_reduced, sym_reduced,
symprec); symprec);
sym_free_symmetry(sym_reduced); sym_free_symmetry(sym_reduced);
sym_reduced = NULL;
if (hall_number > 0) { if (hall_number > 0) {
break; break;
} }
@ -627,6 +640,7 @@ static int search_hall_number(double origin_shift[3],
} }
sym_free_symmetry(conv_symmetry); sym_free_symmetry(conv_symmetry);
conv_symmetry = NULL;
return hall_number; return hall_number;
@ -683,10 +697,10 @@ static int change_basis_monocli(int int_transform_mat[3][3],
{ {
double smallest_lattice[3][3], inv_lattice[3][3], transform_mat[3][3]; double smallest_lattice[3][3], inv_lattice[3][3], transform_mat[3][3];
if (! lat_smallest_lattice_vector_2D(smallest_lattice, if (! del_delaunay_reduce_2D(smallest_lattice,
conv_lattice, conv_lattice,
1, /* unique axis of b */ 1, /* unique axis of b */
symprec)) { symprec)) {
return 0; return 0;
} }
@ -707,7 +721,7 @@ get_initial_conventional_symmetry(const Centering centering,
debug_print("get_initial_conventional_symmetry\n"); debug_print("get_initial_conventional_symmetry\n");
conv_symmetry = NULL; conv_symmetry = NULL;
if (centering == R_CENTER) { if (centering == R_CENTER) {
/* hP for rhombohedral */ /* hP for rhombohedral */
conv_symmetry = get_conventional_symmetry(transform_mat, conv_symmetry = get_conventional_symmetry(transform_mat,
@ -736,7 +750,7 @@ static int match_hall_symbol_db(double origin_shift[3],
SpacegroupType spacegroup_type; SpacegroupType spacegroup_type;
Symmetry * changed_symmetry; Symmetry * changed_symmetry;
double changed_lattice[3][3], inv_lattice[3][3], transform_mat[3][3]; double changed_lattice[3][3], inv_lattice[3][3], transform_mat[3][3];
changed_symmetry = NULL; changed_symmetry = NULL;
spacegroup_type = spgdb_get_spacegroup_type(hall_number); spacegroup_type = spgdb_get_spacegroup_type(hall_number);
@ -757,7 +771,7 @@ static int match_hall_symbol_db(double origin_shift[3],
symmetry, symmetry,
symprec)) {return 1;} symprec)) {return 1;}
break; break;
case ORTHO: case ORTHO:
if (spacegroup_type.number == 48 || if (spacegroup_type.number == 48 ||
spacegroup_type.number == 50 || spacegroup_type.number == 50 ||
@ -766,7 +780,7 @@ static int match_hall_symbol_db(double origin_shift[3],
spacegroup_type.number == 70) { /* uncount origin shift */ spacegroup_type.number == 70) { /* uncount origin shift */
num_hall_types /= 2; num_hall_types /= 2;
} }
if (num_hall_types == 1) { if (num_hall_types == 1) {
if (match_hall_symbol_db_ortho(origin_shift, if (match_hall_symbol_db_ortho(origin_shift,
lattice, lattice,
@ -815,7 +829,8 @@ static int match_hall_symbol_db(double origin_shift[3],
changed_symmetry, changed_symmetry,
2, 2,
symprec); symprec);
sym_free_symmetry(changed_symmetry); sym_free_symmetry(changed_symmetry);
changed_symmetry = NULL;
if (is_found) { if (is_found) {
mat_copy_matrix_d3(lattice, changed_lattice); mat_copy_matrix_d3(lattice, changed_lattice);
return 1; return 1;
@ -843,7 +858,7 @@ static int match_hall_symbol_db(double origin_shift[3],
centering, centering,
symmetry, symmetry,
symprec)) {return 1;} symprec)) {return 1;}
if (hall_number == 501) { /* Try another basis for No.205 */ if (hall_number == 501) { /* Try another basis for No.205 */
mat_multiply_matrix_d3(changed_lattice, mat_multiply_matrix_d3(changed_lattice,
lattice, lattice,
@ -861,13 +876,14 @@ static int match_hall_symbol_db(double origin_shift[3],
changed_symmetry, changed_symmetry,
symprec); symprec);
sym_free_symmetry(changed_symmetry); sym_free_symmetry(changed_symmetry);
changed_symmetry = NULL;
if (is_found) { if (is_found) {
mat_copy_matrix_d3(lattice, changed_lattice); mat_copy_matrix_d3(lattice, changed_lattice);
return 1; return 1;
} }
} }
break; break;
case TRIGO: case TRIGO:
if (centering == R_CENTER) { if (centering == R_CENTER) {
if (hall_number == 433 || if (hall_number == 433 ||
@ -893,6 +909,7 @@ static int match_hall_symbol_db(double origin_shift[3],
changed_symmetry, changed_symmetry,
symprec); symprec);
sym_free_symmetry(changed_symmetry); sym_free_symmetry(changed_symmetry);
changed_symmetry = NULL;
if (is_found) { if (is_found) {
mat_copy_matrix_d3(lattice, changed_lattice); mat_copy_matrix_d3(lattice, changed_lattice);
return 1; return 1;
@ -932,7 +949,7 @@ static int match_hall_symbol_db_monocli(double origin_shift[3],
double changed_lattice[3][3]; double changed_lattice[3][3];
changed_symmetry = NULL; changed_symmetry = NULL;
for (i = 0; i < 18; i++) { for (i = 0; i < 18; i++) {
if (centering == C_FACE) { if (centering == C_FACE) {
changed_centering = change_of_centering_monocli[i]; changed_centering = change_of_centering_monocli[i];
@ -970,6 +987,7 @@ static int match_hall_symbol_db_monocli(double origin_shift[3],
changed_symmetry, changed_symmetry,
symprec); symprec);
sym_free_symmetry(changed_symmetry); sym_free_symmetry(changed_symmetry);
changed_symmetry = NULL;
if (is_found) { if (is_found) {
mat_copy_matrix_d3(lattice, changed_lattice); mat_copy_matrix_d3(lattice, changed_lattice);
return 1; return 1;
@ -996,14 +1014,14 @@ static int match_hall_symbol_db_ortho(double origin_shift[3],
double changed_lattice[3][3]; double changed_lattice[3][3];
changed_symmetry = NULL; changed_symmetry = NULL;
for (i = 0; i < 6; i++) { for (i = 0; i < 6; i++) {
if (centering == C_FACE) { if (centering == C_FACE) {
changed_centering = change_of_centering_ortho[i]; changed_centering = change_of_centering_ortho[i];
} else { } else {
changed_centering = centering; changed_centering = centering;
} }
mat_multiply_matrix_d3(changed_lattice, mat_multiply_matrix_d3(changed_lattice,
lattice, lattice,
change_of_basis_ortho[i]); change_of_basis_ortho[i]);
@ -1034,7 +1052,7 @@ static int match_hall_symbol_db_ortho(double origin_shift[3],
} }
if (norms[0] > norms[1] || norms[1] > norms[2]) {continue;} if (norms[0] > norms[1] || norms[1] > norms[2]) {continue;}
} }
if ((changed_symmetry = get_conventional_symmetry(change_of_basis_ortho[i], if ((changed_symmetry = get_conventional_symmetry(change_of_basis_ortho[i],
PRIMITIVE, PRIMITIVE,
symmetry)) == NULL) { symmetry)) == NULL) {
@ -1048,6 +1066,7 @@ static int match_hall_symbol_db_ortho(double origin_shift[3],
changed_symmetry, changed_symmetry,
symprec); symprec);
sym_free_symmetry(changed_symmetry); sym_free_symmetry(changed_symmetry);
changed_symmetry = NULL;
if (is_found) { if (is_found) {
mat_copy_matrix_d3(lattice, changed_lattice); mat_copy_matrix_d3(lattice, changed_lattice);
return 1; return 1;
@ -1135,7 +1154,7 @@ static Symmetry * get_conventional_symmetry(SPGCONST double transform_mat[3][3],
shift[1][2] = 2. / 3; shift[1][2] = 2. / 3;
multi = 3; multi = 3;
} }
if (centering == FACE) { if (centering == FACE) {
shift[0][0] = 0; shift[0][0] = 0;
shift[0][1] = 0.5; shift[0][1] = 0.5;
@ -1263,7 +1282,7 @@ static Centering get_base_center(SPGCONST int transform_mat[3][3])
/* A center */ /* A center */
for (i = 0; i < 3; i++) { for (i = 0; i < 3; i++) {
if (abs(transform_mat[i][0]) == 1 && if (abs(transform_mat[i][0]) == 1 &&
transform_mat[i][1] == 0 && transform_mat[i][1] == 0 &&
transform_mat[i][2] == 0) { transform_mat[i][2] == 0) {
centering = A_FACE; centering = A_FACE;
@ -1274,7 +1293,7 @@ static Centering get_base_center(SPGCONST int transform_mat[3][3])
/* B center */ /* B center */
for (i = 0; i < 3; i++) { for (i = 0; i < 3; i++) {
if (transform_mat[i][0] == 0 && if (transform_mat[i][0] == 0 &&
abs(transform_mat[i][1]) == 1 && abs(transform_mat[i][1]) == 1 &&
transform_mat[i][2] == 0) { transform_mat[i][2] == 0) {
centering = B_FACE; centering = B_FACE;
goto end; goto end;
@ -1283,7 +1302,7 @@ static Centering get_base_center(SPGCONST int transform_mat[3][3])
/* body center */ /* body center */
if (abs(transform_mat[0][0]) + if (abs(transform_mat[0][0]) +
abs(transform_mat[0][1]) + abs(transform_mat[0][1]) +
abs(transform_mat[0][2]) == 2) { abs(transform_mat[0][2]) == 2) {
centering = BODY; centering = BODY;
goto end; goto end;

20
c/spglib/spg_database.c
View File

@ -531,15 +531,15 @@ static const SpacegroupType spacegroup_types[] = {
{197, "T^3 ", "I 2 2 3 ", "I 2 3 ", "I 2 3 ", "I23 ", " ", BODY, 28 }, /* 491 */ {197, "T^3 ", "I 2 2 3 ", "I 2 3 ", "I 2 3 ", "I23 ", " ", BODY, 28 }, /* 491 */
{198, "T^4 ", "P 2ac 2ab 3 ", "P 2_1 3 ", "P 2_1 3 ", "P2_13 ", " ", PRIMITIVE, 28 }, /* 492 */ {198, "T^4 ", "P 2ac 2ab 3 ", "P 2_1 3 ", "P 2_1 3 ", "P2_13 ", " ", PRIMITIVE, 28 }, /* 492 */
{199, "T^5 ", "I 2b 2c 3 ", "I 2_1 3 ", "I 2_1 3 ", "I2_13 ", " ", BODY, 28 }, /* 493 */ {199, "T^5 ", "I 2b 2c 3 ", "I 2_1 3 ", "I 2_1 3 ", "I2_13 ", " ", BODY, 28 }, /* 493 */
{200, "Th^1 ", "-P 2 2 3 ", "P m 3 ", "P 2/m -3 ", "Pm3 ", " ", PRIMITIVE, 29 }, /* 494 */ {200, "Th^1 ", "-P 2 2 3 ", "P m -3 ", "P 2/m -3 ", "Pm-3 ", " ", PRIMITIVE, 29 }, /* 494 */
{201, "Th^2 ", "P 2 2 3 -1n ", "P n 3 ", "P 2/n -3 ", "Pn3 ", "1 ", PRIMITIVE, 29 }, /* 495 */ {201, "Th^2 ", "P 2 2 3 -1n ", "P n -3 ", "P 2/n -3 ", "Pn-3 ", "1 ", PRIMITIVE, 29 }, /* 495 */
{201, "Th^2 ", "-P 2ab 2bc 3 ", "P n 3 ", "P 2/n -3 ", "Pn3 ", "2 ", PRIMITIVE, 29 }, /* 496 */ {201, "Th^2 ", "-P 2ab 2bc 3 ", "P n -3 ", "P 2/n -3 ", "Pn-3 ", "2 ", PRIMITIVE, 29 }, /* 496 */
{202, "Th^3 ", "-F 2 2 3 ", "F m 3 ", "F 2/m -3 ", "Fm3 ", " ", FACE, 29 }, /* 497 */ {202, "Th^3 ", "-F 2 2 3 ", "F m -3 ", "F 2/m -3 ", "Fm-3 ", " ", FACE, 29 }, /* 497 */
{203, "Th^4 ", "F 2 2 3 -1d ", "F d 3 ", "F 2/d -3 ", "Fd3 ", "1 ", FACE, 29 }, /* 498 */ {203, "Th^4 ", "F 2 2 3 -1d ", "F d -3 ", "F 2/d -3 ", "Fd-3 ", "1 ", FACE, 29 }, /* 498 */
{203, "Th^4 ", "-F 2uv 2vw 3 ", "F d 3 ", "F 2/d -3 ", "Fd3 ", "2 ", FACE, 29 }, /* 499 */ {203, "Th^4 ", "-F 2uv 2vw 3 ", "F d -3 ", "F 2/d -3 ", "Fd-3 ", "2 ", FACE, 29 }, /* 499 */
{204, "Th^5 ", "-I 2 2 3 ", "I m 3 ", "I 2/m -3 ", "Im3 ", " ", BODY, 29 }, /* 500 */ {204, "Th^5 ", "-I 2 2 3 ", "I m -3 ", "I 2/m -3 ", "Im-3 ", " ", BODY, 29 }, /* 500 */
{205, "Th^6 ", "-P 2ac 2ab 3 ", "P a 3 ", "P 2_1/a -3 ", "Pa3 ", " ", PRIMITIVE, 29 }, /* 501 */ {205, "Th^6 ", "-P 2ac 2ab 3 ", "P a -3 ", "P 2_1/a -3 ", "Pa-3 ", " ", PRIMITIVE, 29 }, /* 501 */
{206, "Th^7 ", "-I 2b 2c 3 ", "I a 3 ", "I 2_1/a -3 ", "Ia3 ", " ", BODY, 29 }, /* 502 */ {206, "Th^7 ", "-I 2b 2c 3 ", "I a -3 ", "I 2_1/a -3 ", "Ia-3 ", " ", BODY, 29 }, /* 502 */
{207, "O^1 ", "P 4 2 3 ", "P 4 3 2 ", "P 4 3 2 ", "P432 ", " ", PRIMITIVE, 30 }, /* 503 */ {207, "O^1 ", "P 4 2 3 ", "P 4 3 2 ", "P 4 3 2 ", "P432 ", " ", PRIMITIVE, 30 }, /* 503 */
{208, "O^2 ", "P 4n 2 3 ", "P 4_2 3 2 ", "P 4_2 3 2 ", "P4_232 ", " ", PRIMITIVE, 30 }, /* 504 */ {208, "O^2 ", "P 4n 2 3 ", "P 4_2 3 2 ", "P 4_2 3 2 ", "P4_232 ", " ", PRIMITIVE, 30 }, /* 504 */
{209, "O^3 ", "F 4 2 3 ", "F 4 3 2 ", "F 4 3 2 ", "F432 ", " ", FACE, 30 }, /* 505 */ {209, "O^3 ", "F 4 2 3 ", "F 4 3 2 ", "F 4 3 2 ", "F432 ", " ", FACE, 30 }, /* 505 */
@ -8575,7 +8575,7 @@ SpacegroupType spgdb_get_spacegroup_type(const int hall_number)
spgtype.number = 0; spgtype.number = 0;
if (0 < hall_number || hall_number < 531) { if (0 < hall_number && hall_number < 531) {
spgtype = spacegroup_types[hall_number]; spgtype = spacegroup_types[hall_number];
} else { } else {
spgtype = spacegroup_types[0]; spgtype = spacegroup_types[0];

155
c/spglib/spglib.c
View File

@ -37,10 +37,11 @@
#include <string.h> #include <string.h>
#include "cell.h" #include "cell.h"
#include "debug.h" #include "debug.h"
#include "delaunay.h"
#include "kgrid.h" #include "kgrid.h"
#include "kpoint.h" #include "kpoint.h"
#include "lattice.h"
#include "mathfunc.h" #include "mathfunc.h"
#include "niggli.h"
#include "pointgroup.h" #include "pointgroup.h"
#include "spglib.h" #include "spglib.h"
#include "primitive.h" #include "primitive.h"
@ -467,14 +468,6 @@ int spgat_get_multiplicity(SPGCONST double lattice[3][3],
symprec); symprec);
} }
/* Return 0 if failed */
int spg_get_smallest_lattice(double smallest_lattice[3][3],
SPGCONST double lattice[3][3],
const double symprec)
{
return lat_smallest_lattice_vector(smallest_lattice, lattice, symprec);
}
/* Return 0 if failed */ /* Return 0 if failed */
int spg_get_international(char symbol[11], int spg_get_international(char symbol[11],
SPGCONST double lattice[3][3], SPGCONST double lattice[3][3],
@ -591,6 +584,7 @@ int spg_get_symmetry_from_database(int rotations[192][3][3],
size = symmetry->size; size = symmetry->size;
sym_free_symmetry(symmetry); sym_free_symmetry(symmetry);
symmetry = NULL;
return size; return size;
} }
@ -742,6 +736,24 @@ int spgat_refine_cell(double lattice[3][3],
symprec); symprec);
} }
int spg_delaunay_reduce(double lattice[3][3], const double symprec)
{
int i, j, succeeded;
double red_lattice[3][3];
succeeded = del_delaunay_reduce(red_lattice, lattice, symprec);
if (succeeded) {
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
lattice[i][j] = red_lattice[i][j];
}
}
}
return succeeded;
}
/*---------*/ /*---------*/
/* kpoints */ /* kpoints */
/*---------*/ /*---------*/
@ -832,6 +844,7 @@ int spg_get_grid_points_by_rotations(int rot_grid_points[],
mesh, mesh,
is_shift); is_shift);
mat_free_MatINT(rot); mat_free_MatINT(rot);
rot = NULL;
return 1; return 1;
} }
@ -863,6 +876,7 @@ int spg_get_BZ_grid_points_by_rotations(int rot_grid_points[],
is_shift, is_shift,
bz_map); bz_map);
mat_free_MatINT(rot); mat_free_MatINT(rot);
rot = NULL;
return 1; return 1;
} }
@ -882,64 +896,35 @@ int spg_relocate_BZ_grid_address(int bz_grid_address[][3],
is_shift); is_shift);
} }
/* /\*--------------------*\/ */ /*--------*/
/* /\* tetrahedron method *\/ */ /* Niggli */
/* /\*--------------------*\/ */ /*--------*/
/* void spg_get_neighboring_grid_points(int relative_grid_points[], */ /* Return 0 if failed */
/* const int grid_point, */ int spg_niggli_reduce(double lattice[3][3], const double symprec)
/* SPGCONST int relative_grid_address[][3], */ {
/* const int num_relative_grid_address, */ int i, j, succeeded;
/* const int mesh[3], */ double vals[9];
/* SPGCONST int bz_grid_address[][3], */
/* const int bz_map[]) */ for (i = 0; i < 3; i++) {
/* { */ for (j = 0; j < 3; j++) {
/* thm_get_neighboring_grid_points(relative_grid_points, */ vals[i * 3 + j] = lattice[i][j];
/* grid_point, */ }
/* relative_grid_address, */ }
/* num_relative_grid_address, */
/* mesh, */
/* bz_grid_address, */
/* bz_map); */
/* } */
/* void */ succeeded = niggli_reduce(vals, symprec);
/* spg_get_tetrahedra_relative_grid_address(int relative_grid_address[24][4][3], */
/* SPGCONST double rec_lattice[3][3]) */
/* { */
/* thm_get_relative_grid_address(relative_grid_address, rec_lattice); */
/* } */
/* void */ if (succeeded) {
/* spg_get_all_tetrahedra_relative_grid_address */ for (i = 0; i < 3; i++) {
/* (int relative_grid_address[4][24][4][3]) */ for (j = 0; j < 3; j++) {
/* { */ lattice[i][j] = vals[i * 3 + j];
/* thm_get_all_relative_grid_address(relative_grid_address); */ }
/* } */ }
}
return succeeded;
}
/* double */
/* spg_get_tetrahedra_integration_weight(const double omega, */
/* SPGCONST double tetrahedra_omegas[24][4], */
/* const char function) */
/* { */
/* return thm_get_integration_weight(omega, */
/* tetrahedra_omegas, */
/* function); */
/* } */
/* void */
/* spg_get_tetrahedra_integration_weight_at_omegas */
/* (double integration_weights[], */
/* const int num_omegas, */
/* const double omegas[], */
/* SPGCONST double tetrahedra_omegas[24][4], */
/* const char function) */
/* { */
/* thm_get_integration_weight_at_omegas(integration_weights, */
/* num_omegas, */
/* omegas, */
/* tetrahedra_omegas, */
/* function); */
/* } */
/*=======*/ /*=======*/
@ -974,6 +959,7 @@ static SpglibDataset * get_dataset(SPGCONST double lattice[3][3],
} }
dataset->spacegroup_number = 0; dataset->spacegroup_number = 0;
dataset->hall_number = 0;
strcpy(dataset->international_symbol, ""); strcpy(dataset->international_symbol, "");
strcpy(dataset->hall_symbol, ""); strcpy(dataset->hall_symbol, "");
strcpy(dataset->setting, ""); strcpy(dataset->setting, "");
@ -1000,10 +986,15 @@ static SpglibDataset * get_dataset(SPGCONST double lattice[3][3],
cel_set_cell(cell, lattice, position, types); cel_set_cell(cell, lattice, position, types);
primitive = spa_get_spacegroup(&spacegroup, cell, symprec); if ((primitive = spa_get_spacegroup(&spacegroup, cell, symprec)) == NULL) {
cel_free_cell(cell);
if ((spacegroup.number > 0) && (primitive != NULL)) { cell = NULL;
free(dataset);
dataset = NULL;
return NULL;
}
if (spacegroup.number > 0) {
/* With hall_number > 0, specific choice is searched. */ /* With hall_number > 0, specific choice is searched. */
if (hall_number > 0) { if (hall_number > 0) {
spacegroup_type = spgdb_get_spacegroup_type(hall_number); spacegroup_type = spgdb_get_spacegroup_type(hall_number);
@ -1031,14 +1022,18 @@ static SpglibDataset * get_dataset(SPGCONST double lattice[3][3],
} }
} }
cel_free_cell(cell);
prm_free_primitive(primitive); prm_free_primitive(primitive);
primitive = NULL;
cel_free_cell(cell);
cell = NULL;
return dataset; return dataset;
err: err:
cel_free_cell(cell); cel_free_cell(cell);
cell = NULL;
prm_free_primitive(primitive); prm_free_primitive(primitive);
primitive = NULL;
free(dataset); free(dataset);
dataset = NULL; dataset = NULL;
return NULL; return NULL;
@ -1150,7 +1145,9 @@ static int set_dataset(SpglibDataset * dataset,
} }
cel_free_cell(bravais); cel_free_cell(bravais);
bravais = NULL;
sym_free_symmetry(symmetry); sym_free_symmetry(symmetry);
symmetry = NULL;
dataset->pointgroup_number = spacegroup->pointgroup_number; dataset->pointgroup_number = spacegroup->pointgroup_number;
pointgroup = ptg_get_pointgroup(spacegroup->pointgroup_number); pointgroup = ptg_get_pointgroup(spacegroup->pointgroup_number);
@ -1165,6 +1162,7 @@ static int set_dataset(SpglibDataset * dataset,
} }
if (bravais != NULL) { if (bravais != NULL) {
cel_free_cell(bravais); cel_free_cell(bravais);
bravais = NULL;
} }
if (dataset->equivalent_atoms != NULL) { if (dataset->equivalent_atoms != NULL) {
free(dataset->equivalent_atoms); free(dataset->equivalent_atoms);
@ -1184,6 +1182,7 @@ static int set_dataset(SpglibDataset * dataset,
} }
if (symmetry != NULL) { if (symmetry != NULL) {
sym_free_symmetry(symmetry); sym_free_symmetry(symmetry);
symmetry = NULL;
} }
return 0; return 0;
} }
@ -1268,12 +1267,14 @@ static int get_symmetry_with_collinear_spin(int rotation[][3][3],
0, 0,
symprec)) == NULL) { symprec)) == NULL) {
cel_free_cell(cell); cel_free_cell(cell);
cell = NULL;
goto err; goto err;
} }
if ((sym_nonspin = sym_alloc_symmetry(dataset->n_operations)) == NULL) { if ((sym_nonspin = sym_alloc_symmetry(dataset->n_operations)) == NULL) {
spg_free_dataset(dataset); spg_free_dataset(dataset);
cel_free_cell(cell); cel_free_cell(cell);
cell = NULL;
goto err; goto err;
} }
@ -1289,11 +1290,14 @@ static int get_symmetry_with_collinear_spin(int rotation[][3][3],
spins, spins,
symprec)) == NULL) { symprec)) == NULL) {
sym_free_symmetry(sym_nonspin); sym_free_symmetry(sym_nonspin);
sym_nonspin = NULL;
cel_free_cell(cell); cel_free_cell(cell);
cell = NULL;
goto err; goto err;
} }
sym_free_symmetry(sym_nonspin); sym_free_symmetry(sym_nonspin);
sym_nonspin = NULL;
if (symmetry->size > max_size) { if (symmetry->size > max_size) {
fprintf(stderr, "spglib: Indicated max size(=%d) is less than number ", fprintf(stderr, "spglib: Indicated max size(=%d) is less than number ",
@ -1311,7 +1315,9 @@ static int get_symmetry_with_collinear_spin(int rotation[][3][3],
ret: ret:
sym_free_symmetry(symmetry); sym_free_symmetry(symmetry);
symmetry = NULL;
cel_free_cell(cell); cel_free_cell(cell);
cell = NULL;
return size; return size;
@ -1402,6 +1408,7 @@ static int standardize_primitive(double lattice[3][3],
centering, centering,
symprec); symprec);
cel_free_cell(bravais); cel_free_cell(bravais);
bravais = NULL;
if (primitive == NULL) { if (primitive == NULL) {
goto err; goto err;
@ -1411,6 +1418,7 @@ static int standardize_primitive(double lattice[3][3],
num_prim_atom = primitive->size; num_prim_atom = primitive->size;
cel_free_cell(primitive); cel_free_cell(primitive);
primitive = NULL;
return num_prim_atom; return num_prim_atom;
@ -1501,6 +1509,7 @@ static int get_standardized_cell(double lattice[3][3],
symprec); symprec);
spg_free_dataset(dataset); spg_free_dataset(dataset);
cel_free_cell(cell); cel_free_cell(cell);
cell = NULL;
if (std_cell == NULL) { if (std_cell == NULL) {
goto err; goto err;
@ -1510,6 +1519,7 @@ static int get_standardized_cell(double lattice[3][3],
num_std_atom = std_cell->size; num_std_atom = std_cell->size;
cel_free_cell(std_cell); cel_free_cell(std_cell);
std_cell = NULL;
return num_std_atom; return num_std_atom;
@ -1563,12 +1573,14 @@ static int get_international(char symbol[11],
if ((primitive = spa_get_spacegroup(&spacegroup, cell, symprec)) != NULL) { if ((primitive = spa_get_spacegroup(&spacegroup, cell, symprec)) != NULL) {
prm_free_primitive(primitive); prm_free_primitive(primitive);
primitive = NULL;
if (spacegroup.number > 0) { if (spacegroup.number > 0) {
strcpy(symbol, spacegroup.international_short); strcpy(symbol, spacegroup.international_short);
} }
} }
cel_free_cell(cell); cel_free_cell(cell);
cell = NULL;
return spacegroup.number; return spacegroup.number;
} }
@ -1596,12 +1608,14 @@ static int get_schoenflies(char symbol[10],
if ((primitive = spa_get_spacegroup(&spacegroup, cell, symprec)) != NULL) { if ((primitive = spa_get_spacegroup(&spacegroup, cell, symprec)) != NULL) {
prm_free_primitive(primitive); prm_free_primitive(primitive);
primitive = NULL;
if (spacegroup.number > 0) { if (spacegroup.number > 0) {
strcpy(symbol, spacegroup.schoenflies); strcpy(symbol, spacegroup.schoenflies);
} }
} }
cel_free_cell(cell); cel_free_cell(cell);
cell = NULL;
return spacegroup.number; return spacegroup.number;
} }
@ -1632,6 +1646,7 @@ static int get_symmetry_numerical(int rotation[][3][3],
if ((symmetry = sym_get_operation(cell, symprec)) == NULL) { if ((symmetry = sym_get_operation(cell, symprec)) == NULL) {
cel_free_cell(cell); cel_free_cell(cell);
cell = NULL;
return 0; return 0;
} }
@ -1650,7 +1665,9 @@ static int get_symmetry_numerical(int rotation[][3][3],
ret: ret:
sym_free_symmetry(symmetry); sym_free_symmetry(symmetry);
symmetry = NULL;
cel_free_cell(cell); cel_free_cell(cell);
cell = NULL;
return size; return size;
} }
@ -1688,7 +1705,6 @@ static int get_ir_reciprocal_mesh(int grid_address[][3],
int num_ir, i; int num_ir, i;
MatINT *rotations, *rot_reciprocal; MatINT *rotations, *rot_reciprocal;
dataset = get_dataset(lattice, dataset = get_dataset(lattice,
position, position,
types, types,
@ -1711,7 +1727,9 @@ static int get_ir_reciprocal_mesh(int grid_address[][3],
is_shift, is_shift,
rot_reciprocal); rot_reciprocal);
mat_free_MatINT(rot_reciprocal); mat_free_MatINT(rot_reciprocal);
rot_reciprocal = NULL;
mat_free_MatINT(rotations); mat_free_MatINT(rotations);
rotations = NULL;
spg_free_dataset(dataset); spg_free_dataset(dataset);
return num_ir; return num_ir;
} }
@ -1749,6 +1767,7 @@ static int get_stabilized_reciprocal_mesh(int grid_address[][3],
qpoints); qpoints);
mat_free_MatINT(rot_real); mat_free_MatINT(rot_real);
rot_real = NULL;
return num_ir; return num_ir;
} }

2
c/spglib/spin.c
View File

@ -147,7 +147,9 @@ static Symmetry * get_collinear_operations(SPGCONST Symmetry *sym_nonspin,
} }
mat_free_MatINT(rot); mat_free_MatINT(rot);
rot = NULL;
mat_free_VecDBL(trans); mat_free_VecDBL(trans);
trans = NULL;
return symmetry; return symmetry;
} }

168
c/spglib/symmetry.c
View File

@ -36,18 +36,18 @@
#include <stdio.h> #include <stdio.h>
#include <stdlib.h> #include <stdlib.h>
#include "cell.h" #include "cell.h"
#include "lattice.h" #include "delaunay.h"
#include "mathfunc.h" #include "mathfunc.h"
#include "symmetry.h" #include "symmetry.h"
#include "debug.h" #include "debug.h"
#define NUM_ATOMS_CRITERION_FOR_OPENMP 1000 #define NUM_ATOMS_CRITERION_FOR_OPENMP 1000
#define REDUCE_RATE 0.95 #define ANGLE_REDUCE_RATE 0.95
#define PI 3.14159265358979323846 #define PI 3.14159265358979323846
/* Tolerance of angle between lattice vectors in degrees */ /* Tolerance of angle between lattice vectors in degrees */
/* Negative value invokes converter from symprec. */ /* Negative value invokes converter from symprec. */
static double angle_tolerance = -1.0; static double angle_tolerance = -1.0;
static int relative_axes[][3] = { static int relative_axes[][3] = {
{ 1, 0, 0}, { 1, 0, 0},
@ -88,10 +88,12 @@ static VecDBL * get_translation(SPGCONST int rot[3][3],
const double symprec, const double symprec,
const int is_identity); const int is_identity);
static Symmetry * get_operations(SPGCONST Cell *primitive, static Symmetry * get_operations(SPGCONST Cell *primitive,
const double symprec); const double symprec,
const double angle_symprec);
static Symmetry * reduce_operation(SPGCONST Cell * primitive, static Symmetry * reduce_operation(SPGCONST Cell * primitive,
SPGCONST Symmetry * symmetry, SPGCONST Symmetry * symmetry,
const double symprec); const double symprec,
const double angle_symprec);
static int search_translation_part(int lat_point_atoms[], static int search_translation_part(int lat_point_atoms[],
SPGCONST Cell * cell, SPGCONST Cell * cell,
SPGCONST int rot[3][3], SPGCONST int rot[3][3],
@ -115,10 +117,12 @@ get_space_group_operations(SPGCONST PointSymmetry *lattice_sym,
static void set_axes(int axes[3][3], static void set_axes(int axes[3][3],
const int a1, const int a2, const int a3); const int a1, const int a2, const int a3);
static PointSymmetry get_lattice_symmetry(SPGCONST double cell_lattice[3][3], static PointSymmetry get_lattice_symmetry(SPGCONST double cell_lattice[3][3],
const double symprec); const double symprec,
const double angle_symprec);
static int is_identity_metric(SPGCONST double metric_rotated[3][3], static int is_identity_metric(SPGCONST double metric_rotated[3][3],
SPGCONST double metric_orig[3][3], SPGCONST double metric_orig[3][3],
const double symprec); const double symprec,
const double angle_symprec);
static double get_angle(SPGCONST double metric[3][3], static double get_angle(SPGCONST double metric[3][3],
const int i, const int i,
const int j); const int j);
@ -160,7 +164,7 @@ Symmetry * sym_alloc_symmetry(const int size)
free(symmetry); free(symmetry);
symmetry = NULL; symmetry = NULL;
return NULL; return NULL;
} }
return symmetry; return symmetry;
} }
@ -174,7 +178,6 @@ void sym_free_symmetry(Symmetry *symmetry)
symmetry->trans = NULL; symmetry->trans = NULL;
} }
free(symmetry); free(symmetry);
symmetry = NULL;
} }
/* Return NULL if failed */ /* Return NULL if failed */
@ -184,7 +187,7 @@ Symmetry * sym_get_operation(SPGCONST Cell * primitive,
debug_print("sym_get_operations:\n"); debug_print("sym_get_operations:\n");
return get_operations(primitive, symprec); return get_operations(primitive, symprec, angle_tolerance);
} }
/* Return NULL if failed */ /* Return NULL if failed */
@ -192,7 +195,7 @@ Symmetry * sym_reduce_operation(SPGCONST Cell * primitive,
SPGCONST Symmetry * symmetry, SPGCONST Symmetry * symmetry,
const double symprec) const double symprec)
{ {
return reduce_operation(primitive, symmetry, symprec); return reduce_operation(primitive, symmetry, symprec, angle_tolerance);
} }
/* Return NULL if failed */ /* Return NULL if failed */
@ -208,6 +211,8 @@ VecDBL * sym_get_pure_translation(SPGCONST Cell *cell,
pure_trans = NULL; pure_trans = NULL;
if ((pure_trans = get_translation(identity, cell, symprec, 1)) == NULL) { if ((pure_trans = get_translation(identity, cell, symprec, 1)) == NULL) {
warning_print("spglib: get_translation failed (line %d, %s).\n",
__LINE__, __FILE__);
return NULL; return NULL;
} }
@ -235,7 +240,7 @@ VecDBL * sym_reduce_pure_translation(SPGCONST Cell * cell,
symmetry = NULL; symmetry = NULL;
symmetry_reduced = NULL; symmetry_reduced = NULL;
pure_trans_reduced = NULL; pure_trans_reduced = NULL;
multi = pure_trans->size; multi = pure_trans->size;
if ((symmetry = sym_alloc_symmetry(multi)) == NULL) { if ((symmetry = sym_alloc_symmetry(multi)) == NULL) {
@ -247,24 +252,29 @@ VecDBL * sym_reduce_pure_translation(SPGCONST Cell * cell,
mat_copy_vector_d3(symmetry->trans[i], pure_trans->vec[i]); mat_copy_vector_d3(symmetry->trans[i], pure_trans->vec[i]);
} }
if ((symmetry_reduced = reduce_operation(cell, symmetry, symprec)) == NULL) { if ((symmetry_reduced =
reduce_operation(cell, symmetry, symprec, angle_tolerance)) == NULL) {
sym_free_symmetry(symmetry); sym_free_symmetry(symmetry);
symmetry = NULL;
return NULL; return NULL;
} }
sym_free_symmetry(symmetry); sym_free_symmetry(symmetry);
symmetry = NULL;
multi = symmetry_reduced->size; multi = symmetry_reduced->size;
if ((pure_trans_reduced = mat_alloc_VecDBL(multi)) == NULL) { if ((pure_trans_reduced = mat_alloc_VecDBL(multi)) == NULL) {
sym_free_symmetry(symmetry_reduced); sym_free_symmetry(symmetry_reduced);
symmetry_reduced = NULL;
return NULL; return NULL;
} }
for (i = 0; i < multi; i++) { for (i = 0; i < multi; i++) {
mat_copy_vector_d3(pure_trans_reduced->vec[i], symmetry_reduced->trans[i]); mat_copy_vector_d3(pure_trans_reduced->vec[i], symmetry_reduced->trans[i]);
} }
sym_free_symmetry(symmetry_reduced); sym_free_symmetry(symmetry_reduced);
symmetry_reduced = NULL;
return pure_trans_reduced; return pure_trans_reduced;
} }
@ -289,7 +299,8 @@ double sym_get_angle_tolerance(void)
/* transformed to those of original input cells, if the input cell */ /* transformed to those of original input cells, if the input cell */
/* was not a primitive cell. */ /* was not a primitive cell. */
static Symmetry * get_operations(SPGCONST Cell *primitive, static Symmetry * get_operations(SPGCONST Cell *primitive,
const double symprec) const double symprec,
const double angle_symprec)
{ {
PointSymmetry lattice_sym; PointSymmetry lattice_sym;
Symmetry *symmetry; Symmetry *symmetry;
@ -298,26 +309,27 @@ static Symmetry * get_operations(SPGCONST Cell *primitive,
symmetry = NULL; symmetry = NULL;
lattice_sym = get_lattice_symmetry(primitive->lattice, symprec); lattice_sym = get_lattice_symmetry(primitive->lattice,
symprec,
angle_symprec);
if (lattice_sym.size == 0) { if (lattice_sym.size == 0) {
debug_print("get_lattice_symmetry failed.\n"); return NULL;
goto end;
} }
if ((symmetry = get_space_group_operations(&lattice_sym, if ((symmetry = get_space_group_operations(&lattice_sym,
primitive, primitive,
symprec)) == NULL) { symprec)) == NULL) {
goto end; return NULL;
} }
end:
return symmetry; return symmetry;
} }
/* Return NULL if failed */ /* Return NULL if failed */
static Symmetry * reduce_operation(SPGCONST Cell * primitive, static Symmetry * reduce_operation(SPGCONST Cell * primitive,
SPGCONST Symmetry * symmetry, SPGCONST Symmetry * symmetry,
const double symprec) const double symprec,
const double angle_symprec)
{ {
int i, j, num_sym; int i, j, num_sym;
Symmetry * sym_reduced; Symmetry * sym_reduced;
@ -331,7 +343,12 @@ static Symmetry * reduce_operation(SPGCONST Cell * primitive,
rot = NULL; rot = NULL;
trans = NULL; trans = NULL;
point_symmetry = get_lattice_symmetry(primitive->lattice, symprec); point_symmetry = get_lattice_symmetry(primitive->lattice,
symprec,
angle_symprec);
if (point_symmetry.size == 0) {
return NULL;
}
if ((rot = mat_alloc_MatINT(symmetry->size)) == NULL) { if ((rot = mat_alloc_MatINT(symmetry->size)) == NULL) {
return NULL; return NULL;
@ -339,9 +356,10 @@ static Symmetry * reduce_operation(SPGCONST Cell * primitive,
if ((trans = mat_alloc_VecDBL(symmetry->size)) == NULL) { if ((trans = mat_alloc_VecDBL(symmetry->size)) == NULL) {
mat_free_MatINT(rot); mat_free_MatINT(rot);
rot = NULL;
return NULL; return NULL;
} }
num_sym = 0; num_sym = 0;
for (i = 0; i < point_symmetry.size; i++) { for (i = 0; i < point_symmetry.size; i++) {
for (j = 0; j < symmetry->size; j++) { for (j = 0; j < symmetry->size; j++) {
@ -368,7 +386,9 @@ static Symmetry * reduce_operation(SPGCONST Cell * primitive,
} }
mat_free_MatINT(rot); mat_free_MatINT(rot);
rot = NULL;
mat_free_VecDBL(trans); mat_free_VecDBL(trans);
trans = NULL;
return sym_reduced; return sym_reduced;
} }
@ -497,7 +517,7 @@ static VecDBL * get_translation(SPGCONST int rot[3][3],
ret: ret:
free(is_found); free(is_found);
is_found = NULL; is_found = NULL;
return trans; return trans;
} }
@ -546,7 +566,7 @@ static int is_overlap_all_atoms(const double trans[3],
double pos_rot[3], d_frac[3], d[3]; double pos_rot[3], d_frac[3], d[3];
symprec2 = symprec * symprec; symprec2 = symprec * symprec;
for (i = 0; i < cell->size; i++) { for (i = 0; i < cell->size; i++) {
if (is_identity) { /* Identity matrix is treated as special for speed-up. */ if (is_identity) { /* Identity matrix is treated as special for speed-up. */
for (j = 0; j < 3; j++) { for (j = 0; j < 3; j++) {
@ -600,11 +620,11 @@ static int get_index_with_least_atoms(const Cell *cell)
warning_print("spglib: Memory could not be allocated "); warning_print("spglib: Memory could not be allocated ");
return -1; return -1;
} }
for (i = 0; i < cell->size; i++) { for (i = 0; i < cell->size; i++) {
mapping[i] = 0; mapping[i] = 0;
} }
for (i = 0; i < cell->size; i++) { for (i = 0; i < cell->size; i++) {
for (j = 0; j < cell->size; j++) { for (j = 0; j < cell->size; j++) {
if (cell->types[i] == cell->types[j]) { if (cell->types[i] == cell->types[j]) {
@ -613,7 +633,7 @@ static int get_index_with_least_atoms(const Cell *cell)
} }
} }
} }
min = mapping[0]; min = mapping[0];
min_index = 0; min_index = 0;
for (i = 0; i < cell->size; i++) { for (i = 0; i < cell->size; i++) {
@ -643,7 +663,7 @@ get_space_group_operations(SPGCONST PointSymmetry *lattice_sym,
trans = NULL; trans = NULL;
symmetry = NULL; symmetry = NULL;
if ((trans = (VecDBL**) malloc(sizeof(VecDBL*) * lattice_sym->size)) if ((trans = (VecDBL**) malloc(sizeof(VecDBL*) * lattice_sym->size))
== NULL) { == NULL) {
warning_print("spglib: Memory could not be allocated "); warning_print("spglib: Memory could not be allocated ");
@ -687,6 +707,7 @@ get_space_group_operations(SPGCONST PointSymmetry *lattice_sym,
for (i = 0; i < lattice_sym->size; i++) { for (i = 0; i < lattice_sym->size; i++) {
if (trans[i] != NULL) { if (trans[i] != NULL) {
mat_free_VecDBL(trans[i]); mat_free_VecDBL(trans[i]);
trans[i] = NULL;
} }
} }
free(trans); free(trans);
@ -696,9 +717,11 @@ get_space_group_operations(SPGCONST PointSymmetry *lattice_sym,
} }
static PointSymmetry get_lattice_symmetry(SPGCONST double cell_lattice[3][3], static PointSymmetry get_lattice_symmetry(SPGCONST double cell_lattice[3][3],
const double symprec) const double symprec,
const double angle_symprec)
{ {
int i, j, k, num_sym; int i, j, k, attempt, num_sym;
double angle_tol;
int axes[3][3]; int axes[3][3];
double lattice[3][3], min_lattice[3][3]; double lattice[3][3], min_lattice[3][3];
double metric[3][3], metric_orig[3][3]; double metric[3][3], metric_orig[3][3];
@ -706,54 +729,62 @@ static PointSymmetry get_lattice_symmetry(SPGCONST double cell_lattice[3][3],
debug_print("get_lattice_symmetry:\n"); debug_print("get_lattice_symmetry:\n");
if (! lat_smallest_lattice_vector(min_lattice, lattice_sym.size = 0;
cell_lattice,
symprec)) { if (! del_delaunay_reduce(min_lattice, cell_lattice, symprec)) {
goto err; goto err;
} }
mat_get_metric(metric_orig, min_lattice); mat_get_metric(metric_orig, min_lattice);
angle_tol = angle_symprec;
num_sym = 0; for (attempt = 0; attempt < 100; attempt++) {
for (i = 0; i < 26; i++) { num_sym = 0;
for (j = 0; j < 26; j++) { for (i = 0; i < 26; i++) {
for (k = 0; k < 26; k++) { for (j = 0; j < 26; j++) {
set_axes(axes, i, j, k); for (k = 0; k < 26; k++) {
if (! ((mat_get_determinant_i3(axes) == 1) || set_axes(axes, i, j, k);
(mat_get_determinant_i3(axes) == -1))) { if (! ((mat_get_determinant_i3(axes) == 1) ||
continue; (mat_get_determinant_i3(axes) == -1))) {
} continue;
mat_multiply_matrix_di3(lattice, min_lattice, axes); }
mat_get_metric(metric, lattice); mat_multiply_matrix_di3(lattice, min_lattice, axes);
mat_get_metric(metric, lattice);
if (is_identity_metric(metric, metric_orig, symprec)) {
mat_copy_matrix_i3(lattice_sym.rot[num_sym], axes); if (is_identity_metric(metric, metric_orig, symprec, angle_tol)) {
num_sym++; if (num_sym > 47) {
} angle_tol *= ANGLE_REDUCE_RATE;
warning_print("spglib: Too many lattice symmetries was found.\n");
if (num_sym > 48) { warning_print(" Reduce angle tolerance to %f", angle_tol);
warning_print("spglib: Too many lattice symmetries was found.\n"); warning_print(" (line %d, %s).\n", __LINE__, __FILE__);
warning_print(" Tolerance may be too large "); goto next_attempt;
warning_print("(line %d, %s).\n", __LINE__, __FILE__); }
goto err;
mat_copy_matrix_i3(lattice_sym.rot[num_sym], axes);
num_sym++;
}
} }
} }
} }
if (num_sym < 49 || angle_tol < 0) {
lattice_sym.size = num_sym;
return transform_pointsymmetry(&lattice_sym, cell_lattice, min_lattice);
}
next_attempt:
;
} }
lattice_sym.size = num_sym;
return transform_pointsymmetry(&lattice_sym,
cell_lattice,
min_lattice);
err: err:
lattice_sym.size = 0; debug_print("get_lattice_symmetry failed.\n");
return lattice_sym; return lattice_sym;
} }
static int is_identity_metric(SPGCONST double metric_rotated[3][3], static int is_identity_metric(SPGCONST double metric_rotated[3][3],
SPGCONST double metric_orig[3][3], SPGCONST double metric_orig[3][3],
const double symprec) const double symprec,
const double angle_symprec)
{ {
int i, j, k; int i, j, k;
int elem_sets[3][2] = {{0, 1}, int elem_sets[3][2] = {{0, 1},
@ -769,13 +800,13 @@ static int is_identity_metric(SPGCONST double metric_rotated[3][3],
goto fail; goto fail;
} }
} }
for (i = 0; i < 3; i++) { for (i = 0; i < 3; i++) {
j = elem_sets[i][0]; j = elem_sets[i][0];
k = elem_sets[i][1]; k = elem_sets[i][1];
if (angle_tolerance > 0) { if (angle_symprec > 0) {
if (mat_Dabs(get_angle(metric_orig, j, k) - if (mat_Dabs(get_angle(metric_orig, j, k) -
get_angle(metric_rotated, j, k)) > angle_tolerance) { get_angle(metric_rotated, j, k)) > angle_symprec) {
goto fail; goto fail;
} }
} else { } else {
@ -824,6 +855,8 @@ transform_pointsymmetry(SPGCONST PointSymmetry * lat_sym_orig,
double trans_mat[3][3], inv_mat[3][3], drot[3][3]; double trans_mat[3][3], inv_mat[3][3], drot[3][3];
PointSymmetry lat_sym_new; PointSymmetry lat_sym_new;
lat_sym_new.size = 0;
mat_inverse_matrix_d3(inv_mat, original_lattice, 0); mat_inverse_matrix_d3(inv_mat, original_lattice, 0);
mat_multiply_matrix_d3(trans_mat, inv_mat, new_lattice); mat_multiply_matrix_d3(trans_mat, inv_mat, new_lattice);
@ -856,7 +889,6 @@ transform_pointsymmetry(SPGCONST PointSymmetry * lat_sym_orig,
return lat_sym_new; return lat_sym_new;
err: err:
lat_sym_new.size = 0;
return lat_sym_new; return lat_sym_new;
} }

48
c/spglib_h/lattice.h
View File

@ -1,48 +0,0 @@
/* Copyright (C) 2010 Atsushi Togo */
/* All rights reserved. */
/* This file is part of spglib. */
/* 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 __lattice_H__
#define __lattice_H__
#include "mathfunc.h"
int lat_smallest_lattice_vector(double lattice_new[3][3],
SPGCONST double lattice[3][3],
const double symprec);
int lat_smallest_lattice_vector_2D(double min_lattice[3][3],
SPGCONST double lattice[3][3],
const int unique_axis,
const double symprec);
#endif

7
c/spglib_h/niggli.h
View File

@ -35,6 +35,13 @@
#ifndef __NIGGLI_H__ #ifndef __NIGGLI_H__
#define __NIGGLI_H__ #define __NIGGLI_H__
#define NIGGLI_MAJOR_VERSION 0
#define NIGGLI_MINOR_VERSION 1
#define NIGGLI_MICRO_VERSION 2
int niggli_get_major_version(void);
int niggli_get_minor_version(void);
int niggli_get_micro_version(void);
int niggli_reduce(double *lattice_, const double eps_); int niggli_reduce(double *lattice_, const double eps_);
#endif #endif

3
c/spglib_h/pointgroup.h
View File

@ -36,6 +36,7 @@
#define __pointgroup_H__ #define __pointgroup_H__
#include "mathfunc.h" #include "mathfunc.h"
#include "symmetry.h"
typedef enum { typedef enum {
HOLOHEDRY_NONE, HOLOHEDRY_NONE,
@ -74,4 +75,6 @@ Pointgroup ptg_get_transformation_matrix(int transform_mat[3][3],
SPGCONST int rotations[][3][3], SPGCONST int rotations[][3][3],
const int num_rotations); const int num_rotations);
Pointgroup ptg_get_pointgroup(const int pointgroup_number); Pointgroup ptg_get_pointgroup(const int pointgroup_number);
PointSymmetry ptg_get_pointsymmetry(SPGCONST int rotations[][3][3],
const int num_rotations);
#endif #endif

46
c/spglib_h/spglib.h
View File

@ -245,12 +245,6 @@ int spgat_get_multiplicity(SPGCONST double lattice[3][3],
const double symprec, const double symprec,
const double angle_tolerance); const double angle_tolerance);
/* Considering periodicity of crystal, one of the possible smallest */
/* lattice is searched. The lattice is stored in ``smallest_lattice``. */
int spg_get_smallest_lattice(double smallest_lattice[3][3],
SPGCONST double lattice[3][3],
const double symprec);
/* Space group is found in international table symbol (``symbol``) and */ /* Space group is found in international table symbol (``symbol``) and */
/* number (return value). 0 is returned when it fails. */ /* number (return value). 0 is returned when it fails. */
int spg_get_international(char symbol[11], int spg_get_international(char symbol[11],
@ -321,11 +315,9 @@ int spgat_standardize_cell(double lattice[3][3],
const double symprec, const double symprec,
const double angle_tolerance); const double angle_tolerance);
/************/ /* This is a wrapper of spg_standardize_cell. */
/* Obsolete */ /* A primitive cell is found from an input cell. */
/************/ /* Be careful that ``lattice``, ``position``, and ``types`` are overwritten. */
/* A primitive cell is found from an input cell. Be careful that */
/* ``lattice``, ``position``, and ``types`` are overwritten. */
/* ``num_atom`` is returned as return value. */ /* ``num_atom`` is returned as return value. */
/* When any primitive cell is not found, 0 is returned. */ /* When any primitive cell is not found, 0 is returned. */
int spg_find_primitive(double lattice[3][3], int spg_find_primitive(double lattice[3][3],
@ -341,9 +333,7 @@ int spgat_find_primitive(double lattice[3][3],
const double symprec, const double symprec,
const double angle_tolerance); const double angle_tolerance);
/************/ /* This is a wrapper of spg_standardize_cell. */
/* Obsolete */
/************/
/* Bravais lattice with internal atomic points are returned. */ /* Bravais lattice with internal atomic points are returned. */
/* The arrays are require to have 4 times larger memory space */ /* The arrays are require to have 4 times larger memory space */
/* those of input cell. */ /* those of input cell. */
@ -362,6 +352,9 @@ int spgat_refine_cell(double lattice[3][3],
const double symprec, const double symprec,
const double angle_tolerance); const double angle_tolerance);
/* Delaunay reduction for lattice parameters */
/* ``lattice`` is overwritten when the redution ends succeeded. */
int spg_delaunay_reduce(double lattice[3][3], const double symprec);
/*---------*/ /*---------*/
/* kpoints */ /* kpoints */
@ -475,26 +468,11 @@ void spg_get_neighboring_grid_points(int relative_grid_points[],
SPGCONST int bz_grid_address[][3], SPGCONST int bz_grid_address[][3],
const int bz_map[]); const int bz_map[]);
/* /\*--------------------*\/ */ /*--------*/
/* /\* tetrahedron method *\/ */ /* Niggli */
/* /\*--------------------*\/ */ /*--------*/
/* void */ /* Return 0 if failed */
/* spg_get_tetrahedra_relative_grid_address(int relative_grid_address[24][4][3], */ int spg_niggli_reduce(double lattice[3][3], const double symprec);
/* SPGCONST double rec_lattice[3][3]); */
/* void */
/* spg_get_all_tetrahedra_relative_grid_address */
/* (int relative_grid_address[4][24][4][3]); */
/* double */
/* spg_get_tetrahedra_integration_weight(const double omega, */
/* SPGCONST double tetrahedra_omegas[24][4], */
/* const char function); */
/* void */
/* spg_get_tetrahedra_integration_weight_at_omegas */
/* (double integration_weights[], */
/* const int num_omegas, */
/* const double omegas[], */
/* SPGCONST double tetrahedra_omegas[24][4], */
/* const char function); */
#endif #endif

2
c/spglib_h/version.h
View File

@ -37,7 +37,7 @@
#define SPGLIB_MAJOR_VERSION 1 #define SPGLIB_MAJOR_VERSION 1
#define SPGLIB_MINOR_VERSION 9 #define SPGLIB_MINOR_VERSION 9
#define SPGLIB_MICRO_VERSION 0 #define SPGLIB_MICRO_VERSION 2
#endif #endif

163
phonopy/structure/spglib.py
View File

@ -38,29 +38,46 @@ import numpy as np
def get_version(): def get_version():
return tuple(spg.version()) return tuple(spg.version())
def get_symmetry(bulk, use_magmoms=False, symprec=1e-5, angle_tolerance=-1.0): def get_symmetry(cell, use_magmoms=False, symprec=1e-5, angle_tolerance=-1.0):
""" """This gives crystal symmetry operations from a crystal structure.
Return symmetry operations as hash.
Hash key 'rotations' gives the numpy integer array Args:
of the rotation matrices for scaled positions cell: Crystal structrue given either in Atoms object or tuple.
Hash key 'translations' gives the numpy double array In the case given by a tuple, it has to follow the form below,
of the translation vectors in scaled positions (Lattice parameters in a 3x3 array (see the detail below),
Fractional atomic positions in an Nx3 array,
Integer numbers to distinguish species in a length N array,
(optional) Collinear magnetic moments in a length N array),
where N is the number of atoms.
Lattice parameters are given in the form:
[[a_x, a_y, a_z],
[b_x, b_y, b_z],
[c_x, c_y, c_z]]
use_magmoms:
bool: If True, collinear magnetic polarizatin is considered.
symprec:
float: Symmetry search tolerance in the unit of length.
angle_tolerance:
float: Symmetry search tolerance in the unit of angle deg.
If the value is negative, an internally optimized routine
is used to judge symmetry.
Return:
dictionary: Rotation parts and translation parts.
'rotations': Gives the numpy 'intc' array of the rotation matrices.
'translations': Gives the numpy 'double' array of fractional
translations with respect to a, b, c axes.
""" """
# Atomic positions have to be specified by scaled positions for spglib. lattice, positions, numbers, magmoms = _expand_cell(
positions = np.array(bulk.get_scaled_positions(), dtype='double', order='C') cell, use_magmoms=use_magmoms)
lattice = np.array(bulk.get_cell().T, dtype='double', order='C') multi = 48 * len(positions)
numbers = np.array(bulk.get_atomic_numbers(), dtype='intc')
# Get number of symmetry operations and allocate symmetry operations
# multi = spg.multiplicity(cell, positions, numbers, symprec)
multi = 48 * bulk.get_number_of_atoms()
rotation = np.zeros((multi, 3, 3), dtype='intc') rotation = np.zeros((multi, 3, 3), dtype='intc')
translation = np.zeros((multi, 3), dtype='double') translation = np.zeros((multi, 3), dtype='double')
# Get symmetry operations # Get symmetry operations
if use_magmoms: if use_magmoms:
magmoms = bulk.get_magnetic_moments()
equivalent_atoms = np.zeros(len(magmoms), dtype='intc') equivalent_atoms = np.zeros(len(magmoms), dtype='intc')
num_sym = spg.symmetry_with_collinear_spin(rotation, num_sym = spg.symmetry_with_collinear_spin(rotation,
translation, translation,
@ -90,7 +107,7 @@ def get_symmetry(bulk, use_magmoms=False, symprec=1e-5, angle_tolerance=-1.0):
'translations': np.array(translation[:num_sym], 'translations': np.array(translation[:num_sym],
dtype='double', order='C')} dtype='double', order='C')}
def get_symmetry_dataset(bulk, symprec=1e-5, angle_tolerance=-1.0): def get_symmetry_dataset(cell, symprec=1e-5, angle_tolerance=-1.0):
""" """
number: International space group number number: International space group number
international: International symbol international: International symbol
@ -109,10 +126,8 @@ def get_symmetry_dataset(bulk, symprec=1e-5, angle_tolerance=-1.0):
Standardized unit cell Standardized unit cell
pointgroup_number, pointgroup_symbol: Point group number (see get_pointgroup) pointgroup_number, pointgroup_symbol: Point group number (see get_pointgroup)
""" """
positions = np.array(bulk.get_scaled_positions(), dtype='double', order='C')
lattice = np.array(bulk.get_cell().T, dtype='double', order='C') lattice, positions, numbers, _ = _expand_cell(cell)
numbers = np.array(bulk.get_atomic_numbers(), dtype='intc')
keys = ('number', keys = ('number',
'hall_number', 'hall_number',
'international', 'international',
@ -158,13 +173,13 @@ def get_symmetry_dataset(bulk, symprec=1e-5, angle_tolerance=-1.0):
return dataset return dataset
def get_spacegroup(bulk, symprec=1e-5, angle_tolerance=-1.0, symbol_type=0): def get_spacegroup(cell, symprec=1e-5, angle_tolerance=-1.0, symbol_type=0):
""" """
Return space group in international table symbol and number Return space group in international table symbol and number
as a string. as a string.
""" """
dataset = get_symmetry_dataset(bulk, dataset = get_symmetry_dataset(cell,
symprec=symprec, symprec=symprec,
angle_tolerance=angle_tolerance) angle_tolerance=angle_tolerance)
symbols = spg.spacegroup_type(dataset['hall_number']) symbols = spg.spacegroup_type(dataset['hall_number'])
@ -215,7 +230,7 @@ def get_pointgroup(rotations):
# (symbol, pointgroup_number, transformation_matrix) # (symbol, pointgroup_number, transformation_matrix)
return spg.pointgroup(np.array(rotations, dtype='intc', order='C')) return spg.pointgroup(np.array(rotations, dtype='intc', order='C'))
def standardize_cell(bulk, def standardize_cell(cell,
to_primitive=0, to_primitive=0,
no_idealize=0, no_idealize=0,
symprec=1e-5, symprec=1e-5,
@ -223,16 +238,17 @@ def standardize_cell(bulk,
""" """
Return standardized cell Return standardized cell
""" """
# Atomic positions have to be specified by scaled positions for spglib.
num_atom = bulk.get_number_of_atoms()
lattice = np.array(bulk.get_cell().T, dtype='double', order='C')
pos = np.zeros((num_atom * 4, 3), dtype='double')
pos[:num_atom] = bulk.get_scaled_positions()
lattice, _positions, _numbers, _ = _expand_cell(cell)
# Atomic positions have to be specified by scaled positions for spglib.
num_atom = len(_positions)
positions = np.zeros((num_atom * 4, 3), dtype='double', order='C')
positions[:num_atom] = _positions
numbers = np.zeros(num_atom * 4, dtype='intc') numbers = np.zeros(num_atom * 4, dtype='intc')
numbers[:num_atom] = np.array(bulk.get_atomic_numbers(), dtype='intc') numbers[:num_atom] = _numbers
num_atom_std = spg.standardize_cell(lattice, num_atom_std = spg.standardize_cell(lattice,
pos, positions,
numbers, numbers,
num_atom, num_atom,
to_primitive, to_primitive,
@ -241,45 +257,41 @@ def standardize_cell(bulk,
angle_tolerance) angle_tolerance)
return (np.array(lattice.T, dtype='double', order='C'), return (np.array(lattice.T, dtype='double', order='C'),
np.array(pos[:num_atom_std], dtype='double', order='C'), np.array(positions[:num_atom_std], dtype='double', order='C'),
np.array(numbers[:num_atom_std], dtype='intc')) np.array(numbers[:num_atom_std], dtype='intc'))
def refine_cell(bulk, symprec=1e-5, angle_tolerance=-1.0): def refine_cell(cell, symprec=1e-5, angle_tolerance=-1.0):
""" """
Return refined cell Return refined cell
""" """
# Atomic positions have to be specified by scaled positions for spglib.
num_atom = bulk.get_number_of_atoms()
lattice = np.array(bulk.get_cell().T, dtype='double', order='C')
pos = np.zeros((num_atom * 4, 3), dtype='double')
pos[:num_atom] = bulk.get_scaled_positions()
lattice, _positions, _numbers, _ = _expand_cell(cell)
# Atomic positions have to be specified by scaled positions for spglib.
num_atom = len(_positions)
positions = np.zeros((num_atom * 4, 3), dtype='double', order='C')
positions[:num_atom] = _positions
numbers = np.zeros(num_atom * 4, dtype='intc') numbers = np.zeros(num_atom * 4, dtype='intc')
numbers[:num_atom] = np.array(bulk.get_atomic_numbers(), dtype='intc') numbers[:num_atom] = _numbers
num_atom_std = spg.refine_cell(lattice, num_atom_std = spg.refine_cell(lattice,
pos, positions,
numbers, numbers,
num_atom, num_atom,
symprec, symprec,
angle_tolerance) angle_tolerance)
return (np.array(lattice.T, dtype='double', order='C'), return (np.array(lattice.T, dtype='double', order='C'),
np.array(pos[:num_atom_std], dtype='double', order='C'), np.array(positions[:num_atom_std], dtype='double', order='C'),
np.array(numbers[:num_atom_std], dtype='intc')) np.array(numbers[:num_atom_std], dtype='intc'))
def find_primitive(bulk, symprec=1e-5, angle_tolerance=-1.0): def find_primitive(cell, symprec=1e-5, angle_tolerance=-1.0):
""" """
A primitive cell in the input cell is searched and returned A primitive cell in the input cell is searched and returned
as an object of Atoms class. as an object of Atoms class.
If no primitive cell is found, (None, None, None) is returned. If no primitive cell is found, (None, None, None) is returned.
""" """
# Atomic positions have to be specified by scaled positions for spglib. lattice, positions, numbers, _ = _expand_cell(cell)
positions = np.array(bulk.get_scaled_positions(), dtype='double', order='C')
lattice = np.array(bulk.get_cell().T, dtype='double', order='C')
numbers = np.array(bulk.get_atomic_numbers(), dtype='intc')
# lattice is transposed with respect to the definition of Atoms class
num_atom_prim = spg.primitive(lattice, num_atom_prim = spg.primitive(lattice,
positions, positions,
numbers, numbers,
@ -317,7 +329,7 @@ def get_grid_point_from_address(grid_address, mesh):
def get_ir_reciprocal_mesh(mesh, def get_ir_reciprocal_mesh(mesh,
bulk, cell,
is_shift=np.zeros(3, dtype='intc'), is_shift=np.zeros(3, dtype='intc'),
is_time_reversal=True, is_time_reversal=True,
symprec=1e-5): symprec=1e-5):
@ -327,6 +339,7 @@ def get_ir_reciprocal_mesh(mesh,
is_shift=[0, 0, 0] gives Gamma center mesh. is_shift=[0, 0, 0] gives Gamma center mesh.
""" """
lattice, positions, numbers, _ = _expand_cell(cell)
mapping = np.zeros(np.prod(mesh), dtype='intc') mapping = np.zeros(np.prod(mesh), dtype='intc')
mesh_points = np.zeros((np.prod(mesh), 3), dtype='intc') mesh_points = np.zeros((np.prod(mesh), 3), dtype='intc')
spg.ir_reciprocal_mesh( spg.ir_reciprocal_mesh(
@ -335,9 +348,9 @@ def get_ir_reciprocal_mesh(mesh,
np.array(mesh, dtype='intc'), np.array(mesh, dtype='intc'),
np.array(is_shift, dtype='intc'), np.array(is_shift, dtype='intc'),
is_time_reversal * 1, is_time_reversal * 1,
np.array(bulk.get_cell().T, dtype='double', order='C'), lattice,
np.array(bulk.get_scaled_positions(), dtype='double', order='C'), positions,
np.array(bulk.get_atomic_numbers(), dtype='intc'), numbers,
symprec) symprec)
return mapping, mesh_points return mapping, mesh_points
@ -457,3 +470,49 @@ def get_stabilized_reciprocal_mesh(mesh,
qpoints) qpoints)
return mapping, mesh_points return mapping, mesh_points
def niggli_reduce(lattice, eps=1e-5):
"""Run Niggli reduction
Args:
lattice: Lattice parameters in the form of
[[a_x, a_y, a_z],
[b_x, b_y, b_z],
[c_x, c_y, c_z]]
eps: Tolerance.
Returns:
if the Niggli reduction succeeded:
Reduced lattice parameters are given as a numpy 'double' array:
[[a_x, a_y, a_z],
[b_x, b_y, b_z],
[c_x, c_y, c_z]]
otherwise returns None.
"""
niggli_lattice = np.array(np.transpose(lattice), dtype='double', order='C')
result = spg.niggli_reduce(niggli_lattice, float(eps))
if result == 0:
return None
else:
return np.array(np.transpose(niggli_lattice), dtype='double', order='C')
def _expand_cell(cell, use_magmoms=False):
if isinstance(cell, tuple):
lattice = np.array(np.transpose(cell[0]), dtype='double', order='C')
positions = np.array(cell[1], dtype='double', order='C')
numbers = np.array(cell[2], dtype='intc')
if len(cell) > 3 and use_magmoms:
magmoms = np.array(cell[3], dtype='double')
else:
magmoms = None
else:
lattice = np.array(cell.get_cell().T, dtype='double', order='C')
positions = np.array(cell.get_scaled_positions(),
dtype='double', order='C')
numbers = np.array(cell.get_atomic_numbers(), dtype='intc')
if use_magmoms:
magmoms = np.array(cell.get_magnetic_moments(), dtype='double')
else:
magmoms = None
return (lattice, positions, numbers, magmoms)

2
setup.py
View File

@ -62,10 +62,10 @@ extension_spglib = Extension(
extra_link_args=extra_link_args_spglib, extra_link_args=extra_link_args_spglib,
sources=['c/_spglib.c', sources=['c/_spglib.c',
'c/spglib/cell.c', 'c/spglib/cell.c',
'c/spglib/delaunay.c',
'c/spglib/hall_symbol.c', 'c/spglib/hall_symbol.c',
'c/spglib/kgrid.c', 'c/spglib/kgrid.c',
'c/spglib/kpoint.c', 'c/spglib/kpoint.c',
'c/spglib/lattice.c',
'c/spglib/mathfunc.c', 'c/spglib/mathfunc.c',
'c/spglib/niggli.c', 'c/spglib/niggli.c',
'c/spglib/pointgroup.c', 'c/spglib/pointgroup.c',