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
Marco
"""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=io.read('H2_box.xyz')
atoms.set_cell((3, 3, 3))
atoms.center()
atoms.set_pbc(True)
traj = Trajectory(traj_file, 'w', atoms)
"""Equilibration"""
from ase.calculators.emt import EMT
atoms.set_calculator(EMT())
from ase.optimize import BFGS
dyn = BFGS(atoms)
dyn.attach(traj.write, interval=10)
dyn.run(fmax=1)
"""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)
dyn.run(100)
"""NVT"""
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)
dyn.run(1)
"""NPT"""
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)
dyn.run(1)
Please find the input file below:
16
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
