mirror of https://github.com/phonopy/phonopy.git
Change matrix arrangement to obtain fc2
This commit is contained in:
Atsushi Togo
parent
fca8c3d7b7
commit
837320a0a7
96
c/_phonopy.c
96
c/_phonopy.c
|
|
@ -44,7 +44,6 @@ static PyObject * py_get_dynamical_matrix(PyObject *self, PyObject *args);
|
||||||
static PyObject * py_get_nac_dynamical_matrix(PyObject *self, PyObject *args);
|
static PyObject * py_get_nac_dynamical_matrix(PyObject *self, PyObject *args);
|
||||||
static PyObject * py_get_thermal_properties(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(PyObject *self, PyObject *args);
|
||||||
static PyObject * py_get_rotated_forces(PyObject *self, PyObject *args);
|
|
||||||
|
|
||||||
static double get_free_energy_omega(const double temperature,
|
static double get_free_energy_omega(const double temperature,
|
||||||
const double omega);
|
const double omega);
|
||||||
|
|
@ -62,14 +61,6 @@ static int distribute_fc2(double * fc2,
|
||||||
const int * r,
|
const int * r,
|
||||||
const double * t,
|
const double * t,
|
||||||
const double symprec);
|
const double symprec);
|
||||||
static int get_rotated_forces(double * rotated_forces,
|
|
||||||
const double * pos,
|
|
||||||
const int num_pos,
|
|
||||||
const int atom_number,
|
|
||||||
const double * f,
|
|
||||||
const int * r,
|
|
||||||
const int num_rot,
|
|
||||||
const double symprec);
|
|
||||||
static int nint(const double a);
|
static int nint(const double a);
|
||||||
|
|
||||||
static PyMethodDef functions[] = {
|
static PyMethodDef functions[] = {
|
||||||
|
|
@ -77,7 +68,6 @@ static PyMethodDef functions[] = {
|
||||||
{"nac_dynamical_matrix", py_get_nac_dynamical_matrix, METH_VARARGS, "NAC dynamical matrix"},
|
{"nac_dynamical_matrix", py_get_nac_dynamical_matrix, METH_VARARGS, "NAC dynamical matrix"},
|
||||||
{"thermal_properties", py_get_thermal_properties, METH_VARARGS, "Thermal properties"},
|
{"thermal_properties", py_get_thermal_properties, METH_VARARGS, "Thermal properties"},
|
||||||
{"distribute_fc2", py_distribute_fc2, METH_VARARGS, "Distribute force constants"},
|
{"distribute_fc2", py_distribute_fc2, METH_VARARGS, "Distribute force constants"},
|
||||||
{"rotated_forces", py_get_rotated_forces, METH_VARARGS, "Rotate forces following site-symmetry"},
|
|
||||||
{NULL, NULL, 0, NULL}
|
{NULL, NULL, 0, NULL}
|
||||||
};
|
};
|
||||||
|
|
||||||
|
|
@ -400,92 +390,6 @@ static int distribute_fc2(double * fc2,
|
||||||
return is_found;
|
return is_found;
|
||||||
}
|
}
|
||||||
|
|
||||||
static PyObject * py_get_rotated_forces(PyObject *self, PyObject *args)
|
|
||||||
{
|
|
||||||
PyArrayObject* rotated_forces;
|
|
||||||
PyArrayObject* positions;
|
|
||||||
PyArrayObject* rotations;
|
|
||||||
PyArrayObject* forces;
|
|
||||||
int atom_number;
|
|
||||||
double symprec;
|
|
||||||
|
|
||||||
if (!PyArg_ParseTuple(args, "OOiOOd",
|
|
||||||
&rotated_forces,
|
|
||||||
&positions,
|
|
||||||
&atom_number,
|
|
||||||
&forces,
|
|
||||||
&rotations,
|
|
||||||
&symprec)) {
|
|
||||||
return NULL;
|
|
||||||
}
|
|
||||||
|
|
||||||
const int* r = (int*)rotations->data;
|
|
||||||
const int num_rot = rotations->dimensions[0];
|
|
||||||
const double* pos = (double*)positions->data;
|
|
||||||
const int num_pos = positions->dimensions[0];
|
|
||||||
const double* f = (double*)forces->data;
|
|
||||||
double* rot_forces = (double*)rotated_forces->data;
|
|
||||||
|
|
||||||
get_rotated_forces(rot_forces,
|
|
||||||
pos,
|
|
||||||
num_pos,
|
|
||||||
atom_number,
|
|
||||||
f,
|
|
||||||
r,
|
|
||||||
num_rot,
|
|
||||||
symprec);
|
|
||||||
|
|
||||||
Py_RETURN_NONE;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
static int get_rotated_forces(double * rotated_forces,
|
|
||||||
const double * pos,
|
|
||||||
const int num_pos,
|
|
||||||
const int atom_number,
|
|
||||||
const double * f,
|
|
||||||
const int * r,
|
|
||||||
const int num_rot,
|
|
||||||
const double symprec)
|
|
||||||
{
|
|
||||||
int i, j, k, is_found;
|
|
||||||
double rot_pos[3], diff[3];
|
|
||||||
|
|
||||||
#pragma omp parallel for private(j, k, rot_pos, diff, is_found)
|
|
||||||
for (i = 0; i < num_rot; i++) {
|
|
||||||
for (j = 0; j < 3; j++) {
|
|
||||||
rot_pos[j] = 0;
|
|
||||||
for (k = 0; k < 3; k++) {
|
|
||||||
rot_pos[j] += r[i * 9 + j * 3 + k] * pos[atom_number * 3 + k];
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
for (j = 0; j < num_pos; j++) {
|
|
||||||
is_found = 1;
|
|
||||||
for (k = 0; k < 3; k++) {
|
|
||||||
diff[k] = pos[j * 3 + k] - rot_pos[k];
|
|
||||||
diff[k] -= nint(diff[k]);
|
|
||||||
if (fabs(diff[k]) > symprec) {
|
|
||||||
is_found = 0;
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
if (is_found) {
|
|
||||||
for (k = 0; k < 3; k++) {
|
|
||||||
rotated_forces[i * 3 + k] = f[j * 3 + k];
|
|
||||||
}
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
if (! is_found) {
|
|
||||||
printf("Phonopy encounted symmetry problem (c)\n");
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
return 1;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
static int nint(const double a)
|
static int nint(const double a)
|
||||||
{
|
{
|
||||||
if (a < 0.0)
|
if (a < 0.0)
|
||||||
|
|
|
||||||
|
|
@ -296,114 +296,71 @@ def get_force_constants_disps(force_constants,
|
||||||
symprec = symmetry.get_symmetry_tolerance()
|
symprec = symmetry.get_symmetry_tolerance()
|
||||||
disp_atom_list = np.unique(
|
disp_atom_list = np.unique(
|
||||||
[forces.get_atom_number() for forces in set_of_forces])
|
[forces.get_atom_number() for forces in set_of_forces])
|
||||||
|
lat = supercell.get_cell().T
|
||||||
|
|
||||||
for disp_atom_number in disp_atom_list:
|
for disp_atom_number in disp_atom_list:
|
||||||
site_symmetry = symmetry.get_site_symmetry(disp_atom_number)
|
site_symmetry = symmetry.get_site_symmetry(disp_atom_number)
|
||||||
symmetry_matrices = get_symmetry_matrices(supercell, site_symmetry)
|
positions = supercell.get_scaled_positions()
|
||||||
rot_disps = []
|
pos_center = positions[disp_atom_number].copy()
|
||||||
row_forces = []
|
positions -= pos_center
|
||||||
|
rot_map_syms = get_positions_sent_by_rot_inv(positions,
|
||||||
|
site_symmetry,
|
||||||
|
symprec)
|
||||||
|
site_sym_cart = [similarity_transformation(lat, sym)
|
||||||
|
for sym in site_symmetry]
|
||||||
|
|
||||||
|
rot_disps = []
|
||||||
|
raw_forces = []
|
||||||
|
|
||||||
# Bind several displacements of a displaced atom
|
|
||||||
# with symmetry operations
|
|
||||||
for forces in set_of_forces:
|
for forces in set_of_forces:
|
||||||
if not forces.get_atom_number() == disp_atom_number:
|
if forces.get_atom_number() != disp_atom_number:
|
||||||
continue
|
continue
|
||||||
|
|
||||||
displacement = forces.get_displacement()
|
u = forces.get_displacement()
|
||||||
# Displacement * Rotation (U * A)
|
rot_disps.append([np.dot(sym, u) for sym in site_sym_cart])
|
||||||
rot_disps.append(get_rotated_displacements(displacement,
|
raw_forces.append(forces.get_forces())
|
||||||
symmetry_matrices))
|
|
||||||
row_forces.append(forces.get_forces())
|
|
||||||
|
|
||||||
# Bind forces for several displacements and symmetry
|
rot_disps = np.reshape(rot_disps, (-1, 3))
|
||||||
# operations of a displaced atom
|
inv_displacements = np.linalg.pinv(rot_disps)
|
||||||
rot_disps = np.array(rot_disps).reshape(-1, 9)
|
|
||||||
inv = np.linalg.pinv(rot_disps)
|
|
||||||
|
|
||||||
solve_force_constants_disps(force_constants,
|
for i in range(supercell.get_number_of_atoms()):
|
||||||
supercell,
|
combined_forces = []
|
||||||
disp_atom_number,
|
for forces in raw_forces:
|
||||||
site_symmetry,
|
combined_forces.append(
|
||||||
row_forces,
|
get_rotated_forces(forces[rot_map_syms[:, i]],
|
||||||
inv,
|
site_sym_cart))
|
||||||
symprec)
|
|
||||||
|
combined_forces = np.reshape(combined_forces, (-1, 3))
|
||||||
|
force_constants[disp_atom_number, i] = -np.dot(
|
||||||
|
inv_displacements, combined_forces)
|
||||||
|
|
||||||
return disp_atom_list
|
return disp_atom_list
|
||||||
|
|
||||||
|
def get_positions_sent_by_rot_inv(positions,
|
||||||
def solve_force_constants_disps(force_constants,
|
site_symmetry,
|
||||||
supercell,
|
symprec):
|
||||||
disp_atom_number,
|
rot_map_syms = []
|
||||||
site_symmetry,
|
|
||||||
sets_of_forces,
|
|
||||||
inv_displacements,
|
|
||||||
symprec):
|
|
||||||
for i in range(supercell.get_number_of_atoms()):
|
|
||||||
# Shift positions according to set disp_atom_number is at origin
|
|
||||||
positions = supercell.get_scaled_positions()
|
|
||||||
pos_center = positions[disp_atom_number].copy()
|
|
||||||
positions -= pos_center
|
|
||||||
# Combined forces (F)
|
|
||||||
combined_forces = []
|
|
||||||
try:
|
|
||||||
import phonopy._phonopy as phonoc
|
|
||||||
for forces in sets_of_forces:
|
|
||||||
rotated_forces = np.zeros(len(site_symmetry) * 3,
|
|
||||||
dtype='double')
|
|
||||||
phonoc.rotated_forces(rotated_forces,
|
|
||||||
positions,
|
|
||||||
i,
|
|
||||||
forces,
|
|
||||||
site_symmetry,
|
|
||||||
symprec)
|
|
||||||
combined_forces.append(rotated_forces)
|
|
||||||
except ImportError:
|
|
||||||
for forces in sets_of_forces:
|
|
||||||
combined_forces.append(
|
|
||||||
get_rotated_forces(positions,
|
|
||||||
i,
|
|
||||||
forces,
|
|
||||||
site_symmetry,
|
|
||||||
symprec))
|
|
||||||
|
|
||||||
combined_forces = np.array(combined_forces).reshape(-1, 1)
|
|
||||||
# Force constant (P) = -(U * A)^-1 * (F)
|
|
||||||
force_constants[disp_atom_number,i] = \
|
|
||||||
-np.dot(inv_displacements, combined_forces).reshape(3, 3)
|
|
||||||
|
|
||||||
|
|
||||||
def get_rotated_forces(positions,
|
|
||||||
atom_number,
|
|
||||||
forces,
|
|
||||||
site_symmetry,
|
|
||||||
symprec=1e-5):
|
|
||||||
"""
|
|
||||||
Pack forces on atoms translated by site symmetry
|
|
||||||
|
|
||||||
The relation:
|
|
||||||
R [ F(r) ] = F(R.r)
|
|
||||||
where R is a rotation cetering at displaced atom.
|
|
||||||
(This is not the transformation of a function,
|
|
||||||
but just the rotation of force vector at r.)
|
|
||||||
"""
|
|
||||||
rot_forces = []
|
|
||||||
for sym in site_symmetry:
|
for sym in site_symmetry:
|
||||||
rot_pos = np.dot(positions[atom_number], sym.T)
|
rot_map = np.zeros(len(positions), dtype='intc')
|
||||||
|
rot_pos = np.dot(positions, sym.T)
|
||||||
is_found = False
|
is_found = False
|
||||||
for i, p in enumerate(positions):
|
for i, rot_pos_i in enumerate(rot_pos):
|
||||||
diff = p - rot_pos
|
diff = positions - rot_pos_i
|
||||||
if (abs(diff - diff.round()) < symprec).all():
|
j = np.nonzero((
|
||||||
rot_forces.append(forces[i])
|
np.abs(diff - np.rint(diff)) < symprec).all(axis=1))[0]
|
||||||
is_found = True
|
rot_map[j] = i
|
||||||
break
|
|
||||||
|
|
||||||
if not is_found:
|
rot_map_syms.append(rot_map)
|
||||||
print "Phonopy encontered symmetry problem"
|
|
||||||
|
return np.intc(rot_map_syms)
|
||||||
|
|
||||||
|
def get_rotated_forces(forces_syms, site_sym_cart):
|
||||||
|
rot_forces = []
|
||||||
|
for forces, sym_cart in zip(forces_syms, site_sym_cart):
|
||||||
|
rot_forces.append(np.dot(forces, sym_cart.T))
|
||||||
|
|
||||||
return rot_forces
|
return rot_forces
|
||||||
|
|
||||||
|
|
||||||
def set_tensor_symmetry(force_constants, supercell, symmetry):
|
def set_tensor_symmetry(force_constants, supercell, symmetry):
|
||||||
"""
|
"""
|
||||||
Full force constants are symmetrized using crystal symmetry.
|
Full force constants are symmetrized using crystal symmetry.
|
||||||
|
|
@ -537,68 +494,7 @@ def force_constants_log(force_constants):
|
||||||
i + 1, j + 1, v[0], v[1], v[2])
|
i + 1, j + 1, v[0], v[1], v[2])
|
||||||
|
|
||||||
|
|
||||||
# Shared functions to calculate force constant
|
|
||||||
def get_symmetry_matrices(cell, site_symmetry):
|
|
||||||
"""
|
|
||||||
Transformation of 2nd order force constant
|
|
||||||
|
|
||||||
In the phonopy implementation (Cartesian coords.)
|
|
||||||
|
|
||||||
(R.F)^T = -(R.U)^T Psi' --> F^T = -U^T.R^T.Psi'.R
|
|
||||||
Therefore,
|
|
||||||
Psi = R^T.Psi'.R --> Psi' = R.Psi.R^T
|
|
||||||
|
|
||||||
The symmetrical relation between Psi and Psi' can be represented
|
|
||||||
by a 9x9 matrix. What we want is transformation matrix A defined
|
|
||||||
by
|
|
||||||
|
|
||||||
P' = A.P
|
|
||||||
|
|
||||||
where P' and P are the 9x1 matrices and A is the 9x9 matrices.
|
|
||||||
"""
|
|
||||||
matrices = []
|
|
||||||
for reduced_rot in site_symmetry:
|
|
||||||
mat = []
|
|
||||||
rot = similarity_transformation(cell.get_cell().T, reduced_rot)
|
|
||||||
for i in range(3):
|
|
||||||
for j in range(3):
|
|
||||||
psi = []
|
|
||||||
for k in range(3):
|
|
||||||
for l in range(3):
|
|
||||||
psi.append(rot[i, k] * rot[j, l])
|
|
||||||
mat.append(psi)
|
|
||||||
matrices.append(mat)
|
|
||||||
return np.array(matrices)
|
|
||||||
|
|
||||||
def similarity_transformation(rot, mat):
|
def similarity_transformation(rot, mat):
|
||||||
""" R x M x R^-1 """
|
""" R x M x R^-1 """
|
||||||
return np.dot(rot, np.dot(mat, np.linalg.inv(rot)))
|
return np.dot(rot, np.dot(mat, np.linalg.inv(rot)))
|
||||||
|
|
||||||
def get_rotated_displacements(displacement, symmetry_matrices):
|
|
||||||
"""
|
|
||||||
U x A
|
|
||||||
[123456789]
|
|
||||||
[2 ]
|
|
||||||
[3 ]
|
|
||||||
[ d_x 0 0 d_y 0 0 d_z 0 0 ] [4 ]
|
|
||||||
U [ 0 d_x 0 0 d_y 0 0 d_z 0 ] [5 A ]
|
|
||||||
[ 0 0 d_x 0 0 d_y 0 0 d_z ] [6 ]
|
|
||||||
[7 ]
|
|
||||||
[8 ]
|
|
||||||
[9 ]
|
|
||||||
"""
|
|
||||||
rot_disps = []
|
|
||||||
for sym in symmetry_matrices:
|
|
||||||
rot_disps.append(np.dot(expand_displacement(displacement), sym))
|
|
||||||
return rot_disps
|
|
||||||
|
|
||||||
def expand_displacement(displacement):
|
|
||||||
"""
|
|
||||||
[ d_x 0 0 d_y 0 0 d_z 0 0 ]
|
|
||||||
[ 0 d_x 0 0 d_y 0 0 d_z 0 ]
|
|
||||||
[ 0 0 d_x 0 0 d_y 0 0 d_z ]
|
|
||||||
"""
|
|
||||||
d = displacement
|
|
||||||
disp = np.hstack((np.eye(3)*d[0], np.eye(3)*d[1], np.eye(3)*d[2]))
|
|
||||||
return disp
|
|
||||||
|
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue