1# Copyright (C) 2011 Atsushi Togo 2# All rights reserved. 3# 4# This file is part of phonopy. 5# 6# Redistribution and use in source and binary forms, with or without 7# modification, are permitted provided that the following conditions 8# are met: 9# 10# * Redistributions of source code must retain the above copyright 11# notice, this list of conditions and the following disclaimer. 12# 13# * Redistributions in binary form must reproduce the above copyright 14# notice, this list of conditions and the following disclaimer in 15# the documentation and/or other materials provided with the 16# distribution. 17# 18# * Neither the name of the phonopy project nor the names of its 19# contributors may be used to endorse or promote products derived 20# from this software without specific prior written permission. 21# 22# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 23# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 24# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS 25# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE 26# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, 27# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, 28# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 29# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER 30# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 31# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN 32# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 33# POSSIBILITY OF SUCH DAMAGE. 34 35import warnings 36import numpy as np 37from phonopy.units import Kb, THzToEv, EvTokJmol 38 39 40def mode_cv(temp, freqs): # freqs (eV) 41 x = freqs / Kb / temp 42 expVal = np.exp(x) 43 return Kb * x ** 2 * expVal / (expVal - 1.0) ** 2 44 45 46def mode_F(temp, freqs): 47 return (Kb * temp * np.log(1.0 - np.exp((- freqs) / (Kb * temp))) 48 + freqs / 2) 49 50 51def mode_S(temp, freqs): 52 val = freqs / (2 * Kb * temp) 53 return (1 / (2 * temp) * freqs * np.cosh(val) / np.sinh(val) 54 - Kb * np.log(2 * np.sinh(val))) 55 56 57def mode_ZPE(temp, freqs): 58 return freqs / 2 59 60 61def mode_zero(temp, freqs): 62 return np.zeros_like(freqs) 63 64 65class ThermalPropertiesBase(object): 66 def __init__(self, 67 mesh, 68 is_projection=False, 69 band_indices=None, 70 cutoff_frequency=None, 71 pretend_real=False): 72 self._is_projection = is_projection 73 self._band_indices = None 74 75 if cutoff_frequency is None or cutoff_frequency < 0: 76 self._cutoff_frequency = 0.0 77 else: 78 self._cutoff_frequency = cutoff_frequency 79 80 if band_indices is not None: 81 bi = np.hstack(band_indices).astype('intc') 82 self._band_indices = bi 83 self._frequencies = np.array(mesh.frequencies[:, bi], 84 dtype='double', order='C') 85 if mesh.eigenvectors is not None: 86 self._eigenvectors = np.array(mesh.eigenvectors[:, :, bi], 87 dtype='double', order='C') 88 else: 89 self._frequencies = mesh.frequencies 90 self._eigenvectors = mesh.eigenvectors 91 92 if pretend_real: 93 self._frequencies = abs(self._frequencies) 94 self._frequencies = np.array(self._frequencies, 95 dtype='double', order='C') * THzToEv 96 self._weights = mesh.weights 97 self._num_modes = self._frequencies.shape[1] * self._weights.sum() 98 self._num_integrated_modes = np.sum( 99 self._weights * (self._frequencies > 100 self._cutoff_frequency).sum(axis=1)) 101 102 def run_free_energy(self, t): 103 if t > 0: 104 free_energy = self._calculate_thermal_property(mode_F, t) 105 else: 106 free_energy = self._calculate_thermal_property(mode_ZPE, None) 107 return free_energy / np.sum(self._weights) * EvTokJmol 108 109 def run_heat_capacity(self, t): 110 if t > 0: 111 cv = self._calculate_thermal_property(mode_cv, t) 112 else: 113 cv = self._calculate_thermal_property(mode_zero, None) 114 return cv / np.sum(self._weights) * EvTokJmol 115 116 def run_entropy(self, t): 117 if t > 0: 118 entropy = self._calculate_thermal_property(mode_S, t) 119 else: 120 entropy = self._calculate_thermal_property(mode_zero, None) 121 return entropy / np.sum(self._weights) * EvTokJmol 122 123 def _calculate_thermal_property(self, func, t): 124 if not self._is_projection: 125 t_property = 0.0 126 for freqs, w in zip(self._frequencies, self._weights): 127 cond = freqs > self._cutoff_frequency 128 t_property += np.sum(func(t, freqs[cond])) * w 129 return t_property 130 else: 131 t_property = np.zeros(len(self._frequencies[0]), dtype='double') 132 for freqs, eigvecs2, w in zip(self._frequencies, 133 np.abs(self._eigenvectors) ** 2, 134 self._weights): 135 cond = freqs > self._cutoff_frequency 136 t_property += np.dot(eigvecs2[:, cond], 137 func(t, freqs[cond])) * w 138 return t_property 139 140 141class ThermalProperties(ThermalPropertiesBase): 142 def __init__(self, 143 mesh, 144 is_projection=False, 145 band_indices=None, 146 cutoff_frequency=None, 147 pretend_real=False): 148 ThermalPropertiesBase.__init__(self, 149 mesh, 150 is_projection=is_projection, 151 band_indices=band_indices, 152 cutoff_frequency=cutoff_frequency, 153 pretend_real=pretend_real) 154 self._thermal_properties = None 155 self._temperatures = None 156 self._high_T_entropy = None 157 self._zero_point_energy = None 158 self._projected_thermal_properties = None 159 160 self._set_high_T_entropy_and_zero_point_energy() 161 162 @property 163 def temperatures(self): 164 return self._temperatures 165 166 def get_temperatures(self): 167 warnings.warn("ThermalProperties.get_temperatures is deprecated." 168 "Use temperatures attribute.", 169 DeprecationWarning) 170 return self.temperatures 171 172 @temperatures.setter 173 def temperatures(self, temperatures): 174 t_array = np.array(temperatures, dtype='double') 175 self._temperatures = np.array( 176 np.extract(np.invert(t_array < 0), t_array), dtype='double') 177 178 def set_temperatures(self, temperatures): 179 warnings.warn("ThermalProperties.set_temperatures is deprecated." 180 "Use temperatures attribute.", 181 DeprecationWarning) 182 self.temperatures = temperatures 183 184 @property 185 def thermal_properties(self): 186 return self._thermal_properties 187 188 def get_thermal_properties(self): 189 warnings.warn("ThermalProperties.get_thermal_properties is deprecated." 190 "Use thermal_properties attribute.", 191 DeprecationWarning) 192 return self.thermal_properties 193 194 @property 195 def zero_point_energy(self): 196 return self._zero_point_energy 197 198 def get_zero_point_energy(self): 199 warnings.warn("ThermalProperties.get_zero_point_energy is deprecated." 200 "Use zero_point_energy attribute.", 201 DeprecationWarning) 202 return self.zero_point_energy 203 204 @property 205 def high_T_entropy(self): 206 return self._high_T_entropy 207 208 def get_high_T_entropy(self): 209 warnings.warn("ThermalProperties.get_high_T_entropy is deprecated." 210 "Use high_T_entropy attribute.", 211 DeprecationWarning) 212 return self.high_T_entropy 213 214 @property 215 def number_of_integrated_modes(self): 216 """Number of phonon modes used for integration on sampling mesh""" 217 return self._num_integrated_modes 218 219 def get_number_of_integrated_modes(self): 220 warnings.warn("ThermalProperties.get_number_of_integrated_modes is " 221 "deprecated. Use number_of_integrated_modes attribute.", 222 DeprecationWarning) 223 return self.number_of_integrated_modes 224 225 @property 226 def number_of_modes(self): 227 """Number of phonon modes on sampling mesh""" 228 return self._num_modes 229 230 def get_number_of_modes(self): 231 warnings.warn("ThermalProperties.get_number_of_modes is " 232 "deprecated. Use number_of_modes attribute.", 233 DeprecationWarning) 234 return self.number_of_modes 235 236 def set_temperature_range(self, t_min=None, t_max=None, t_step=None): 237 if t_min is None: 238 _t_min = 10 239 elif t_min < 0: 240 _t_min = 0 241 else: 242 _t_min = t_min 243 244 if t_max is None: 245 _t_max = 1000 246 elif t_max > _t_min: 247 _t_max = t_max 248 else: 249 _t_max = _t_min 250 251 if t_step is None: 252 _t_step = 10 253 elif t_step > 0: 254 _t_step = t_step 255 else: 256 _t_step = 10 257 258 self._temperatures = np.arange(_t_min, _t_max + _t_step / 2.0, _t_step, 259 dtype='double') 260 261 def plot(self, plt): 262 temps, fe, entropy, cv = self._thermal_properties 263 264 plt.plot(temps, fe, 'r-') 265 plt.plot(temps, entropy, 'b-') 266 plt.plot(temps, cv, 'g-') 267 plt.legend(('Free energy [kJ/mol]', 'Entropy [J/K/mol]', 268 r'C$_\mathrm{V}$ [J/K/mol]'), 269 loc='best') 270 plt.grid(True) 271 plt.xlabel('Temperature [K]') 272 273 def run(self, t_step=None, t_max=None, t_min=None, lang='C'): 274 import warnings 275 if (t_step is not None or t_max is not None or t_min is not None): 276 warnings.warn("keywords for this method are depreciated. " 277 "Use \'set_temperature_range\' or " 278 "\'set_temperature_range\' method instead.", 279 DeprecationWarning) 280 self.set_temperature_range(t_min=t_min, t_max=t_max, t_step=t_step) 281 282 if lang == 'C': 283 import phonopy._phonopy as phonoc 284 self._run_c_thermal_properties() 285 else: 286 self._run_py_thermal_properties() 287 288 if self._is_projection: 289 fe = [] 290 entropy = [] 291 cv = [] 292 for t in self._temperatures: 293 fe.append(self.run_free_energy(t)) 294 entropy.append(self.run_entropy(t) * 1000,) 295 cv.append(self.run_heat_capacity(t) * 1000) 296 297 self._projected_thermal_properties = [ 298 self._temperatures, 299 np.array(fe, dtype='double'), 300 np.array(entropy, dtype='double'), 301 np.array(cv, dtype='double')] 302 303 def write_yaml(self, filename='thermal_properties.yaml', volume=None): 304 lines = self._get_tp_yaml_lines(volume=volume) 305 if self._is_projection: 306 lines += self._get_projected_tp_yaml_lines() 307 with open(filename, 'w') as f: 308 f.write("\n".join(lines)) 309 310 def _run_c_thermal_properties(self): 311 import phonopy._phonopy as phonoc 312 313 props = np.zeros((len(self._temperatures), 3), 314 dtype='double', order='C') 315 phonoc.thermal_properties(props, 316 self._temperatures, 317 self._frequencies, 318 self._weights, 319 self._cutoff_frequency) 320 # for f, w in zip(self._frequencies, self._weights): 321 # phonoc.thermal_properties( 322 # props, 323 # self._temperatures, 324 # np.array(f, dtype='double', order='C')[None, :], 325 # np.array([w], dtype='intc'), 326 # cutoff_frequency) 327 328 props /= np.sum(self._weights) 329 fe = props[:, 0] * EvTokJmol + self._zero_point_energy 330 entropy = props[:, 1] * EvTokJmol * 1000 331 cv = props[:, 2] * EvTokJmol * 1000 332 self._thermal_properties = [self._temperatures, fe, entropy, cv] 333 334 def _run_py_thermal_properties(self): 335 fe = [] 336 entropy = [] 337 cv = [] 338 for t in self._temperatures: 339 props = self._get_py_thermal_properties(t) 340 fe.append(props[0]) 341 entropy.append(props[1] * 1000) 342 cv.append(props[2] * 1000) 343 self._thermal_properties = [ 344 self._temperatures, 345 np.array(fe, dtype='double'), 346 np.array(entropy, dtype='double'), 347 np.array(cv, dtype='double')] 348 349 def _get_tp_yaml_lines(self, volume=None): 350 lines = [] 351 lines.append("# Thermal properties / unit cell (natom)") 352 lines.append("") 353 lines.append("unit:") 354 lines.append(" temperature: K") 355 lines.append(" free_energy: kJ/mol") 356 lines.append(" entropy: J/K/mol") 357 lines.append(" heat_capacity: J/K/mol") 358 lines.append("") 359 lines.append("natom: %-5d" % (self._frequencies[0].shape[0] // 3)) 360 if volume is not None: 361 lines.append("volume: %-20.10f" % volume) 362 lines.append("cutoff_frequency: %8.3f" % self._cutoff_frequency) 363 lines.append("num_modes: %d" % self._num_modes) 364 lines.append("num_integrated_modes: %d" % self._num_integrated_modes) 365 if self._band_indices is not None: 366 bi = self._band_indices + 1 367 lines.append("band_index: [ " + ("%d, " * (len(bi) - 1)) % 368 tuple(bi[:-1]) + ("%d ]" % bi[-1])) 369 lines.append("") 370 lines.append("zero_point_energy: %15.7f" % self._zero_point_energy) 371 lines.append("high_T_entropy: %15.7f" % 372 (self._high_T_entropy * 1000)) 373 lines.append("") 374 lines.append("thermal_properties:") 375 temperatures, fe, entropy, cv = self._thermal_properties 376 for i, t in enumerate(temperatures): 377 lines.append("- temperature: %15.7f" % t) 378 lines.append(" free_energy: %15.7f" % fe[i]) 379 lines.append(" entropy: %15.7f" % entropy[i]) 380 # Sometimes 'nan' of C_V is returned at low temperature. 381 if np.isnan(cv[i]): 382 lines.append(" heat_capacity: %15.7f" % 0) 383 else: 384 lines.append(" heat_capacity: %15.7f" % cv[i]) 385 lines.append(" energy: %15.7f" % 386 (fe[i] + entropy[i] * t / 1000)) 387 lines.append("") 388 return lines 389 390 def _get_projected_tp_yaml_lines(self): 391 lines = [] 392 lines.append("projected_thermal_properties:") 393 temperatures, fe, entropy, cv = self._projected_thermal_properties 394 for i, t in enumerate(temperatures): 395 lines.append("- temperature: %13.7f" % t) 396 line = " free_energy: [ " 397 line += ", ".join(["%13.7f" % x for x in fe[i]]) 398 line += " ] # %13.7f" % np.sum(fe[i]) 399 lines.append(line) 400 line = " entropy: [ " 401 line += ", ".join(["%13.7f" % x for x in entropy[i]]) 402 line += " ] # %13.7f" % np.sum(entropy[i]) 403 lines.append(line) 404 # Sometimes 'nan' of C_V is returned at low temperature. 405 line = " heat_capacity: [ " 406 sum_cv = 0.0 407 for j, cv_i in enumerate(cv[i]): 408 if np.isnan(cv_i): 409 line += "%13.7f" % 0 410 else: 411 sum_cv += cv_i 412 line += "%13.7f" % cv_i 413 if j < len(cv[i]) - 1: 414 line += ", " 415 else: 416 line += " ]" 417 line += " # %13.7f" % sum_cv 418 lines.append(line) 419 energy = fe[i] + entropy[i] * t / 1000 420 line = " energy: [ " 421 line += ", ".join(["%13.7f" % x for x in energy]) 422 line += " ] # %13.7f" % np.sum(energy) 423 lines.append(line) 424 return lines 425 426 def _get_py_thermal_properties(self, t): 427 return (self.run_free_energy(t), 428 self.run_entropy(t), 429 self.run_heat_capacity(t)) 430 431 def _set_high_T_entropy_and_zero_point_energy(self): 432 zp_energy = 0.0 433 entropy = 0.0 434 for freqs, w in zip(self._frequencies, self._weights): 435 positive_fs = np.extract(freqs > 0.0, freqs) 436 entropy -= np.sum(np.log(positive_fs)) * w 437 zp_energy += np.sum(positive_fs) * w / 2 438 self._high_T_entropy = entropy * Kb / np.sum(self._weights) * EvTokJmol 439 self._zero_point_energy = zp_energy / np.sum(self._weights) * EvTokJmol 440