I am testing ASE with a basic H2 MD simulation. I am reading the input file from a xyz file (generated with Packmol), I set pbc True . Then, I relax the system with BFGS and move to NVT and NPT equilibration. The cell just explodes.

I have tried to set initial velocities to a Maxwell distribution (300K) and have a VelocityVerlet run after BSGF but it also explodes. I have also tried reducing taut to 0.01 but I still cannot get a meaningful result.

Can you comment on how to improve my simulation? Why is the volume exploding? Shall I impose the external pressure before?

Thank you


"""Demonstrates molecular dynamics with constant energy."""

import ase
from ase import units

import numpy as np
from ase import Atoms
from ase.io.trajectory import Trajectory, TrajectoryReader

from ase import io

atoms.set_cell((3, 3, 3))

traj = Trajectory(traj_file, 'w', atoms)

from ase.calculators.emt import EMT 

from ase.optimize import BFGS
dyn = BFGS(atoms)

dyn.attach(traj.write, interval=10)

"""Guess velocities"""
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
from ase.md.verlet import VelocityVerlet

# Set the momenta corresponding to T=300K
MaxwellBoltzmannDistribution(atoms, temperature_K=300)

# We want to run MD with constant energy using the VelocityVerlet algorithm.
dyn = VelocityVerlet(atoms, 5 * units.fs)  # 5 fs time step.

dyn.attach(traj.write, interval=10)

from ase.md.nvtberendsen import NVTBerendsen

temperature = 300  # K
timestep = 1.0  # fs
total_time = 100  # ps
nsteps = int(total_time / timestep)

# Set up the integrator
dyn = VelocityVerlet(atoms, timestep)
dyn = NVTBerendsen(atoms, timestep, temperature_K=temperature, taut=0.01)
for step in range(nsteps):
    dyn.attach(traj.write, interval=10)

from ase.md.npt import NPT

pressure = 1.0  # bar
dyn = NPT(atoms, timestep, temperature_K=temperature, externalstress=pressure)
for step in range(nsteps):
    dyn.attach(traj.write, interval=10)

Please find the input file below:

    Energy:      -2.2240932
H          3.21205        0.01401       -0.00776
H          2.99965        0.66478        0.17258
H         -0.25392        0.17187       -0.08020
H          0.41341       -0.00174        0.07928
H         -0.29417        2.69561        3.24792
H         -0.34925        2.57862        2.55199
H          2.98091        0.28598        3.19910
H          3.09647       -0.25178        2.75376
H          3.26389        3.34259       -0.30528
H          3.35353        3.37297        0.39616
H          0.59907        2.89804        0.03153
H         -0.07742        2.95948       -0.16759
H          0.29545       -0.06128        3.14226
H         -0.24225       -0.21404        2.70812
H          2.37813        2.85750        2.82622
H          2.85583        3.19932        3.22112

2 Answers 2


Without further looking into the details, I can tell you straight-up why you cannot get any meaningful result. You are using the EMT calculator, which is a semi-empirical potential parameterized for transition metals. Non-metal elements like H, C, O are supported in the ASE EMT calculator only for testing purposes.

See here: the parameters for non-metal elements are completely made up so they are only for fun!

  • $\begingroup$ I see, thanks. That makes very much sense :) I will have a look into a better calculator. Any suggestion? $\endgroup$ Commented Jul 25, 2023 at 10:04

This was actually addressed on the ASE Element chat. The primary problem in this calculation is actually that packmol does not consider periodic boundary conditions. A simple solution involves simply running packmol and adding a buffer around the cell of empty space which is filled during equilibration. The EMT potential is also problematic, but the resulting cell as presented will include overlapping atoms which will give any potential issues.


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .