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

Minor changes of command options and fixes of ofrce fits

This commit is contained in:
Atsushi Togo 2014-06-16 17:39:15 +09:00
parent a396f294fe
commit 9ce31451ca
6 changed files with 164 additions and 139 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

29
anharmonic/phonon3/__init__.py
View File

@ -163,7 +163,7 @@ class Phono3py:
forces_fc2,
displacement_dataset=None,
is_permutation_symmetry=False,
is_translational_symmetry=False):
translational_symmetry_type=0):
if displacement_dataset is None:
disp_dataset = self._displacement_dataset
else:
@ -176,14 +176,16 @@ class Phono3py:
disp_dataset)
if is_permutation_symmetry:
set_permutation_symmetry(self._fc2)
if is_translational_symmetry:
set_translational_invariance(self._fc2)
if translational_symmetry_type:
set_translational_invariance(
self._fc2,
translational_symmetry_type=translational_symmetry_type)
def produce_fc3(self,
forces_fc3,
displacement_dataset=None,
cutoff_distance=None, # set fc3 zero
is_translational_symmetry=False,
translational_symmetry_type=0,
is_permutation_symmetry=False):
if displacement_dataset is None:
disp_dataset = self._displacement_dataset
@ -195,8 +197,10 @@ class Phono3py:
fc2 = get_fc2(self._supercell, self._symmetry, disp_dataset)
if is_permutation_symmetry:
set_permutation_symmetry(fc2)
if is_translational_symmetry:
set_translational_invariance(fc2)
if translational_symmetry_type:
set_translational_invariance(
fc2,
translational_symmetry_type=translational_symmetry_type)
count = len(disp_dataset['first_atoms'])
for disp1 in disp_dataset['first_atoms']:
@ -208,7 +212,7 @@ class Phono3py:
disp_dataset,
fc2,
self._symmetry,
is_translational_symmetry=is_translational_symmetry,
translational_symmetry_type=translational_symmetry_type,
is_permutation_symmetry=is_permutation_symmetry,
verbose=self._log_level)
@ -235,11 +239,16 @@ class Phono3py:
if self._fc3 is not None:
set_permutation_symmetry_fc3(self._fc3)
def set_translational_invariance(self):
def set_translational_invariance(self,
translational_symmetry_type=1):
if self._fc2 is not None:
set_translational_invariance(self._fc2)
set_translational_invariance(
self._fc2,
translational_symmetry_type=translational_symmetry_type)
if self._fc3 is not None:
set_translational_invariance_fc3(self._fc3)
set_translational_invariance_fc3(
self._fc3,
translational_symmetry_type=translational_symmetry_type)
def get_interaction_strength(self):
return self._interaction

114
anharmonic/phonon3/fc3.py
View File

