LAMMPS FFS harness
#!/usr/bin/python
# -*- coding: utf-8 -*-
# Copyright (c) 2013 Kai Kratzer, Universität Stuttgart, ICP,
# Allmandring 3, 70569 Stuttgart, Germany; all rights
# reserved unless otherwise stated.
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston,
# MA 02111-1307 USA
########################################################################
# Harness script which uses LAMMPS as python lib, calculates
# a trajectory and returns the metadata
########################################################################
import sys
# lammps
from lammps import lammps
def save_snapshot(dumpfile):
# save here the simulation snapshot with the relevant variables and coordinates
# tricking: we have to know the timestep when reading in again
lmp.command('reset_timestep 1')
dumpcommand = 'dump cfp all custom 1 ' + dumpfile + ' id x y z vx vy vz'
lmp.command(dumpcommand)
# tricking: run 0 writes out dump but does not integrate
lmp.command('run 0')
def load_snapshot(dumpfile):
# read snapshot example
dumpcommand = 'read_dump ' + dumpfile + ' 1 x y z vx vy vz box yes replace yes'
lmp.command(dumpcommand)
# ready for new simulation run
lmp.command('reset_timestep 0')
def calc_rc():
# calculate and return reaction coordinate
def init_A():
# initialize and equilibrate system in state A
# create LAMMPS instance
lmp = lammps()
# parse command line arguments (alternatively, argparse can be used, but this is dynamic)
optval = {}
for eli in range(len(sys.argv)):
if (eli+1) % 2 == 0:
optval[re.sub('-','', sys.argv[eli],1)] = sys.argv[eli+1]
# extract variables
metafifo = optval['metadata_fifoname']
seed = int(optval['seed'])
act_lambda = int(optval['act_lambda'])
tmpdir = optval['tmpdir']
outdump = optval['back_fifoname']
clientname = optval['clientname']
timestamp = optval['timestamp']
if 'max_steps' in optval:
max_steps = int(optval['max_steps'])
else:
max_steps = 0
try:
uuid = optval['uuid']
except:
uuid = ''
harnessdir = re.sub('job_script','',sys.argv[0])
if act_lambda == 0:
escape_flux = True
rc_max = float(optval['B'])
else:
escape_flux = False
# TODO: Set the integrator, seed, etc, e.g.
lmp.command('fix 2 all langevin 1.0 1.0 2.0 ' + str(seed))
# init the system
init_A()
if escape_flux:
if rc < rc_min:
comefromok = True
else:
comefromok = False
# decide to load the initial config (then we are in A) or a dump.
# This is not necessary if an initial dump exists
if not "initial_config.dat" in str(sys.argv) and 'in_fifoname' in optval:
# loading snapshot from fifo
indump = optval['in_fifoname']
# this will set the last coordinates, velocities and the box size
load_snapshot(indump)
rc = calc_rc()
else:
# Not loading dump file, using initial state with equilibration.
eqc = 0
init_A()
rc = calc_rc()
while rc >= rc_next:
print "Equilibrating again"
init_A()
rc = calc_rc()
eqc += 1
max_rc = rc
all_steps = 0
in_progress = True
step_abort = False
# READY TO GO: This is an example for a script which checks the reaction coordinate itself
while in_progress:
# integrate
lmp.command('run ' + str(steps))
all_steps += steps
if rc > max_rc:
max_rc = rc
rc = calc_rc()
if escape_flux:
if rc >= rc_next and comefromok:
#print "### RC info for abort escape ", rc, rc_min, rc_next
#time.sleep(5)
in_progress = False
elif rc < rc_min and not comefromok:
comefromok = True
elif rc >= rc_max:
init_A(datfile, seed, dt)
rc = calc_rc()
comefromok = True
else:
if rc >= rc_next or rc <= rc_min:
in_progress = False
if max_steps > 0:
if all_steps >= max_steps:
print "Max steps reached!"
step_abort = True
in_progress = False
# build metadata
metadata = "\"rc\": %.14f, \"steps\": %i, \"time\": %.14f, \"max_lam\": %.14f" % (rc, all_steps, all_steps*dt,max_rc)
if step_abort:
metadata += ", \"step_abort\": True"
# save snapshot to fifo
save_snapshot(outdump)
# write to output metadata
fmeta=open(metafifo,'w')
fmeta.write("{" + metadata + "}\n" )
fmeta.close()
# destroy lammps instance
lmp.close()