From 6bb3cd46e7c88c5a73b88e7b2573c8158fb0d839 Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Fri, 2 Mar 2018 15:41:35 +0900 Subject: [PATCH] Make an option to choose either full fc or compact fc --- c/_phonopy.c | 9 ++++-- c/harmonic/dynmat.c | 3 +- c/harmonic_h/dynmat.h | 1 + phonopy/api_phonopy.py | 17 +++++------ phonopy/harmonic/dynamical_matrix.py | 2 ++ phonopy/harmonic/dynmat_to_fc.py | 43 +++++++++++++++++----------- 6 files changed, 48 insertions(+), 27 deletions(-) diff --git a/c/_phonopy.c b/c/_phonopy.c index 31f9a4b9..cb00e95e 100644 --- a/c/_phonopy.c +++ b/c/_phonopy.c @@ -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); diff --git a/c/harmonic/dynmat.c b/c/harmonic/dynmat.c index 956e5cac..99782b99 100644 --- a/c/harmonic/dynmat.c +++ b/c/harmonic/dynmat.c @@ -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; } } diff --git a/c/harmonic_h/dynmat.h b/c/harmonic_h/dynmat.h index 42d99710..8dae1139 100644 --- a/c/harmonic_h/dynmat.h +++ b/c/harmonic_h/dynmat.h @@ -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); diff --git a/phonopy/api_phonopy.py b/phonopy/api_phonopy.py index 2f38df9c..356d84b5 100644 --- a/phonopy/api_phonopy.py +++ b/phonopy/api_phonopy.py @@ -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() diff --git a/phonopy/harmonic/dynamical_matrix.py b/phonopy/harmonic/dynamical_matrix.py index 614e0e47..fbff77ea 100644 --- a/phonopy/harmonic/dynamical_matrix.py +++ b/phonopy/harmonic/dynamical_matrix.py @@ -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()) diff --git a/phonopy/harmonic/dynmat_to_fc.py b/phonopy/harmonic/dynmat_to_fc.py index cab72531..24e8836d 100644 --- a/phonopy/harmonic/dynmat_to_fc.py +++ b/phonopy/harmonic/dynmat_to_fc.py @@ -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,