@ -9,7 +9,7 @@ def get_fc3(supercell,
disp_dataset,
fc2,
symmetry,
is_translational_symmetry=False,
translational_symmetry_type=0,
is_permutation_symmetry=False,
verbose=False):
num_atom = supercell.get_number_of_atoms()
@ -30,7 +30,7 @@ def get_fc3(supercell,
disp_dataset,
fc2,
symmetry,
is_translational_symmetry,
translational_symmetry_type,
is_permutation_symmetry,
verbose)
@ -63,19 +63,21 @@ def get_fc3(supercell,
disp_dataset,
symmetry,
verbose=verbose)
if is_translational_symmetry:
if translational_symmetry_type:
if verbose:
print "Imposing translational invariance symmetry to fc3"
set_translational_invariance_fc3(fc3)
set_translational_invariance_fc3(
fc3, translational_symmetry_type=translational_symmetry_type)
if is_permutation_symmetry:
if verbose:
print "Imposing index permulation symmetry to fc3"
set_permutation_symmetry_fc3(fc3)
else:
if is_translational_symmetry:
if translational_symmetry_type:
if verbose:
print "Imposing translational invariance symmetry to fc3"
set_translational_invariance_fc3_per_index(fc3)
set_translational_invariance_fc3_per_index(
fc3, translational_symmetry_type=translational_symmetry_type)
if is_permutation_symmetry:
if verbose:
print "Imposing index permulation symmetry to fc3"
@ -185,51 +187,49 @@ def set_permutation_symmetry_fc3_elem(fc3, a, b, c, divisor=6):
fc3[c, b, a, k, j, i]) / divisor
return tensor3
def set_translational_invariance_fc3(fc3):
def set_translational_invariance_fc3(fc3,
translational_symmetry_type=1):
for i in range(3):
set_translational_invariance_fc3_per_index(fc3, index=i)
# set_translational_invariance_fc3_per_index_weighted(fc3, index=i)
set_translational_invariance_fc3_per_index(
fc3,
index=i,
translational_symmetry_type=translational_symmetry_type)
def set_translational_invariance_fc3_per_index(fc3, index=0):
def set_translational_invariance_fc3_per_index(fc3,
index=0,
translational_symmetry_type=1):
for i in range(fc3.shape[(1 + index) % 3]):
for j in range(fc3.shape[(2 + index) % 3]):
for k in range(fc3.shape[3]):
for l in range(fc3.shape[4]):
for m in range(fc3.shape[5]):
if index == 0:
fc3[:, i, j, k, l, m] -= np.sum(
fc3[:, i, j, k, l, m]) / fc3.shape[0]
elif index == 1:
fc3[j, :, i, k, l, m] -= np.sum(
fc3[j, :, i, k, l, m]) / fc3.shape[1]
elif index == 2:
fc3[i, j, :, k, l, m] -= np.sum(
fc3[i, j, :, k, l, m]) / fc3.shape[2]
def set_translational_invariance_fc3_per_index_weighted(fc3, index=0):
for i in range(fc3.shape[(1 + index) % 3]):
for j in range(fc3.shape[(2 + index) % 3]):
for k in range(fc3.shape[3]):
for l in range(fc3.shape[4]):
for m in range(fc3.shape[5]):
if index == 0:
fc_abs = np.abs(fc3[:, i, j, k, l, m])
fc_sum = np.sum(fc3[:, i, j, k, l, m])
fc_abs_sum = np.sum(fc_abs)
fc3[:, i, j, k, l, m] -= (
fc_sum / fc_abs_sum * fc_abs)
elif index == 1:
fc_abs = np.abs(fc3[i, :, j, k, l, m])
fc_sum = np.sum(fc3[i, :, j, k, l, m])
fc_abs_sum = np.sum(fc_abs)
fc3[i, :, j, k, l, m] -= (
fc_sum / fc_abs_sum * fc_abs)
elif index == 2:
fc_abs = np.abs(fc3[i, j, :, k, l, m])
fc_sum = np.sum(fc3[i, j, :, k, l, m])
fc_abs_sum = np.sum(fc_abs)
fc3[i, j, :, k, l, m] -= (
fc_sum / fc_abs_sum * fc_abs)
for k, l, m in list(np.ndindex(3, 3, 3)):
if translational_symmetry_type == 2: # Type 2
if index == 0:
fc_abs = np.abs(fc3[:, i, j, k, l, m])
fc_sum = np.sum(fc3[:, i, j, k, l, m])
fc_abs_sum = np.sum(fc_abs)
fc3[:, i, j, k, l, m] -= (
fc_sum / fc_abs_sum * fc_abs)
elif index == 1:
fc_abs = np.abs(fc3[i, :, j, k, l, m])
fc_sum = np.sum(fc3[i, :, j, k, l, m])
fc_abs_sum = np.sum(fc_abs)
fc3[i, :, j, k, l, m] -= (
fc_sum / fc_abs_sum * fc_abs)
elif index == 2:
fc_abs = np.abs(fc3[i, j, :, k, l, m])
fc_sum = np.sum(fc3[i, j, :, k, l, m])
fc_abs_sum = np.sum(fc_abs)
fc3[i, j, :, k, l, m] -= (
fc_sum / fc_abs_sum * fc_abs)
else: # Type 1
if index == 0:
fc3[:, i, j, k, l, m] -= np.sum(
fc3[:, i, j, k, l, m]) / fc3.shape[0]
elif index == 1:
fc3[j, :, i, k, l, m] -= np.sum(
fc3[j, :, i, k, l, m]) / fc3.shape[1]
elif index == 2:
fc3[i, j, :, k, l, m] -= np.sum(
fc3[i, j, :, k, l, m]) / fc3.shape[2]
def third_rank_tensor_rotation(rot_cart, tensor):
rot_tensor = np.zeros((3,3,3), dtype='double')
@ -278,14 +278,14 @@ def get_delta_fc2(dataset_second_atoms,
fc2,
supercell,
reduced_site_sym,
is_translational_symmetry,
translational_symmetry_type,
is_permutation_symmetry,
symprec):
disp_fc2 = get_constrained_fc2(supercell,
dataset_second_atoms,
atom1,
reduced_site_sym,
is_translational_symmetry,
translational_symmetry_type,
is_permutation_symmetry,
symprec)
return disp_fc2 - fc2
@ -294,7 +294,7 @@ def get_constrained_fc2(supercell,
dataset_second_atoms,
atom1,
reduced_site_sym,
is_translational_symmetry,
translational_symmetry_type,
is_permutation_symmetry,
symprec):
"""
@ -346,8 +346,10 @@ def get_constrained_fc2(supercell,
dtype='double'),
symprec)
if is_translational_symmetry:
set_translational_invariance(fc2)
if translational_symmetry_type:
set_translational_invariance(
fc2,
translational_symmetry_type=translational_symmetry_type)
if is_permutation_symmetry:
set_permutation_symmetry(fc2)
@ -467,7 +469,7 @@ def _get_fc3_least_atoms(fc3,
disp_dataset,
fc2,
symmetry,
is_translational_symmetry,
translational_symmetry_type,
is_permutation_symmetry,
verbose):
symprec = symmetry.get_symmetry_tolerance()
@ -480,7 +482,7 @@ def _get_fc3_least_atoms(fc3,
fc2,
first_atom_num,
symmetry.get_site_symmetry(first_atom_num),
is_translational_symmetry,
translational_symmetry_type,
is_permutation_symmetry,
symprec,
verbose)
@ -491,7 +493,7 @@ def _get_fc3_one_atom(fc3,
fc2,
first_atom_num,
site_symmetry,
is_translational_symmetry,
translational_symmetry_type,
is_permutation_symmetry,
symprec,
verbose):
@ -515,7 +517,7 @@ def _get_fc3_one_atom(fc3,
fc2,
supercell,
reduced_site_sym,
is_translational_symmetry,
translational_symmetry_type,
is_permutation_symmetry,
symprec))

45
phonopy/cui/settings.py
View File

@ -59,7 +59,7 @@ class Settings:
self._is_plusminus_displacement = 'auto'
self._is_symmetry = True
self._is_tetrahedron_method = False
self._is_time_symmetry = True
self._is_time_reversal_symmetry = True
self._is_translational_symmetry = False
self._is_trigonal_displacement = False
self._fc_decimals = None
@ -270,12 +270,16 @@ class Settings:
def get_temperature_step(self):
return self._tstep
def set_time_symmetry(self, time_symmetry=True):
self._is_time_symmetry = time_symmetry
def set_time_reversal_symmetry(self, time_reversal_symmetry=True):
self._is_time_reversal_symmetry = time_reversal_symmetry
def get_time_symmetry(self):
return self._is_time_symmetry
def get_time_reversal_symmetry(self):
return self._is_time_reversal_symmetry
# Translational symmetry type
# 0: No imposition
# 1: Simple sum, sum(fc) / N
# 2: Weighted sum, sum(fc) * abs(fc) / sum(abs(fc))
def set_tsym_type(self, tsym_type):
self._tsym_type = tsym_type
@ -1003,19 +1007,19 @@ class PhonopySettings(Settings):
def set_mesh(self,
mesh,
mesh_shift=[0.,0.,0.],
is_time_symmetry=True,
is_time_reversal_symmetry=True,
is_mesh_symmetry=True,
is_gamma_center=False):
self._mesh = mesh
self._mesh_shift = mesh_shift
self._is_time_symmetry = is_time_symmetry
self._is_time_reversal_symmetry = is_time_reversal_symmetry
self._is_mesh_symmetry = is_mesh_symmetry
self._is_gamma_center = is_gamma_center
def get_mesh(self):
return (self._mesh,
self._mesh_shift,
self._is_time_symmetry,
self._is_time_reversal_symmetry,
self._is_mesh_symmetry,
self._is_gamma_center)
@ -1232,9 +1236,9 @@ class PhonopyConfParser(ConfParser):
self.setting_error("MP_SHIFT is incorrectly set.")
self.set_parameter('mp_shift', vals[:3])
if conf_key == 'time_symmetry':
if confs['time_symmetry'] == '.false.':
self.set_parameter('is_time_symmetry', False)
if conf_key == 'time_reversal_symmetry':
if confs['time_reversal_symmetry'] == '.false.':
self.set_parameter('is_time_reversal_symmetry', False)
if conf_key == 'gamma_center':
if confs['gamma_center'] == '.true.':
@ -1419,10 +1423,10 @@ class PhonopyConfParser(ConfParser):
else:
shift = [0.,0.,0.]
time_symmetry = True
if params.has_key('is_time_symmetry'):
if not params['is_time_symmetry']:
time_symmetry = False
time_reversal_symmetry = True
if params.has_key('is_time_reversal_symmetry'):
if not params['is_time_reversal_symmetry']:
time_reversal_symmetry = False
mesh_symmetry = True
if params.has_key('is_mesh_symmetry'):
@ -1434,11 +1438,12 @@ class PhonopyConfParser(ConfParser):
if params['is_gamma_center']:
gamma_center = True
self._settings.set_mesh(params['mesh_numbers'],
mesh_shift=shift,
is_time_symmetry=time_symmetry,
is_mesh_symmetry=mesh_symmetry,
is_gamma_center=gamma_center)
self._settings.set_mesh(
params['mesh_numbers'],
mesh_shift=shift,
is_time_reversal_symmetry=time_reversal_symmetry,
is_mesh_symmetry=mesh_symmetry,
is_gamma_center=gamma_center)
# band mode
if params.has_key('band_paths'):

68
phonopy/harmonic/force_constants.py
View File

@ -438,50 +438,44 @@ def set_tensor_symmetry(force_constants,
# Take average and set to new force cosntants
force_constants[i, j] = tmp_fc / len(rotations)
def set_translational_invariance(force_constants):
def set_translational_invariance(force_constants,
translational_symmetry_type=1):
"""
Translational invariance is imposed. This is quite simple
Translational invariance is imposed. The type1 is quite simple
implementation, which is just taking sum of the force constants in
an axis and an atom index. The sum has to be zero due to the
translational invariance. If the sum is not zero, this error is
uniformly subtracted from force constants.
"""
set_translational_invariance_per_index(force_constants, index=0)
set_translational_invariance_per_index(force_constants, index=1)
# set_translational_invariance_per_index_weighted(force_constants, index=0)
# set_translational_invariance_per_index_weighted(force_constants, index=1)
for i in range(2):
set_translational_invariance_per_index(
force_constants,
index=i,
translational_symmetry_type=translational_symmetry_type)
def set_translational_invariance_per_index(force_constants, index=0):
if index == 0:
for i in range(force_constants.shape[1]):
for j in range(force_constants.shape[2]):
for k in range(force_constants.shape[3]):
force_constants[:, i, j, k] -= np.sum(
force_constants[:, i, j, k]) / force_constants.shape[0]
elif index == 1:
for i in range(force_constants.shape[0]):
for j in range(force_constants.shape[2]):
for k in range(force_constants.shape[3]):
force_constants[i, :, j, k] -= np.sum(
force_constants[i, :, j, k]) / force_constants.shape[1]
def set_translational_invariance_per_index_weighted(force_constants, index=0):
if index == 0:
for i in range(force_constants.shape[1]):
for j in range(force_constants.shape[2]):
for k in range(force_constants.shape[3]):
fc_abs = np.abs(force_constants[:, i, j, k])
fc_sum = np.sum(force_constants[:, i, j, k])
fc_abs_sum = np.sum(fc_abs)
force_constants[:, i, j, k] -= fc_sum / fc_abs_sum * fc_abs
elif index == 1:
for i in range(force_constants.shape[0]):
for j in range(force_constants.shape[2]):
for k in range(force_constants.shape[3]):
fc_abs = np.abs(force_constants[i, :, j, k])
fc_sum = np.sum(force_constants[i, :, j, k])
fc_abs_sum = np.sum(fc_abs)
force_constants[i, :, j, k] -= fc_sum / fc_abs_sum * fc_abs
def set_translational_invariance_per_index(fc2,
index=0,
translational_symmetry_type=1):
for i in range(fc2.shape[1 - index]):
for j, k in list(np.ndindex(3, 3)):
if translational_symmetry_type == 2: # Type 2
if index == 0:
fc_abs = np.abs(fc2[:, i, j, k])
fc_sum = np.sum(fc2[:, i, j, k])
fc_abs_sum = np.sum(fc_abs)
fc2[:, i, j, k] -= fc_sum / fc_abs_sum * fc_abs
else:
fc_abs = np.abs(fc2[i, :, j, k])
fc_sum = np.sum(fc2[i, :, j, k])
fc_abs_sum = np.sum(fc_abs)
fc2[i, :, j, k] -= fc_sum / fc_abs_sum * fc_abs
else: # Type 1
if index == 0:
fc2[:, i, j, k] -= np.sum(
fc2[:, i, j, k]) / fc2.shape[0]
else:
fc2[i, :, j, k] -= np.sum(
fc2[i, :, j, k]) / fc2.shape[1]
def set_permutation_symmetry(force_constants):
"""

