Temperature dependent force constants calculation using pypolymlp and symfc

Temperature dependent force constants calculation using pypolymlp and symfc#

This is an experimental feature, and its usage may change occasionally.

With the --pypolymlp option, phonopy can interface with the polynomial machine learning potential (MLP) code, pypolymlp, to perform training and evaluation tasks of MLPs. This feature aims to reduce the computational cost of anharmonic force constant calculations by using MLPs as an intermediary layer, efficiently representing atomic interactions. The example is found at example/KCl-SSCHA.

The training process involves using a dataset consisting of supercell displacements, forces, and energies. The trained MLPs are then employed to compute forces for supercells with specific displacements.

For further details on combining phonopy calculations with pypolymlp, refer to A. Togo and A. Seko, J. Chem. Phys. 160, 211001 (2024) [doi] [arxiv].

Using the polynomial MLPs, stochastic self-consistent harmonic approximation (SSCHA) calculation is performed in the following sections. By this, temperature dependent force constants are calculated within SSCHA. About SSCHA, please refer the paper by L. Monacelli et al., J. Phys.: Condens. Matter 33 363001 (2021) and A. van Roekeghem A, et al., Comput. Phys. Commun. 263 107945 (2021). Technically, the computational procedure introduced here is equivalent to the approach of the latter paper.

Citation of pypolymlp#

“Tutorial: Systematic development of polynomial machine learning potentials for elemental and alloy systems”, A. Seko, J. Appl. Phys. 133, 011101 (2023) [doi].

@article{pypolymlp,
  author = {Seko, Atsuto},
  title = "{"Tutorial: Systematic development of polynomial machine learning potentials for elemental and alloy systems"}",
  journal = {J. Appl. Phys.},
  volume = {133},
  number = {1},
  pages = {011101},
  year = {2023},
  month = {01},
}

Citation of symfc#

“Projector-based efficient estimation of force constants”, A. Seko and A. Togo, Phys. Rev. B, 110, 214302 (2024) [doi] [arxiv].

@article{PhysRevB.110.214302,
  title = {Projector-based efficient estimation of force constants},
  author = {Seko, Atsuto and Togo, Atsushi},
  journal = {Phys. Rev. B},
  volume = {110},
  issue = {21},
  pages = {214302},
  numpages = {18},
  year = {2024},
  month = {Dec},
}

Requirements#

  • pypolymlp >= 0.10.0

    For linux (x86-64), a compiled package of pypolymlp can be installed via conda-forge (recommended). Otherwise, pypolymlp can be installed from source-code.

How to calculate#

Workflow#

  1. Generate random displacements in supercells. Use –rd option.

  2. Calculate corresponding forces and energies in supercells. Use of VASP interface is recommended for –sp option is supported.

  3. Prepare dataset composed of displacements, forces, and energies in supercells. The dataset must be stored in a phonopy-yaml-like file, e.g., phonopy_params.yaml. Use -f and –sp option simultaneously.

  4. Develop MLPs. By default, 90 and 10 percents of the dataset are used for the training and test, respectively. At this step polymlp.yaml is saved.

  5. Generate random displacements in supercells

  6. Evaluate MLPs for forces of the supercells generated in step 5.

  7. Calculate force constants from displacement-force dataset from steps 5 and 6.

  8. Temperature dependent force constants calculation

The steps 4-7 are executed in running phonopy with --pypolymlp option.

Steps 1-3: Dataset preparation#

For the training, the following supercell data are required in the phonopy setting to use pypolymlp:

  • Displacements

  • Forces

  • Total energies

These data must be stored in phonopy.yaml-like file.

The supercells with displacements are generated by

% phonopy --pa auto --rd 1000 -c POSCAR-unitcell --dim 2 2 2 --amin 0.03 --amax 1.5
        _
  _ __ | |__   ___  _ __   ___   _ __  _   _
 | '_ \| '_ \ / _ \| '_ \ / _ \ | '_ \| | | |
 | |_) | | | | (_) | | | | (_) || |_) | |_| |
 | .__/|_| |_|\___/|_| |_|\___(_) .__/ \__, |
 |_|                            |_|    |___/
                                      2.31.1

