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

Merge branch 'ExpHP-distribute-fc2-all-with-perm' into develop

This commit is contained in:
Atsushi Togo 2017-08-14 13:41:36 +09:00
parent e9944690fa f217a2e500
commit 095c2ec57a
2 changed files with 234 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

159
c/_phonopy.c
View File

@ -52,6 +52,7 @@ static PyObject * py_get_derivative_dynmat(PyObject *self, PyObject *args);
static PyObject * py_get_thermal_properties(PyObject *self, PyObject *args);
static PyObject * py_distribute_fc2(PyObject *self, PyObject *args);
static PyObject * py_distribute_fc2_all(PyObject *self, PyObject *args);
static PyObject * py_distribute_fc2_with_mappings(PyObject *self, PyObject *args);
static PyObject * py_compute_permutation(PyObject *self, PyObject *args);
static int distribute_fc2(double *fc2,
@ -65,6 +66,16 @@ static int distribute_fc2(double *fc2,
const double t[3],
const double symprec);
static void distribute_fc2_with_mappings(double (*fc2)[3][3],
const int * atom_list,
const int len_atom_list,
PHPYCONST double (*r_carts)[3][3],
const int * permutations,
const int * map_atoms,
const int * map_syms,
const int num_rot,
const int num_pos);
static int compute_permutation(int * rot_atom,
PHPYCONST double lat[3][3],
PHPYCONST double (*pos)[3],
@ -128,6 +139,8 @@ static PyMethodDef _phonopy_methods[] = {
{"distribute_fc2", py_distribute_fc2, METH_VARARGS, "Distribute force constants"},
{"distribute_fc2_all", py_distribute_fc2_all, METH_VARARGS,
"Distribute force constants for all atoms in atom_list"},
{"distribute_fc2_with_mappings", py_distribute_fc2_with_mappings, METH_VARARGS,
"Distribute force constants for all atoms in atom_list using precomputed symmetry mappings."},
{"compute_permutation", py_compute_permutation, METH_VARARGS,
"Compute indices of original points in a set of rotated points."},
{"neighboring_grid_points", py_thm_neighboring_grid_points,
@ -220,6 +233,8 @@ static PyObject * py_compute_permutation(PyObject *self, PyObject *args)
double (*rot_pos)[3];
int num_pos;
int is_found;
if (!PyArg_ParseTuple(args, "OOOOd",
&permutation,
&lattice,
@ -235,12 +250,12 @@ static PyObject * py_compute_permutation(PyObject *self, PyObject *args)
rot_pos = (double(*)[3])PyArray_DATA(permuted_positions);
num_pos = PyArray_DIMS(positions)[0];
int is_found = compute_permutation(rot_atoms,
lat,
pos,
rot_pos,
num_pos,
symprec);
is_found = compute_permutation(rot_atoms,
lat,
pos,
rot_pos,
num_pos,
symprec);
return Py_BuildValue("i", is_found);
}
@ -686,12 +701,12 @@ static int compute_permutation(int * rot_atom,
rot_atom[i] = -1;
}
// optimization: Iterate primarily by pos instead of rot_pos.
// (find where 0 belongs in rot_atom, then where 1 belongs, etc.)
// Then track the first unassigned index.
//
// This works best if the permutation is close to the identity.
// (more specifically, if the max value of 'rot_atom[i] - i' is small)
/* optimization: Iterate primarily by pos instead of rot_pos. */
/* (find where 0 belongs in rot_atom, then where 1 belongs, etc.) */
/* Then track the first unassigned index. */
/* */
/* This works best if the permutation is close to the identity. */
/* (more specifically, if the max value of 'rot_atom[i] - i' is small) */
search_start = 0;
for (i = 0; i < num_pos; i++) {
while (rot_atom[search_start] >= 0) {
@ -878,6 +893,73 @@ static PyObject * py_distribute_fc2_all(PyObject *self, PyObject *args)
Py_RETURN_NONE;
}
static PyObject * py_distribute_fc2_with_mappings(PyObject *self, PyObject *args)
{
PyArrayObject* force_constants_py;
PyArrayObject* permutations_py;
PyArrayObject* map_atoms_py;
PyArrayObject* map_syms_py;
PyArrayObject* atom_list_py;
PyArrayObject* rotations_cart_py;
double (*r_carts)[3][3];
double (*fc2)[3][3];
int *permutations;
int *map_atoms;
int *map_syms;
int *atom_list;
int num_pos, num_rot, len_atom_list;
if (!PyArg_ParseTuple(args, "OOOOOO",
&force_constants_py,
&atom_list_py,
&rotations_cart_py,
&permutations_py,
&map_atoms_py,
&map_syms_py)) {
return NULL;
}
fc2 = (double(*)[3][3])PyArray_DATA(force_constants_py);
atom_list = (int*)PyArray_DATA(atom_list_py);
len_atom_list = PyArray_DIMS(atom_list_py)[0];
permutations = (int*)PyArray_DATA(permutations_py);
map_atoms = (int*)PyArray_DATA(map_atoms_py);
map_syms = (int*)PyArray_DATA(map_syms_py);
r_carts = (double(*)[3][3])PyArray_DATA(rotations_cart_py);
num_rot = PyArray_DIMS(permutations_py)[0];
num_pos = PyArray_DIMS(permutations_py)[1];
if (PyArray_NDIM(map_atoms_py) != 1 || PyArray_DIMS(map_atoms_py)[0] != num_pos)
{
PyErr_SetString(PyExc_ValueError, "wrong shape for map_atoms");
return NULL;
}
if (PyArray_NDIM(map_syms_py) != 1 || PyArray_DIMS(map_syms_py)[0] != num_pos)
{
PyErr_SetString(PyExc_ValueError, "wrong shape for map_syms");
return NULL;
}
if (PyArray_DIMS(rotations_cart_py)[0] != num_rot)
{
PyErr_SetString(PyExc_ValueError, "permutations and rotations are different length");
return NULL;
}
distribute_fc2_with_mappings(fc2,
atom_list,
len_atom_list,
r_carts,
permutations,
map_atoms,
map_syms,
num_rot,
num_pos);
Py_RETURN_NONE;
}
static PyObject *py_thm_neighboring_grid_points(PyObject *self, PyObject *args)
{
PyArrayObject* relative_grid_points_py;
@ -1272,6 +1354,59 @@ static int distribute_fc2(double *fc2,
return is_found;
}
/* Distributes all force constants using precomputed data about symmetry mappings. */
static void distribute_fc2_with_mappings(double (*fc2)[3][3], /* shape[num_pos][num_pos] */
const int * atom_list,
const int len_atom_list,
PHPYCONST double (*r_carts)[3][3], /* shape[num_rot] */
const int * permutations, /* shape[num_rot][num_pos] */
const int * map_atoms, /* shape [num_pos] */
const int * map_syms, /* shape [num_pos] */
const int num_rot,
const int num_pos)
{
int i, j, k, l, m;
int atom_todo, atom_done, atom_other;
int sym_index;
double (*fc2_done)[3];
double (*fc2_todo)[3];
double (*r_cart)[3];
const int * permutation;
for (i = 0; i < len_atom_list; i++) {
/* look up how this atom maps into the done list. */
atom_todo = atom_list[i];
atom_done = map_atoms[atom_todo];
sym_index = map_syms[atom_todo];
/* skip the atoms in the done list, */
/* which are easily identified because they map to themselves. */
if (atom_todo == atom_done) {
continue;
}
/* look up information about the rotation */
r_cart = r_carts[sym_index];
permutation = &permutations[sym_index * num_pos]; /* shape[num_pos] */
/* distribute terms from atom_done to atom_todo */
for (atom_other = 0; atom_other < num_pos; atom_other++) {
fc2_done = fc2[atom_done * num_pos + permutation[atom_other]];
fc2_todo = fc2[atom_todo * num_pos + atom_other];
for (j = 0; j < 3; j++) {
for (k = 0; k < 3; k++) {
for (l = 0; l < 3; l++) {
for (m = 0; m < 3; m++) {
/* P' = R^-1 P R */
fc2_todo[j][k] += r_cart[l][j] * r_cart[m][k] * fc2_done[l][m];
}
}
}
}
}
}
}
static int check_overlap(PHPYCONST double (*pos)[3],
const int num_pos,
const double pos_orig[3],

91
phonopy/harmonic/force_constants.py
View File

@ -145,11 +145,33 @@ def distribute_force_constants(force_constants,
atom_list,
atom_list_done,
lattice, # column vectors
positions,
rotations,
trans,
positions, # scaled (fractional)
rotations, # scaled (fractional)
trans, # scaled (fractional)
symprec):
if True:
permutations = _compute_all_sg_permutations(positions,
rotations,
trans,
lattice,
symprec)
map_atoms, map_syms = _get_sym_mappings_from_permutations(
permutations, atom_list_done)
rots_cartesian = np.array([similarity_transformation(lattice, r)
for r in rotations],
dtype='double', order='C')
import phonopy._phonopy as phonoc
phonoc.distribute_fc2_with_mappings(force_constants,
np.array(atom_list, dtype='intc'),
rots_cartesian,
permutations,
np.array(map_atoms, dtype='intc'),
np.array(map_syms, dtype='intc'))
elif True:
rots_cartesian = np.array([similarity_transformation(lattice, r)
for r in rotations],
dtype='double', order='C')
@ -807,6 +829,25 @@ def _get_atom_mapping_by_symmetry(atom_list_done,
print("or something wrong (e.g. crystal structure does not match).")
raise ValueError
# Compute a permutation for every space group operation.
# See '_compute_permutation_for_rotation' for more info.
#
# Output has shape (num_rot, num_pos)
def _compute_all_sg_permutations(positions, # scaled positions
rotations, # scaled
translations, # scaled
lattice, # column vectors
symprec):
out = [] # Finally the shape is fixed as (num_sym, num_pos_of_supercell).
for (sym, t) in zip(rotations, translations):
rotated_positions = np.dot(positions, sym.T) + t
out.append(_compute_permutation_for_rotation(positions,
rotated_positions,
lattice,
symprec))
return np.array(out, dtype='intc', order='C')
# Get the overall permutation such that
#
# positions_a[perm[i]] == positions_b[i] (modulo the lattice)
@ -831,7 +872,7 @@ def _compute_permutation_for_rotation(positions_a, # scaled positions
def sort_by_lattice_distance(fracs):
carts = np.dot(fracs - np.rint(fracs), lattice.T)
perm = np.argsort(np.sum(carts**2, axis=1))
sorted_fracs = fracs[perm]
sorted_fracs = np.array(fracs[perm], dtype='double', order='C')
return perm, sorted_fracs
(perm_a, sorted_a) = sort_by_lattice_distance(positions_a)
@ -894,3 +935,45 @@ def _compute_permutation_c(positions_a, # scaled positions
permutation_error()
return permutation
# This can be thought of as computing 'map_atom_disp' and 'map_sym' for all atoms,
# except done using permutations instead of by computing overlaps.
#
# Input:
# * permutations, shape [num_rot][num_pos]
# * atom_list_done
#
# Output:
# * map_atoms, shape [num_pos].
# Maps each atom in the full structure to its equivalent atom in atom_list_done.
# (each entry will be an integer found in atom_list_done)
#
# * map_syms, shape [num_pos].
# For each atom, provides the index of a rotation that maps it into atom_list_done.
# (there might be more than one such rotation, but only one will be returned)
# (each entry will be an integer 0 <= i < num_rot)
def _get_sym_mappings_from_permutations(permutations,
atom_list_done):
assert permutations.ndim == 2
num_pos = permutations.shape[1]
# filled with -1
map_atoms = np.zeros((num_pos,), dtype='intc') - 1
map_syms = np.zeros((num_pos,), dtype='intc') - 1
atom_list_done = set(atom_list_done)
for atom_todo in range(num_pos):
for (sym_index, permutation) in enumerate(permutations):
if permutation[atom_todo] in atom_list_done:
map_atoms[atom_todo] = permutation[atom_todo]
map_syms[atom_todo] = sym_index
break
else:
print("Input forces are not enough to calculate force constants,")
print("or something wrong (e.g. crystal structure does not match).")
raise ValueError
assert set(map_atoms) & set(atom_list_done) == set(map_atoms)
assert -1 not in map_atoms
assert -1 not in map_syms
return map_atoms, map_syms