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

JDOS with phonon occupation numbers

This commit is contained in:
Atsushi Togo 2014-08-19 16:22:47 +09:00
parent 91cbfb6908
commit 9b31b58f1c
4 changed files with 83 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
anharmonic/file_IO.py
View File

@ -757,12 +757,42 @@ def write_joint_dos(gp,
frequencies,
jdos,
sigma=None,
temperatures=None,
filename=None,
is_nosym=False):
if temperatures is None:
_write_joint_dos_at_t(gp,
mesh,
frequencies,
jdos,
sigma=sigma,
temperature=None,
filename=filename,
is_nosym=is_nosym)
else:
for jdos_at_t, t in zip(jdos, temperatures):
_write_joint_dos_at_t(gp,
mesh,
frequencies,
jdos_at_t,
sigma=sigma,
temperature=t,
filename=filename,
is_nosym=is_nosym)
def _write_joint_dos_at_t(gp,
mesh,
frequencies,
jdos,
sigma=None,
temperature=None,
filename=None,
is_nosym=False):
jdos_filename = "jdos-m%d%d%d-g%d" % (mesh[0], mesh[1], mesh[2], gp)
if sigma is not None:
jdos_filename += ("-s%f" % sigma).rstrip('0').rstrip('\.')
if temperature is not None:
jdos_filename += ("-t%f" % temperature).rstrip('0').rstrip('\.')
if is_nosym:
jdos_filename += ".nosym"
if filename is not None:

7
anharmonic/phonon3/__init__.py
View File

@ -672,8 +672,10 @@ class Phono3pyJointDos:
nac_params=None,
nac_q_direction=None,
sigmas=[],
cutoff_frequency=1e-4,
frequency_step=None,
num_frequency_points=None,
temperatures=None,
frequency_factor_to_THz=VaspToTHz,
frequency_scale_factor=None,
is_nosym=False,
@ -687,8 +689,10 @@ class Phono3pyJointDos:
self._nac_params = nac_params
self._nac_q_direction = nac_q_direction
self._sigmas = sigmas
self._cutoff_frequency = cutoff_frequency
self._frequency_step = frequency_step
self._num_frequency_points = num_frequency_points
self._temperatures = temperatures
self._frequency_factor_to_THz = frequency_factor_to_THz
self._frequency_scale_factor = frequency_scale_factor
self._is_nosym = is_nosym
@ -703,8 +707,10 @@ class Phono3pyJointDos:
self._fc2,
nac_params=self._nac_params,
nac_q_direction=self._nac_q_direction,
cutoff_frequency=self._cutoff_frequency,
frequency_step=self._frequency_step,
num_frequency_points=self._num_frequency_points,
temperatures=self._temperatures,
frequency_factor_to_THz=self._frequency_factor_to_THz,
frequency_scale_factor=self._frequency_scale_factor,
is_nosym=self._is_nosym,
@ -747,6 +753,7 @@ class Phono3pyJointDos:
self._jdos.get_frequency_points(),
self._jdos.get_joint_dos(),
sigma=sigma,
temperatures=self._temperatures,
filename=self._filename,
is_nosym=self._is_nosym)

43
anharmonic/phonon3/joint_dos.py
View File