33
scripts/force-fit
View File

@ -46,7 +46,7 @@ from phonopy.file_IO import parse_FORCE_SETS
from anharmonic.force_fit.fc2 import FC2Fit
from anharmonic.force_fit.fc3 import FC3Fit
from anharmonic.force_fit.fc4 import FC4Fit
from anharmonic.file_IO import parse_disp_fc4_yaml, parse_FORCES_FOURTH, parse_disp_fc3_yaml, parse_FORCES_THIRD, parse_FORCES_SECOND
from anharmonic.file_IO import parse_disp_fc4_yaml, parse_FORCES_FOURTH, parse_disp_fc3_yaml, parse_FORCES_FC3, parse_disp_fc2_yaml, parse_FORCES_FC2
from anharmonic.file_IO import write_fc3_dat, write_fc4_dat, write_fc4_to_hdf5, write_fc3_to_hdf5, write_fc2_to_hdf5
from anharmonic.phonon3.fc3 import show_drift_fc3
from anharmonic.phonon4.fc4 import show_drift_fc4
@ -136,9 +136,11 @@ if options.fc2:
disp_dataset = parse_FORCE_SETS("FORCE_SETS")
else:
file_exists("disp_fc3.yaml")
file_exists("FORCES_SECOND")
disp_dataset = parse_disp_fc3_yaml()
parse_FORCES_THIRD(disp_dataset)
file_exists("FORCES_FC2")
disp_dataset = parse_disp_fc2_yaml()
forces_fc2 = parse_FORCES_FC2(disp_dataset)
for forces, disp1 in zip(forces_fc2, disp_dataset['first_atoms']):
disp1['forces'] = forces
if options.coef_invariants is not None:
print "Adjustment parameter: %e" % options.coef_invariants
if options.pinv_cutoff is not None:
@ -153,21 +155,28 @@ if options.fc2:
fc2fit.run()
fc2 = fc2fit.get_fc2()
print "Writing fc2..."
write_fc2_to_hdf5(fc2, 'fc2-fit.hdf5')
write_fc2_to_hdf5(fc2, 'fc2.fit.hdf5')
if options.fc3:
file_exists("disp_fc3.yaml")
file_exists("FORCES_SECOND")
file_exists("FORCES_THIRD")
file_exists("FORCES_FC3")
disp_dataset = parse_disp_fc3_yaml()
parse_FORCES_THIRD(disp_dataset)
forces_fc3 = parse_FORCES_FC3(disp_dataset)
for forces, disp1 in zip(forces_fc3, disp_dataset['first_atoms']):
disp1['forces'] = forces
count = len(disp_dataset['first_atoms'])
for disp1 in disp_dataset['first_atoms']:
for disp2 in disp1['second_atoms']:
disp2['forces'] = forces_fc3[count]
count += 1
fc3fit = FC3Fit(supercell, disp_dataset, symmetry)
fc3fit.run()
fc3 = fc3fit.get_fc3()
print "Calculating drift fc3..."
show_drift_fc3(fc3)
print "Writing fc3..."
write_fc3_to_hdf5(fc3, 'fc3-fit.hdf5')
write_fc3_to_hdf5(fc3, 'fc3.fit.hdf5')
if options.fc4:
file_exists("disp_fc4.yaml")
@ -182,15 +191,15 @@ if options.fc4:
print "Calculating drift fc4..."
show_drift_fc4(fc4)
print "Writing fc4..."
write_fc4_to_hdf5(fc4, 'fc4-fit.hdf5')
write_fc4_to_hdf5(fc4, 'fc4.fit.hdf5')
fc3 = fc4fit.get_fc3()
print "Calculating drift of fc3..."
show_drift_fc3(fc3)
print "Writing fc3..."
write_fc3_to_hdf5(fc3, 'fc3-fit.hdf5')
write_fc3_to_hdf5(fc3, 'fc3.fit.hdf5')
fc2 = fc4fit.get_fc2()
show_drift_force_constants(fc2, name="fc2")
print "Writing fc2..."
write_fc2_to_hdf5(fc2, 'fc2-fit.hdf5')
write_fc2_to_hdf5(fc2, 'fc2.fit.hdf5')

