1# Copyright (C) 2020 Atsushi Togo 2# All rights reserved. 3# 4# This file is part of phono3py. 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 35"""SSCHA calculation 36 37Formulae implemented are based on these papers: 38 39 Ref. 1: Bianco et al. https://doi.org/10.1103/PhysRevB.96.014111 40 Ref. 2: Aseginolaza et al. https://doi.org/10.1103/PhysRevB.100.214307 41 42""" 43 44import numpy as np 45from phonopy.units import VaspToTHz 46from phonopy.harmonic.dynmat_to_fc import DynmatToForceConstants 47from phono3py.phonon.func import mode_length 48 49 50class SupercellPhonon(object): 51 """Supercell phonon class 52 53 Dynamical matrix is created for supercell atoms and solved in real. 54 All phonons at commensurate points are folded to those at Gamma point. 55 Phonon eigenvectors are represented in real type. 56 57 Attributes 58 ---------- 59 eigenvalues : ndarray 60 Phonon eigenvalues of supercell dynamical matrix. 61 shape=(3 * num_satom, ), dtype='double', order='C' 62 eigenvectors : ndarray 63 Phonon eigenvectors of supercell dynamical matrix. 64 shape=(3 * num_satom, 3 * num_satom), dtype='double', order='C' 65 frequencies : ndarray 66 Phonon frequencies of supercell dynamical matrix. Frequency conversion 67 factor to THz is multiplied. 68 shape=(3 * num_satom, ), dtype='double', order='C' 69 force_constants : ndarray 70 Supercell force constants. The array shape is different from 71 the input force constants. 72 shape=(3 * num_satom, 3 * num_satom), dtype='double', order='C' 73 supercell : PhonopyAtoms or its derived class 74 Supercell. 75 76 """ 77 78 def __init__(self, 79 supercell, 80 force_constants, 81 frequency_factor_to_THz=VaspToTHz): 82 """ 83 84 Parameters 85 ---------- 86 supercell : PhonopyAtoms 87 Supercell. 88 force_constants : array_like 89 Second order force constants. 90 shape=(num_satom, num_satom, 3, 3), dtype='double', order='C' 91 frequency_factor_to_THz : float 92 Frequency conversion factor to THz. 93 94 """ 95 96 self._supercell = supercell 97 N = len(supercell) 98 _fc2 = np.array(np.transpose(force_constants, axes=[0, 2, 1, 3]), 99 dtype='double', order='C') 100 _fc2 = _fc2.reshape((3 * N, 3 * N)) 101 _fc2 = np.array(_fc2, dtype='double', order='C') 102 inv_sqrt_masses = 1.0 / np.repeat(np.sqrt(supercell.masses), 3) 103 dynmat = np.array(inv_sqrt_masses * (inv_sqrt_masses * _fc2).T, 104 dtype='double', order='C') 105 eigvals, eigvecs = np.linalg.eigh(dynmat) 106 freqs = np.sqrt(np.abs(eigvals)) * np.sign(eigvals) 107 freqs *= frequency_factor_to_THz 108 self._eigenvalues = np.array(eigvals, dtype='double', order='C') 109 self._eigenvectors = np.array(eigvecs, dtype='double', order='C') 110 self._frequencies = np.array(freqs, dtype='double', order='C') 111 self._force_constants = _fc2 112 113 @property 114 def eigenvalues(self): 115 return self._eigenvalues 116 117 @property 118 def eigenvectors(self): 119 return self._eigenvectors 120 121 @property 122 def frequencies(self): 123 return self._frequencies 124 125 @property 126 def force_constants(self): 127 return self._force_constants 128 129 @property 130 def supercell(self): 131 return self._supercell 132 133 134class DispCorrMatrix(object): 135 """Calculate displacement correlation matrix <u_a u_b> 136 137 Attributes 138 ---------- 139 psi_matrix : ndarray 140 Displacement-displacement correlation matrix at temperature. 141 Physical unit is [Angstrom^2]. 142 shape=(3 * num_satom, 3 * num_satom), dtype='double', order='C' 143 upsilon_matrix : ndarray 144 Inverse displacement-displacement correlation matrix at temperature. 145 Physical unit is [1/Angstrom^2]. 146 shape=(3 * num_satom, 3 * num_satom), dtype='double', order='C' 147 supercell_phonon : SupercellPhonon 148 Supercell phonon object. Phonons at Gamma point, where 149 eigenvectors are not complex type but real type. 150 151 """ 152 153 def __init__(self, 154 supercell_phonon, 155 cutoff_frequency=1e-5): 156 """ 157 158 Parameters 159 ---------- 160 supercell_phonon : SupercellPhonon 161 Supercell phonon object. Phonons at Gamma point, where 162 eigenvectors are not complex type but real type. 163 cutoff_frequency : float 164 Phonons are ignored if they have frequencies less than this value. 165 166 """ 167 168 self._supercell_phonon = supercell_phonon 169 self._cutoff_frequency = cutoff_frequency 170 self._psi_matrix = None 171 self._upsilon_matrix = None 172 self._determinant = None 173 174 def run(self, T): 175 freqs = self._supercell_phonon.frequencies 176 eigvecs = self._supercell_phonon.eigenvectors 177 sqrt_masses = np.repeat( 178 np.sqrt(self._supercell_phonon.supercell.masses), 3) 179 inv_sqrt_masses = np.repeat( 180 1.0 / np.sqrt(self._supercell_phonon.supercell.masses), 3) 181 182 # ignore zero and imaginary frequency modes 183 condition = freqs > self._cutoff_frequency 184 _freqs = np.where(condition, freqs, 1) 185 _a = mode_length(_freqs, T) 186 a2 = np.where(condition, _a ** 2, 0) 187 a2_inv = np.where(condition, 1 / _a ** 2, 0) 188 189 matrix = np.dot(a2 * eigvecs, eigvecs.T) 190 self._psi_matrix = np.array( 191 inv_sqrt_masses * (inv_sqrt_masses * matrix).T, 192 dtype='double', order='C') 193 194 matrix = np.dot(a2_inv * eigvecs, eigvecs.T) 195 self._upsilon_matrix = np.array( 196 sqrt_masses * (sqrt_masses * matrix).T, 197 dtype='double', order='C') 198 199 self._determinant = np.prod(2 * np.pi * np.extract(condition, a2)) 200 201 @property 202 def upsilon_matrix(self): 203 return self._upsilon_matrix 204 205 @property 206 def psi_matrix(self): 207 return self._psi_matrix 208 209 @property 210 def supercell_phonon(self): 211 return self._supercell_phonon 212 213 @property 214 def determinant(self): 215 return self._determinant 216 217 218class DispCorrMatrixMesh(object): 219 """Calculate upsilon and psi matrix 220 221 This calculation is similar to the transformation from 222 dynamical matrices to force constants. Instead of creating 223 dynamcial matrices from eigenvalues and eigenvectors, 224 1/a**2 or a**2 and eigenvectors are used, where a is mode length. 225 226 psi_matrix : ndarray 227 Displacement-displacement correlation matrix at temperature. 228 Physical unit is [Angstrom^2]. 229 shape=(3 * num_satom, 3 * num_satom), dtype='double', order='C' 230 upsilon_matrix : ndarray 231 Inverse displacement-displacement correlation matrix at temperature. 232 Physical unit is [1/Angstrom^2]. 233 shape=(3 * num_satom, 3 * num_satom), dtype='double', order='C' 234 commensurate_points : ndarray 235 Commensurate q-points of transformation matrix from primitive cell 236 to supercell. 237 shape=(det(transformation_matrix), 3), dtype='double', order='C'. 238 239 """ 240 241 def __init__(self, 242 primitive, 243 supercell, 244 cutoff_frequency=1e-5): 245 self._d2f = DynmatToForceConstants( 246 primitive, supercell, is_full_fc=True) 247 self._masses = supercell.masses 248 self._cutoff_frequency = cutoff_frequency 249 250 self._psi_matrix = None 251 self._upsilon_matrix = None 252 253 @property 254 def commensurate_points(self): 255 return self._d2f.commensurate_points 256 257 def run(self, frequencies, eigenvectors, T): 258 """ 259 260 Parameters 261 ---------- 262 frequencies : ndarray 263 Supercell phonon frequencies in THz (without 2pi). 264 shape=(grid_point, band), dtype='double', order='C' 265 eigenvectors : ndarray 266 Supercell phonon eigenvectors. 267 shape=(grid_point, band, band), dtype='double', order='C' 268 269 """ 270 271 condition = frequencies > self._cutoff_frequency 272 _freqs = np.where(condition, frequencies, 1) 273 _a = mode_length(_freqs, T) 274 a2 = np.where(condition, _a ** 2, 0) 275 a2_inv = np.where(condition, 1 / _a ** 2, 0) 276 N = len(self._masses) 277 shape = (N * 3, N * 3) 278 279 self._d2f.create_dynamical_matrices(a2_inv, eigenvectors) 280 self._d2f.run() 281 matrix = self._d2f.force_constants 282 matrix = np.transpose(matrix, axes=[0, 2, 1, 3]).reshape(shape) 283 self._upsilon_matrix = np.array(matrix, dtype='double', order='C') 284 285 self._d2f.create_dynamical_matrices(a2, eigenvectors) 286 self._d2f.run() 287 matrix = self._d2f.force_constants 288 for i, m_i in enumerate(self._masses): 289 for j, m_j in enumerate(self._masses): 290 matrix[i, j] /= m_i * m_j 291 matrix = np.transpose(matrix, axes=[0, 2, 1, 3]).reshape(shape) 292 self._psi_matrix = np.array(matrix, dtype='double', order='C') 293 294 @property 295 def upsilon_matrix(self): 296 return self._upsilon_matrix 297 298 @property 299 def psi_matrix(self): 300 return self._psi_matrix 301 302 303class SecondOrderFC(object): 304 r"""SSCHA second order force constants by ensemble average 305 306 This class is made just for the test of the ensemble average in 307 Ref. 1, and will not be used for usual fc2 calculation. 308 309 \Phi_{ab} = - \sum_{p} \Upsilon_{ap} \left< u^p f_b 310 \right>_{tilde{\rho}_{\mathfrak{R},\Phi}} 311 312 Attributes 313 ---------- 314 displacements : ndarray 315 shape=(snap_shots, 3 * num_satoms), dtype='double', order='C' 316 forces : ndarray 317 shape=(snap_shots, 3 * num_satoms), dtype='double', order='C' 318 fc2 : ndarray 319 shape=(num_satom, num_satom, 3, 3) 320 321 """ 322 323 def __init__(self, 324 displacements, 325 forces, 326 supercell_phonon, 327 cutoff_frequency=1e-5, 328 log_level=0): 329 """ 330 331 Parameters 332 ---------- 333 displacements : ndarray 334 shape=(snap_shots, num_satom, 3), dtype='double', order='C' 335 forces : ndarray 336 shape=(snap_shots, num_satom, 3), dtype='double', order='C' 337 supercell_phonon : SupercellPhonon 338 Supercell phonon object. Phonons at Gamma point, where 339 eigenvectors are not complex type but real type. 340 cutoff_frequency : float 341 Phonons are ignored if they have frequencies less than this value. 342 343 """ 344 345 assert (displacements.shape == forces.shape) 346 shape = displacements.shape 347 u = np.array(displacements.reshape(shape[0], -1), 348 dtype='double', order='C') 349 f = np.array(forces.reshape(shape[0], -1), 350 dtype='double', order='C') 351 self._displacements = u 352 self._forces = f 353 self._uu = DispCorrMatrix( 354 supercell_phonon, cutoff_frequency=cutoff_frequency) 355 self._force_constants = supercell_phonon.force_constants 356 self._cutoff_frequency = cutoff_frequency 357 self._log_level = log_level 358 359 @property 360 def displacements(self): 361 return self._displacements 362 363 @property 364 def forces(self): 365 return self._forces 366 367 @property 368 def fc2(self): 369 return self._fc2 370 371 def run(self, T=300.0): 372 """Calculate fc2 stochastically 373 374 As displacement correlation matrix <u_a u_b>^-1, two choices exist. 375 376 1. Upsilon matrix 377 self._uu.run(T) 378 Y = self._uu.upsilon_matrix 379 2. From displacements used for the force calculations 380 Y = np.linalg.pinv(np.dot(u.T, u) / u.shape[0]) 381 382 They should give close results, otherwise self-consistency 383 is not reached well. The later seems better agreement with 384 least square fit to fc2 under not very good self-consistency 385 condition. 386 387 """ 388 389 u = self._displacements 390 f = self._forces 391 Y = np.linalg.pinv(np.dot(u.T, u) / u.shape[0]) 392 u_inv = np.dot(u, Y) 393 394 if self._log_level: 395 nelems = np.prod(u.shape) 396 # print("sum u_inv:", u_inv.sum(axis=0) / u.shape[0]) 397 print("sum all u_inv:", u_inv.sum() / nelems) 398 print("rms u_inv:", np.sqrt((u_inv ** 2).sum() / nelems)) 399 print("rms u:", np.sqrt((u ** 2).sum() / nelems)) 400 print("rms forces:", np.sqrt((self._forces ** 2).sum() / nelems)) 401 # print("drift forces:", 402 # self._forces.sum(axis=0) / self._forces.shape[0]) 403 404 fc2 = - np.dot(u_inv.T, f) / f.shape[0] 405 N = Y.shape[0] // 3 406 self._fc2 = np.array( 407 np.transpose(fc2.reshape(N, 3, N, 3), axes=[0, 2, 1, 3]), 408 dtype='double', order='C') 409 410 411class ThirdOrderFC(object): 412 r"""SSCHA third order force constants 413 414 Eq. 45a in Ref.1 (See top docstring of this file) 415 416 \Phi_{abc} = - \sum_{pq} \Upsilon_{ap} \Upsilon_{bq} 417 \left< u^p u^q \mathfrak{f}_c \right>_{tilde{\rho}_{\mathfrak{R},\Phi}} 418 419 \mathfrak{f}_i = f_i - \left[ 420 \left< f_i \right>_{\tilde{\rho}_{\mathfrak{R},\Phi}} 421 - \sum_j \Phi_{ij}u^j \right] 422 423 Attributes 424 ---------- 425 displacements : ndarray 426 shape=(snap_shots, 3 * num_satoms), dtype='double', order='C' 427 forces : ndarray 428 shape=(snap_shots, 3 * num_satoms), dtype='double', order='C' 429 fc3 : ndarray 430 shape=(num_satom, num_satom, num_satom, 3, 3, 3) 431 432 """ 433 434 def __init__(self, 435 displacements, 436 forces, 437 supercell_phonon, 438 cutoff_frequency=1e-5, 439 log_level=0): 440 """ 441 442 Parameters 443 ---------- 444 displacements : ndarray 445 shape=(snap_shots, num_satom, 3), dtype='double', order='C' 446 forces : ndarray 447 shape=(snap_shots, num_satom, 3), dtype='double', order='C' 448 upsilon_matrix : DispCorrMatrix 449 Displacement-displacement correlation matrix class instance. 450 cutoff_frequency : float 451 Phonons are ignored if they have frequencies less than this value. 452 453 """ 454 455 assert (displacements.shape == forces.shape) 456 shape = displacements.shape 457 u = np.array(displacements.reshape(shape[0], -1), 458 dtype='double', order='C') 459 f = np.array(forces.reshape(shape[0], -1), 460 dtype='double', order='C') 461 self._displacements = u 462 self._forces = f 463 self._uu = DispCorrMatrix( 464 supercell_phonon, cutoff_frequency=cutoff_frequency) 465 self._force_constants = supercell_phonon.force_constants 466 self._cutoff_frequency = cutoff_frequency 467 self._log_level = log_level 468 469 self._drift = u.sum(axis=0) / u.shape[0] 470 self._fmat = None 471 self._fc3 = None 472 473 @property 474 def displacements(self): 475 return self._displacements 476 477 @property 478 def forces(self): 479 return self._forces 480 481 @property 482 def fc3(self): 483 return self._fc3 484 485 @property 486 def ff(self): 487 return self._fmat 488 489 def run(self, T=300.0): 490 if self._fmat is None: 491 self._fmat = self._run_fmat() 492 493 fc3 = self._run_fc3_ave(T) 494 N = fc3.shape[0] // 3 495 fc3 = fc3.reshape((N, 3, N, 3, N, 3)) 496 self._fc3 = np.array( 497 np.transpose(fc3, axes=[0, 2, 4, 1, 3, 5]), 498 dtype='double', order='C') 499 500 def _run_fmat(self): 501 f = self._forces 502 u = self._displacements 503 fc2 = self._force_constants 504 return f - f.sum(axis=0) / f.shape[0] + np.dot(u, fc2) 505 506 def _run_fc3_ave(self, T): 507 # self._uu.run(T) 508 # Y = self._uu.upsilon_matrix 509 f = self._fmat 510 u = self._displacements 511 Y = np.linalg.pinv(np.dot(u.T, u) / u.shape[0]) 512 513 # This is faster than multiplying Y after ansemble average at least 514 # in python implementation. 515 u_inv = np.dot(u, Y) 516 517 if self._log_level: 518 N = np.prod(u.shape) 519 print("rms u_inv:", np.sqrt((u_inv ** 2).sum() / N)) 520 print("rms u:", np.sqrt((u ** 2).sum() / N)) 521 print("rms forces:", np.sqrt((self._forces ** 2).sum() / N)) 522 print("rms f:", np.sqrt((f ** 2).sum() / N)) 523 524 return - np.einsum('li,lj,lk->ijk', u_inv, u_inv, f) / f.shape[0] 525