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

A little code cleanup

This commit is contained in:
Atsushi Togo 2014-06-17 18:32:56 +09:00
parent 9ce31451ca
commit db6d0b0c27
4 changed files with 50 additions and 102 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

33
anharmonic/force_fit/fc3.py
View File

@ -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

2
anharmonic/phonon3/displacement_fc3.py
View File

@ -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,

3
anharmonic/phonon3/fc3.py
View File

@ -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)

112
phonopy/harmonic/force_constants.py
View File

@ -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)
map_atom_disp, map_sym = get_atom_mapping_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)
_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):