Compiled with OpenMP support (max 10 threads).
Python version 3.12.6
Spglib version 2.5.0

Crystal structure was read from "POSCAR-unitcell".
Unit of length: angstrom
Displacements creation mode
  Number of supercells with random displacements: 1000
  Min displacement distance: 0.03
  Max displacement distance: 1.5
Settings:
  Supercell: [2 2 2]
  Primitive matrix (Auto):
    [0.  0.5 0.5]
    [0.5 0.  0.5]
    [0.5 0.5 0. ]
Spacegroup: Fm-3m (225)
Number of symmetry operations in supercell: 1536
Use -v option to watch primitive cell, unit cell, and supercell structures.

"phonopy_disp.yaml" and supercells have been created.

Summary of calculation was written in "phonopy_disp.yaml".
                 _
   ___ _ __   __| |
  / _ \ '_ \ / _` |
 |  __/ | | | (_| |
  \___|_| |_|\__,_|

For the generated supercells, forces and energies are calculated. Here it is assumed to use the VASP code. Once the calculations are complete, the data (forces and energies) can be extracted using the following command:

% phonopy-init --sp -f vasprun_xmls/vasprun-{001..120}.xml

This command extracts the necessary data and stores it in the phonopy_params.yaml file. For more details, refer to the description of the –sp option. Currently, supercell energy extraction from calculator outputs is only supported when using the VASP interface.

Steps 4: Develop MLPs#

Having phonopy_params.yaml, phonopy is executed with --pypolymlp option,

% phonopy phonopy_mlpsscha_params_KCl-120.yaml.xz --pypolymlp --mlp-params="ntrain=100, ntest=20" -v
        _
  _ __ | |__   ___  _ __   ___   _ __  _   _
 | '_ \| '_ \ / _ \| '_ \ / _ \ | '_ \| | | |
 | |_) | | | | (_) | | | | (_) || |_) | |_| |
 | .__/|_| |_|\___/|_| |_|\___(_) .__/ \__, |
 |_|                            |_|    |___/
                                 4.4.1.dev12

-------------------------[time 2026-07-20 11:48:04]-------------------------
Rust backend (phonors) using rayon (10 threads).
Running in phonopy.load mode.
Python version 3.13.11
Spglib version 2.7.0

Crystal structure was read from "phonopy_mlpsscha_params_KCl-120.yaml.xz".
Unit of length: angstrom
Settings:
  Supercell: [2 2 2]
  Primitive matrix:
    [0.  0.5 0.5]
    [0.5 0.  0.5]
    [0.5 0.5 0. ]
Spacegroup: Fm-3m (225)
Number of symmetry operations in supercell: 1536
(Crystal structures are omitted here.)

NAC parameters were read from "phonopy_mlpsscha_params_KCl-120.yaml.xz".
(Dielectric constant and Born effective charges are omitted here.)

Displacement-force dataset was read from "phonopy_mlpsscha_params_KCl-120.yaml.xz".
----------------------------- pypolymlp start ------------------------------
Pypolymlp version 0.18.9.post40
Pypolymlp is a generator of polynomial machine learning potentials.
Please cite the paper: A. Seko, J. Appl. Phys. 133, 011101 (2023).
Pypolymlp is developed at https://github.com/sekocha/pypolymlp.
Parameters:
  cutoff: 8.0
  model_type: 3
  max_p: 2
  gtinv_order: 3
  gtinv_maxl: (8, 8)
  gaussian_params1: (1.0, 1.0, 1)
  gaussian_params2: (0.0, 7.0, 10)
  ntrain: 100
  ntest: 20
Developing MLPs by pypolymlp...
n_features: 8283
Minimum memory required for Cholesky solver in GB: 1.1
Memory required for allocating X additionally.
----- Dataset: data1 -----
Structures: 100 / 100
 Matrix shape (X): (19900,8283)
 Required memory for X: 1.3187 (GB)
 Estimated peak memory allocation (X.T @ X, X): 2.42 (GB)
Compute X.T @ X and X.T @ y
Regression: Use standard ridge solver.
Regression: cholesky decomposition.
- alpha: 0.001
  Compute X.T @ X + alpha @ I
  Solve linear equation
