Zero stress tensors in NPT simulations

[lnnbig]

I use the ANI-1xnr potential, outlined in a recent publication (Nat. Chem., 2024. https://doi.org/10.1038/s41557-023-01427-3), alongside TorchANI, to explore O-H fluids at 1000 K and 1 GPa. The simulation box comprises 54 oxygen and 122 hydrogen atoms, with pressure and temperature regulated using the NPTBerendsen ensemble.

Throughout the MD simulation, I observed a clear adjustment in volume as the system approached equilibrium. However, despite this, the stress tensors remain essentially at zero. I use "get_volume" and "get_stress" functions in my Python script to inspect the volume and pressure.

I would greatly appreciate it if you could review these files and provide insights into any potential errors I may have overlooked. Thank you very much.

I've attached both the input and output files: test_ANI_1xnr.tar.gz. Below is a brief overview of the contents of the "load_from_neurochem.py" file.

###############################################################################

To begin with, let's first import the modules we will use:

import os
import torch
import torchani
import ase
import importlib_metadata

import numpy as np
import time
from ase.md.langevin import Langevin
from ase.io.trajectory import Trajectory
from ase import units
from ase.optimize.fire import FIRE as QuasiNewton
from ase.md.nvtberendsen import NVTBerendsen
from ase.md.nptberendsen import NPTBerendsen
from ase.md.npt import NPT
from ase.md import MDLogger
from ase.io import read, write
from ase.optimize import BFGS, LBFGS
from ase.io.vasp import write_vasp_xdatcar
from ase.build import make_supercell
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
from ase.cell import Cell
from ase.io.vasp import read_vasp
from ase.io.vasp import write_vasp
from ase.io.xyz import write_xyz
###############################################################################
try:
path = os.path.dirname(os.path.realpath(file))
except NameError:
path = os.getcwd()
const_file = os.path.join(path, './model/ani-1xnr/rHCNO-5.2R_32-3.5A_a8-4.params') # noqa: E501
consts = torchani.neurochem.Constants(const_file)
aev_computer = torchani.AEVComputer(**consts)
sae_file = os.path.join(path, './model/ani-1xnr/sae_linfit.dat') # noqa: E501
energy_shifter = torchani.neurochem.load_sae(sae_file)

###############################################################################

Now let's read a whole ensemble of models.

model_prefix = os.path.join(path, './model/ani-1xnr/train') # noqa: E501
ensemble = torchani.neurochem.load_model_ensemble(consts.species, model_prefix, 8) # noqa: E501

Or alternatively a single model.

model_dir = os.path.join(path, './model/ani-1xnr/train0/networks') # noqa: E501
model = torchani.neurochem.load_model(consts.species, model_dir)

###############################################################################
nnp2 = torchani.nn.Sequential(aev_computer, model, energy_shifter)
calculator2 = torchani.ase.Calculator(consts.species, nnp2)

###############################################################################
water_unitcell = read_vasp('POSCAR_H2O')
###############################################################################

Run Optimization

water_unitcell.set_calculator(calculator2)
start_time = time.time()
dyn = LBFGS(water_unitcell, trajectory='0.traj')
dyn.run(fmax=0.5)
write('temp_cell.xyz',water_unitcell)

Run MD

supercell = make_supercell(read('temp_cell.xyz'),[[1, 0, 0], [0, 1, 0], [0, 0, 1]])
supercell.set_calculator(calculator2)
MaxwellBoltzmannDistribution(supercell, temperature_K=2000)
#dyn = NVTBerendsen(supercell, 0.5 * units.fs, 1000.0, taut=1.01000units.fs, trajectory='1.traj')
dyn = NPTBerendsen(supercell,timestep=1units.fs,temperature_K=1000,taut=100units.fs,pressure_au=10000units.bar,
taup=1000units.fs,compressibility_au=4.6e-5/units.bar,trajectory='1.traj')
#dyn = NPT(supercell, timestep=1 * units.fs, temperature_K=1000, ttime=25 * units.fs, externalstress=0.06, pfactor= 30 * units.fs, mask=(1, 1, 1),trajectory='1.traj')

dyn.attach(MDLogger(dyn, supercell, 'md.log', header=True, stress=True, peratom=True, mode="w"), 1)

def printenergy(a=supercell): # store a reference to atoms in the definition.
"""Function to print the potential, kinetic and total energy."""
epot = a.get_potential_energy() / len(a)
ekin = a.get_kinetic_energy() / len(a)
vol = a.get_volume()
print('Energy per atom: Epot = %.3feV Ekin = %.3feV (T=%3.0fK) Etot = %.3feV Vol = %.3f' % (epot, ekin, ekin / (1.5 * units.kB), epot + ekin, vol))

dyn.attach(printenergy, interval=10)
printenergy()

dyn.run(100) # Do 1ps of MD
p=supercell.get_stress(include_ideal_gas=True)
print('Stress', p)

[WardLT]

I found a problem with the ASE calculator in torchani where it did not clear out results from the previous calculation when evaluating a new structure, which could explain the weird behavior here. It lead to the same stress being reported continually with on-the-fly for NPTBerendsen dynamics for me.

Solution: #650

It kind of works like this:

1. The MD timestep of NPTBerendsen ends with a force calculation, which results in the "forces" and "energy" of the calculator being updated but the previous stress calculation being kept. The new calculation also updates the atoms object used by ASE to determine whether a new calculation is needed.
2. The next MD step starts by requesting a stress calculation. ASE sees that the structure has not updated since the last invocation of the calculator, and retrieves the stress. That stress is whatever has been in stress since the last time stress was computed because the stress field was not cleared in 1.
3. The NPT dynamics keeps propagating with the same stress for all time. In your case, that might be the zero stress that the NPT run started with.