Merge pull request #191 from phonopy/fix-unit-conversion

Fix unit conversion

.pre-commit-config.yaml

@@ -19,7 +19,7 @@ repos:
     - "--ignore=E203,W503"
 
 - repo: https://github.com/psf/black
-  rev: 23.11.0
+  rev: 23.10.1
   hooks:
   - id: black
     args:

phono3py/api_phono3py.py

@@ -34,6 +34,7 @@
 # POSSIBILITY OF SUCH DAMAGE.
 
 import warnings
+from typing import Optional
 
 import numpy as np
 from phonopy.exception import ForceCalculatorRequiredError
@@ -155,7 +156,7 @@ def __init__(
         store_dense_gp_map=True,
         store_dense_svecs=True,
         symprec=1e-5,
-        calculator=None,
+        calculator: Optional[str] = None,
         log_level=0,
         lapack_zheev_uplo=None,
     ):
@@ -254,7 +255,7 @@ def __init__(
         self._store_dense_svecs = store_dense_svecs
 
         self._cutoff_frequency = cutoff_frequency
-        self._calculator = calculator
+        self._calculator: Optional[str] = calculator
         self._log_level = log_level
 
         # Create supercell and primitive cell
@@ -399,7 +400,7 @@ def get_version(self):
         return self.version
 
     @property
-    def calculator(self):
+    def calculator(self) -> Optional[str]:
         """Return calculator interface name.
 
         str

phono3py/cui/create_force_constants.py

@@ -33,6 +33,7 @@
 # ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 # POSSIBILITY OF SUCH DAMAGE.
 
+import copy
 import os
 import sys
 from typing import Optional
@@ -78,7 +79,12 @@ def create_phono3py_force_constants(
     output_filename: Optional[str] = None,
     log_level=1,
 ):
-    """Read or calculate force constants."""
+    """Read or calculate force constants.
+
+    This function is for the 'phonopy' command only and not for the
+    'phonopy-load' command.
+
+    """
     # Only for build-in fc calculator.
     # These are not applied to external fc calculators.
     symmetrize_fc3r = settings.is_symmetrize_fc3_r or settings.fc_symmetry
@@ -205,19 +211,43 @@ def parse_forces(
     fc_type=None,
     log_level=0,
 ):
-    """Read displacements and forces."""
-    filename_read_from = None
+    """Read displacements and forces.
+
+    Physical units of displacements and forces are converted following the
+    calculator name. The calculator name may be given as the user input or found
+    in phono3py-yaml file. When dumping to phono3py-yaml file, it is assumed
+    that displacements and forces are written in the default units (A and eV/A)
+    without writing calculator name in it.
+
+    """
+    filename_read_from: Optional[str] = None
+
+    calculator = phono3py.calculator
+    # Get dataset from ph3py_yaml. dataset can be None.
+    # physical_units can be overwritten if calculator is found in ph3py_yaml.
+    dataset = _extract_dataset_from_ph3py_yaml(ph3py_yaml, fc_type)
+    if dataset and ph3py_yaml.calculator:
+        calculator = ph3py_yaml.calculator
+
+    physical_units = get_default_physical_units(calculator)
 
     if fc_type == "phonon_fc2":
         natom = len(phono3py.phonon_supercell)
     else:
         natom = len(phono3py.supercell)
 
-    # Get dataset from ph3py_yaml. dataset can be None.
-    dataset = _extract_dataset_from_ph3py_yaml(ph3py_yaml, fc_type)
     if dataset:
         filename_read_from = phono3py_yaml_filename
 
+        # Units of displacements and forces are converted. If forces don't
+        # exist, the convesion will not be performed for forces.
+        if calculator is not None:
+            _convert_unit_in_dataset(
+                dataset,
+                distance_to_A=physical_units["distance_to_A"],
+                force_to_eVperA=physical_units["force_to_eVperA"],
+            )
+
     # Try to read FORCES_FC* if type-2 and return dataset.
     # None is returned unless type-2.
     # can emit FileNotFoundError.
@@ -230,6 +260,15 @@ def parse_forces(
             filename_read_from = force_filename
             dataset = _dataset
 
+            # Units of displacements and forces are converted.
+            if calculator is not None:
+                _convert_unit_in_dataset(
+                    dataset,
+                    distance_to_A=physical_units["distance_to_A"],
+                    force_to_eVperA=physical_units["force_to_eVperA"],
+                )
+
+    # No dataset is found. "disp_fc*.yaml" is read. But this is deprecated.
     if dataset is None:
         if disp_filename is None:
             msg = (
@@ -242,6 +281,20 @@ def parse_forces(
         dataset = _read_disp_fc_yaml(disp_filename, fc_type)
         filename_read_from = disp_filename
 
+        # No forces should exist. Therefore only unit of displacements is
+        # converted.
+        if calculator is None:
+            if log_level:
+                print(
+                    f'Displacements are read from "{disp_filename}", but '
+                    " the unit has not converted."
+                )
+        else:
+            _convert_unit_in_dataset(
+                dataset,
+                distance_to_A=physical_units["distance_to_A"],
+            )
+
     if "natom" in dataset and dataset["natom"] != natom:
         msg = (
             "Number of atoms in supercell is not consistent with "
@@ -273,12 +326,18 @@ def parse_forces(
         else:
             parse_FORCES_FC3(dataset, filename=force_filename)
 
+        # Unit of displacements is already converted.
+        # Therefore, only unit of forces is converted.
+        if calculator is not None:
+            _convert_unit_in_dataset(
+                dataset,
+                force_to_eVperA=physical_units["force_to_eVperA"],
+            )
+
         if log_level:
             print('Sets of supercell forces were read from "%s".' % force_filename)
             sys.stdout.flush()
 
-    _convert_unit_in_dataset(dataset, phono3py.calculator)
-
     return dataset
 
 
@@ -567,11 +626,16 @@ def _create_phono3py_phonon_fc2(
     )
 
 
-def _convert_unit_in_dataset(dataset, calculator):
-    physical_units = get_default_physical_units(calculator)
-    force_to_eVperA = physical_units["force_to_eVperA"]
-    distance_to_A = physical_units["distance_to_A"]
+def _convert_unit_in_dataset(
+    dataset: dict,
+    distance_to_A: Optional[float] = None,
+    force_to_eVperA: Optional[float] = None,
+) -> None:
+    """Convert physical units of displacements and forces in dataset.
+
+    dataset is overwritten.
 
+    """
     if "first_atoms" in dataset:
         for d1 in dataset["first_atoms"]:
             if distance_to_A is not None:
@@ -608,8 +672,8 @@ def _extract_dataset_from_ph3py_yaml(ph3py_yaml: Optional[Phono3pyYaml], fc_type
     dataset = None
     if fc_type == "phonon_fc2":
         if ph3py_yaml and ph3py_yaml.phonon_dataset is not None:
-            dataset = ph3py_yaml.phonon_dataset
+            dataset = copy.deepcopy(ph3py_yaml.phonon_dataset)
     else:
         if ph3py_yaml and ph3py_yaml.dataset is not None:
-            dataset = ph3py_yaml.dataset
+            dataset = copy.deepcopy(ph3py_yaml.dataset)
     return dataset