scnc
/
phono3py
mirror of https://github.com/phonopy/phono3py.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #372 from phonopy/physical-unit 

Migration of the physical units system
This commit is contained in:
Atsushi Togo 2025-04-29 19:56:36 +09:00 committed by GitHub
parent c171367edd 7c0b4a2d0e
commit 6264dc5e35
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
29 changed files with 171 additions and 100 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
c/phonoc_const.h
View File

@ -36,6 +36,5 @@
#define __phonoc_const_H__
#define M_2PI 6.283185307179586
#define KB 8.6173382568083159E-05
#endif

3
phono3py/api_isotope.py
View File

@ -35,7 +35,6 @@
# POSSIBILITY OF SUCH DAMAGE.
import numpy as np
from phonopy.units import VaspToTHz
from phono3py.other.isotope import Isotope
@ -50,7 +49,7 @@ class Phono3pyIsotope:
mass_variances=None, # length of list is num_atom.
band_indices=None,
sigmas=None,
frequency_factor_to_THz=VaspToTHz,
frequency_factor_to_THz=None,
use_grg=False,
symprec=1e-5,
cutoff_frequency=None,

9
phono3py/api_jointdos.py
View File

@ -35,9 +35,9 @@
# POSSIBILITY OF SUCH DAMAGE.
import numpy as np
from phonopy.physical_units import get_physical_units
from phonopy.structure.cells import Primitive, Supercell
from phonopy.structure.symmetry import Symmetry
from phonopy.units import VaspToTHz
from phono3py.file_IO import write_joint_dos
from phono3py.phonon.grid import BZGrid
@ -65,7 +65,7 @@ class Phono3pyJointDos:
num_frequency_points=None,
num_points_in_batch=None,
temperatures=None,
frequency_factor_to_THz=VaspToTHz,
frequency_factor_to_THz=None,
frequency_scale_factor=None,
use_grg=False,
SNF_coordinates="reciprocal",
@ -87,7 +87,10 @@ class Phono3pyJointDos:
else:
self._sigmas = sigmas
self._cutoff_frequency = cutoff_frequency
self._frequency_factor_to_THz = frequency_factor_to_THz
if frequency_factor_to_THz is None:
self._frequency_factor_to_THz = get_physical_units().defaultToTHz
else:
self._frequency_factor_to_THz = frequency_factor_to_THz
self._frequency_scale_factor = frequency_scale_factor
self._is_mesh_symmetry = is_mesh_symmetry
self._is_symmetry = is_symmetry

9
phono3py/api_phono3py.py
View File

