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.