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)