@ -58,6 +58,7 @@ from phonopy.interface.mlp import PhonopyMLP
from phonopy.interface.pypolymlp import (
PypolymlpParams,
)
from phonopy.physical_units import get_physical_units
from phonopy.structure.atoms import PhonopyAtoms
from phonopy.structure.cells import (
Primitive,
@ -69,7 +70,6 @@ from phonopy.structure.cells import (
shape_supercell_matrix,
)
from phonopy.structure.symmetry import Symmetry
from phonopy.units import VaspToTHz
from phono3py.conductivity.init_direct_solution import get_thermal_conductivity_LBTE
from phono3py.conductivity.init_rta import get_thermal_conductivity_RTA
@ -149,7 +149,7 @@ class Phono3py:
primitive_matrix=None,
phonon_supercell_matrix=None,
cutoff_frequency=1e-4,
frequency_factor_to_THz=VaspToTHz,
frequency_factor_to_THz=None,
is_symmetry=True,
is_mesh_symmetry=True,
use_grg=False,
@ -218,7 +218,10 @@ class Phono3py:
"""
self._symprec = symprec
self._frequency_factor_to_THz = frequency_factor_to_THz
if frequency_factor_to_THz is None:
self._frequency_factor_to_THz = get_physical_units().defaultToTHz
else:
self._frequency_factor_to_THz = frequency_factor_to_THz
self._is_symmetry = is_symmetry
self._is_mesh_symmetry = is_mesh_symmetry
self._use_grg = use_grg

29
phono3py/conductivity/base.py
View File

@ -42,7 +42,7 @@ from typing import List, Optional, Tuple, Union
import numpy as np
from phonopy.phonon.group_velocity import GroupVelocity
from phonopy.phonon.thermal_properties import mode_cv
from phonopy.units import EV, Angstrom, Kb, THz, THzToEv
from phonopy.physical_units import get_physical_units
from phono3py.other.isotope import Isotope
from phono3py.phonon.grid import get_grid_points_by_rotations, get_ir_grid_points
@ -50,9 +50,17 @@ from phono3py.phonon3.collision_matrix import CollisionMatrix
from phono3py.phonon3.imag_self_energy import ImagSelfEnergy
from phono3py.phonon3.interaction import Interaction
unit_to_WmK = (
(THz * Angstrom) ** 2 / (Angstrom**3) * EV / THz / (2 * np.pi)
) # 2pi comes from definition of lifetime.
def get_unit_to_WmK() -> float:
"""Return conversion factor to WmK."""
unit_to_WmK = (
(get_physical_units().THz * get_physical_units().Angstrom) ** 2
/ (get_physical_units().Angstrom ** 3)
* get_physical_units().EV
/ get_physical_units().THz
/ (2 * np.pi)
) # 2pi comes from definition of lifetime.
return unit_to_WmK
class HeatCapacityMixIn:
@ -96,15 +104,18 @@ class HeatCapacityMixIn:
"""
grid_point = self._grid_points[i_gp]
freqs = self._frequencies[grid_point][self._pp.band_indices] * THzToEv
cutoff = self._pp.cutoff_frequency * THzToEv
freqs = (
self._frequencies[grid_point][self._pp.band_indices]
* get_physical_units().THzToEv
)
cutoff = self._pp.cutoff_frequency * get_physical_units().THzToEv
cv = np.zeros((len(self._temperatures), len(freqs)), dtype="double")
# x=freq/T has to be small enough to avoid overflow of exp(x).
# x < 100 is the hard-corded criterion.
# Otherwise just set 0.
for i, f in enumerate(freqs):
if f > cutoff:
condition = f < 100 * self._temperatures * Kb
condition = f < 100 * self._temperatures * get_physical_units().Kb
cv[:, i] = np.where(
condition,
mode_cv(np.where(condition, self._temperatures, 10000), f),
@ -379,7 +390,7 @@ class ConductivityBase(ABC):
self._gamma_iso: Optional[np.ndarray] = None
volume = self._pp.primitive.volume
self._conversion_factor = unit_to_WmK / volume
self._conversion_factor = get_unit_to_WmK() / volume
self._averaged_pp_interaction = None
@ -903,7 +914,7 @@ class ConductivityBase(ABC):
for ll in range(num_band):
g_boundary[ll] = (
np.linalg.norm(gv[i_gp, ll])
* Angstrom
* get_physical_units().Angstrom
* 1e6
/ (4 * np.pi * self._boundary_mfp)
)

21
phono3py/conductivity/direct_solution_base.py
View File

@ -42,7 +42,7 @@ from typing import Optional
import numpy as np
from phonopy.phonon.degeneracy import degenerate_sets
from phonopy.units import Kb, THzToEv
from phonopy.physical_units import get_physical_units
from phono3py.conductivity.base import ConductivityBase
from phono3py.conductivity.utils import select_colmat_solver
@ -791,11 +791,18 @@ class ConductivityLBTEBase(ConductivityBase):
t = self._temperatures[i_temp]
sinh = np.where(
freqs > self._pp.cutoff_frequency,
np.sinh(freqs * THzToEv / (2 * Kb * t)),
np.sinh(
freqs * get_physical_units().THzToEv / (2 * get_physical_units().Kb * t)
),
-1.0,
)
inv_sinh = np.where(sinh > 0, 1.0 / sinh, 0)
freqs_sinh = freqs * THzToEv * inv_sinh / (4 * Kb * t**2)
freqs_sinh = (
freqs
* get_physical_units().THzToEv
* inv_sinh
/ (4 * get_physical_units().Kb * t**2)
)
for i, f in enumerate(freqs_sinh):
X[i] *= weights[i]
@ -1080,7 +1087,9 @@ class ConductivityLBTEBase(ConductivityBase):
mode_kappa[i_sigma, i_temp, i, j, k] = sum_k[vxf]
t = self._temperatures[i_temp]
mode_kappa[i_sigma, i_temp] *= self._conversion_factor * Kb * t**2
mode_kappa[i_sigma, i_temp] *= (
self._conversion_factor * get_physical_units().Kb * t**2
)
def _set_mode_kappa_Chaput(self, mode_kappa, i_sigma, i_temp, weights):
"""Calculate mode kappa by the way in Laurent Chaput's PRL paper.
@ -1124,7 +1133,7 @@ class ConductivityLBTEBase(ConductivityBase):
vals = vals.reshape(num_ir_grid_points, num_band)
mode_kappa[i_sigma, i_temp, :, :, i] += vals
factor = self._conversion_factor * Kb * t**2
factor = self._conversion_factor * get_physical_units().Kb * t**2
mode_kappa[i_sigma, i_temp] *= factor
def _set_mode_kappa_from_mfp(self, weights, rotations_cartesian, i_sigma, i_temp):
@ -1151,7 +1160,7 @@ class ConductivityLBTEBase(ConductivityBase):
if cv < 1e-10:
continue
self._mfp[i_sigma, i_temp, i, j] = (
-2 * t * np.sqrt(Kb / cv) * f / (2 * np.pi)
-2 * t * np.sqrt(get_physical_units().Kb / cv) * f / (2 * np.pi)
)
def _show_log(self, i):

4
phono3py/conductivity/init_direct_solution.py
View File

@ -37,7 +37,7 @@
import sys
from typing import Optional, Union
from phono3py.conductivity.base import unit_to_WmK
from phono3py.conductivity.base import get_unit_to_WmK
from phono3py.conductivity.direct_solution import ConductivityLBTE
from phono3py.conductivity.direct_solution_base import ConductivityLBTEBase
from phono3py.conductivity.utils import (
@ -438,7 +438,7 @@ class ConductivityLBTEWriter:
weight=weights,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
kappa_unit_conversion=unit_to_WmK / volume,
kappa_unit_conversion=get_unit_to_WmK() / volume,
compression=compression,
filename=filename,
verbose=log_level,

8
phono3py/conductivity/init_rta.py
View File

@ -38,7 +38,7 @@ from typing import Optional, Union, cast
import numpy as np
from phono3py.conductivity.base import unit_to_WmK
from phono3py.conductivity.base import get_unit_to_WmK
from phono3py.conductivity.kubo_rta import ConductivityKuboRTA
from phono3py.conductivity.rta import ConductivityRTA
from phono3py.conductivity.rta_base import ConductivityRTABase
@ -362,7 +362,7 @@ class ConductivityRTAWriter:
grid_point=gp,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
kappa_unit_conversion=unit_to_WmK / volume,
kappa_unit_conversion=get_unit_to_WmK() / volume,
compression=compression,
filename=filename,
verbose=verbose,
@ -417,7 +417,7 @@ class ConductivityRTAWriter:
band_index=bi,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
kappa_unit_conversion=unit_to_WmK / volume,
kappa_unit_conversion=get_unit_to_WmK() / volume,
compression=compression,
filename=filename,
verbose=verbose,
@ -544,7 +544,7 @@ class ConductivityRTAWriter:
weight=weights,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
kappa_unit_conversion=unit_to_WmK / volume,
kappa_unit_conversion=get_unit_to_WmK() / volume,
compression=compression,
filename=filename,
verbose=log_level,

8
phono3py/conductivity/kubo_base.py
View File

@ -35,7 +35,7 @@
# POSSIBILITY OF SUCH DAMAGE.
import numpy as np
from phonopy.units import Kb, THzToEv
from phonopy.physical_units import get_physical_units
from phono3py.phonon.group_velocity_matrix import GroupVelocityMatrix
from phono3py.phonon.heat_capacity_matrix import mode_cv_matrix
@ -71,11 +71,11 @@ class ConductivityKuboMixIn:
"""
irgp = self._grid_points[i_gp]
freqs = self._frequencies[irgp] * THzToEv
cutoff = self._pp.cutoff_frequency * THzToEv
freqs = self._frequencies[irgp] * get_physical_units().THzToEv
cutoff = self._pp.cutoff_frequency * get_physical_units().THzToEv
for i_temp, temp in enumerate(self._temperatures):
if (freqs / (temp * Kb) > 100).any():
if (freqs / (temp * get_physical_units().Kb) > 100).any():
continue
cvm = mode_cv_matrix(temp, freqs, cutoff=cutoff)
self._cv_mat[i_temp, i_data] = cvm[self._pp.band_indices, :]

13
phono3py/conductivity/wigner_base.py
View File

@ -38,7 +38,7 @@ import textwrap
import numpy as np
from phonopy.phonon.degeneracy import degenerate_sets
from phonopy.units import EV, Angstrom, Hbar, THz
from phonopy.physical_units import get_physical_units
from phono3py.conductivity.base import HeatCapacityMixIn
from phono3py.phonon.grid import get_grid_points_by_rotations
@ -215,8 +215,9 @@ class ConductivityWignerMixIn(HeatCapacityMixIn):
def get_conversion_factor_WTE(volume):
"""Return conversion factor of thermal conductivity."""
return (
(THz * Angstrom) ** 2 # ----> group velocity
* EV # ----> specific heat is in eV/
* Hbar # ----> transform lorentzian_div_hbar from eV^-1 to s
/ (volume * Angstrom**3)
) # ----> unit cell volume
(get_physical_units().THz * get_physical_units().Angstrom)
** 2 # --> group velocity
* get_physical_units().EV # --> specific heat is in eV/
* get_physical_units().Hbar # --> transform lorentzian_div_hbar from eV^-1 to s
/ (volume * get_physical_units().Angstrom ** 3)
) # --> unit cell volume

3
phono3py/conductivity/wigner_direct_solution.py
View File

@ -35,7 +35,7 @@
# POSSIBILITY OF SUCH DAMAGE.
import numpy as np
from phonopy.units import THzToEv
from phonopy.physical_units import get_physical_units
from phono3py.conductivity.direct_solution import (
ConductivityLBTEBase,
@ -282,6 +282,7 @@ class ConductivityWignerLBTE(ConductivityWignerMixIn, ConductivityLBTEBase):
N = self._num_sampling_grid_points
num_band = len(self._pp.primitive) * 3
# num_ir_grid_points = len(self._ir_grid_points)
THzToEv = get_physical_units().THzToEv
for i, gp in enumerate(self._ir_grid_points):
# linewidths at qpoint i, sigma i_sigma, and temperature i_temp
g = self._get_main_diagonal(i, i_sigma, i_temp) * 2.0 # linewidth (FWHM)

3
phono3py/conductivity/wigner_rta.py
View File

@ -35,7 +35,7 @@
# POSSIBILITY OF SUCH DAMAGE.
import numpy as np
from phonopy.units import THzToEv
from phonopy.physical_units import get_physical_units
from phono3py.conductivity.rta_base import ConductivityRTABase
from phono3py.conductivity.wigner_base import (
@ -119,6 +119,7 @@ class ConductivityWignerRTA(ConductivityWignerMixIn, ConductivityRTABase):
"""
num_band = len(self._pp.primitive) * 3
THzToEv = get_physical_units().THzToEv
for i, _ in enumerate(self._grid_points):
cv = self._cv[:, i, :]
gp = self._grid_points[i]

4
phono3py/cui/create_force_constants.py
View File

@ -50,7 +50,7 @@ from phonopy.harmonic.force_constants import (
symmetrize_compact_force_constants,
symmetrize_force_constants,
)
from phonopy.interface.calculator import get_default_physical_units
from phonopy.interface.calculator import get_calculator_physical_units
from phonopy.interface.fc_calculator import fc_calculator_names
from phonopy.interface.pypolymlp import PypolymlpParams, parse_mlp_params
from phonopy.interface.symfc import parse_symfc_options
@ -254,7 +254,7 @@ def parse_forces(
if dataset:
filename_read_from = phono3py_yaml_filename
physical_units = get_default_physical_units(calculator)
physical_units = get_calculator_physical_units(calculator)
# Forces are not yet found in dataset. Then try to read from FORCES_FC3 or
# FORCES_FC2.

8
phono3py/cui/load.py
View File

@ -43,10 +43,10 @@ from typing import Optional, Union
import numpy as np
import phonopy.cui.load_helper as load_helper
from phonopy.harmonic.force_constants import show_drift_force_constants
from phonopy.interface.calculator import get_default_physical_units
from phonopy.interface.calculator import get_calculator_physical_units
from phonopy.physical_units import get_physical_units
from phonopy.structure.atoms import PhonopyAtoms
from phonopy.structure.cells import determinant
from phonopy.units import VaspToTHz
from phono3py import Phono3py
from phono3py.cui.create_force_constants import (
@ -310,12 +310,12 @@ def load(
_nac_params = None
# Convert distance unit of unit cell to Angstrom
physical_units = get_default_physical_units(_calculator)
physical_units = get_calculator_physical_units(_calculator)
factor_to_A = physical_units["distance_to_A"]
cell.cell = cell.cell * factor_to_A
if factor is None:
_factor = VaspToTHz
_factor = get_physical_units().defaultToTHz
else:
_factor = factor
ph3py = Phono3py(

16
phono3py/cui/phono3py_script.py
View File

@ -58,10 +58,10 @@ from phonopy.cui.phonopy_script import (
from phonopy.exception import ForceCalculatorRequiredError
from phonopy.file_IO import is_file_phonopy_yaml
from phonopy.harmonic.force_constants import show_drift_force_constants
from phonopy.interface.calculator import get_default_physical_units
from phonopy.interface.calculator import get_calculator_physical_units
from phonopy.phonon.band_structure import get_band_qpoints
from phonopy.physical_units import get_physical_units
from phonopy.structure.cells import isclose as cells_isclose
from phonopy.units import Bohr, Hartree, VaspToTHz
from phono3py import Phono3py, Phono3pyIsotope, Phono3pyJointDos
from phono3py.cui.create_force_constants import (
@ -159,7 +159,7 @@ def finalize_phono3py(
else:
yaml_filename = filename
_physical_units = get_default_physical_units(phono3py.calculator)
_physical_units = get_calculator_physical_units(phono3py.calculator)
ph3py_yaml = Phono3pyYaml(
configuration=confs_dict,
@ -390,7 +390,7 @@ def get_default_values(settings):
temperatures = settings.temperatures # For others
if settings.frequency_conversion_factor is None:
frequency_factor_to_THz = VaspToTHz
frequency_factor_to_THz = get_physical_units().defaultToTHz
else:
frequency_factor_to_THz = settings.frequency_conversion_factor
@ -475,7 +475,7 @@ def init_phono3py(
settings, cell_info, interface_mode, symprec, log_level
) -> tuple[Phono3py, dict]:
"""Initialize phono3py and update settings by default values."""
physical_units = get_default_physical_units(interface_mode)
physical_units = get_calculator_physical_units(interface_mode)
distance_to_A = physical_units["distance_to_A"]
# Change unit of lattice parameters to angstrom
@ -566,7 +566,7 @@ def create_supercells_with_displacements(
cell_info["phonopy_yaml"],
unitcell_filename,
log_level,
nac_factor=Hartree * Bohr,
nac_factor=get_physical_units().Hartree * get_physical_units().Bohr,
load_phonopy_yaml=load_phono3py_yaml,
)
@ -730,7 +730,7 @@ def run_gruneisen_then_exit(phono3py, settings, output_filename, log_level):
nac_params=phono3py.nac_params,
nac_q_direction=settings.nac_q_direction,
ion_clamped=settings.ion_clamped,
factor=VaspToTHz,
factor=get_physical_units().defaultToTHz,
symprec=phono3py.symmetry.tolerance,
output_filename=output_filename,
log_level=log_level,
@ -1149,7 +1149,7 @@ def main(**argparse_control):
cell_info["phonopy_yaml"],
unitcell_filename,
log_level,
nac_factor=Hartree * Bohr,
nac_factor=get_physical_units().Hartree * get_physical_units().Bohr,
load_phonopy_yaml=load_phono3py_yaml,
)

9
phono3py/other/isotope.py
View File

@ -41,12 +41,12 @@ from typing import Optional, Union
import numpy as np
from phonopy.harmonic.dynamical_matrix import get_dynamical_matrix
from phonopy.phonon.tetrahedron_mesh import get_tetrahedra_frequencies
from phonopy.physical_units import get_physical_units
from phonopy.structure.atoms import PhonopyAtoms
from phonopy.structure.atoms import isotope_data as phonopy_isotope_data
from phonopy.structure.cells import Primitive
from phonopy.structure.symmetry import Symmetry
from phonopy.structure.tetrahedron_method import TetrahedronMethod
from phonopy.units import VaspToTHz
from phono3py.other.tetrahedron_method import (
get_integration_weights,
@ -100,7 +100,7 @@ class Isotope:
band_indices=None,
sigma=None,
bz_grid=None,
frequency_factor_to_THz=VaspToTHz,
frequency_factor_to_THz=None,
use_grg=False,
symprec=1e-5,
cutoff_frequency=None,
@ -122,7 +122,10 @@ class Isotope:
self._cutoff_frequency = 0
else:
self._cutoff_frequency = cutoff_frequency
self._frequency_factor_to_THz = frequency_factor_to_THz
if frequency_factor_to_THz is None:
self._frequency_factor_to_THz = get_physical_units().defaultToTHz
else:
self._frequency_factor_to_THz = frequency_factor_to_THz
self._lapack_zheev_uplo = lapack_zheev_uplo
self._nac_q_direction = None

14
phono3py/phonon/func.py
View File

@ -35,7 +35,7 @@
# POSSIBILITY OF SUCH DAMAGE.
import numpy as np
from phonopy.units import AMU, EV, Angstrom, Hbar, Kb, THz, THzToEv
from phonopy.physical_units import get_physical_units
def gaussian(x, sigma):
@ -62,7 +62,9 @@ def bose_einstein(x, T):
Temperature in K
"""
return 1.0 / (np.exp(THzToEv * x / (Kb * T)) - 1)
return 1.0 / (
np.exp(get_physical_units().THzToEv * x / (get_physical_units().Kb * T)) - 1
)
def sigma_squared(x, T):
@ -96,7 +98,13 @@ def sigma_squared(x, T):
n = bose_einstein(x, T)
# factor=1.0107576777968994
factor = Hbar * EV / (2 * np.pi * THz) / AMU / Angstrom**2
factor = (
get_physical_units().Hbar
* get_physical_units().EV
/ (2 * np.pi * get_physical_units().THz)
/ get_physical_units().AMU
/ get_physical_units().Angstrom ** 2
)
#########################
np.seterr(**old_settings)

3
phono3py/phonon/group_velocity_matrix.py
View File

@ -37,7 +37,6 @@
import numpy as np
from phonopy.phonon.degeneracy import degenerate_sets
from phonopy.phonon.group_velocity import GroupVelocity
from phonopy.units import VaspToTHz
from phonopy.utils import similarity_transformation
@ -57,7 +56,7 @@ class GroupVelocityMatrix(GroupVelocity):
dynamical_matrix,
q_length=None,
symmetry=None,
frequency_factor_to_THz=VaspToTHz,
frequency_factor_to_THz=None,
cutoff_frequency=1e-4,
):
"""Init method.

3
phono3py/phonon/heat_capacity_matrix.py
View File

@ -36,7 +36,7 @@
import numpy as np
from phonopy.phonon.thermal_properties import mode_cv
from phonopy.units import Kb
from phonopy.physical_units import get_physical_units
def mode_cv_matrix(temp, freqs, cutoff=1e-4):
@ -67,6 +67,7 @@ def mode_cv_matrix(temp, freqs, cutoff=1e-4):
shape=(num_band, num_band), dtype='double', order='C'.
"""
Kb = get_physical_units().Kb
x = freqs / Kb / temp
shape = (len(freqs), len(freqs))
cvm = np.zeros(shape, dtype="double", order="C")

13
phono3py/phonon/solver.py
View File

@ -35,8 +35,8 @@
# POSSIBILITY OF SUCH DAMAGE.
import numpy as np
from phonopy.physical_units import get_physical_units
from phonopy.structure.cells import sparse_to_dense_svecs
from phonopy.units import VaspToTHz
def run_phonon_solver_c(
@ -47,7 +47,7 @@ def run_phonon_solver_c(
grid_points,
grid_address,
QDinv,
frequency_conversion_factor=VaspToTHz,
frequency_conversion_factor=None,
nac_q_direction=None, # in reduced coordinates
lapack_zheev_uplo="L",
):
@ -84,6 +84,11 @@ def run_phonon_solver_c(
import phono3py._phono3py as phono3c
import phono3py._phononcalc as phononcalc
if frequency_conversion_factor is None:
_frequency_conversion_factor = get_physical_units().defaultToTHz
else:
_frequency_conversion_factor = frequency_conversion_factor
(
svecs,
multi,
@ -154,7 +159,7 @@ def run_phonon_solver_c(
masses,
fc_p2s,
fc_s2p,
frequency_conversion_factor,
_frequency_conversion_factor,
born,
dielectric,
rec_lattice,
@ -179,7 +184,7 @@ def run_phonon_solver_c(
frequencies[phonon_undone] = (
np.sign(frequencies[phonon_undone])
* np.sqrt(np.abs(frequencies[phonon_undone]))
* frequency_conversion_factor
* _frequency_conversion_factor
)

9
phono3py/phonon/velocity_operator.py
View File

@ -36,7 +36,7 @@
import numpy as np
from phonopy.phonon.group_velocity import GroupVelocity
from phonopy.units import VaspToTHz
from phonopy.physical_units import get_physical_units
class VelocityOperator(GroupVelocity):
@ -47,7 +47,7 @@ class VelocityOperator(GroupVelocity):
dynamical_matrix,
q_length=None,
symmetry=None,
frequency_factor_to_THz=VaspToTHz,
frequency_factor_to_THz=None,
cutoff_frequency=1e-4,
):
"""Init method.
@ -75,7 +75,10 @@ class VelocityOperator(GroupVelocity):
if self._q_length is None:
self._q_length = 5e-6
self._symmetry = symmetry
self._factor = frequency_factor_to_THz
if frequency_factor_to_THz is None:
self._factor = get_physical_units().defaultToTHz
else:
self._factor = frequency_factor_to_THz
self._cutoff_frequency = cutoff_frequency
self._directions = np.array(

8
phono3py/phonon3/collision_matrix.py
View File

@ -37,7 +37,7 @@
from typing import Optional
import numpy as np
from phonopy.units import Kb, THzToEv
from phonopy.physical_units import get_physical_units
from phono3py.phonon3.imag_self_energy import ImagSelfEnergy
from phono3py.phonon3.interaction import Interaction
@ -311,7 +311,11 @@ class CollisionMatrix(ImagSelfEnergy):
freqs = self._frequencies[gp]
sinh = np.where(
freqs > self._cutoff_frequency,
np.sinh(freqs * THzToEv / (2 * Kb * self._temperature)),
np.sinh(
freqs
* get_physical_units().THzToEv
/ (2 * get_physical_units().Kb * self._temperature)
),
-1.0,
)
inv_sinh = np.where(sinh > 0, 1.0 / sinh, 0)

9
phono3py/phonon3/gruneisen.py
View File

@ -38,10 +38,10 @@ import sys
import numpy as np
from phonopy.harmonic.dynamical_matrix import get_dynamical_matrix
from phonopy.physical_units import get_physical_units
from phonopy.structure.atoms import PhonopyAtoms
from phonopy.structure.cells import Primitive
from phonopy.structure.grid_points import get_qpoints
from phonopy.units import VaspToTHz
def run_gruneisen_parameters(
@ -133,7 +133,7 @@ class Gruneisen:
nac_params=None,
nac_q_direction=None,
ion_clamped=False,
factor=VaspToTHz,
factor=None,
symprec=1e-5,
):
"""Init method."""
@ -142,7 +142,10 @@ class Gruneisen:
self._scell = supercell
self._pcell = primitive
self._ion_clamped = ion_clamped
self._factor = factor
if factor is None:
self._factor = get_physical_units().defaultToTHz
else:
self._factor = factor
self._symprec = symprec
self._dm = get_dynamical_matrix(
self._fc2,

8
phono3py/phonon3/imag_self_energy.py
View File

@ -40,7 +40,7 @@ from typing import List, Optional
import numpy as np
from phonopy.phonon.degeneracy import degenerate_sets
from phonopy.units import EV, Hbar, THz
from phonopy.physical_units import get_physical_units
from phono3py.file_IO import (
write_gamma_detail_to_hdf5,
@ -97,7 +97,11 @@ class ImagSelfEnergy:
# Unit to THz of Gamma
self._unit_conversion = (
18 * np.pi / (Hbar * EV) ** 2 / (2 * np.pi * THz) ** 2 * EV**2
18
* np.pi
/ (get_physical_units().Hbar * get_physical_units().EV) ** 2
/ (2 * np.pi * get_physical_units().THz) ** 2
* get_physical_units().EV ** 2
)
def run(self):

21
phono3py/phonon3/interaction.py
View File

@ -41,9 +41,9 @@ from typing import Literal, Optional, Union
import numpy as np
from phonopy.harmonic.dynamical_matrix import DynamicalMatrix, get_dynamical_matrix
from phonopy.physical_units import get_physical_units
from phonopy.structure.cells import Primitive, compute_all_sg_permutations
from phonopy.structure.symmetry import Symmetry
from phonopy.units import AMU, EV, Angstrom, Hbar, THz, VaspToTHz
from phono3py.phonon.grid import (
BZGrid,
@ -97,7 +97,7 @@ class Interaction:
fc3: Optional[np.ndarray] = None,
band_indices: Optional[Union[np.ndarray, Sequence]] = None,
constant_averaged_interaction: Optional[float] = None,
frequency_factor_to_THz: float = VaspToTHz,
frequency_factor_to_THz: Optional[float] = None,
frequency_scale_factor: Optional[float] = None,
unit_conversion: Optional[float] = None,
is_mesh_symmetry: bool = True,
@ -115,7 +115,10 @@ class Interaction:
self._band_indices = None
self._set_band_indices(band_indices)
self._constant_averaged_interaction = constant_averaged_interaction
self._frequency_factor_to_THz = frequency_factor_to_THz
if frequency_factor_to_THz is None:
self._frequency_factor_to_THz = get_physical_units().defaultToTHz
else:
self._frequency_factor_to_THz = frequency_factor_to_THz
self._frequency_scale_factor = frequency_scale_factor
if fc3 is not None:
@ -125,15 +128,15 @@ class Interaction:
if unit_conversion is None:
num_grid = np.prod(self.mesh_numbers)
self._unit_conversion = (
(Hbar * EV) ** 3
(get_physical_units().Hbar * get_physical_units().EV) ** 3
/ 36
/ 8
* EV**2
/ Angstrom**6
/ (2 * np.pi * THz) ** 3
/ AMU**3
* get_physical_units().EV ** 2
/ get_physical_units().Angstrom ** 6
/ (2 * np.pi * get_physical_units().THz) ** 3
/ get_physical_units().AMU ** 3
/ num_grid
/ EV**2
/ get_physical_units().EV ** 2
)
else:
self._unit_conversion = unit_conversion

9
phono3py/phonon3/joint_dos.py
View File

@ -39,8 +39,8 @@ import warnings
import numpy as np
from phonopy.harmonic.dynamical_matrix import DynamicalMatrix, get_dynamical_matrix
from phonopy.physical_units import get_physical_units
from phonopy.structure.cells import Primitive, Supercell
from phonopy.units import VaspToTHz
from phono3py.phonon.func import bose_einstein
from phono3py.phonon.grid import BZGrid, get_grid_point_from_address
@ -66,7 +66,7 @@ class JointDos:
sigma=None,
sigma_cutoff=None,
cutoff_frequency=None,
frequency_factor_to_THz=VaspToTHz,
frequency_factor_to_THz=None,
frequency_scale_factor=1.0,
is_mesh_symmetry=True,
symprec=1e-5,
@ -90,7 +90,10 @@ class JointDos:
self._cutoff_frequency = 0
else:
self._cutoff_frequency = cutoff_frequency
self._frequency_factor_to_THz = frequency_factor_to_THz
if frequency_factor_to_THz is None:
self._frequency_factor_to_THz = get_physical_units().defaultToTHz
else:
self._frequency_factor_to_THz = frequency_factor_to_THz
self._frequency_scale_factor = frequency_scale_factor
self._is_mesh_symmetry = is_mesh_symmetry
self._symprec = symprec

9
phono3py/phonon3/real_self_energy.py
View File

@ -38,7 +38,7 @@ import sys
import numpy as np
from phonopy.phonon.degeneracy import degenerate_sets
from phonopy.units import EV, Hbar, THz
from phonopy.physical_units import get_physical_units
from phono3py.file_IO import (
write_real_self_energy_at_grid_point,
@ -115,7 +115,12 @@ class RealSelfEnergy:
self._real_self_energies = None
# Unit to THz of Delta
self._unit_conversion = 18 / (Hbar * EV) ** 2 / (2 * np.pi * THz) ** 2 * EV**2
self._unit_conversion = (
18
/ (get_physical_units().Hbar * get_physical_units().EV) ** 2
/ (2 * np.pi * get_physical_units().THz) ** 2
* get_physical_units().EV ** 2
)
def run(self):
"""Calculate real-part of self-energies."""

11
phono3py/sscha/sscha.py
View File

@ -43,7 +43,7 @@ Formulae implemented are based on these papers:
import numpy as np
from phonopy.harmonic.dynmat_to_fc import DynmatToForceConstants
from phonopy.units import VaspToTHz
from phonopy.physical_units import get_physical_units
from phono3py.phonon.func import sigma_squared
@ -97,7 +97,7 @@ class SupercellPhonon:
"""
def __init__(self, supercell, force_constants, frequency_factor_to_THz=VaspToTHz):
def __init__(self, supercell, force_constants, frequency_factor_to_THz=None):
"""Init method.
Parameters
@ -107,7 +107,7 @@ class SupercellPhonon:
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_factor_to_THz : float, optional
Frequency conversion factor to THz.
"""
@ -124,7 +124,10 @@ class SupercellPhonon:
)
eigvals, eigvecs = np.linalg.eigh(dynmat)
freqs = np.sqrt(np.abs(eigvals)) * np.sign(eigvals)
freqs *= frequency_factor_to_THz
if frequency_factor_to_THz is None:
freqs *= get_physical_units().defaultToTHz
else:
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")

6
test/phonon/test_velocity_operator.py
View File

@ -2,7 +2,7 @@
import numpy as np
from phonopy import Phonopy
from phonopy.units import THzToCm, VaspToTHz
from phonopy.physical_units import get_physical_units
from phono3py.phonon.velocity_operator import VelocityOperator
@ -194,7 +194,7 @@ def test_gv_operator_nacl(ph_nacl: Phonopy):
eigvals, eigvecs = np.linalg.eigh(dm)
np.testing.assert_allclose(
eigvals * VaspToTHz * THzToCm,
eigvals * get_physical_units().defaultToTHz * get_physical_units().THzToCm,
eigvals_NaCl_Ref,
atol=0.00001,
rtol=0.00001,
@ -406,7 +406,7 @@ def test_gv_operator_si(ph_si: Phonopy):
eigvals, eigvecs = np.linalg.eigh(dm)
np.testing.assert_allclose(
eigvals * VaspToTHz * THzToCm,
eigvals * get_physical_units().defaultToTHz * get_physical_units().THzToCm,
eigvals_si_Ref,
atol=0.00001,
rtol=0.00001,