Fix handling datasets with random disps and update docs

doc/command-options.md

@@ -174,6 +174,25 @@ created from `FORCES_FC2` and `phono3py_disp.yaml` instead of `FORCES_FC3` and
 % phono3py --cfs --dim-fc2="x x x"
 ```
 
+### `--sp` or `--save-params`
+
+Instead of `FORCES_FC3`, `phono3py_params.yaml` is generated. This option must
+be used with `--cf3`, and optionally with `--cf2`. If the force calculator
+supports reading energy of supercell, those are written into
+`phono3y_params.yaml`. These energies are necessary for using `--pypolymlp`
+option.
+
+```bash
+% phono3py --cf3 disp-{00001..00755}/vasprun.xml --sp
+```
+
+When using with `--cf2`, `--cf3` has to be specified simultaneously as below,
+
+```bash
+% phono3py --cf3 disp-{00001..00755}/vasprun.xml --cf2 disp_fc2-{00001..00002}/vasprun.xml --sp
+```
+
+
 ## Supercell, primitive cell, masses, magnetic moments
 
 (dim_option)=

doc/random-displacements.md

@@ -16,7 +16,7 @@ presented with using symfc that can be installed via pip or conda easily.
 
 The option `--rd NUM` is used instead of `-d` in generating displacements as follows:
 
-```
+```bash
 % phono3py --rd 100 --dim 2 2 2 --pa auto -c POSCAR-unitcell
 ```
 
@@ -25,3 +25,12 @@ must be specified, and the initial guess may be from around the number of
 supecells generated for the systematic displacements by `-d`. In the case of the
 `NaCl-rd` example, 146 supercells are generated with `-d`, so similar
 number `--rd 100` was chosen here.
+
+If random directional displacements for fc2 are expected, `--rd-fc2` and
+`--dim-fc2` have to be specified:
+
+```bash
+% phono3py --rd 100 --dim 2 2 2 --rd-fc2 2 --dim-fc2 4 4 4 --pa auto -c POSCAR-unitcell
+```
+
+where `--dim` is necessary but `--rd` is not.

example/NaCl-rd/README.md

@@ -6,7 +6,7 @@ This example utilizes an NaCl calculation result from A. Togo and A. Seko, J.
 Chem. Phys. 160, 211001 (2024). Supercells of 2x2x2 and 4x4x4 conventional unit
 cells are chosen for the third-order force constants (fc3) and second-order
 force constants (fc2), respectively. Displacement-force datasets consisting of
-100 supercells for fc3 and 4 supercells for fc2 are extracted and stored in
+100 supercells for fc3 and 2 supercells for fc2 are extracted and stored in
 `phono3py_params_NaCl.yaml.xz`. Random directional displacements of a constant
 0.03 Angstrom are used.
 
@@ -26,29 +26,30 @@ Lattice thermal conductivity (LTC) is calculated with the following command:
 % phono3py-load phono3py_params_NaCl.yaml.xz --br --ts 300 --mesh 50
 ```
 
-By this, LTC is obtained around 7.8 W/m-k.
+By this, LTC is obtained around 8.3 W/m-k.
 
 
 ## How to use pypolymlp
 
-The polynomial machine learning potential (poly-MLP) by pypolymlp can be used to
+When supercell energies are included in `phono3py_params.yaml` like file, the
+polynomial machine learning potential (poly-MLP) by pypolymlp can be used to
 calculate fc3 by the following command:
 
-```
+```bash
 % phono3py-load phono3py_params_NaCl.yaml.xz --pypolymlp --symfc --rd 400 -v
 ```
 
 the procedure below is performed:
 
 1. Poly-MLPs are computed from the displacement-force dataset for fc3. This is
-   activated by the `--pypolymlp` option.
+   activated by `--pypolymlp` option.
 2. 800=400+400 supercells for random directional displacements are generated,
    where 400+400 means 400 supercells with random displacements (u) and 400
