mirror of https://github.com/phonopy/phono3py.git
SSCHA code was updated.
This commit is contained in:
Atsushi Togo
parent
d37fe70e6b
commit
3801539196
|
|
@ -36,8 +36,8 @@
|
|||
|
||||
Formulae implemented are based on these papers:
|
||||
|
||||
Bianco et al. https://doi.org/10.1103/PhysRevB.96.014111
|
||||
Aseginolaza et al. https://doi.org/10.1103/PhysRevB.100.214307
|
||||
Ref. 1: Bianco et al. https://doi.org/10.1103/PhysRevB.96.014111
|
||||
Ref. 2: Aseginolaza et al. https://doi.org/10.1103/PhysRevB.100.214307
|
||||
|
||||
"""
|
||||
|
||||
|
|
@ -48,7 +48,52 @@ from phono3py.phonon.func import mode_length
|
|||
|
||||
|
||||
class SupercellPhonon(object):
|
||||
def __init__(self, supercell, force_constants, factor=VaspToTHz):
|
||||
"""Supercell phonon class
|
||||
|
||||
Dynamical matrix is created for supercell atoms and solved in real.
|
||||
All phonons at commensurate points are folded to those at Gamma point.
|
||||
Phonon eigenvectors are represented in real type.
|
||||
|
||||
Attributes
|
||||
----------
|
||||
eigenvalues : ndarray
|
||||
Phonon eigenvalues of supercell dynamical matrix.
|
||||
shape=(3 * num_satom, ), dtype='double', order='C'
|
||||
eigenvectors : ndarray
|
||||
Phonon eigenvectors of supercell dynamical matrix.
|
||||
shape=(3 * num_satom, 3 * num_satom), dtype='double', order='C'
|
||||
frequencies : ndarray
|
||||
Phonon frequencies of supercell dynamical matrix. Frequency conversion
|
||||
factor to THz is multiplied.
|
||||
shape=(3 * num_satom, ), dtype='double', order='C'
|
||||
dynamical_matrix : ndarray
|
||||
Numerical matrix corresponding to supercell dynamical matrix given
|
||||
at Gamma point. Therefore no phase factor is multiplied, i.e.,
|
||||
force constants weighted by atomic masses, 1/sqrt(M_a, M_b).
|
||||
shape=(3 * num_satom, 3 * num_satom), dtype='double', order='C'
|
||||
supercell : PhonopyAtoms or its derived class
|
||||
Supercell.
|
||||
|
||||
"""
|
||||
|
||||
def __init__(self,
|
||||
supercell,
|
||||
force_constants,
|
||||
frequency_factor_to_THz=VaspToTHz):
|
||||
"""
|
||||
|
||||
Parameters
|
||||
----------
|
||||
supercell : PhonopyAtoms
|
||||
Supercell.
|
||||
force_constants : array_like
|
||||
Second order force constants.
|
||||
shape=(num_satom, num_satom, 3, 3), dtype='double', order='C'
|
||||
frequency_factor_to_THz : float
|
||||
Frequency conversion factor to THz.
|
||||
|
||||
"""
|
||||
|
||||
self._supercell = supercell
|
||||
_fc2 = np.swapaxes(force_constants, 1, 2)
|
||||
_fc2 = _fc2.reshape(-1, np.prod(_fc2.shape[-2:]))
|
||||
|
|
@ -56,7 +101,8 @@ class SupercellPhonon(object):
|
|||
dynmat = np.array(_fc2 / np.sqrt(np.outer(masses, masses)),
|
||||
dtype='double', order='C')
|
||||
eigvals, eigvecs = np.linalg.eigh(dynmat)
|
||||
freqs = np.sqrt(np.abs(eigvals)) * np.sign(eigvals) * factor
|
||||
freqs = np.sqrt(np.abs(eigvals)) * np.sign(eigvals)
|
||||
freqs *= frequency_factor_to_THz
|
||||
self._eigenvalues = np.array(eigvals, dtype='double', order='C')
|
||||
self._eigenvectors = np.array(eigvecs, dtype='double', order='C')
|
||||
self._frequencies = np.array(freqs, dtype='double', order='C')
|
||||
|
|
@ -84,11 +130,30 @@ class SupercellPhonon(object):
|
|||
|
||||
|
||||
class DispCorrMatrix(object):
|
||||
"""Calculate gamma matrix"""
|
||||
"""Calculate Upsilon matrix
|
||||
|
||||
Attributes
|
||||
----------
|
||||
upsilon_matrix : ndarray
|
||||
Displacement-displacement correlation matrix at temperature.
|
||||
Physical unit is [1/Angstrom^2].
|
||||
shape=(3 * num_satom, 3 * num_satom), dtype='double', order='C'
|
||||
|
||||
"""
|
||||
|
||||
def __init__(self, supercell_phonon):
|
||||
"""
|
||||
|
||||
Parameters
|
||||
----------
|
||||
supercell_phonon : SupercellPhonon
|
||||
Supercell phonon object. Phonons at Gamma point, where
|
||||
eigenvectors are not complex type but real type.
|
||||
|
||||
"""
|
||||
|
||||
self._supercell_phonon = supercell_phonon
|
||||
self._gamma_matrix = None
|
||||
self._upsilon_matrix = None
|
||||
|
||||
def run(self, T):
|
||||
freqs = self._supercell_phonon.frequencies
|
||||
|
|
@ -96,11 +161,13 @@ class DispCorrMatrix(object):
|
|||
a = mode_length(freqs, T)
|
||||
masses = np.repeat(self._supercell_phonon.supercell.masses, 3)
|
||||
gamma = np.dot(eigvecs, np.dot(np.diag(1.0 / a ** 2), eigvecs.T))
|
||||
self._gamma_matrix = gamma * np.sqrt(np.outer(masses, masses))
|
||||
self._upsilon_matrix = np.array(
|
||||
gamma * np.sqrt(np.outer(masses, masses)),
|
||||
dtype='double', order='C')
|
||||
|
||||
@property
|
||||
def gamma_matrix(self):
|
||||
return self._gamma_matrix
|
||||
def upsilon_matrix(self):
|
||||
return self._upsilon_matrix
|
||||
|
||||
|
||||
class DispCorrMatrixMesh(object):
|
||||
|
|
@ -121,7 +188,7 @@ class DispCorrMatrixMesh(object):
|
|||
def commensurate_points(self):
|
||||
return self._d2f.commensurate_points
|
||||
|
||||
def create_gamma_matrix(self, frequencies, eigenvectors, T):
|
||||
def create_upsilon_matrix(self, frequencies, eigenvectors, T):
|
||||
"""
|
||||
|
||||
Parameters
|
||||
|
|
@ -142,5 +209,68 @@ class DispCorrMatrixMesh(object):
|
|||
self._d2f.run()
|
||||
|
||||
@property
|
||||
def gamma_matrix(self):
|
||||
def upsilon_matrix(self):
|
||||
return self._d2f.force_constants
|
||||
|
||||
|
||||
class ThirdOrderFC(object):
|
||||
r"""SSCHA third order force constants
|
||||
|
||||
Eq. 45a in Ref.1 (See top docstring of this file)
|
||||
|
||||
\Phi_{abc} = - \sum_{pq} \Upsilon_{ap} \Upsilon_{bq}
|
||||
\left< u^p u^q \mathfrak{f}_c \right>_{tilde{\rho}_{\mathfrak{R},\Phi}}
|
||||
|
||||
\mathfrak{f}_i = f_i - \left[
|
||||
\left< f_i \right>_{\tilde{\rho}_{\mathfrak{R},\Phi}}
|
||||
- \sum_j \Phi_{ij}u^j \right]
|
||||
|
||||
Attributes
|
||||
----------
|
||||
displacements : ndarray
|
||||
shape=(3 * num_satoms, snap_shots), dtype='double', order='C'
|
||||
forces : ndarray
|
||||
shape=(3 * num_satoms, snap_shots), dtype='double', order='C'
|
||||
fc3 : ndarray
|
||||
shape=(num_satom, num_satom, num_satom, 3, 3, 3)
|
||||
|
||||
"""
|
||||
|
||||
def __init__(self, displacements, forces, upsilon_matrix):
|
||||
"""
|
||||
|
||||
Parameters
|
||||
----------
|
||||
upsilon_matrix : ndarray
|
||||
Displacement-displacement correlation matrix at a temperature
|
||||
point.
|
||||
shape=(3 * num_satom, 3 * num_satom), dtype='double', order='C'
|
||||
displacements : ndarray
|
||||
shape=(snap_shots, num_satom, 3), dtype='double', order='C'
|
||||
forces : ndarray
|
||||
shape=(snap_shots, num_satom, 3), dtype='double', order='C'
|
||||
|
||||
"""
|
||||
|
||||
self._upsilon_matrix = upsilon_matrix
|
||||
assert (displacements.shape == forces.shape)
|
||||
shape = displacements.shape
|
||||
self._displacements = np.array(displacements.reshape(-1, shape[0]),
|
||||
dtype='double', order='C')
|
||||
self._forces = np.array(forces.reshape(-1, shape[0]),
|
||||
dtype='double', order='C')
|
||||
|
||||
@property
|
||||
def displacements(self):
|
||||
return self._displacements
|
||||
|
||||
@property
|
||||
def displacements(self):
|
||||
return self._forces
|
||||
|
||||
@property
|
||||
def fc3(self):
|
||||
return self._fc3
|
||||
|
||||
def run(self):
|
||||
pass
|
||||
|
|
|
|||
|
|
@ -1,5 +1,6 @@
|
|||
import os
|
||||
import pytest
|
||||
import phonopy
|
||||
import phono3py
|
||||
|
||||
current_dir = os.path.dirname(os.path.abspath(__file__))
|
||||
|
|
@ -18,3 +19,10 @@ def si_pbesol():
|
|||
def si_pbesol_111():
|
||||
yaml_filename = os.path.join(current_dir, "phono3py_params_Si111.yaml")
|
||||
return phono3py.load(yaml_filename, log_level=1)
|
||||
|
||||
|
||||
@pytest.fixture(scope='session')
|
||||
def si_pbesol_iterha_111():
|
||||
yaml_filename = os.path.join(current_dir,
|
||||
"phonopy_params-Si111-iterha.yaml.xz")
|
||||
return phonopy.load(yaml_filename, log_level=1, produce_fc=False)
|
||||
|
|
|
|||
Binary file not shown.
|
|
@ -1,15 +1,16 @@
|
|||
import pytest
|
||||
import numpy as np
|
||||
from phono3py.sscha.sscha import (
|
||||
DispCorrMatrix, DispCorrMatrixMesh, SupercellPhonon)
|
||||
DispCorrMatrix, DispCorrMatrixMesh, SupercellPhonon, ThirdOrderFC)
|
||||
from phonopy.phonon.qpoints import QpointsPhonon
|
||||
|
||||
|
||||
si_pbesol_gamma0_0 = [[3.849187e+02, 0, 0],
|
||||
[0, 3.849187e+02, 0],
|
||||
[0, 0, 3.849187e+02]]
|
||||
si_pbesol_gamma1_34 = [[1.886404, -1.549705, -1.126055],
|
||||
[-1.549705, 1.886404, -1.126055],
|
||||
[-1.126055, -1.126055, -0.006187]]
|
||||
si_pbesol_upsilon0_0 = [[3.849187e+02, 0, 0],
|
||||
[0, 3.849187e+02, 0],
|
||||
[0, 0, 3.849187e+02]]
|
||||
si_pbesol_upsilon1_34 = [[1.886404, -1.549705, -1.126055],
|
||||
[-1.549705, 1.886404, -1.126055],
|
||||
[-1.126055, -1.126055, -0.006187]]
|
||||
si_pbesol_111_freqs = [
|
||||
0.00000, 0.00000, 0.00000, 4.02839, 4.02839, 4.02839,
|
||||
4.02839, 4.02839, 4.02839, 12.13724, 12.13724, 12.13724,
|
||||
|
|
@ -23,7 +24,7 @@ def get_supercell_phonon(ph3):
|
|||
fc2 = ph3.dynamical_matrix.force_constants
|
||||
supercell = ph3.phonon_supercell
|
||||
factor = ph3.unit_conversion_factor
|
||||
return SupercellPhonon(supercell, fc2, factor=factor)
|
||||
return SupercellPhonon(supercell, fc2, frequency_factor_to_THz=factor)
|
||||
|
||||
|
||||
def test_SupercellPhonon(si_pbesol_111):
|
||||
|
|
@ -32,30 +33,48 @@ def test_SupercellPhonon(si_pbesol_111):
|
|||
si_pbesol_111_freqs, sph.frequencies, atol=1e-4)
|
||||
|
||||
|
||||
def test_gamma_matrix_mesh(si_pbesol):
|
||||
def test_upsilon_matrix_mesh(si_pbesol):
|
||||
si_pbesol.mesh_numbers = [9, 9, 9]
|
||||
si_pbesol.init_phph_interaction()
|
||||
dynmat = si_pbesol.dynamical_matrix
|
||||
gmat = DispCorrMatrixMesh(dynmat.primitive, dynmat.supercell)
|
||||
qpoints_phonon = QpointsPhonon(gmat.commensurate_points,
|
||||
upmat = DispCorrMatrixMesh(dynmat.primitive, dynmat.supercell)
|
||||
qpoints_phonon = QpointsPhonon(upmat.commensurate_points,
|
||||
dynmat,
|
||||
with_eigenvectors=True)
|
||||
freqs = qpoints_phonon.frequencies
|
||||
eigvecs = qpoints_phonon.eigenvectors
|
||||
gmat.create_gamma_matrix(freqs, eigvecs, 300.0)
|
||||
gmat.run()
|
||||
upmat.create_upsilon_matrix(freqs, eigvecs, 300.0)
|
||||
upmat.run()
|
||||
np.testing.assert_allclose(
|
||||
si_pbesol_gamma0_0, gmat.gamma_matrix[0, 0], atol=1e-4)
|
||||
si_pbesol_upsilon0_0, upmat.upsilon_matrix[0, 0], atol=1e-4)
|
||||
np.testing.assert_allclose(
|
||||
si_pbesol_gamma1_34, gmat.gamma_matrix[1, 34], atol=1e-4)
|
||||
si_pbesol_upsilon1_34, upmat.upsilon_matrix[1, 34], atol=1e-4)
|
||||
|
||||
|
||||
def test_gamma_matrix(si_pbesol):
|
||||
def test_upsilon_matrix(si_pbesol):
|
||||
supercell_phonon = get_supercell_phonon(si_pbesol)
|
||||
gmat = DispCorrMatrix(supercell_phonon)
|
||||
gmat.run(300.0)
|
||||
upmat = DispCorrMatrix(supercell_phonon)
|
||||
upmat.run(300.0)
|
||||
np.testing.assert_allclose(
|
||||
si_pbesol_gamma0_0, gmat.gamma_matrix[0:3, 0:3], atol=1e-4)
|
||||
si_pbesol_upsilon0_0, upmat.upsilon_matrix[0:3, 0:3], atol=1e-4)
|
||||
np.testing.assert_allclose(
|
||||
si_pbesol_gamma1_34, gmat.gamma_matrix[1 * 3: 2 * 3, 34 * 3: 35 * 3],
|
||||
si_pbesol_upsilon1_34, upmat.upsilon_matrix[1 * 3: 2 * 3, 34 * 3: 35 * 3],
|
||||
atol=1e-4)
|
||||
|
||||
|
||||
def test_fc3(si_pbesol_iterha_111):
|
||||
try:
|
||||
import alm
|
||||
except ModuleNotFoundError:
|
||||
pytest.skip("Skip this test because ALM module was not found.")
|
||||
|
||||
ph = si_pbesol_iterha_111
|
||||
ph.produce_force_constants(calculate_full_force_constants=True,
|
||||
fc_calculator='alm')
|
||||
supercell_phonon = SupercellPhonon(
|
||||
ph.supercell, ph.force_constants,
|
||||
frequency_factor_to_THz=ph.unit_conversion_factor)
|
||||
upmat = DispCorrMatrix(supercell_phonon)
|
||||
upmat.run(300.0)
|
||||
|
||||
ThirdOrderFC(ph.displacements, ph.forces, upmat.upsilon_matrix)
|
||||
|
|
|
|||
Loading…
Reference in New Issue