Hello everyone,
I’m currently working on a Molecular Dynamics simulation using the ASE (Atomic Simulation Environment) package in Python. However, due to a time limit imposed on the HPC system I’m using, the simulation is interrupted before completion. Now, I’m seeking guidance on how to effectively restart these simulations from where they left off.
Here’s a brief overview of what I’m working with:
from ase.io import read, write
import ase.units as units
from ase.optimize import BFGS
from ase import Atoms
from ase.io.trajectory import Trajectory
from ase.md.langevin import Langevin
from ase.md import velocitydistribution
import numpy as np
import torchani
import glob, os, os.path
atoms = read("box_water_ani1ccx_nvt_pbc.pdb")
vol = ((25.0))
atoms.set_cell((vol, vol, vol))
atoms.set_pbc(True)
calculator = torchani.models.ANI1ccx().ase()
atoms.set_calculator(calculator)
def printenergy_equil(a=atoms):
"""Function to print the potential, kinetic and total energy."""
epot = a.get_potential_energy() / len(a)
ekin = a.get_kinetic_energy() / len(a)
outfile = "energy_equil.out"
f = open(outfile, "a")
f.write('Energy per atom: Epot = %.3f eV Ekin = %.3f eV (T=%3.0f K) '
'Etot = %.3f eV \n' % (epot, ekin, ekin / (1.5 * units.kB), epot + ekin))
#print('Energy per atom: Epot = %.3f eV Ekin = %.3f eV (T=%3.0f K) '
# 'Etot = %.3f eV' % (epot, ekin, ekin / (1.5 * units.kB), epot + ekin))
def printenergy_prod(a=atoms):
"""Function to print the potential, kinetic and total energy."""
epot = a.get_potential_energy() / len(a)
ekin = a.get_kinetic_energy() / len(a)
outfile = "energy_prod.out"
f = open(outfile, "a")
f.write('Energy per atom: Epot = %.3f eV Ekin = %.3f eV (T=%3.0f K) '
'Etot = %.3f eV \n' % (epot, ekin, ekin / (1.5 * units.kB), epot + ekin))
def printxyz_equil(a=atoms):
checkpoint_out = "equil.pdb"
write(checkpoint_out, atoms, append=True)
def printxyz_prod(a=atoms):
"""Function to print out xyz coordinates"""
checkpoint_out = "prod.pdb"
write(checkpoint_out, atoms, append=True)
#os.chdir('..')
dyn = Langevin(atoms, 0.1 * units.fs, 300.15 * units.kB, 0.002)
velocitydistribution.MaxwellBoltzmannDistribution(atoms, 300.15 * units.kB)
print("Begin minimizing...")
opt = BFGS(atoms)
opt.run(fmax=1.0)
#print(len(atoms))
tag = 'water_box'
traj = Trajectory(tag + '.traj', 'w', atoms)
dyn.attach(traj.write, interval=1000)
dyn.attach(printenergy_equil, interval=100)
dyn.attach(printxyz_equil, interval=1000)
dyn.run(50000)
print("equilibration complete!")
dyn1 = Langevin(atoms, 0.1 * units.fs, units.kB * 300.15, 0.002)
tag = 'water_box_prod'
traj = Trajectory(tag + '.traj', 'w', atoms)
dyn1.attach(traj.write, interval=1000)
dyn1.attach(printenergy_prod, interval=100)
dyn1.attach(printxyz_prod, interval=1000)
dyn1.run(100000)
Specifically, I want to ensure that I can seamlessly continue the simulation without losing any data or compromising the integrity of the results.
If anyone has experience or insights into restarting ASE Molecular Dynamics simulations, especially using the provided code snippet, I would greatly appreciate your guidance. Any tips, suggestions, or pointers in the right direction would be incredibly helpful.