diff --git a/anharmonic/force_fit/fc3.py b/anharmonic/force_fit/fc3.py index 51554863..5e44114f 100644 --- a/anharmonic/force_fit/fc3.py +++ b/anharmonic/force_fit/fc3.py @@ -79,8 +79,12 @@ class FC3Fit: rot_map_syms) fc = self._solve(rot_disps, rot_forces) fc2 = fc[:, 1:4, :].reshape((self._num_atom, 3, 3)) + fc2_2 = fc[:, 5:8, :].reshape((self._num_atom, 3, 3)) fc3 = fc[:, 7:16, :].reshape((self._num_atom, 3, 3, 3)) - self._fc3[first_atom_num, second_atom_num] = fc3 + fc3_21 = fc[:, 16:25, :].reshape((self._num_atom, 3, 3, 3)) + for i, j in list(np.ndindex(3, 3)): + self._fc3[first_atom_num, second_atom_num, :, i, j, :] = ( + fc3[:, i, j, :] + fc3_21[:, j, i, :]) / 2 def _solve(self, rot_disps, rot_forces): fc = [] @@ -104,10 +108,10 @@ class FC3Fit: for map_sym, rot_atom_num, sym in zip( rot_map_syms, rot_atom_map, site_syms_cart): for forces in sets_of_forces_u1[rot_atom_num]: - force_matrix_atom.append(np.dot(sym, - forces[map_sym[i]])) + force_matrix_atom.append( + np.dot(sym, forces[map_sym[i]])) force_matrix.append(force_matrix_atom) - return np.array(force_matrix, dtype='double') + return np.array(force_matrix, dtype='double', order='C') def _create_displacement_matrix(self, disp_pairs, @@ -116,6 +120,7 @@ class FC3Fit: rot_disp1s = [] rot_disp2s = [] rot_pair12 = [] + rot_pair21 = [] rot_pair11 = [] rot_pair22 = [] @@ -127,18 +132,19 @@ class FC3Fit: Su2 = np.dot(ssym_c, u2) rot_disp1s.append(Su1) rot_disp2s.append(Su2) - rot_pair12.append(np.outer(Su1, Su2).flatten()) - rot_pair11.append(np.outer(Su1, Su1).flatten()) - rot_pair22.append(np.outer(Su2, Su2).flatten()) + rot_pair12.append(np.outer(Su1, Su2).flatten() / 2) + rot_pair21.append(np.outer(Su2, Su1).flatten() / 2) + rot_pair11.append(np.outer(Su1, Su1).flatten() / 2) + rot_pair22.append(np.outer(Su2, Su2).flatten() / 2) ones = np.ones(len(rot_disp1s)).reshape((-1, 1)) return np.hstack((ones, rot_disp1s, rot_disp2s, - rot_pair12, rot_pair11, rot_pair22)) + rot_pair12, rot_pair21, rot_pair11, rot_pair22)) def _collect_disp_pairs_and_forces(self, dataset_1st): second_atom_nums = [x['number'] for x in dataset_1st['second_atoms']] - unique_second_atom_nums = np.unique(np.array(second_atom_nums)) + unique_second_atom_nums = np.unique(second_atom_nums) set_of_disps = [] sets_of_forces = [] @@ -181,10 +187,7 @@ class FC3Fit: disp_pairs = [] for second_atom_num in range(self._num_atom): - if set_of_disps[second_atom_num] is not None: - disp_pairs.append( - [[disp1, d] for d in set_of_disps[second_atom_num]]) - else: + if set_of_disps[second_atom_num] is None: (sets_of_forces_atom2, set_of_disps_atom2) = self._copy_dataset_2nd( second_atom_num, @@ -193,10 +196,12 @@ class FC3Fit: unique_second_atom_nums, reduced_site_sym, rot_map_syms) - disp_pairs.append( [[disp1, d] for d in set_of_disps_atom2]) sets_of_forces[second_atom_num] = sets_of_forces_atom2 + else: + disp_pairs.append( + [[disp1, d] for d in set_of_disps[second_atom_num]]) return disp_pairs, sets_of_forces diff --git a/anharmonic/phonon3/displacement_fc3.py b/anharmonic/phonon3/displacement_fc3.py index 01c1e5d8..1745e65c 100644 --- a/anharmonic/phonon3/displacement_fc3.py +++ b/anharmonic/phonon3/displacement_fc3.py @@ -159,7 +159,7 @@ def get_reduced_site_symmetry(site_sym, direction, symprec=1e-5): for rot in site_sym: if (abs(direction - np.dot(direction, rot.T)) < symprec).all(): reduced_site_sym.append(rot) - return np.intc(reduced_site_sym) + return np.array(reduced_site_sym, dtype='intc') def get_bond_symmetry(site_symmetry, positions, diff --git a/anharmonic/phonon3/fc3.py b/anharmonic/phonon3/fc3.py index 2dfc058f..aa81bacb 100644 --- a/anharmonic/phonon3/fc3.py +++ b/anharmonic/phonon3/fc3.py @@ -341,7 +341,8 @@ def get_constrained_fc2(supercell, atom_list_done, lattice, positions, - np.intc(reduced_site_sym).copy(), + np.array(reduced_site_sym, + dtype='intc', order='C'), np.zeros((len(reduced_site_sym), 3), dtype='double'), symprec) diff --git a/phonopy/harmonic/force_constants.py b/phonopy/harmonic/force_constants.py index 08882b63..6e29256c 100644 --- a/phonopy/harmonic/force_constants.py +++ b/phonopy/harmonic/force_constants.py @@ -154,68 +154,41 @@ def distribute_force_constants(force_constants, positions, rotations, trans, - symprec, - is_symmetrize=False): + symprec): for atom_disp in atom_list: if atom_disp in atom_list_done: continue - if is_symmetrize: - map_atom_disps, map_syms = get_all_atom_mappings_by_symmetry( - atom_list_done, - atom_disp, - rotations, - trans, - positions, - symprec) - - for i, pos_i in enumerate(positions): - for map_atom_disp, map_sym in zip(map_atom_disps, map_syms): - # L R L^-1 - rot_cartesian = np.array(similarity_transformation( - lattice, rotations[map_sym]), dtype='double') - distribute_fc2_part(force_constants, - positions, - atom_disp, - map_atom_disp, - i, - rot_cartesian, - rotations[map_sym], - trans[map_sym], - symprec) - - force_constants[atom_disp, i] /= len(map_atom_disps) - else: - map_atom_disp, map_sym = get_atom_mapping_by_symmetry( - atom_list_done, - atom_disp, - rotations, - trans, - positions, - symprec) - - # L R L^-1 - rot_cartesian = np.array(similarity_transformation( - lattice, rotations[map_sym]), dtype='double') - - _distribute_fc2_part(force_constants, - positions, - atom_disp, - map_atom_disp, - rot_cartesian, - rotations[map_sym], - trans[map_sym], - symprec) + map_atom_disp, map_sym = get_atom_mapping_by_symmetry( + atom_list_done, + atom_disp, + rotations, + trans, + positions, + symprec) + + _distribute_fc2_part(force_constants, + positions, + atom_disp, + map_atom_disp, + lattice, + rotations[map_sym], + trans[map_sym], + symprec) def _distribute_fc2_part(force_constants, positions, atom_disp, map_atom_disp, - rot_cartesian, + lattice, r, t, symprec): + # L R L^-1 + rot_cartesian = np.array( + similarity_transformation(lattice, r), dtype='double', order='C') + try: import phonopy._phonopy as phonoc phonoc.distribute_fc2(force_constants, @@ -223,8 +196,8 @@ def _distribute_fc2_part(force_constants, atom_disp, map_atom_disp, rot_cartesian, - r, - t, + np.array(r, dtype='intc', order='C'), + np.array(t, dtype='double'), symprec) except ImportError: for i, pos_i in enumerate(positions): @@ -244,8 +217,7 @@ def _distribute_fc2_part(force_constants, # R^-1 P R (inverse transformation) force_constants[atom_disp, i] += similarity_transformation( rot_cartesian.T, - force_constants[map_atom_disp, - rot_atom]) + force_constants[map_atom_disp, rot_atom]) def get_atom_mapping_by_symmetry(atom_list_done, @@ -269,34 +241,6 @@ def get_atom_mapping_by_symmetry(atom_list_done, print 'or something wrong (e.g. crystal structure does not match).' raise ValueError -def get_all_atom_mappings_by_symmetry(atom_list_done, - atom_number, - rotations, - translations, - positions, - symprec=1e-5): - """ - Find mappings from an atom to atoms in the atom list done. - """ - - map_atoms = [] - map_syms = [] - for i, (r, t) in enumerate(zip(rotations, translations)): - rot_pos = np.dot(positions[atom_number], r.T) + t - for j in atom_list_done: - diff = positions[j] - rot_pos - if (abs(diff - np.rint(diff)) < symprec).all(): - map_atoms.append(j) - map_syms.append(i) - break - - if len(map_atoms) == 0: - print 'Input forces are not enough to calculate force constants,' - print 'or something wrong (e.g. crystal structure does not match).' - raise ValueError - - return map_atoms, map_syms - def _get_force_constants_disps(force_constants, supercell, dataset, @@ -359,9 +303,7 @@ def solve_force_constants(force_constants, force_constants[disp_atom_number, i] = -np.dot( inv_displacements, combined_forces) -def get_positions_sent_by_rot_inv(positions, - site_symmetry, - symprec): +def get_positions_sent_by_rot_inv(positions, site_symmetry, symprec): rot_map_syms = [] for sym in site_symmetry: rot_map = np.zeros(len(positions), dtype='intc') @@ -375,7 +317,7 @@ def get_positions_sent_by_rot_inv(positions, rot_map_syms.append(rot_map) - return np.array(rot_map_syms, dtype='intc') + return np.array(rot_map_syms, dtype='intc', order='C') def get_rotated_displacement(displacements, site_sym_cart): rot_disps = [] @@ -514,7 +456,7 @@ def rotational_invariance(force_constants, mat = np.zeros((3, 3), dtype='double') for s in range(supercell.get_number_of_atoms()): vecs = np.array(get_equivalent_smallest_vectors( - s, p, supercell, primitive.get_cell(), symprec)) + s, p, supercell, primitive.get_cell(), symprec)) m = len(vecs) v = np.dot(vecs[:,:].sum(axis=0) / m, primitive.get_cell()) for j in range(3):