Make an option to choose either full fc or compact fc

This commit is contained in:
Atsushi Togo 2018-03-02 15:41:35 +09:00
parent c67d0a543b
commit 6bb3cd46e7
6 changed files with 48 additions and 27 deletions

View File

@ -242,6 +242,7 @@ static PyObject * py_transform_dynmat_to_fc(PyObject *self, PyObject *args)
PyArrayObject* py_multiplicities;
PyArrayObject* py_masses;
PyArrayObject* py_s2pp_map;
PyArrayObject* py_fc_index_map;
double* fc;
double* dm;
@ -250,17 +251,19 @@ static PyObject * py_transform_dynmat_to_fc(PyObject *self, PyObject *args)
double* masses;
int* multiplicities;
int* s2pp_map;
int* fc_index_map;
int num_patom;
int num_satom;
if (!PyArg_ParseTuple(args, "OOOOOOO",
if (!PyArg_ParseTuple(args, "OOOOOOOO",
&py_force_constants,
&py_dynamical_matrices,
&py_commensurate_points,
&py_shortest_vectors,
&py_multiplicities,
&py_masses,
&py_s2pp_map)) {
&py_s2pp_map,
&py_fc_index_map)) {
return NULL;
}
@ -271,6 +274,7 @@ static PyObject * py_transform_dynmat_to_fc(PyObject *self, PyObject *args)
masses = (double*)PyArray_DATA(py_masses);
multiplicities = (int*)PyArray_DATA(py_multiplicities);
s2pp_map = (int*)PyArray_DATA(py_s2pp_map);
fc_index_map = (int*)PyArray_DATA(py_fc_index_map);
num_patom = PyArray_DIMS(py_multiplicities)[1];
num_satom = PyArray_DIMS(py_multiplicities)[0];
@ -281,6 +285,7 @@ static PyObject * py_transform_dynmat_to_fc(PyObject *self, PyObject *args)
multiplicities,
masses,
s2pp_map,
fc_index_map,
num_patom,
num_satom);

View File