- alpha: 0.01
  Compute X.T @ X + alpha @ I
  Solve linear equation
- alpha: 0.1
  Compute X.T @ X + alpha @ I
  Solve linear equation
- alpha: 1.0
  Compute X.T @ X + alpha @ I
  Solve linear equation
- alpha: 10.0
  Compute X.T @ X + alpha @ I
  Solve linear equation
----- Dataset: data2 -----
Structures: 20 / 20
 Matrix shape (X): (3980,8283)
 Required memory for X: 0.26373 (GB)
 Estimated peak memory allocation (X.T @ X, X): 1.36 (GB)
Compute X.T @ X and X.T @ y
Regression: model selection ...
- alpha = 1.000e-03 : rmse (train, test) = 0.02432 0.23669
- alpha = 1.000e-02 : rmse (train, test) = 0.03613 0.16766
- alpha = 1.000e-01 : rmse (train, test) = 0.07193 0.22140
- alpha = 1.000e+00 : rmse (train, test) = 0.11563 0.26042
- alpha = 1.000e+01 : rmse (train, test) = 0.19375 0.31767
prediction: train-data1
  rmse_energy:      26.14364 (meV/atom)
  rmse_force:        9.98259 (eV/ang)
  mae_energy:        9.09094 (meV/atom)
  mae_force:         0.60389 (eV/ang)
prediction: test-data2
  rmse_energy:      40.23820 (meV/atom)
  rmse_force:        2.15935 (eV/ang)
  mae_energy:       17.99147 (meV/atom)
  mae_force:         0.41367 (eV/ang)
MLPs were written into "polymlp.yaml"
------------------------------ pypolymlp end -------------------------------
Use --rd or -d option for running phonon calculations with pypolymlp.

Summary of calculation was written in "phonopy.yaml".
-------------------------[time 2026-07-20 11:49:09]-------------------------
                 _
   ___ _ __   __| |
  / _ \ '_ \ / _` |
 |  __/ | | | (_| |
  \___|_| |_|\__,_|

Information about the development of MLPs using pypolymlp is provided between the pypolymlp start and pypolymlp end sections. The -v option is specified here to show the prediction errors of the developed MLPs. Those reported for test-data2 are the errors for the test dataset, which was not used for the training, and therefore indicate how accurately the MLPs predict energies and forces of unseen supercells. The polynomial MLPs are saved in the polymlp.yaml file. This file is automatically searched in subsequent phonopy executions with the --pypolymlp option and reused.

Step 5-8: Temperature dependent force constants calculation#

After the MLPs are developed, random displacements are generated with a displacement distance of 0.01 Angstrom. The forces for these supercells are then evaluated using pypolymlp. Both the generated displacements and the corresponding forces are stored in the phonopy_mlp_eval_dataset.yaml file. The calculated force constants may be referred as the harmonic force constants.

After the last step, the polymlp.yaml file exists in the current directory. This file is read automatically in the next calculation with the --pypolymlp option. If the developed MLPs can predict well forces at relatively large displacements, temperature dependent force constants are calculated with the --sscha NUMBER_OF_ITERATIONS option.

% phonopy phonopy_mlpsscha_params_KCl-120.yaml.xz --pypolymlp --sscha 10 --rd-temperature 300 --rd 1000
        _
  _ __ | |__   ___  _ __   ___   _ __  _   _
 | '_ \| '_ \ / _ \| '_ \ / _ \ | '_ \| | | |
 | |_) | | | | (_) | | | | (_) || |_) | |_| |
 | .__/|_| |_|\___/|_| |_|\___(_) .__/ \__, |
 |_|                            |_|    |___/
                                 4.4.1.dev12

-------------------------[time 2026-07-20 11:20:18]-------------------------
Rust backend (phonors) using rayon (10 threads).
Running in phonopy.load mode.
Python version 3.13.11
Spglib version 2.7.0

Crystal structure was read from "phonopy_mlpsscha_params_KCl-120.yaml.xz".
Unit of length: angstrom
Pypolymlp displacements creation mode
  Random displacements
  Number of supercells with random displacements: 1000
  Temperatuere to generate random displacements: 300.0