-   supercells with opposite displacement vectors (-u). This is activated by the
-   `--rd 400` option. The default displacement distance is 0.001 Angstrom in the
+   supercells with opposite displacement vectors (-u). This is activated by
+   `--rd 400` option. The default displacement distance is 0.001 Angstrom in
    `--pypolymlp` mode. Since random displacements are generated `--symfc` has to
    be specified for fc3. In this example, random displacements are used for fc2,
-   too, `--symfc` is applied to both of fc3 and fc2. Without the `--rd` option,
+   too, `--symfc` is applied to both of fc3 and fc2. Without `--rd` option,
    systematic displacements are generated, for which the option `--fc-calc "symfc|"`
    has to be specified instead of `--symfc` (equivalent to `--fc-calc "symfc|symfc")`).
 3. Forces on atoms in these 800 supercells are calculated using poly-MLP.
@@ -58,9 +59,22 @@ the procedure below is performed:
 The `fc3.hdf5` and `fc2.hdf5` are obtained. Using these force constants, LTC is
 calculated by
 
-```
+```bash
 % phono3py-load phono3py_params_NaCl.yaml.xz --br --ts 300 --mesh 50
 ```
 
-and the LTC value of around 7.8 W/m-k is obtained. This LTC value is equivalent
-to above, but of course, it can be different.
+and the LTC value of around 8.2 W/m-k is obtained.
+
+## Generating phono3py_params.yaml from vasprun.xml's
+
+`phono3py_params.yaml` is generated from
+
+```bash
+% phono3py phono3py_disp.yaml --cf3 NaCl-vasprun/vasprun-{00001..00100}.xml --cf2 NaCl-vasprun/vasprun-ph0000{1,2}.xml --sp
+```
+
+This command reads electronic energies of supercells from `vasprun.xml`s and
+writes them into `phono3py_params.yaml`, too. Here, `phono3py_disp.yaml` is not
+included in this example, but `phono3py_params_NaCl.yaml.xz` can be used to run
+this example since corresponding information of displacements is included in
+this file, too.

phono3py/cui/create_supercells.py

@@ -35,8 +35,10 @@
 # POSSIBILITY OF SUCH DAMAGE.
 
 from phonopy.interface.calculator import write_supercells_with_displacements
+from phonopy.structure.cells import print_cell
 
 from phono3py import Phono3py
+from phono3py.cui.show_log import print_supercell_matrix
 from phono3py.interface.calculator import (
     get_additional_info_to_write_fc2_supercells,
     get_additional_info_to_write_supercells,
@@ -92,6 +94,16 @@ def create_phono3py_supercells(
     if log_level:
         print("")
         print('Unit cell was read from "%s".' % optional_structure_info[0])
+        print("-" * 32 + " unit cell " + "-" * 33)  # 32 + 11 + 33 = 76
+        print_cell(phono3py.unitcell)
+        print("-" * 76)
+        print_supercell_matrix(
+            phono3py.supercell_matrix, phono3py.phonon_supercell_matrix
+        )
+        if phono3py.primitive_matrix is not None:
+            print("Primitive matrix:")
+            for v in phono3py.primitive_matrix:
+                print("  %s" % v)
         print("Displacement distance: %s" % distance)
 
     ids = []
@@ -125,6 +137,7 @@ def create_phono3py_supercells(
             print("Number of displacement supercell files created: %d" % num_disp_files)
 
     if phono3py.phonon_supercell_matrix is not None:
+        num_disps = len(phono3py.phonon_supercells_with_displacements)
         additional_info = get_additional_info_to_write_fc2_supercells(
             interface_mode, phono3py.phonon_supercell_matrix
         )

phono3py/cui/load.py

@@ -588,9 +588,9 @@ def _get_dataset_phonon_dataset_or_fc2(
         )
     elif (
         forces_fc2_filename is not None or pathlib.Path("FORCES_FC2").exists()
-    ) and ph3py.phonon_supercell_matrix:
+    ) and ph3py.phonon_supercell_matrix is not None:
         if forces_fc2_filename is None:
-            force_filename = forces_fc2_filename
+            force_filename = "FORCES_FC2"
         else:
             force_filename = forces_fc2_filename
         phonon_dataset = _get_dataset_for_fc2(

phono3py/cui/phono3py_script.py

@@ -36,8 +36,10 @@
 
 from __future__ import annotations
 
+import argparse
 import datetime
 import sys
+from typing import Optional
 
 import numpy as np
 from phonopy.cui.collect_cell_info import collect_cell_info
@@ -75,7 +77,7 @@ from phono3py.cui.load import (
     set_dataset_and_force_constants,
 )
 from phono3py.cui.phono3py_argparse import get_parser
-from phono3py.cui.settings import Phono3pyConfParser
+from phono3py.cui.settings import Phono3pyConfParser, Phono3pySettings
 from phono3py.cui.show_log import (
     show_general_settings,
     show_phono3py_cells,
@@ -212,7 +214,7 @@ def get_run_mode(settings):
     return run_mode
 
 
-def start_phono3py(**argparse_control):
+def start_phono3py(**argparse_control) -> tuple[argparse.Namespace, int]:
     """Parse arguments and set some basic parameters."""
     parser, deprecated = get_parser(**argparse_control)
     args = parser.parse_args()
@@ -310,7 +312,9 @@ def get_input_output_filenames_from_args(args):
     return input_filename, output_filename
 
 
-def get_cell_info(settings, cell_filename, log_level):
+def get_cell_info(
+    settings: Phono3pySettings, cell_filename: str, log_level: int
+) -> dict:
     """Return calculator interface and crystal structure information."""
     cell_info = collect_cell_info(
         supercell_matrix=settings.supercell_matrix,
@@ -513,6 +517,46 @@ def grid_addresses_to_grid_points(grid_addresses, bz_grid):
     return bz_grid.grg2bzg[grid_points]
 
 
+def create_supercells(
+    settings: Phono3pySettings,
+    cell_info: dict,
+    confs_dict: dict,
+    interface_mode: Optional[str],
+    symprec: float,
+    log_level: int,
+):
+    """Create supercells and write displacements."""
+    if (
+        settings.create_displacements
+        or settings.random_displacements
+        or settings.random_displacements_fc2
+    ):
+        phono3py = create_phono3py_supercells(
+            cell_info,
+            settings,
+            symprec,
+            interface_mode=interface_mode,
+            log_level=log_level,
+        )
+
+        if log_level:
+            if phono3py.supercell.magnetic_moments is None:
+                print("Spacegroup: %s" % phono3py.symmetry.get_international_table())
+            else:
+                print(
+                    "Number of symmetry operations in supercell: %d"
+                    % len(phono3py.symmetry.symmetry_operations["rotations"])
+                )
+
+        finalize_phono3py(
+            phono3py,
+            confs_dict,
+            log_level,
+            write_displacements=True,
+            filename="phono3py_disp.yaml",
+        )
+
+
 def store_force_constants(
     phono3py: Phono3py,
     settings,
@@ -853,7 +897,7 @@ def main(**argparse_control):
     # warnings.simplefilter("error")
     load_phono3py_yaml = argparse_control.get("load_phono3py_yaml", False)
 
-    if "args" in argparse_control:  # For pytest
+    if "args" in argparse_control:  # This is for pytest.
         args = argparse_control["args"]
         log_level = args.log_level
     else:
@@ -926,29 +970,9 @@ def main(**argparse_control):
     ######################################################
     # Create supercells with displacements and then exit #
     ######################################################
-    if settings.create_displacements:
-        phono3py = create_phono3py_supercells(
-            cell_info,
-            settings,
-            symprec,
-            interface_mode=interface_mode,
-            log_level=log_level,
-        )
-
-        if phono3py.supercell.magnetic_moments is None:
-            print("Spacegroup: %s" % phono3py.symmetry.get_international_table())
-        else:
-            print(
-                "Number of symmetry operations in supercell: %d"
-                % len(phono3py.symmetry.symmetry_operations["rotations"])
-            )
-
-        finalize_phono3py(
-            phono3py,
-            confs_dict,
-            log_level,
-            write_displacements=True,
-            filename="phono3py_disp.yaml",
+    if not settings.use_pypolymlp:
+        create_supercells(
+            settings, cell_info, confs_dict, interface_mode, symprec, log_level
         )
 
     #######################

phono3py/cui/show_log.py

@@ -35,6 +35,8 @@
 # POSSIBILITY OF SUCH DAMAGE.
 
 import sys
+from collections.abc import Sequence
+from typing import Optional, Union
 
 import numpy as np
 from phonopy.structure.cells import print_cell
@@ -68,6 +70,23 @@ def show_general_settings(
         print("Calculator interface: %s" % phono3py.calculator)
     print('Crystal structure was read from "%s".' % cell_filename)
 
+    print_supercell_matrix(supercell_matrix, phonon_supercell_matrix)
+
+    if is_primitive_axes_auto:
+        print("Primitive matrix (Auto):")
+        for v in primitive_matrix:
+            print("  %s" % v)
+    elif primitive_matrix is not None:
+        print("Primitive matrix:")
+        for v in primitive_matrix:
+            print("  %s" % v)
+
+
+def print_supercell_matrix(
+    supercell_matrix: Union[Sequence, np.ndarray],
+    phonon_supercell_matrix: Optional[Union[Sequence, np.ndarray]] = None,
+):
+    """Print supercell matrix."""
     if (np.diag(np.diag(supercell_matrix)) - supercell_matrix).any():
         print("Supercell matrix (dim):")
         for v in supercell_matrix:
@@ -81,14 +100,6 @@ def show_general_settings(
                 print("  %s" % v)
         else:
             print("Phonon supercell (dim-fc2): %s" % np.diag(phonon_supercell_matrix))
-    if is_primitive_axes_auto:
-        print("Primitive matrix (Auto):")
-        for v in primitive_matrix:
-            print("  %s" % v)
-    elif primitive_matrix is not None:
-        print("Primitive matrix:")
-        for v in primitive_matrix:
-            print("  %s" % v)
 
 
 def show_phono3py_cells(phono3py: Phono3py):

phono3py/file_IO.py

@@ -200,16 +200,26 @@ def write_FORCES_FC2(disp_dataset, forces_fc2=None, fp=None, filename="FORCES_FC
     else:
         w = fp
 
-    for i, disp1 in enumerate(disp_dataset["first_atoms"]):
-        w.write("# File: %-5d\n" % (i + 1))
-        w.write("# %-5d " % (disp1["number"] + 1))
-        w.write("%20.16f %20.16f %20.16f\n" % tuple(disp1["displacement"]))
-        if "forces" in disp1 and forces_fc2 is None:
-            force_set = disp1["forces"]
+    if "first_atoms" in disp_dataset:
+        for i, disp1 in enumerate(disp_dataset["first_atoms"]):
+            w.write("# File: %-5d\n" % (i + 1))
+            w.write("# %-5d " % (disp1["number"] + 1))
+            w.write("%20.16f %20.16f %20.16f\n" % tuple(disp1["displacement"]))
+            if "forces" in disp1 and forces_fc2 is None:
+                force_set = disp1["forces"]
+            else:
+                force_set = forces_fc2[i]
+            for forces in force_set:
+                w.write("%15.10f %15.10f %15.10f\n" % tuple(forces))
+    else:
+        if "forces" in disp_dataset:
+            write_FORCE_SETS(disp_dataset, filename="FORCES_FC2")
         else:
-            force_set = forces_fc2[i]
-        for forces in force_set:
-            w.write("%15.10f %15.10f %15.10f\n" % tuple(forces))
+            if forces_fc2 is None:
+                raise RuntimeError("No forces are found.")
+            dataset = disp_dataset.copy()
+            dataset["forces"] = forces_fc2
+            write_FORCE_SETS(dataset, filename="FORCES_FC2")
 
     if fp is None:
         w.close()