@ -15,8 +15,10 @@ class JointDos:
nac_params=None,
nac_q_direction=None,
sigma=None,
cutoff_frequency=None,
frequency_step=None,
num_frequency_points=None,
temperatures=None,
frequency_factor_to_THz=VaspToTHz,
frequency_scale_factor=1.0,
is_nosym=False,
@ -36,8 +38,13 @@ class JointDos:
self._sigma = None
self.set_sigma(sigma)
if cutoff_frequency is None:
self._cutoff_frequency = 0
else:
self._cutoff_frequency = cutoff_frequency
self._frequency_step = frequency_step
self._num_frequency_points = num_frequency_points
self._temperatures = temperatures
self._frequency_factor_to_THz = frequency_factor_to_THz
self._frequency_scale_factor = frequency_scale_factor
self._is_nosym = is_nosym
@ -142,7 +149,17 @@ class JointDos:
num_freq_points = len(self._frequency_points)
num_mesh = np.prod(self._mesh)
jdos = np.zeros((num_freq_points, 2), dtype='double')
if self._temperatures is None:
jdos = np.zeros((num_freq_points, 2), dtype='double')
else:
num_temps = len(self._temperatures)
jdos = np.zeros((num_temps, num_freq_points, 2), dtype='double')
occ_phonons = []
for t in self._temperatures:
freqs = self._frequencies[self._triplets_at_q[:, 1:]]
occ_phonons.append(np.where(freqs > self._cutoff_frequency,
occupation(freqs, t), 0))
for i, freq_point in enumerate(self._frequency_points):
g = get_triplets_integration_weights(
@ -151,14 +168,22 @@ class JointDos:
self._sigma,
is_collision_matrix=True,
neighboring_phonons=(i == 0))
jdos[i, 0] = np.sum(
np.tensordot(g[0, :, 0], self._weights_at_q, axes=(0, 0)))
gx = g[2] - g[0]
jdos[i, 1] = np.sum(
np.tensordot(gx[:, 0], self._weights_at_q, axes=(0, 0)))
if self._log_level > 1:
print "%4d %f %e %e" % ((i + 1, freq_point,) +
tuple(jdos[i] / num_mesh))
if self._temperatures is None:
jdos[i, 0] = np.sum(
np.tensordot(g[0, :, 0], self._weights_at_q, axes=(0, 0)))
gx = g[2] - g[0]
jdos[i, 1] = np.sum(
np.tensordot(gx[:, 0], self._weights_at_q, axes=(0, 0)))
else:
for j, n in enumerate(occ_phonons):
for k, l in list(np.ndindex(g.shape[3:])):
jdos[j, i, 0] += np.dot(
(n[:, 0, k] + n[:, 1, l] + 1) *
g[0, :, 0, k, l], self._weights_at_q)
jdos[j, i, 1] += np.dot((n[:, 0, k] - n[:, 1, l]) *
g[1, :, 0, k, l],
self._weights_at_q)
self._joint_dos = jdos / num_mesh

19
scripts/phono3py
View File

@ -545,11 +545,14 @@ else:
if settings.get_is_tetrahedron_method():
sigmas = [None] + sigmas
if settings.get_temperatures() is None:
t_max=settings.get_max_temperature()
t_min=settings.get_min_temperature()
t_step=settings.get_temperature_step()
temperature_points = [0.0, 300.0] # For spectra
temperatures = np.arange(t_min, t_max + float(t_step) / 10, t_step)
if options.is_joint_dos:
temperature_points = None
else:
t_max=settings.get_max_temperature()
t_min=settings.get_min_temperature()
t_step=settings.get_temperature_step()
temperature_points = [0.0, 300.0] # For spectra
temperatures = np.arange(t_min, t_max + float(t_step) / 10, t_step)
else:
temperature_points = settings.get_temperatures() # For spectra
temperatures = settings.get_temperatures() # For others
@ -921,9 +924,9 @@ if log_level:
print " %.1f" % temperatures[-1]
else:
print ("%.1f " * len(temperatures)) % tuple(temperatures)
elif settings.get_is_isotope() or options.is_joint_dos:
elif settings.get_is_isotope():
pass
else:
elif temperature_points is not None:
print "Temperatures:",
print ("%.1f " * len(temperature_points)) % tuple(temperature_points)
if settings.get_scattering_event_class() is not None:
@ -964,8 +967,10 @@ if options.is_joint_dos:
nac_params=nac_params,
nac_q_direction=nac_q_direction,
sigmas=sigmas,
cutoff_frequency=cutoff_frequency,
frequency_step=frequency_step,
num_frequency_points=num_frequency_points,
temperatures=temperature_points,
frequency_factor_to_THz=frequency_factor_to_THz,
frequency_scale_factor=frequency_scale_factor,
is_nosym=options.is_nosym,