14
scripts/phono3py
View File

@ -553,6 +553,12 @@ if settings.get_cutoff_frequency() is None:
cutoff_frequency = 1e-2
else:
cutoff_frequency = settings.get_cutoff_frequency()
if settings.get_is_translational_symmetry():
tsym_type = 1
elif settings.get_tsym_type() > 0:
tsym_type = settings.get_tsym_type()
else:
tsym_type = 0
phono3py = Phono3py(
unitcell,
@ -601,7 +607,7 @@ if log_level:
print "----------------------------- General settings ------------------------------"
print "Imposing translational symmetry to fc3 and fc2:",
print settings.get_is_translational_symmetry()
print (tsym_type > 0)
print "Imposing symmetry of index exchange to fc2:",
print options.is_symmetrize_fc2
print "Imposing symmetry of index exchange to fc3 in real space:",
@ -687,7 +693,7 @@ else:
forces_fc3,
displacement_dataset=disp_dataset,
cutoff_distance=settings.get_cutoff_fc3_distance(),
is_translational_symmetry=settings.get_is_translational_symmetry(),
translational_symmetry_type=tsym_type,
is_permutation_symmetry=options.is_symmetrize_fc3_r)
if output_filename is None:
filename = 'fc3.hdf5'
@ -741,7 +747,7 @@ else:
phono3py.produce_fc2(
forces_fc3,
displacement_dataset=disp_dataset,
is_translational_symmetry=settings.is_translational_symmetry(),
translational_symmetry_type=tsym_type,
is_permutation_symmetry=options.is_symmetrize_fc2)
else:
if input_filename is None:
@ -759,7 +765,7 @@ else:
phono3py.produce_fc2(
forces_fc2,
displacement_dataset=disp_dataset,
is_translational_symmetry=settings.get_is_translational_symmetry(),
translational_symmetry_type=tsym_type,
is_permutation_symmetry=options.is_symmetrize_fc2)
if output_filename is None: