diff --git a/phono3py/phonon/velocity_operator.py b/phono3py/phonon/velocity_operator.py new file mode 100644 index 00000000..840c0765 --- /dev/null +++ b/phono3py/phonon/velocity_operator.py @@ -0,0 +1,233 @@ +"""Velocity operator of Simoncelli, Marzari, and Mauri.""" +# Copyright (C) 2013 Atsushi Togo +# All rights reserved. +# +# This file is part of phonopy. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions +# are met: +# +# * Redistributions of source code must retain the above copyright +# notice, this list of conditions and the following disclaimer. +# +# * Redistributions in binary form must reproduce the above copyright +# notice, this list of conditions and the following disclaimer in +# the documentation and/or other materials provided with the +# distribution. +# +# * Neither the name of the phonopy project nor the names of its +# contributors may be used to endorse or promote products derived +# from this software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS +# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE +# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, +# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, +# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN +# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. + +import numpy as np +from phonopy.phonon.group_velocity import GroupVelocity +from phonopy.units import VaspToTHz + + +class VelocityOperator(GroupVelocity): + """Class to calculate velocity operator of phonons.""" + + def __init__( + self, + dynamical_matrix, + q_length=None, + symmetry=None, + frequency_factor_to_THz=VaspToTHz, + cutoff_frequency=1e-4, + ): + """Init method. + + dynamical_matrix : DynamicalMatrix or DynamicalMatrixNAC + Dynamical matrix class instance. + q_length : float + This is used such as D(q + q_length) - D(q - q_length) for + calculating finite difference of dynamical matrix. + Default is None, which gives 1e-5. + symmetry : Symmetry + This is used to symmetrize group velocity at each q-points. + Default is None, which means no symmetrization. + frequency_factor_to_THz : float + Unit conversion factor to convert to THz. Default is VaspToTHz. + cutoff_frequency : float + Group velocity is set zero if phonon frequency is below this value. + + """ + self._dynmat = dynamical_matrix + primitive = dynamical_matrix.primitive + self._reciprocal_lattice_inv = primitive.cell + self._reciprocal_lattice = np.linalg.inv(self._reciprocal_lattice_inv) + self._q_length = q_length + if self._q_length is None: + self._q_length = 5e-6 + self._symmetry = symmetry + self._factor = frequency_factor_to_THz + self._cutoff_frequency = cutoff_frequency + + self._directions = np.array( + [ + [ + 1, + 2, + 3, + ], + # this is a random direction, not used and left here for historical + # reasons. + [1, 0, 0], # x + [0, 1, 0], # y + [0, 0, 1], + ], + dtype="double", + ) # z + self._directions[0] /= np.linalg.norm( + self._directions[0] + ) # normalize the random direction + + self._q_points = None + self._velocity_operators = None + self._perturbation = None + + def run(self, q_points, perturbation=None): + """Velocity operators are computed at q-points. + + q_points : Array-like + List of q-points such as [[0, 0, 0], [0.1, 0.2, 0.3], ...]. + perturbation : Array-like + Direction in fractional coordinates of reciprocal space. + + """ + self._q_points = q_points + self._perturbation = perturbation + if perturbation is None: + # Give a random direction to break symmetry + self._directions[0] = np.array([1, 0, 0]) # normalized later + else: + self._directions[0] = np.dot(self._reciprocal_lattice, perturbation) + self._directions[0] /= np.linalg.norm(self._directions[0]) + + gv_operator = [ + self._calculate_velocity_operator_at_q(q) for q in self._q_points + ] + self._velocity_operators = np.array(gv_operator, dtype="complex", order="C") + + @property + def velocity_operators(self): + """Return velocity operators.""" + return self._velocity_operators + + def _calculate_velocity_operator_at_q(self, q): + + self._dynmat.set_dynamical_matrix(q) + dm = self._dynmat.get_dynamical_matrix() + eigvals, eigvecs = np.linalg.eigh(dm) + eigvals = eigvals.real + freqs = np.sqrt(abs(eigvals)) * np.sign(eigvals) * self._factor + nat3 = len(freqs) + gv_operator = np.zeros((nat3, nat3, 3), dtype="complex", order="C") + # + # compute the finite differences along several directions + # ddms is a list of FD derivatives of the dynamical matrix + # computed along several directions + ddms = self._get_dsqrtD_FD(np.array(q)) + # + # ddms[0] cointains the FD derivative in the direction in which the velocity + # operator is diagonalized + for id_dir in range(0, 3): + gv_operator[:, :, id_dir] = ( + np.matmul(eigvecs.T.conj(), np.matmul(ddms[id_dir + 1], eigvecs)) + * self._factor + ) + # + # enforce the velocity operator to be Hermitian + gv_operator[:, :, id_dir] = ( + gv_operator[:, :, id_dir].T.conj() + gv_operator[:, :, id_dir] + ) / 2.0 + + return gv_operator + + def _get_dsqrtD_FD(self, q): + """Compute finite difference of sqrt of dynamical matrices.""" + # + ddm = [] + # _q_length is a float specifying the modulus of the q-vector displacement + # used to compute the finite differences + for dqc in self._directions * self._q_length: + dq = np.dot(self._reciprocal_lattice_inv, dqc) + ddm.append( + self._delta_sqrt_dynamical_matrix(q, dq, self._dynmat) + / self._q_length + / 2 + ) + return np.array(ddm) + + def _sqrt_dynamical_matrix(self, flag_gamma, dm): + # returns the sqrt of the dynamical matrix in the cartesian basis + eigvals, eigvecs = np.linalg.eigh(dm) + eigvals = eigvals.real + # will give NaN in case of negative frequencies + freqs = np.sign(eigvals) * np.sqrt(abs(eigvals)) + # print('raw=',freqs) + # + # enforce eigenvalues and eigenvectors to zero for acoustic modes at Gamma point + if flag_gamma: + freqs[0] = 0.0 + eigvecs[:, 0] = 0.0 + freqs[1] = 0.0 + eigvecs[:, 1] = 0.0 + freqs[2] = 0.0 + eigvecs[:, 2] = 0.0 + # + if any(f < 0.0 for f in freqs): + print("ERROR: negative frequency=", freqs) + # this is a diagonal matrix containing the frequencies on the diagonal + # we take the element-wise sqrt of a diagonal matrix, in eigenmodes basis + omega_matrix = np.sqrt(np.matmul(eigvecs.T.conj(), np.matmul(dm, eigvecs))) + # now we go back to the Cartesian basis + sqrt_dm = np.matmul(eigvecs, np.matmul(omega_matrix, eigvecs.T.conj())) + # + return sqrt_dm + + def _delta_sqrt_dynamical_matrix(self, q, delta_q, dynmat): + # + flag_gamma = False + if np.linalg.norm(q) < np.linalg.norm(delta_q): + flag_gamma = True + + if ( + (self._dynmat.is_nac()) + and (self._dynmat.nac_method == "gonze") + and flag_gamma + ): + dynmat.set_dynamical_matrix( + q - delta_q, q_direction=(q - delta_q) / np.linalg.norm(q - delta_q) + ) + dm1 = dynmat.get_dynamical_matrix() + sqrt_dm1 = self._sqrt_dynamical_matrix(flag_gamma, dm1) + dynmat.set_dynamical_matrix( + q + delta_q, q_direction=(q + delta_q) / np.linalg.norm(q + delta_q) + ) + dm2 = dynmat.get_dynamical_matrix() + sqrt_dm2 = self._sqrt_dynamical_matrix(flag_gamma, dm2) + else: + dynmat.set_dynamical_matrix(q - delta_q) + dm1 = dynmat.get_dynamical_matrix() + sqrt_dm1 = self._sqrt_dynamical_matrix(flag_gamma, dm1) + dynmat.set_dynamical_matrix(q + delta_q) + dm2 = dynmat.get_dynamical_matrix() + sqrt_dm2 = self._sqrt_dynamical_matrix(flag_gamma, dm2) + # + # + return sqrt_dm2 - sqrt_dm1 diff --git a/test/BORN_NaCl b/test/BORN_NaCl new file mode 100644 index 00000000..af42902d --- /dev/null +++ b/test/BORN_NaCl @@ -0,0 +1,4 @@ +14.400 +2.43533967 0 0 0 2.43533967 0 0 0 2.43533967 +1.086875 0 0 0 1.086875 0 0 0 1.086875 +-1.086875 0 0 0 -1.086875 0 0 0 -1.086875 diff --git a/test/FORCE_SETS_NaCl b/test/FORCE_SETS_NaCl new file mode 100644 index 00000000..eb9f19a4 --- /dev/null +++ b/test/FORCE_SETS_NaCl @@ -0,0 +1,136 @@ +64 +2 + +1 + 0.0100000000000000 0.0000000000000000 0.0000000000000000 + -0.0180619400 0.0000000000 0.0000000000 + 0.0030240400 0.0000000000 0.0000000000 + -0.0000434000 0.0000000000 0.0000000000 + 0.0001066900 0.0000000000 0.0000000000 + -0.0000434000 0.0000000000 0.0000000000 + 0.0001066900 0.0000000000 0.0000000000 + -0.0001589900 0.0000000000 0.0000000000 + 0.0001344200 0.0000000000 0.0000000000 + -0.0005716800 -0.0000021200 -0.0000021200 + 0.0003817600 0.0000002400 0.0000002400 + -0.0005716800 0.0000021200 -0.0000021200 + 0.0003817600 -0.0000002400 0.0000002400 + -0.0005716800 -0.0000021200 0.0000021200 + 0.0003817600 0.0000002400 -0.0000002400 + -0.0005716800 0.0000021200 0.0000021200 + 0.0003817600 -0.0000002400 -0.0000002400 + 0.0004845700 0.0000000000 0.0019354900 + 0.0004850100 0.0000000000 -0.0019214300 + -0.0000183800 0.0000000000 -0.0000178800 + -0.0000173700 0.0000000000 0.0000170400 + 0.0004845700 0.0000000000 -0.0019354900 + 0.0004850100 0.0000000000 0.0019214300 + -0.0000183800 0.0000000000 0.0000178800 + -0.0000173700 0.0000000000 -0.0000170400 + 0.0004845700 0.0019354900 0.0000000000 + 0.0004850100 -0.0019214300 0.0000000000 + 0.0004845700 -0.0019354900 0.0000000000 + 0.0004850100 0.0019214300 0.0000000000 + -0.0000183800 -0.0000178800 0.0000000000 + -0.0000173700 0.0000170400 0.0000000000 + -0.0000183800 0.0000178800 0.0000000000 + -0.0000173700 -0.0000170400 0.0000000000 + -0.0000810800 -0.0004282000 -0.0004282000 + -0.0000835100 0.0004258800 0.0004258800 + -0.0000810800 0.0004282000 -0.0004282000 + -0.0000835100 -0.0004258800 0.0004258800 + -0.0000810800 -0.0004282000 0.0004282000 + -0.0000835100 0.0004258800 -0.0004258800 + -0.0000810800 0.0004282000 0.0004282000 + -0.0000835100 -0.0004258800 -0.0004258800 + 0.0050201700 0.0000000000 0.0000000000 + 0.0046140400 0.0000000000 0.0000000000 + 0.0000408300 0.0000000000 0.0000000000 + 0.0000373700 0.0000000000 0.0000000000 + 0.0000408300 0.0000000000 0.0000000000 + 0.0000373700 0.0000000000 0.0000000000 + -0.0000105300 0.0000000000 0.0000000000 + -0.0000118000 0.0000000000 0.0000000000 + 0.0017148800 -0.0000110200 0.0000000000 + -0.0009600300 -0.0000013600 0.0000000000 + 0.0017148800 0.0000110200 0.0000000000 + -0.0009600300 0.0000013600 0.0000000000 + 0.0001936100 0.0000011800 0.0000000000 + -0.0002359000 0.0000000900 0.0000000000 + 0.0001936100 -0.0000011800 0.0000000000 + -0.0002359000 -0.0000000900 0.0000000000 + 0.0017148800 0.0000000000 -0.0000110200 + -0.0009600300 0.0000000000 -0.0000013600 + 0.0001936100 0.0000000000 0.0000011800 + -0.0002359000 0.0000000000 0.0000000900 + 0.0017148800 0.0000000000 0.0000110200 + -0.0009600300 0.0000000000 0.0000013600 + 0.0001936100 0.0000000000 -0.0000011800 + -0.0002359000 0.0000000000 -0.0000000900 + +33 + 0.0100000000000000 0.0000000000000000 0.0000000000000000 + -0.0000721800 -0.0004247300 -0.0004247300 + -0.0000667600 0.0004272300 0.0004272300 + -0.0000721800 0.0004247300 -0.0004247300 + -0.0000667600 -0.0004272300 0.0004272300 + -0.0000721800 -0.0004247300 0.0004247300 + -0.0000667600 0.0004272300 -0.0004272300 + -0.0000721800 0.0004247300 0.0004247300 + -0.0000667600 -0.0004272300 -0.0004272300 + 0.0046345100 0.0000000000 0.0000000000 + 0.0050314800 0.0000000000 0.0000000000 + 0.0000499300 0.0000000000 0.0000000000 + 0.0000510900 0.0000000000 0.0000000000 + 0.0000499300 0.0000000000 0.0000000000 + 0.0000510900 0.0000000000 0.0000000000 + -0.0000030400 0.0000000000 0.0000000000 + 0.0000012000 0.0000000000 0.0000000000 + 0.0017234100 0.0000077800 0.0000000000 + -0.0009445700 -0.0000004700 0.0000000000 + 0.0017234100 -0.0000077800 0.0000000000 + -0.0009445700 0.0000004700 0.0000000000 + 0.0002042400 0.0000019500 0.0000000000 + -0.0002202400 -0.0000015500 0.0000000000 + 0.0002042400 -0.0000019500 0.0000000000 + -0.0002202400 0.0000015500 0.0000000000 + 0.0017234100 0.0000000000 0.0000077800 + -0.0009445700 0.0000000000 -0.0000004700 + 0.0002042400 0.0000000000 0.0000019500 + -0.0002202400 0.0000000000 -0.0000015500 + 0.0017234100 0.0000000000 -0.0000077800 + -0.0009445700 0.0000000000 0.0000004700 + 0.0002042400 0.0000000000 -0.0000019500 + -0.0002202400 0.0000000000 0.0000015500 + -0.0236731700 0.0000000000 0.0000000000 + 0.0013201600 0.0000000000 0.0000000000 + -0.0005159800 0.0000000000 0.0000000000 + 0.0006463500 0.0000000000 0.0000000000 + -0.0005159800 0.0000000000 0.0000000000 + 0.0006463500 0.0000000000 0.0000000000 + -0.0001691500 0.0000000000 0.0000000000 + 0.0001742700 0.0000000000 0.0000000000 + -0.0008016000 0.0000039800 0.0000039800 + 0.0005170500 -0.0000010700 -0.0000010700 + -0.0008016000 -0.0000039800 0.0000039800 + 0.0005170500 0.0000010700 -0.0000010700 + -0.0008016000 0.0000039800 -0.0000039800 + 0.0005170500 -0.0000010700 0.0000010700 + -0.0008016000 -0.0000039800 -0.0000039800 + 0.0005170500 0.0000010700 0.0000010700 + 0.0013674700 0.0000000000 0.0020506700 + 0.0013743700 0.0000000000 -0.0020695000 + -0.0000146100 0.0000000000 0.0001925900 + -0.0000113300 0.0000000000 -0.0001898400 + 0.0013674700 0.0000000000 -0.0020506700 + 0.0013743700 0.0000000000 0.0020695000 + -0.0000146100 0.0000000000 -0.0001925900 + -0.0000113300 0.0000000000 0.0001898400 + 0.0013674700 0.0020506700 0.0000000000 + 0.0013743700 -0.0020695000 0.0000000000 + 0.0013674700 -0.0020506700 0.0000000000 + 0.0013743700 0.0020695000 0.0000000000 + -0.0000146100 0.0001925900 0.0000000000 + -0.0000113300 -0.0001898400 0.0000000000 + -0.0000146100 -0.0001925900 0.0000000000 + -0.0000113300 0.0001898400 0.0000000000 diff --git a/test/conftest.py b/test/conftest.py index baa7eb68..0ef0aa89 100644 --- a/test/conftest.py +++ b/test/conftest.py @@ -3,6 +3,7 @@ import os import phonopy import pytest +from phonopy import Phonopy from phonopy.interface.phonopy_yaml import read_cell_yaml import phono3py @@ -253,3 +254,19 @@ def aln_lda(request): store_dense_svecs=enable_v2, log_level=1, ) + + +@pytest.fixture(scope="session") +def ph_nacl() -> Phonopy: + """Return Phonopy class instance of NaCl 2x2x2.""" + yaml_filename = os.path.join(current_dir, "phonopy_disp_NaCl.yaml") + force_sets_filename = os.path.join(current_dir, "FORCE_SETS_NaCl") + born_filename = os.path.join(current_dir, "BORN_NaCl") + return phonopy.load( + yaml_filename, + force_sets_filename=force_sets_filename, + born_filename=born_filename, + is_compact_fc=False, + log_level=1, + produce_fc=True, + ) diff --git a/test/phonon/test_velocity_operator.py b/test/phonon/test_velocity_operator.py new file mode 100644 index 00000000..00173342 --- /dev/null +++ b/test/phonon/test_velocity_operator.py @@ -0,0 +1,363 @@ +"""Tests for velocity operator.""" +import numpy as np +from phonopy import Phonopy + +from phono3py.phonon.velocity_operator import VelocityOperator + + +def test_gv_operator_nacl(ph_nacl: Phonopy): + """Test of GroupVelocity by NaCl.""" + gv_operator_square_modulus_ref = {} + # direction 0 + gv_operator_square_modulus_ref[0] = np.array( + [ + [ + 2.42970100e02, + -6.20236119e00, + -1.82999046e-01, + -2.09295390e-01, + 4.07782967e01, + 2.10142118e-01, + ], + [ + -6.20236119e00, + 2.73261377e02, + -3.71620328e02, + 8.44686890e-01, + -1.61985084e00, + -5.01748557e00, + ], + [ + -1.82999046e-01, + -3.71620328e02, + 1.03489135e03, + -2.58482318e01, + 1.88975807e-02, + 1.33902037e01, + ], + [ + -2.09295390e-01, + 8.44686890e-01, + -2.58482318e01, + 1.30173952e00, + -1.56714339e-01, + 3.74138211e00, + ], + [ + 4.07782967e01, + -1.61985084e00, + 1.88975807e-02, + -1.56714339e-01, + 1.07994337e01, + -2.56681660e00, + ], + [ + 2.10142118e-01, + -5.01748557e00, + 1.33902037e01, + 3.74138211e00, + -2.56681660e00, + 5.34677654e02, + ], + ] + ) + + # direction 1 + gv_operator_square_modulus_ref[1] = np.array( + [ + [190.26961913, 24.42080327, 18.66103712, 3.499692, 17.12151154, 7.82170714], + [ + 24.42080327, + 753.95642101, + 449.00523399, + 56.90864201, + -3.00513131, + 178.25703584, + ], + [ + 18.66103712, + 449.00523399, + 282.90229148, + 124.96202594, + -3.1864197, + 46.43438272, + ], + [ + 3.499692, + 56.90864201, + 124.96202594, + 628.61784617, + -13.41223828, + -240.42545741, + ], + [ + 17.12151154, + -3.00513131, + -3.1864197, + -13.41223828, + 4.49707039, + 2.33350832, + ], + [ + 7.82170714, + 178.25703584, + 46.43438272, + -240.42545741, + 2.33350832, + 394.18935848, + ], + ] + ) + # direction 2 + gv_operator_square_modulus_ref[2] = np.array( + [ + [ + 3.21847047e02, + -2.87032091e01, + -2.45808366e01, + -1.57578343e00, + -1.32856444e02, + -1.34389001e01, + ], + [ + -2.87032091e01, + 5.54329513e01, + 1.16917174e02, + -3.18740464e00, + 8.04928014e00, + 6.19917390e01, + ], + [ + -2.45808366e01, + 1.16917174e02, + 2.55127746e02, + 5.66979440e-01, + 9.79228317e00, + 7.96414618e01, + ], + [ + -1.57578343e00, + -3.18740464e00, + 5.66979440e-01, + 1.61522304e01, + 1.98704126e01, + -1.14635470e02, + ], + [ + -1.32856444e02, + 8.04928014e00, + 9.79228317e00, + 1.98704126e01, + 8.60482934e02, + -1.68196151e01, + ], + [ + -1.34389001e01, + 6.19917390e01, + 7.96414618e01, + -1.14635470e02, + -1.68196151e01, + 8.58940749e02, + ], + ] + ) + + gv_operator = VelocityOperator( + ph_nacl.dynamical_matrix, symmetry=ph_nacl.primitive_symmetry + ) + # we chose an 'ugly' q-point because we want to avoid degeneracies. + # degeneracies are tested in phono3py + gv_operator.run([[0.1, 0.22, 0.33]]) + square_modulus = np.zeros((6, 6, 3)) + for direction in range(2, 3): + square_modulus[:, :, direction] = np.matmul( + gv_operator.velocity_operators[0][:, :, direction], + gv_operator.velocity_operators[0][:, :, direction].conjugate().T, + ).real + np.testing.assert_allclose( + square_modulus[:, :, direction].ravel(), + gv_operator_square_modulus_ref[direction].ravel(), + atol=1e-4, + ) + + +def test_gv_operator_si(ph_si: Phonopy): + """Test of GroupVelocity by Si.""" + gv_operator_square_modulus_ref = {} + # direction 0 + gv_operator_square_modulus_ref[0] = np.array( + [ + [ + 2.73572808e03, + -4.06210738e00, + 1.33571471e03, + -1.11264018e00, + -2.72825378e02, + 3.23442510e-01, + ], + [ + -4.06210738e00, + 2.72065300e03, + -3.20535769e01, + -4.41527949e02, + -1.29661972e01, + 5.12752451e02, + ], + [ + 1.33571471e03, + -3.20535769e01, + 4.16058788e03, + -1.19524446e00, + -1.26829836e01, + -1.15732958e00, + ], + [ + -1.11264018e00, + -4.41527949e02, + -1.19524446e00, + 1.45914531e03, + 1.90180987e00, + -8.45069100e02, + ], + [ + -2.72825378e02, + -1.29661972e01, + -1.26829836e01, + 1.90180987e00, + 1.78584698e03, + -2.25355370e-01, + ], + [ + 3.23442510e-01, + 5.12752451e02, + -1.15732958e00, + -8.45069100e02, + -2.25355370e-01, + 7.25820137e02, + ], + ] + ) + + # direction 1 + gv_operator_square_modulus_ref[1] = np.array( + [ + [ + 1.95698890e03, + -8.55308045e00, + -6.97133428e02, + 8.21442883e-01, + -1.08587694e03, + -9.83370396e-02, + ], + [ + -8.55308045e00, + 2.20579768e03, + 4.90281444e00, + -2.83807879e02, + 3.01470115e00, + -7.45274905e02, + ], + [ + -6.97133428e02, + 4.90281444e00, + 2.96863421e03, + -2.26085394e-01, + -7.66976692e00, + 2.51720167e-01, + ], + [ + 8.21442883e-01, + -2.83807879e02, + -2.26085394e-01, + 4.24208287e02, + 6.40518671e-03, + 4.96247028e02, + ], + [ + -1.08587694e03, + 3.01470115e00, + -7.66976692e00, + 6.40518671e-03, + 8.39180146e02, + -8.39817213e-01, + ], + [ + -9.83370396e-02, + -7.45274905e02, + 2.51720167e-01, + 4.96247028e02, + -8.39817213e-01, + 7.45082308e02, + ], + ] + ) + + # direction 2 + gv_operator_square_modulus_ref[2] = np.array( + [ + [ + 1.53807537e03, + -3.89302418e00, + -1.75424484e02, + -1.65219640e00, + 1.31051315e02, + 2.03331851e-02, + ], + [ + -3.89302418e00, + 2.57583012e03, + 1.04592533e01, + -9.42612420e02, + -1.42765014e01, + 3.33706787e01, + ], + [ + -1.75424484e02, + 1.04592533e01, + 2.92924900e03, + -8.53691987e-02, + -2.44376137e02, + -1.00559293e00, + ], + [ + -1.65219640e00, + -9.42612420e02, + -8.53691987e-02, + 4.18562364e02, + -3.54778535e-01, + -1.90266774e02, + ], + [ + 1.31051315e02, + -1.42765014e01, + -2.44376137e02, + -3.54778535e-01, + 4.36890621e01, + 6.82908650e-01, + ], + [ + 2.03331851e-02, + 3.33706787e01, + -1.00559293e00, + -1.90266774e02, + 6.82908650e-01, + 1.51072731e03, + ], + ] + ) + + gv_operator = VelocityOperator( + ph_si.dynamical_matrix, symmetry=ph_si.primitive_symmetry + ) + gv_operator.run([[0.1, 0.22, 0.33]]) + square_modulus = np.zeros((6, 6, 3)) + for direction in range(2, 3): + square_modulus[:, :, direction] = np.matmul( + gv_operator.velocity_operators[0][:, :, direction], + gv_operator.velocity_operators[0][:, :, direction].conj().T, + ) + np.testing.assert_allclose( + square_modulus[:, :, direction].ravel(), + gv_operator_square_modulus_ref[direction].ravel(), + atol=1e-5, + ) diff --git a/test/phonopy_disp_NaCl.yaml b/test/phonopy_disp_NaCl.yaml new file mode 100644 index 00000000..2f7d82ec --- /dev/null +++ b/test/phonopy_disp_NaCl.yaml @@ -0,0 +1,356 @@ +phonopy: + version: 2.7.0 + frequency_unit_conversion_factor: 15.633302 + symmetry_tolerance: 1.00000e-05 + configuration: + cell_filename: "POSCAR-unitcell" + create_displacements: ".true." + primitive_axes: "auto" + dim: "2 2 2" + +physical_unit: + atomic_mass: "AMU" + length: "angstrom" + force_constants: "eV/angstrom^2" + +space_group: + type: "Fm-3m" + number: 225 + Hall_symbol: "-F 4 2 3" + +primitive_matrix: +- [ 0.000000000000000, 0.500000000000000, 0.500000000000000 ] +- [ 0.500000000000000, 0.000000000000000, 0.500000000000000 ] +- [ 0.500000000000000, 0.500000000000000, 0.000000000000000 ] + +supercell_matrix: +- [ 2, 0, 0 ] +- [ 0, 2, 0 ] +- [ 0, 0, 2 ] + +primitive_cell: + lattice: + - [ 0.000000000000000, 2.845150738087836, 2.845150738087836 ] # a + - [ 2.845150738087836, 0.000000000000000, 2.845150738087836 ] # b + - [ 2.845150738087836, 2.845150738087836, 0.000000000000000 ] # c + points: + - symbol: Na # 1 + coordinates: [ 0.000000000000000, 0.000000000000000, 0.000000000000000 ] + mass: 22.989769 + - symbol: Cl # 2 + coordinates: [ 0.500000000000000, 0.500000000000000, 0.500000000000000 ] + mass: 35.453000 + reciprocal_lattice: # without 2pi + - [ -0.175737613233118, 0.175737613233118, 0.175737613233118 ] # a* + - [ 0.175737613233118, -0.175737613233118, 0.175737613233118 ] # b* + - [ 0.175737613233118, 0.175737613233118, -0.175737613233118 ] # c* + +unit_cell: + lattice: + - [ 5.690301476175671, 0.000000000000000, 0.000000000000000 ] # a + - [ 0.000000000000000, 5.690301476175671, 0.000000000000000 ] # b + - [ 0.000000000000000, 0.000000000000000, 5.690301476175671 ] # c + points: + - symbol: Na # 1 + coordinates: [ 0.000000000000000, 0.000000000000000, 0.000000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Na # 2 + coordinates: [ 0.000000000000000, 0.500000000000000, 0.500000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Na # 3 + coordinates: [ 0.500000000000000, 0.000000000000000, 0.500000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Na # 4 + coordinates: [ 0.500000000000000, 0.500000000000000, 0.000000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Cl # 5 + coordinates: [ 0.500000000000000, 0.500000000000000, 0.500000000000000 ] + mass: 35.453000 + reduced_to: 5 + - symbol: Cl # 6 + coordinates: [ 0.500000000000000, 0.000000000000000, 0.000000000000000 ] + mass: 35.453000 + reduced_to: 5 + - symbol: Cl # 7 + coordinates: [ 0.000000000000000, 0.500000000000000, 0.000000000000000 ] + mass: 35.453000 + reduced_to: 5 + - symbol: Cl # 8 + coordinates: [ 0.000000000000000, 0.000000000000000, 0.500000000000000 ] + mass: 35.453000 + reduced_to: 5 + +supercell: + lattice: + - [ 11.380602952351342, 0.000000000000000, 0.000000000000000 ] # a + - [ 0.000000000000000, 11.380602952351342, 0.000000000000000 ] # b + - [ 0.000000000000000, 0.000000000000000, 11.380602952351342 ] # c + points: + - symbol: Na # 1 + coordinates: [ 0.000000000000000, 0.000000000000000, 0.000000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Na # 2 + coordinates: [ 0.500000000000000, 0.000000000000000, 0.000000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Na # 3 + coordinates: [ 0.000000000000000, 0.500000000000000, 0.000000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Na # 4 + coordinates: [ 0.500000000000000, 0.500000000000000, 0.000000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Na # 5 + coordinates: [ 0.000000000000000, 0.000000000000000, 0.500000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Na # 6 + coordinates: [ 0.500000000000000, 0.000000000000000, 0.500000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Na # 7 + coordinates: [ 0.000000000000000, 0.500000000000000, 0.500000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Na # 8 + coordinates: [ 0.500000000000000, 0.500000000000000, 0.500000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Na # 9 + coordinates: [ 0.000000000000000, 0.250000000000000, 0.250000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Na # 10 + coordinates: [ 0.500000000000000, 0.250000000000000, 0.250000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Na # 11 + coordinates: [ 0.000000000000000, 0.750000000000000, 0.250000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Na # 12 + coordinates: [ 0.500000000000000, 0.750000000000000, 0.250000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Na # 13 + coordinates: [ 0.000000000000000, 0.250000000000000, 0.750000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Na # 14 + coordinates: [ 0.500000000000000, 0.250000000000000, 0.750000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Na # 15 + coordinates: [ 0.000000000000000, 0.750000000000000, 0.750000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Na # 16 + coordinates: [ 0.500000000000000, 0.750000000000000, 0.750000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Na # 17 + coordinates: [ 0.250000000000000, 0.000000000000000, 0.250000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Na # 18 + coordinates: [ 0.750000000000000, 0.000000000000000, 0.250000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Na # 19 + coordinates: [ 0.250000000000000, 0.500000000000000, 0.250000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Na # 20 + coordinates: [ 0.750000000000000, 0.500000000000000, 0.250000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Na # 21 + coordinates: [ 0.250000000000000, 0.000000000000000, 0.750000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Na # 22 + coordinates: [ 0.750000000000000, 0.000000000000000, 0.750000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Na # 23 + coordinates: [ 0.250000000000000, 0.500000000000000, 0.750000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Na # 24 + coordinates: [ 0.750000000000000, 0.500000000000000, 0.750000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Na # 25 + coordinates: [ 0.250000000000000, 0.250000000000000, 0.000000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Na # 26 + coordinates: [ 0.750000000000000, 0.250000000000000, 0.000000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Na # 27 + coordinates: [ 0.250000000000000, 0.750000000000000, 0.000000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Na # 28 + coordinates: [ 0.750000000000000, 0.750000000000000, 0.000000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Na # 29 + coordinates: [ 0.250000000000000, 0.250000000000000, 0.500000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Na # 30 + coordinates: [ 0.750000000000000, 0.250000000000000, 0.500000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Na # 31 + coordinates: [ 0.250000000000000, 0.750000000000000, 0.500000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Na # 32 + coordinates: [ 0.750000000000000, 0.750000000000000, 0.500000000000000 ] + mass: 22.989769 + reduced_to: 1 + - symbol: Cl # 33 + coordinates: [ 0.250000000000000, 0.250000000000000, 0.250000000000000 ] + mass: 35.453000 + reduced_to: 33 + - symbol: Cl # 34 + coordinates: [ 0.750000000000000, 0.250000000000000, 0.250000000000000 ] + mass: 35.453000 + reduced_to: 33 + - symbol: Cl # 35 + coordinates: [ 0.250000000000000, 0.750000000000000, 0.250000000000000 ] + mass: 35.453000 + reduced_to: 33 + - symbol: Cl # 36 + coordinates: [ 0.750000000000000, 0.750000000000000, 0.250000000000000 ] + mass: 35.453000 + reduced_to: 33 + - symbol: Cl # 37 + coordinates: [ 0.250000000000000, 0.250000000000000, 0.750000000000000 ] + mass: 35.453000 + reduced_to: 33 + - symbol: Cl # 38 + coordinates: [ 0.750000000000000, 0.250000000000000, 0.750000000000000 ] + mass: 35.453000 + reduced_to: 33 + - symbol: Cl # 39 + coordinates: [ 0.250000000000000, 0.750000000000000, 0.750000000000000 ] + mass: 35.453000 + reduced_to: 33 + - symbol: Cl # 40 + coordinates: [ 0.750000000000000, 0.750000000000000, 0.750000000000000 ] + mass: 35.453000 + reduced_to: 33 + - symbol: Cl # 41 + coordinates: [ 0.250000000000000, 0.000000000000000, 0.000000000000000 ] + mass: 35.453000 + reduced_to: 33 + - symbol: Cl # 42 + coordinates: [ 0.750000000000000, 0.000000000000000, 0.000000000000000 ] + mass: 35.453000 + reduced_to: 33 + - symbol: Cl # 43 + coordinates: [ 0.250000000000000, 0.500000000000000, 0.000000000000000 ] + mass: 35.453000 + reduced_to: 33 + - symbol: Cl # 44 + coordinates: [ 0.750000000000000, 0.500000000000000, 0.000000000000000 ] + mass: 35.453000 + reduced_to: 33 + - symbol: Cl # 45 + coordinates: [ 0.250000000000000, 0.000000000000000, 0.500000000000000 ] + mass: 35.453000 + reduced_to: 33 + - symbol: Cl # 46 + coordinates: [ 0.750000000000000, 0.000000000000000, 0.500000000000000 ] + mass: 35.453000 + reduced_to: 33 + - symbol: Cl # 47 + coordinates: [ 0.250000000000000, 0.500000000000000, 0.500000000000000 ] + mass: 35.453000 + reduced_to: 33 + - symbol: Cl # 48 + coordinates: [ 0.750000000000000, 0.500000000000000, 0.500000000000000 ] + mass: 35.453000 + reduced_to: 33 + - symbol: Cl # 49 + coordinates: [ 0.000000000000000, 0.250000000000000, 0.000000000000000 ] + mass: 35.453000 + reduced_to: 33 + - symbol: Cl # 50 + coordinates: [ 0.500000000000000, 0.250000000000000, 0.000000000000000 ] + mass: 35.453000 + reduced_to: 33 + - symbol: Cl # 51 + coordinates: [ 0.000000000000000, 0.750000000000000, 0.000000000000000 ] + mass: 35.453000 + reduced_to: 33 + - symbol: Cl # 52 + coordinates: [ 0.500000000000000, 0.750000000000000, 0.000000000000000 ] + mass: 35.453000 + reduced_to: 33 + - symbol: Cl # 53 + coordinates: [ 0.000000000000000, 0.250000000000000, 0.500000000000000 ] + mass: 35.453000 + reduced_to: 33 + - symbol: Cl # 54 + coordinates: [ 0.500000000000000, 0.250000000000000, 0.500000000000000 ] + mass: 35.453000 + reduced_to: 33 + - symbol: Cl # 55 + coordinates: [ 0.000000000000000, 0.750000000000000, 0.500000000000000 ] + mass: 35.453000 + reduced_to: 33 + - symbol: Cl # 56 + coordinates: [ 0.500000000000000, 0.750000000000000, 0.500000000000000 ] + mass: 35.453000 + reduced_to: 33 + - symbol: Cl # 57 + coordinates: [ 0.000000000000000, 0.000000000000000, 0.250000000000000 ] + mass: 35.453000 + reduced_to: 33 + - symbol: Cl # 58 + coordinates: [ 0.500000000000000, 0.000000000000000, 0.250000000000000 ] + mass: 35.453000 + reduced_to: 33 + - symbol: Cl # 59 + coordinates: [ 0.000000000000000, 0.500000000000000, 0.250000000000000 ] + mass: 35.453000 + reduced_to: 33 + - symbol: Cl # 60 + coordinates: [ 0.500000000000000, 0.500000000000000, 0.250000000000000 ] + mass: 35.453000 + reduced_to: 33 + - symbol: Cl # 61 + coordinates: [ 0.000000000000000, 0.000000000000000, 0.750000000000000 ] + mass: 35.453000 + reduced_to: 33 + - symbol: Cl # 62 + coordinates: [ 0.500000000000000, 0.000000000000000, 0.750000000000000 ] + mass: 35.453000 + reduced_to: 33 + - symbol: Cl # 63 + coordinates: [ 0.000000000000000, 0.500000000000000, 0.750000000000000 ] + mass: 35.453000 + reduced_to: 33 + - symbol: Cl # 64 + coordinates: [ 0.500000000000000, 0.500000000000000, 0.750000000000000 ] + mass: 35.453000 + reduced_to: 33 + +displacements: +- atom: 1 + displacement: + [ 0.0100000000000000, 0.0000000000000000, 0.0000000000000000 ] +- atom: 33 + displacement: + [ 0.0100000000000000, 0.0000000000000000, 0.0000000000000000 ]