diff --git a/phono3py/sscha/sscha.py b/phono3py/sscha/sscha.py index 8db4310e..a3267de0 100644 --- a/phono3py/sscha/sscha.py +++ b/phono3py/sscha/sscha.py @@ -172,11 +172,10 @@ class SupercellPhonon(object): def __init__(self, supercell, force_constants, factor=VaspToTHz): self._supercell = supercell _fc2 = np.swapaxes(force_constants, 1, 2) - self._force_constants = np.array( - _fc2.reshape(-1, np.prod(_fc2.shape[-2:])), - dtype='double', order='C') + _fc2 = _fc2.reshape(-1, np.prod(_fc2.shape[-2:])) masses = np.repeat(supercell.masses, 3) - dynmat = self._force_constants / np.sqrt(np.outer(masses, masses)) + 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 self._eigenvalues = np.array(eigvals, dtype='double', order='C') @@ -204,10 +203,6 @@ class SupercellPhonon(object): def supercell(self): return self._supercell - @property - def force_constants(self): - return self._force_constants - class DispCorrMatrixMesh(object): """Calculate gamma matrix diff --git a/test/conftest.py b/test/conftest.py index a588bc86..09d4a90a 100644 --- a/test/conftest.py +++ b/test/conftest.py @@ -12,3 +12,9 @@ def si_pbesol(): return phono3py.load(yaml_filename, forces_fc3_filename=forces_fc3_filename, log_level=1) + + +@pytest.fixture(scope='session') +def si_pbesol_111(): + yaml_filename = os.path.join(current_dir, "phono3py_params_Si111.yaml") + return phono3py.load(yaml_filename, log_level=1) diff --git a/test/sscha/test_sscha.py b/test/sscha/test_sscha.py index 02f8c7f8..7a77921b 100644 --- a/test/sscha/test_sscha.py +++ b/test/sscha/test_sscha.py @@ -1,6 +1,6 @@ import numpy as np from phono3py.sscha.sscha import ( - DispCorrMatrix, DispCorrMatrixMesh, SupercellPhonon) + DispCorrMatrix, DispCorrMatrixMesh, SupercellPhonon, LambdaTensor) from phonopy.phonon.qpoints import QpointsPhonon @@ -10,6 +10,29 @@ si_pbesol_gamma0_0 = [[3.849187e+02, 0, 0], si_pbesol_gamma1_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, + 12.13724, 12.13724, 12.13724, 13.71746, 13.71746, 13.71746, + 13.71746, 13.71746, 13.71746, 15.24974, 15.24974, 15.24974] + + +def get_supercell_phonon(ph3): + ph3.mesh_numbers = [1, 1, 1] + ph3.init_phph_interaction() + fc2 = ph3.dynamical_matrix.force_constants + supercell = ph3.phonon_supercell + factor = ph3.unit_conversion_factor + return SupercellPhonon(supercell, fc2, factor=factor) + + +# def test_run_Lambda(si_pbesol_111): +# lt = LambdaTensor() + +def test_SupercellPhonon(si_pbesol_111): + sph = get_supercell_phonon(si_pbesol_111) + np.testing.assert_allclose( + si_pbesol_111_freqs, sph.frequencies, atol=1e-4) def test_gamma_matrix_mesh(si_pbesol): @@ -24,28 +47,18 @@ def test_gamma_matrix_mesh(si_pbesol): eigvecs = qpoints_phonon.eigenvectors gmat.create_gamma_matrix(freqs, eigvecs, 300.0) gmat.run() - np.testing.assert_allclose(si_pbesol_gamma0_0, - gmat.gamma_matrix[0, 0], - atol=1e-4) - np.testing.assert_allclose(si_pbesol_gamma1_34, - gmat.gamma_matrix[1, 34], - atol=1e-4) - print(gmat.gamma_matrix[1, 2]) + np.testing.assert_allclose( + si_pbesol_gamma0_0, gmat.gamma_matrix[0, 0], atol=1e-4) + np.testing.assert_allclose( + si_pbesol_gamma1_34, gmat.gamma_matrix[1, 34], atol=1e-4) def test_gamma_matrix(si_pbesol): - si_pbesol.mesh_numbers = [1, 1, 1] - si_pbesol.init_phph_interaction() - fc2 = si_pbesol.dynamical_matrix.force_constants - supercell = si_pbesol.phonon_supercell - factor = si_pbesol.unit_conversion_factor - supercell_phonon = SupercellPhonon(supercell, fc2, factor=factor) + supercell_phonon = get_supercell_phonon(si_pbesol) gmat = DispCorrMatrix(supercell_phonon) gmat.run(300.0) - np.testing.assert_allclose(si_pbesol_gamma0_0, - gmat.gamma_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], - atol=1e-4) - print(gmat.gamma_matrix[1 * 3: 2 * 3, 34 * 3: 35 * 3]) + np.testing.assert_allclose( + si_pbesol_gamma0_0, gmat.gamma_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], + atol=1e-4)