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