Add anisotropic NpT with MTK equations

Kohei Shinohararequested to merge
lan496/ase:aniso-mtk-pr into master

Summary

This MR adds anisotropic NPT by MTK equations, which uses Nose-Hoover-like chains for both thermostats and barostats. Both cell volume and shape fluctuate by this integrator. The barostat for the cell fluctuations is controlled by damping parameter pdamp (in time unit).

This MR is the successor of !3550 (merged), and my implementation note for the MTK equations is in the same PDF for the previous Nose-Hoover chain implementation, https://github.com/lan496/lan496.github.io/blob/main/notes/nose_hoover_chain/main.pdf.

Testing

Conserved quantity

In ase/test/md/test_nose_hoover_chain.py, I have checked the new anisotropic NPT integrator conserves the following energy-like quantity.

image

I tried to generalize MTK equations to allow fluctuating only specified axes. However, I failed to modify the energy-like quantity to be conserved during MD. Thus, I would like to proceed with this MR without considering mask keyword for now.

External pressure

First, I performed NPT with this new integrator by changing external pressures pressure_au from 0 to 100 GPa. The figure below shows the volumes during NPT trajectories, and the equilibrium volumes decrease by applying the external pressure as expected.

image

Barostat damping parameter (pdamp) dependency

Next, I performed NPT by changing pdamp parameter with pressure_au=10 * GPa and tdamp=100*timestep. The figure below shows the instantaneous stress tensors during NPT trajectories. Here I define the stress tensors by pxx = -atoms.get_stress()[0, 0] and so on. All components of the diagonal stress (pxx, pyy, pzz) are close to the external pressure (10 GPa) as expected. Also, the stress tensors are well-controlled in wide range of pdamp from 100*timestep to 10000*timestep.

image

I checked the distributions of instantaneous pressure with Q-Q plot for these trajectories. If the distribution is normal, the ordered points should line up straight. All cases distribute like normal distributions within 2σ regions.

image

Scripts

Scripts for the testings 
from typing import Union, IO
import weakref
from concurrent.futures import ThreadPoolExecutor, as_completed

import numpy as np
import asap3
from ase.parallel import world
from ase.md.md import MolecularDynamics
from ase.md.logger import MDLogger
from ase import Atoms
import ase.units
import ase.build
from ase.calculators.emt import EMT
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution, Stationary

from ase.md.nose_hoover_chain import NoseHooverChainNVT, IsotropicMTKNPT

class CustomMDLogger(MDLogger):

    def __init__(
        self,
        dyn: MolecularDynamics,
        atoms: Atoms,
        logfile: Union[IO, str],
    ):
        self.dyn: MolecularDynamics = weakref.proxy(dyn)
        self.atoms = atoms

        self.logfile = self.openfile(file=logfile, mode='w', comm=world)

        self.hdr = "step,time(ps),epot(eV/atom),ekin(eV/atom),etot(eV/atom),temperature(K),volume(Ang3/atom),pxx(GPa),pyy(GPa),pzz(GPa),pyz(GPa),pzx(GPa),pxy(GPa)"
        self.fmt = "%d" + ",%.8f" * 12 + "\n"

        self.logfile.write(self.hdr + "\n")

    def __call__(self):
        num_atoms = len(self.atoms)
        epot_eV_per_atom = self.atoms.get_potential_energy() / num_atoms
        ekin_eV_per_atom = self.atoms.get_kinetic_energy() / num_atoms
        etot_eV_per_atom = epot_eV_per_atom + ekin_eV_per_atom
        temperature_K = self.atoms.get_temperature()
        volume_Ang3_per_atom = self.atoms.get_volume() / num_atoms

        stress_GPa = self.atoms.get_stress(include_ideal_gas=True, voigt=False) / ase.units.GPa
        pxx_GPa = -stress_GPa[0, 0]
        pyy_GPa = -stress_GPa[1, 1]
        pzz_GPa = -stress_GPa[2, 2]
        pyz_GPa = -stress_GPa[1, 2]
        pzx_GPa = -stress_GPa[2, 0]
        pxy_GPa = -stress_GPa[0, 1]

        steps = self.dyn.nsteps
        time_ps = self.dyn.get_time() / (1000 * ase.units.fs)

        dat = (
            steps,
            time_ps,
            epot_eV_per_atom,
            ekin_eV_per_atom,
            etot_eV_per_atom,
            temperature_K,
            volume_Ang3_per_atom,
            pxx_GPa,
            pyy_GPa,
            pzz_GPa,
            pyz_GPa,
            pzx_GPa,
            pxy_GPa,
        )

        self.logfile.write(self.fmt % dat)
        self.logfile.flush()

def run_aniso_npt(pdamp_timestep: int, pressure_GPa: int):
    atoms = ase.build.bulk("Cu", crystalstructure='hcp', a=2.53, c=4.0).repeat((4, 4, 4))
    atoms.calc = asap3.EMT()

    temperature_K = 300
    rng = np.random.default_rng(0)
    MaxwellBoltzmannDistribution(atoms, temperature_K=temperature_K, force_temp=True, rng=rng)
    Stationary(atoms)

    timestep = 1.0 * ase.units.fs
    md = MTKNPT(
        atoms,
        timestep=timestep,
        temperature_K=temperature_K,
        pressure_au=float(pressure_GPa) * ase.units.GPa,
        tdamp=100 * timestep,
        pdamp=pdamp_timestep * timestep,
    )
    logfile = f'aniso-npt_pressure-{pressure_GPa}_pdamp-{pdamp_timestep}.log'
    md.attach(CustomMDLogger(md, atoms, logfile), interval=100)
    md.run(100_000)

Checklist

Merge request reports