Replaying a Lammps Trajectory

EDIT: Turns out I just didn't google the correct words, LAMMPS does have one. It's the command rerun. Oops.

LAMMPS, the molecular dynamics engine, doesn't have a replay function allowing one to go back and calculate properties on a trajectory, especially forces and energies. This came up in my work and I got started on writing some new LAMMPS code to accomplish this. I wanted to be able to calculate forces with a different potential on a trajectory. Anyway, I realized you can accomplish this just by using the LAMMPS looping variables. No new code necessary. It can be done by using the read_dump command to load particular trajectory frame. A sample input script is below

	#the 3 is the number of frames in the trajectory to replay
	variable i loop 3
	#the stride is the time between frames
	variable stride equal 1
	#set-up run here, including reading data/restart file
	...
	...
	...
	...
	#Paranoid. Rebuild neighbor list every step.
	#Haven't figured out if read_dump forces a neighbor list rebuild.
	neigh_modify delay 0 every 1
	#now replay trajectory
	label replay
	#calculate the timestep
	variable t equal "(v_i-1)*v_stride"
	#read the file with read_dump
	#Reading a dcd requires the USER-MOLFILE package to be compiled
	read_dump traj.dcd $t x y z format molfile box yes reader dcd
	#Using XYZ reader instead
	#NOTE: If you use xyz instead, set stride to 1 since lammps doesn't read the timestep
	#read_dump traj.xyz $t x y z box no format xyz
	#set the timestep
	reset_timestep $t
	#calculate forces
	run 0
	#next frame
	next i
	jump SELF replay

view raw gistfile1.txt hosted with ❤ by GitHub