@ -328,6 +328,7 @@ void dym_transform_dynmat_to_fc(double *fc,
const int *multiplicities,
const double *masses,
const int *s2pp_map,
const int *fc_index_map,
const int num_patom,
const int num_satom)
{
@ -361,7 +362,7 @@ void dym_transform_dynmat_to_fc(double *fc,
for (m = 0; m < 3; m++) {
adrs = k * num_patom * num_patom * 18 + i * num_patom * 18 +
l * num_patom * 6 + s2pp_map[j] * 6 + m * 2;
fc[i * num_satom * 9 + j * 9 + l * 3 + m] +=
fc[fc_index_map[i] * num_satom * 9 + j * 9 + l * 3 + m] +=
(dm[adrs] * cos_phase - dm[adrs + 1] * sin_phase) * coef;
}
}

View File

@ -88,6 +88,7 @@ void dym_transform_dynmat_to_fc(double *fc,
const int *multiplicities,
const double *masses,
const int *s2pp_map,
const int *fc_index_map,
const int num_patom,
const int num_satom);

View File

@ -349,7 +349,7 @@ class Phonopy(object):
self._supercell)
self.set_displacement_dataset(displacement_dataset)
def produce_force_constants(self,
def produce_force_constants(self,
forces=None,
calculate_full_force_constants=True,
computation_algorithm="svd"):
@ -361,16 +361,17 @@ class Phonopy(object):
if 'forces' not in disp:
return False
if calculate_full_force_constants:
# if calculate_full_force_constants:
if True:
self._run_force_constants_from_forces(
decimals=self._force_constants_decimals,
computation_algorithm=computation_algorithm)
else:
p2s_map = self._primitive.get_primitive_to_supercell_map()
self._run_force_constants_from_forces(
distributed_atom_list=p2s_map,
decimals=self._force_constants_decimals,
computation_algorithm=computation_algorithm)
# else:
# p2s_map = self._primitive.get_primitive_to_supercell_map()
# self._run_force_constants_from_forces(
# distributed_atom_list=p2s_map,
# decimals=self._force_constants_decimals,
# computation_algorithm=computation_algorithm)
self._set_dynamical_matrix()

View File

@ -443,8 +443,10 @@ class DynamicalMatrixNAC(DynamicalMatrix):
self._dynamical_matrix += dm_dd
def _set_Gonze_force_constants(self):
fc_shape = self._force_constants.shape
d2f = DynmatToForceConstants(self._pcell,
self._scell,
is_full_fc=(fc_shape[0] == fc_shape[1]),
symprec=self._symprec)
dynmat = []
num_q = len(d2f.get_commensurate_points())

View File

@ -54,6 +54,7 @@ class DynmatToForceConstants(object):
supercell,
frequencies=None,
eigenvectors=None,
is_full_fc=False,
symprec=1e-5):
self._primitive = primitive
self._supercell = supercell
@ -65,13 +66,13 @@ class DynmatToForceConstants(object):
self._dynmat = None
n_s = self._supercell.get_number_of_atoms()
n_p = self._primitive.get_number_of_atoms()
# Full fc
# self._force_constants = np.zeros((n_p, n_s, 3, 3),
# dtype='double', order='C')
# Compact fc
self._force_constants = np.zeros((n_p, n_s, 3, 3),
dtype='double', order='C')
itemsize = self._force_constants.itemsize
if is_full_fc:
fc_shape = (n_s, n_s, 3, 3)
else:
fc_shape = (n_p, n_s, 3, 3)
self._fc = np.zeros(fc_shape, dtype='double', order='C')
itemsize = self._fc.itemsize
self._dtype_complex = ("c%d" % (itemsize * 2))
if frequencies is not None and eigenvectors is not None:
@ -83,7 +84,7 @@ class DynmatToForceConstants(object):
# self._distribute_force_constants()
def get_force_constants(self):
return self._force_constants
return self._fc
def get_commensurate_points(self):
return self._commensurate_points
@ -114,6 +115,9 @@ class DynmatToForceConstants(object):
except ImportError:
self._py_inverse_transformation()
if self._fc.shape[0] == self._fc.shape[1]:
self._distribute_force_constants()
def _c_inverse_transformation(self):
import phonopy._phonopy as phonoc
@ -121,20 +125,26 @@ class DynmatToForceConstants(object):
p2p = self._primitive.get_primitive_to_primitive_map()
s2pp = np.array([p2p[i] for i in s2p], dtype='intc')
phonoc.transform_dynmat_to_fc(self._force_constants,
if self._fc.shape[0] == self._fc.shape[1]:
fc_index_map = self._primitive.get_primitive_to_supercell_map()
else:
fc_index_map = np.arange(self._fc.shape[0], dtype='intc')
phonoc.transform_dynmat_to_fc(self._fc,
self._dynmat.view(dtype='double'),
self._commensurate_points,
self._shortest_vectors,
self._multiplicity,
self._primitive.get_masses(),
s2pp)
s2pp,
fc_index_map)
def _py_inverse_transformation(self):
s2p = self._primitive.get_supercell_to_primitive_map()
p2s = self._primitive.get_primitive_to_supercell_map()
p2p = self._primitive.get_primitive_to_primitive_map()
fc = self._force_constants
fc = self._fc
m = self._primitive.get_masses()
N = (self._supercell.get_number_of_atoms() /
self._primitive.get_number_of_atoms())
@ -142,10 +152,11 @@ class DynmatToForceConstants(object):
for p_i, s_i in enumerate(p2s):
for s_j, p_j in enumerate([p2p[i] for i in s2p]):
coef = np.sqrt(m[p_i] * m[p_j]) / N
# Full fc
# fc[s_i, s_j] = self._sum_q(p_i, s_j, p_j) * coef
# Compact fc
fc[p_i, s_j] = self._sum_q(p_i, s_j, p_j) * coef
fc_elem = self._sum_q(p_i, s_j, p_j) * coef
if fc.shape[0] == fc.shape[1]:
fc[s_i, s_j] = fc_elem
else:
fc[p_i, s_j] = fc_elem
def _sum_q(self, p_i, s_j, p_j):
multi = self._multiplicity[s_j, p_i]
@ -169,7 +180,7 @@ class DynmatToForceConstants(object):
dtype='double', order='C')
rotations = np.array([np.eye(3, dtype='intc')] * len(trans),
dtype='intc', order='C')
distribute_force_constants(self._force_constants,
distribute_force_constants(self._fc,
range(self._supercell.get_number_of_atoms()),
p2s,
lattice,