Settings:
  Supercell: [2 2 2]
  Primitive matrix:
    [0.  0.5 0.5]
    [0.5 0.  0.5]
    [0.5 0.5 0. ]
Spacegroup: Fm-3m (225)
Number of symmetry operations in supercell: 1536
Use -v option to watch primitive cell, unit cell, and supercell structures.

NAC parameters were read from "phonopy_mlpsscha_params_KCl-120.yaml.xz".
Displacement-force dataset was read from "phonopy_mlpsscha_params_KCl-120.yaml.xz".
----------------------------- pypolymlp start ------------------------------
Pypolymlp version 0.18.9.post40
Pypolymlp is a generator of polynomial machine learning potentials.
Please cite the paper: A. Seko, J. Appl. Phys. 133, 011101 (2023).
Pypolymlp is developed at https://github.com/sekocha/pypolymlp.
Load MLPs from "polymlp.yaml".
------------------------------ pypolymlp end -------------------------------
Generate random displacements
  Displacement distance: 0.01
  Plus-minus displacements: auto
Evaluate forces in 1000 supercells by pypolymlp
Dataset generated using MLPs was written in "phonopy_mlp_eval_dataset.yaml".
-------------------------------- Symfc start -------------------------------
Symfc version 1.7.1 (https://github.com/symfc/symfc)
Citation: A. Seko and A. Togo, Phys. Rev. B, 110, 214302 (2024)
Computing [2] order force constants.
Increase log-level to watch detailed symfc log.
--------------------------------- Symfc end --------------------------------
Max drift of force constants: -0.00000000 (yy) -0.00000000 (yy)
Max drift after symmetrization by symfc projector: -0.00000000 (yy) -0.00000000 (yy)

------------------------------- SSCHA start --------------------------------
Use provided force constants.

[ SSCHA iteration 1 / 10 ]
Generate 1000 supercells with displacements at 300.0 K
  [0.006, 0.081] ****
  [0.081, 0.156] *****************
  [0.156, 0.231] ***************************
  [0.231, 0.306] **************************
  [0.306, 0.381] ****************
  [0.381, 0.456] *******
  [0.456, 0.532] **
  [0.532, 0.607] *
  [0.607, 0.682]
  [0.682, 0.757]
Evaluate MLP to obtain forces using pypolymlp
Calculate force constants using symfc
SSCHA free energy: -98.107 +/- 0.089 meV
SSCHA force constants are written into "phonopy_sscha_fc_1.yaml.xz".

[ SSCHA iteration 2 / 10 ]
Generate 1000 supercells with displacements at 300.0 K
  [0.004, 0.083] ****
  [0.083, 0.162] *******************
  [0.162, 0.241] *****************************
  [0.241, 0.320] *************************
  [0.320, 0.399] ***************
  [0.399, 0.478] ******
  [0.478, 0.557] **
  [0.557, 0.636]
  [0.636, 0.715]
  [0.715, 0.793]
Evaluate MLP to obtain forces using pypolymlp
Calculate force constants using symfc
SSCHA free energy: -98.193 +/- 0.069 meV
SSCHA force constants are written into "phonopy_sscha_fc_2.yaml.xz".

(Iterations 3 to 9 are omitted here.)

[ SSCHA iteration 10 / 10 ]
Generate 1000 supercells with displacements at 300.0 K
  [0.005, 0.079] ***
  [0.079, 0.152] ****************
  [0.152, 0.225] **************************
  [0.225, 0.298] **************************
  [0.298, 0.371] *****************
  [0.371, 0.445] ********
  [0.445, 0.518] ***
  [0.518, 0.591] *
  [0.591, 0.664]
  [0.664, 0.738]
Evaluate MLP to obtain forces using pypolymlp
Calculate force constants using symfc
SSCHA free energy: -98.153 +/- 0.070 meV
SSCHA force constants are written into "phonopy_sscha_fc_10.yaml.xz".

-------------------------------- SSCHA end ---------------------------------
----------------------------------------------------------------------------
 No run mode was specified, so no phonon calculation was performed.
 Specify one of the following to calculate phonons.
 - Mesh sampling (MESH, --mesh)
 - Q-points (QPOINTS, --qpoints)
 - Band structure (BAND, --band)
 - Animation (ANIME, --anime)
 - Modulation (MODULATION, --modulation)
 - Characters of Irreps (IRREPS, --irreps)
----------------------------------------------------------------------------

Summary of calculation was written in "phonopy.yaml".
-------------------------[time 2026-07-20 11:23:39]-------------------------
                 _
   ___ _ __   __| |
  / _ \ '_ \ / _` |
 |  __/ | | | (_| |
  \___|_| |_|\__,_|

The final force constants are stored in files named phonopy_sscha_fc_NUM.yaml.xz, where NUM represents the integer corresponding to the iteration step. By performing a sufficient number of SSCHA iterations and utilizing a sufficiently large set of supercells with random displacements at a given temperature, the SSCHA force constants can be reliably determined.

SSCHA free energy#

The SSCHA free energy printed at each iteration is defined for the force constants \(\Phi\) by

\[\mathcal{F}_\Phi = \tilde{F}_\Phi - \langle \tilde{V}_\Phi \rangle_{\tilde{\rho}_\Phi} + \langle V \rangle_{\tilde{\rho}_\Phi},\]

where \(\tilde{F}_\Phi\) and \(\langle \tilde{V}_\Phi \rangle_{\tilde{\rho}_\Phi}\) are the harmonic Helmholtz free energy and potential energy of \(\Phi\), respectively, and \(\langle V \rangle_{\tilde{\rho}_\Phi}\) is the potential energy. The averages are taken over the harmonic density matrix \(\tilde{\rho}_\Phi\) at the temperature, which is sampled by the supercells with random displacements. \(\langle V \rangle_{\tilde{\rho}_\Phi}\) is obtained from the supercell energies evaluated by the MLPs relative to the energy of the supercell without displacements, and the harmonic potential energy is evaluated from the displacements as

\[\langle \tilde{V}_\Phi \rangle_{\tilde{\rho}_\Phi} = \frac{1}{2} \sum_{l\kappa j, l'\kappa' j'} \Phi_{l\kappa j, l'\kappa' j'} \langle u_{l\kappa j} u_{l'\kappa' j'} \rangle_{\tilde{\rho}_\Phi}.\]

These terms are given per primitive cell. The notation and the description used here are those of equations (B1) and (B3) in appendix B of A. Togo et al., J. Phys.: Condens. Matter 34, 365401 (2022) [doi], where it is also shown that evaluating \(\langle \tilde{V}_\Phi \rangle_{\tilde{\rho}_\Phi}\) from the displacements gives a more stable measure of the convergence than evaluating it from the phonon frequencies and eigenvectors.

The convergence of these force constants is monitored through the SSCHA free energy printed at each iteration. The value after +/- is its statistical error, which originates from the random sampling of the supercells and is reduced by increasing the number of supercells given by the --rd option. The iterations are regarded as converged when the free energies of the successive iterations scatter within this error, because the remaining variation is then indistinguishable from the sampling noise. In the log shown above, the free energies of the ten iterations scatter by 0.081 meV, which is comparable to the statistical error of 0.07 meV, and hence they are converged.

As an illustration, the phonon band structures corresponding to the SSCHA force constants at the iteration steps can be overlaid:

% phonopy phonopy_mlp_eval_dataset.yaml --band auto --band-points 101; mv band.yaml band-0.yaml
% for i in {1..10}; do phonopy phonopy_sscha_fc_$i.yaml.xz --band auto --band-points 101; mv band.yaml band-$i.yaml; done
% phonopy-bandplot band-{0..10}.yaml --legend
_images/sscha-bands.png

band-0.yaml is the harmonic band structure, and the others are those of the SSCHA iterations. The SSCHA band structures lie on top of each other, whereas the harmonic one deviates from them, which shows that the SSCHA force constants are converged and are distinct from the harmonic force constants.

SSCHA free energy of saved force constants#

The SSCHA free energy is printed at every iteration, but it can also be evaluated afterwards from a saved phonopy_sscha_fc_NUM.yaml.xz file. Since this file contains force constants only, the anharmonic part is re-evaluated by sampling displacements from the canonical ensemble of these force constants and by evaluating the supercell energies using the MLPs in polymlp.yaml.

import phonopy
from phonopy.sscha.core import MLPSSCHA

ph = phonopy.load("phonopy_sscha_fc_10.yaml.xz", log_level=0)
ph.load_mlp("polymlp.yaml")

sscha = MLPSSCHA(ph, ph.mlp, temperature=300.0, number_of_snapshots=1000)
sscha.sample_supercells()
sscha.calculate_free_energy()

print(
    f"SSCHA free energy: {sscha.free_energy * 1000:.3f} "
    f"+/- {sscha.free_energy_error * 1000:.3f} meV"
)

This does not reproduce the value printed at the corresponding iteration exactly. In the iteration, the displacements are sampled from the force constants of the previous iteration, whereas here they are sampled from the force constants stored in the file. The two values therefore agree only within the statistical error.

Parameters for developing MLPs#

A few parameters can be specified using the --mlp-params option for the development of MLPs. The parameters are provided as a string, e.g.,

% phonopy phonopy_params.yaml --pypolymlp --mlp-params="ntrain=80, ntest=20"

Parameters are separated by commas for configuration. A brief explanation of the available parameters can be found in the docstring of PypolymlpParams that is found by

In [1]: from phonopy.interface.pypolymlp import PypolymlpParams

In [2]: help(PypolymlpParams)

ntrain and ntest are implemented in phonopy, while the remaining parameters are directly passed to pypolymlp. Optimizing pypolymlp parameters can be difficult, both in terms of achieving accuracy and managing the computational resources required. The current default parameters are likely suitable for systems up to ternary compounds. For binary systems, the calculations can generally be run on standard laptop computers, but for ternary systems, around 40 GB of memory or more may be necessary.

For parameter adjustments, it is recommended to consult the pypolymlp documentation and review the relevant research papers.

ntrain and ntest#

This method provides a straightforward dataset split: the first ntrain supercells from the list are used for training, while the last ntest supercells are reserved for testing.

Convergence with respect to dataset size#

In general, increasing the amount of data improves the accuracy of representing force constants. Therefore, it is recommended to check the convergence of the target property with respect to the number of supercells in the training dataset. The SSCHA free energy is convenient to monitor, since it is a single number at each temperature and is printed at every SSCHA iteration.

For example, by preparing an initial set with 100 supercell data, calculations can then be performed by varying the size of the training dataset while keeping the test dataset unchanged as follows:

% for n in 20 40 60 80 100; do
    mkdir ntrain-$n
    cd ntrain-$n
    phonopy ../phonopy_params.yaml --pypolymlp --mlp-params="ntrain=$n, ntest=20" \
            --sscha 10 --rd-temperature 300 --rd 1000 | tee log
    cd ..
  done
% grep "SSCHA free energy" ntrain-*/log

Each calculation is performed in a separate directory. This is necessary because an existing polymlp.yaml in the current directory is loaded and reused, by which the ntrain setting would be silently ignored.

The SSCHA free energies are compared against the size of the training dataset to observe the convergence. Since the SSCHA free energy fluctuates by the stochastic sampling of random displacements, the comparison should be made using the same number of supercells (--rd) and the values at the later iterations. If it has not converged, an additional set of supercell data (e.g., forces and energies in the next 100 supercells) will be computed and included. With this procedure in mind, it may be convenient to generate a sufficiently large number of supercells with random displacements in advance, such as 1000 supercells, before starting the temperature dependent force constants calculation with pypolymlp.

Converting phonopy.pmlp to polymlp.yaml#

In older versions, polynomial MLPs were stored in phonopy.pmlp. This file can be converted to polymlp.yaml using the following Python snippet.

from pypolymlp.mlp_dev.pypolymlp import Pypolymlp
polymlp = Pypolymlp()
polymlp.convert_to_yaml(filename_txt="phonopy.pmlp", filename_yaml="polymlp.yaml”)