run_ibverbs.py

Tuesday, 16 August 2011 07:10
run_ibverbs.py Version:3.1
Python submit script for NAMD with Sun Grid Engin (SGE) This Python script contains a template for consecutive NAMD runs. The first run has to be set up manually and all the resulting simulation files have to be put in the folder "run000" The script then resumes the first run from dir "run000" and moves all new files into "run001" then it resumes the simulation from "run001" and so on. There are some checks to make sure that the simulation terminated with the correct number of steps. If something fails, the script reruns the current run for X times and then terminates. I polished the script to be PEP8 conform (except the template part) hopfully I didn't introduce any typos or bugs.

GNU/GPL 2011-08-16 English Linux 24.51 KB 139
#!/usr/bin/python
# -*- coding: utf-8 -*-
#
#------------------------------------------------------------------------------
#run_ibverbs.py v3.1, Copyright Bjoern Olausson
#------------------------------------------------------------------------------
#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 2 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.
#
#To view the license visit
#http://www.gnu.org/licenses/old-licenses/gpl-2.0.html
#or write to
#Free Software Foundation, Inc.
#51 Franklin St, Fifth Floor
#Boston, MA
#02110-1301 USA
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------
#
#################################################
##QSUB Options
#################################################
#$ -S /usr/bin/python
 
# name of the job
#$ -N some_name
 
# Send Mail on b)Beginning, e)End or a)Aborted
#$ -M 
 This e-mail address is being protected from spambots. You need JavaScript enabled to view it
 
#$ -m ea
 
# combine error and output
#$ -j yes
 
# which queue to submit to: rack[0-4].q
#$ -q all.q
 
# parallel environment and number of procs
#$ -pe charm 168
 
# execute job from current working directory
#$ -cwd
 
# use current environment vars in the job (-V)
# -V
 
# Set environental vars for the run:
# -v MV2_ENABLE_AFFINITY=0
 
# Wait for job ID to complete
# -hold_jid ID
#################################################
#################################################
import os
import subprocess
import shutil
import string
import sys
import glob
import linecache
import re
 
# Save NAMD2 configuration to file NAMDconf
NAMDconf = "some_name.conf"
# Name of the simulation related files namd produces
# (.coor, .dcd, .xst, .xsc, .vel)
NAMDoutputname = "some_name"
# List of Additional files which should be moved by the function clean.
add_files = ['FFTW*']
# List of files to be deleted.
del_files = []
 
# Since I stumbled over some broken DCD files but according to the NAMD log,
# NAMD terminated successfully with full step count and the extended restart
# file also containt the right nuber of steps I added a file size check which
# compare the previous run with the current run.
# Set this to "False" if you make changes to anything between runs:
CHK_SIZE = False
 
# Frequency of writing out the coordinates to the DCD file
NAMDfreq = 500
# Number of steps to simulate 100000 Steps = 0.2ns (at 2fs Timestep)
# 500000 Steps = 1ns
NAMDsteps = 500000
# Number of simulations to perform
RUNS = 100
#  Restart the simulatin n times when it crases befor we terminate the job
MAX_RESTARTS = 3
# Dir to store initial simulation.
# This tool needs the following files to be in RUNDIR,
# (RUNDIR = RUN_PREFIX + RUN):
# ${NAMDoutputname}.restart.coor
# ${NAMDoutputname}.restart.vel
# ${NAMDoutputname}.restart.xsc
RUN_PREFIX = "run"
# Has to be a three digit number
# When restarting a simulation set this to the last dir.
RUN = "000"
 
# Get the some environmental variables
TMPDIR = os.environ.get("TMPDIR")
NCPU = os.environ.get("NSLOTS")
 
# Save the "machines" files
try:
    machin_file = open(TMPDIR + '/machines', 'r').read()
except Exception:
    print >>sys.stderr, "Faild to read %s/machines" % (TMPDIR)
    raise Exception("Faild to read %s/machines" % (TMPDIR))
 
#NAMD2 flavour to load
NAMD2 = ['/cvos/shared/apps/namd2/intel/64/2.8/ibverbs/namd2']
NAMD2opts = ['++nodelist',
    TMPDIR + '/machines',
    '+p' + str(NCPU),
    '+setcpuaffinity']
 
#CHARMRUN
CHARMRUN = ['/cvos/shared/apps/namd2/intel/64/2.8/ibverbs/charmrun']
CHARMopts = ['++remote-shell', 'ssh', '++scalable-start']
 
# Intel C Compiler lib dir
ICC_DIR = "/cvos/shared/apps/intel/Compiler/11.1/046/lib/intel64"
os.environ['LD_LIBRARY_PATH'] = ICC_DIR
 
#cvos-launcher
CVOSLAUNCHER = ['/cvos/shared/apps/sge/6.2/cvos/cvos-launcher']
 
#reduce DCD (True/False)
REDUCE = False
CATDCD = "/cvos/shared/apps/catdcd/4.0/catdcd"
STRIDE = "100"
 
# Assamble the directory
RUNDIR = RUN_PREFIX + str(RUN)
 
 
# This is the template for the NAMD configuration file.
# Change it to fit your needs.
# All strings starting with $ can be substituted $$ results
# in a normal $. # ${subs}test substitutes only the part between {}
def mkconf(NAMDfirsttimestep, RunDir):
    # Disable the caolvars restart file for the first run since
    # we do not provide on via the firstrun script.
    if RunDir == RUN_PREFIX + "000":
        ColvarRestartFile = "#colvarsInput"
    else:
        ColvarRestartFile = "colvarsInput"
 
    NAMDconf_form = '''
#############################################################
## JOB DESCRIPTION ##
#############################################################
 
# Production run
 
#############################################################
## ADJUSTABLE PARAMETERS ##
#############################################################
 
structure		step4_equilibration.xplor.psf
coordinates		step4_equilibration.pdb
set temp		294.15
set outputname		$NAMDoutputnameT
firsttimestep		$NAMDfirsttimestepT
 
# Continuing a job from the restart files
set inputname		$NAMDoutputnameT
 
binCoordinates		${RunDirT}/$$inputname.restart.coor
# remove the "temperature" entry if you use this!
binVelocities		${RunDirT}/$$inputname.restart.vel
extendedSystem		${RunDirT}/$$inputname.restart.xsc
 
 
#
# IMD Settings (can view sim in VMD)
#
 
# Specifies whether or not to listen for an IMD connection (def:off)
IMDon								off
 
# This is a free port number on the machine that node 0 is running on. This
# number will have to be entered into VMD.
IMDport							3000
 
# This allows coordinates to be sent less often, which may increase NAMD
# performance or be necessary due to a slow network.
IMDfreq							2
 
#  If no, NAMD will proceed with calculations whether a connection is present
# or not. If yes, NAMD will pause at startup until a connection is made, and
# pause when the connection is lost.
IMDwait							no
 
# If yes, NAMD will ignore any steering forces generated by VMD to allow a
# simulation to be monitored without the possibility of perturbing it.
IMDignore						yes
 
#############################################################
## SIMULATION PARAMETERS ##
#############################################################
 
# Input
paraTypeCharmm		on
parameters					par_all27_prot_lipid.prm
 
#
# GBIS parameters
#
 
# Turns on GBIS method in NAMD (def: off)
GBIS								off
 
# Defines the dielectric of the solvent, usually 78.5 or 80 (def: 78.5).
solventDielectric			80
 
# This offset shrinks the intrinsic radius of atoms (used only for
# calculating Born radius) to give better agreement with Poisson Boltzmann
# calculations. Most users should not change this parameter (def:0.09).
intrinsicRadiusOffset	0.09
 
# An ion concentration of 0 M represents distilled water. Increasing the ion
# concentration increases the electrostatic screening (def:0.2).
ionConcentration			0.003
 
# Parameter for calculating Born radii (def:0.3).
GBISDelta						1.0
 
# Parameter for calculating Born radii (def:0.8).
GBISBeta						0.8
 
# Parameter for calculating Born radii (def:4.85).
GBISGamma					4.85
 
# Cutoff used for calculating Born radius. Only atoms within this
# cutoff de-screen an atom. Though alphaCutoff can bet set to be larger
# or shorter than cutoff, since atom forces are more sensitive to changes in
# position than changes in Born radius, users should generally set alphaCutoff
# to be shorter than cutoff (def:15).
alphaCutoff					14
 
#
# Force-Field Parameters
#
 
# This parameter speciﬁes which pairs of bonded atoms should be excluded from
# non-bonded interactions.
exclude							scaled1-4
 
# Scaling factor for 1-4 interactions (def:1.0)
1-4scaling						1.0
 
# Local interaction distance common to both electrostatic and van der Waals
# calcu-lations
cutoff								12.0
 
# If switching is turned on, then smoothing functions are applied to both the
# electrostatics and van der Waals forces
switching						on
 
# Distance (≤ cutoff) at which to activate switching/splitting function for
# electrostatic and van der Waals calculations
switchdist						10.0
 
# Distance between pairs for inclusion in pair lists (≥ cutoff). A value of at
# least one greater than cutoff is recommended.
pairlistdist						16
 
# Regenerate pair lists x times per cycle (def:2)
pairlistsPerCycle			2
 
#
# Integrator Parameters
#
 
# The timestep size to use when integrating each step of the simulation. The
# value is speciﬁed in femtoseconds. (def:1.0)
timestep							2
 
# Controls if and how ShakeH is used (needed for 2fs steps)
# (values:all,water,none) (def:none)
rigidBonds						all
 
# This parameter speciﬁes how often short-range nonbonded interactions should
# be calculated. Setting nonbondedFreq between 1 and fullElectFrequency allows
# triple timestepping where, for example, one could evaluate bonded forces
# every 1 fs, short-range nonbonded forces every 2 fs, and long-range
# electrostatics every 4 fs.
nonbondedFreq			1
 
# Number of timesteps between full electrostatic evaluations (Positive integer
# factor of stepspercycle) (def:nonbondedFreq)
fullElectFrequency		1
 
# Number of timesteps in each cycle. Each cycle represents the number of
# timesteps between atom reassignments (def:20)
stepspercycle				10
 
#
# Constant Temperature Control
#
 
# Do langevin dynamics
langevin							on
 
# Damping coefficient (gamma) of 5/ps
langevinDamping			1
 
# Temperature to which atoms aﬀected by Langevin dynamics will be adjusted.
langevinTemp				$$temp
 
# If langevinDamping is set then setting langevinHydrogen to off will turn oﬀ
# Langevin dynamics for hydrogen atoms. This parameter has no eﬀect if Langevin
# coupling coeﬃcients are read from a PDB ﬁle.
langevinHydrogen		off
 
#
# Constant Pressure Control (variable volume)
#
 
# Is required in conjunction with rigidBonds (SHAKE).
useGroupPressure		yes
 
# Anisotropic cell ﬂuctuations
useFlexibleCell				yes
 
# When enabled, NAMD keeps the dimension of the unit cell in the x-y plane
# constant, while allowing ﬂuctuations along the z axis. This is not currently
# implemented in Berendsen’s method.
useConstantArea			yes
 
# When enabled, NAMD keeps the ratio of the unit cell in the x-y plane
# constant, while allowing ﬂuctuations along all axes. The useFlexibleCell
# option is required for this option.
useConstantRatio			no
 
# Use Langevin piston pressure control?
LangevinPiston				on
 
# Speciﬁes target pressure for Langevin piston method. A typical value would be
# 1.01325 bar, atmospheric pressure at sea level.
LangevinPistonTarget	1.01325
 
# Barostat oscillation time scale for Langevin piston method. Having a longer
# period results in more averaging over pressure measurements and hence slower
# ﬂuctuations in the cell volume. A reasonable choice for the piston period
# would be 200 fs.
LangevinPistonPeriod	200
 
# Speciﬁes barostat damping time scale for Langevin piston method
# (<= LangevinPistonPeriod)
LangevinPistonDecay	100
 
# Values: positive decimal Speciﬁes barostat noise temperature for Langevin
# piston method (should be ~temperature)
LangevinPistonTemp	$$temp
 
# Speciﬁes surface tension target. Must be used with useFlexibleCell and
# periodic boundary conditions (dyn/cm). The pressure speciﬁed in
# LangevinPistonTarget becomes the pressure along the z axis, and surface
# tension is applied in the x-y plane.
#SurfaceTensionTarget	30
 
#
# Periodic Boundary Conditions
#
cellBasisVector1			98.	0.	0.
cellBasisVector2			0.	98.	0.
cellBasisVector3			0.	0.	98.
cellOrigin						0.	0.	0.
 
# Wrap all coordinates around periodic boundaries
wrapAll							on
 
#
# PME (for full-system periodic electrostatics)
#
# Use particle mesh Ewald for electrostatics
PME									yes
 
# Maximum space between grid points (calculates PMEGridSize with small factors
# (2,3,5) in mind)
PMEGridSpacing			1.0
 
#PMEGridSizeX				96
#PMEGridSizeY				96
#PMEGridSizeZ				60
 
#############################################################
## CONSTRAINTS ##
#############################################################
 
# Use VMD to set the BETA vaule to the desired constraint
# (make sure all other BETA values are 0)
# set cons [atomselect top "protein"]
# $$cons set beta 10
# set sel [atomselect top all]
# $$sel writepdb k.pdb
 
# Specifies whether or not harmonic constraints are active.
# If it is set to off, then no harmonic constraints are computed. If it is set
# to on, then harmonic constraints are calculated using the values specified
# by the parameters consref, conskfile, conskcol, and consexp.
constraints						off
 
# PDB file to use for reference positions for harmonic constraints. Each atom
# that has an active constraint will be constrained about the position
# specified in this file.
consref							k.pdb
 
# PDB file to use for force constants for harmonic constraints.
conskfile							k.pdb
 
# Column of the PDB file to use for the harmonic constraint force constant.
# This parameter may specify any of the floating point fields of the PDB file,
# either X, Y, Z, occupancy, or beta-coupling (temperature-coupling).
# Regardless of which column is used, a value of 0 indicates that the atom
# should not be constrained. Otherwise, the value specified is used as the
# force constant for that atoms restraining potential.
conskcol							B
 
# The harmonic constraint energy function is multiplied by this parameter,
# making it possible to gradually turn off constraints during equilibration.
# This parameter is used only if constraints is set to on.
constraintScaling			1.0
 
 
#############################################################
## ABF ##
#############################################################
# Disable the collective variables module since it does not work with
# minimization
colvars							off
 
# Configuration file for the collective variables
colvarsConfig				colvars.in
 
# When continuing a previous simulation run, this file contains the current
# state of all collective variables and their biasing methods.
colvarsInput					$$inputname.restart.colvars.state	;
 
 
#############################################################
## Accelerated MD ##
#############################################################
# Specifies if accelerated MD is active. (def: off)
accelMD							off
 
# XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
accelMDdihe					on
 
# Specifies the threshold energy E in the aMD equations.
accelMDE						2000
 
# Specifies the acceleration a factor in the aMD equations.
accelMDalpha				160
 
# When accelMDdual is on, aMD switches to the dual boost mode. Two independent
# boost potentials will be applied: one to the dihedral potential that is
# controlled by the parameters accelMDE and accelMDalpha, and a second to the
# (Total - Dihedral) potential that is controlled by the accelMDTE and
# accelMDTalpha parameters described below. (def: off)
accelMDdual					off
 
# Specifies the threshold energy E used in the calculation of boost energy for
# the (Total - Dihedral) potential. This option is only available when
# accelMDdual is turned on.
accelMDTE						REAL
 
# Specifies the acceleration factor a used in the calculation of boost energy
# for the (Total - Dihedral) potential. This option is only available when
# accelMDdual is turned on.
accelMDTalpha				REAL
 
# Accelerated MD will only be performed when the current step is equal to or
# higher than accelMDFirstStep, and equal to or lower than accelMDLastStep.
# Otherwise regular MD will be performed. (def: 0)
accelMDFirstStep			0
 
# Accelerated MD will only be performed when the current step is equal to or
# higher than accelMDFirstStep, and equal to or lower than accelMDLastStep.
# Otherwise regular MD will be performed. Note that the accelMDLastStep
# parameter only has an effect when it is positive. When accelMDLastStep is set
# to zero (the default), aMD is "open-ended" and will be performed till the end
# of the simulation. (def: 0)
accelMDLastStep			0
 
# An aMD output line will be printed to the log file at the frequency specified
# by accelMDOutFreq. The aMD output will contain the boost potential (dV) at
# the current timestep, the average boost potential (dV AVG) since the last aMD
# output, and various potential energy values at the current timestep.
# The boost potential dV can be used to reconstruct the ensemble average
# described earlier. (def: 1)
accelMDOutFreq			100
 
 
#############################################################
## EXTRA PARAMETERS ##
#############################################################
# twoAwayX has the unique advantage of also improving the scalability of PME.
#twoAwayX					no
 
# Similar options twoAwayY and twoAwayZ also exist, and may be used
# in combination, but this greatly increases the number of compute objects.
#twoAwayY					no
#twoAwayZ					no
#margin							2.5
 
# Output
outputName					$$outputname
 
restartfreq						$NAMDfreqT
dcdfreq							$NAMDfreqT
xstFreq							$NAMDfreqT
outputEnergies				100
outputPressure				100
 
#############################################################
## EXECUTION SCRIPT ##
#############################################################
 
# 500000 = 1ns @ 2fs timestep
run									$NAMDstepsT
'''
 
    NAMDconf_template = string.Template(NAMDconf_form)
    # Define all substitutable strings in the template
    NAMDconf_file = NAMDconf_template.substitute(RunDirT=RunDir,
        NAMDoutputnameT=NAMDoutputname,
        NAMDstepsT=NAMDsteps,
        NAMDfirsttimestepT=NAMDfirsttimestep,
        NAMDfreqT=NAMDfreq,
        ColResFil=ColvarRestartFile)
    # Write out the config file
    open(NAMDconf, "w").write(NAMDconf_file)
 
NAMDerr = NAMDconf.replace('.conf', '.err.out')
NAMDout = NAMDconf.replace('.conf', '.out')
CATerr = "catdcd.err"
CATout = "catdcd.out"
 
# Just check for RUNDIR
if os.path.exists(RUNDIR) == False:
    MSG = "There is nothing to resume. Please run the very first simulation\
        manually and put the resulting files ( %s.restart.coor,\
        %s.restart.vel, %s.restart.xsc) in %s" % (NAMDoutputname,
            NAMDoutputname,
            NAMDoutputname,
            RUNDIR)
    raise Exception(MSG)
 
 
# This function moves all files generated by the current run to a new folder
def clean(DIR):
    os.mkdir(DIR, 0755)
    shutil.move(NAMDout, DIR)
    shutil.move(NAMDerr, DIR)
    shutil.move(NAMDconf, DIR)
    if REDUCE:
        shutil.move(CATerr, DIR)
        shutil.move(CATout, DIR)
    files = glob.glob(NAMDoutputname + "*")
    # We want to keep the jobfiles from SGE
    JOBFILES = re.compile('.+\.p?o[0-9]+')
    for element in files:
        if JOBFILES.match(element):
            shutil.copy(element, DIR)
        else:
            shutil.move(element, DIR)
    if len(add_files) > 0:
        for element in add_files:
            #If "*" is found in a parameter, glob it.
                if '*' in element:
                    for file in glob.glob(element):
                        try:
                            shutil.move(file, DIR)
                        except IOError:
                            MSG = "Faild to move %s to %s." % (file, DIR)
                            print >>sys.stderr, MSG
                else:
                        try:
                            shutil.move(element, DIR)
                        except IOError:
                            MSG = "Faild to move %s to %s." % (file, DIR)
                            print >>sys.stderr, MSG
    if len(del_files) > 0:
        for element in del_files:
            #If "*" is found in a parameter, glob it.
                if '*' in element:
                    for file in glob.glob(element):
                        try:
                            os.remove(file)
                        except OSError:
                            MSG = "Faild to delete %s to %s." % (file, DIR)
                            print >>sys.stderr, MSG
                else:
                        try:
                            os.remove(element)
                        except OSError:
                            MSG = "Faild to delete %s to %s." % (file, DIR)
                            print >>sys.stderr, MSG
 
 
# This function gets the last timestep from the .restart.xsc file
def GetLastTimeStep(DIR):
    if DIR == ".":
        XSC = NAMDoutputname + '.restart.xsc'
    else:
        XSC = DIR + "/" + NAMDoutputname + '.restart.xsc'
    try:
        open(XSC, "r")
    except IOError:
            raise Exception(XSC + " does not exist!")
    LastTimeStep = linecache.getline(XSC, 3).split()[0]
    linecache.clearcache()
    RestartTimestep = int(LastTimeStep)
    return LastTimeStep, RestartTimestep
 
 
# Just checking if the current simulations terminated successfully by comparing
# the current timestep with the expected timestep.
# If the current Steps are < the expected steps, we just rerun the simulation.
def check(DIR):
    CurTimeStep = int(GetLastTimeStep(".")[0])
    PreTimeStrp = int(GetLastTimeStep(DIR)[0])
    AllSteps = PreTimeStrp + NAMDsteps
    if CurTimeStep < AllSteps:
        MSG = \
        "CurrentTimeStep (%s) does not match the expected timesteps (%s)." \
        % (CurTimeStep, AllSteps)
        print >>sys.stderr, MSG
        check_ts = False
    else:
        check_ts = True
 
    # Since I stumbled over some broken DCD files but according to the
    # NAMD log, NAMD terminated successfully with full step count and
    # the extended restart file also containt the right nuber of steps,
    # I added these sanity checks:
    if CHK_SIZE:
        DCDsizeOLD = os.path.getsize(DIR + "/" + NAMDoutputname + ".dcd")
        DCDsizeNEW = os.path.getsize(NAMDoutputname + ".dcd")
        if DCDsizeNEW >= DCDsizeOLD:
            check_size = True
        else:
            MSG = \
            "DCD file-size (%i) does not match the previous runs size (%i)." \
            % (DCDsizeNEW, DCDsizeOLD)
            print >>sys.stderr, MSG
            check_size = False
    else:
        check_size = True
 
    if check_size and check_ts:
        return True
    else:
        return False
 
#INITIALrun = int(''.join(c for c in RUNDIR if c.isdigit()))
INITIALrun = int(RUN)
 
 
def RunNAMD():
    run_counter = 1
    OLDDIR = RUNDIR
    mkconf(GetLastTimeStep(RUNDIR)[1], RUNDIR)
    RESTARTS = 0
    while run_counter <= RUNS:
        NEWDIR = RUN_PREFIX + str(INITIALrun + run_counter).rjust(3, "0")
        # Check if the machines file is accessible
        if os.path.exists(TMPDIR):
            try:
                machin_file_check = open(TMPDIR + '/machines', 'r').read()
            except IOError:
                machin_file_regen = open(TMPDIR + '/machines', 'w')
                machin_file_regen.write(machin_file)
                machin_file_regen.close()
        else:
            os.mkdir(TMPDIR)
            machin_file_regen = open(TMPDIR + '/machines', 'w')
            machin_file_regen.write(machin_file)
            machin_file_regen.close()
 
        cmd = CHARMRUN + CHARMopts + NAMD2 + NAMD2opts
        cmd.append(NAMDconf)
        try:
            namd = subprocess.Popen(cmd, stderr=open(NAMDerr, "w"),
                stdout=open(NAMDout, "w"))
        except OSError:
            MSG = "Faild to either run cvos-launcher, mpirun_rsh or NAMD2"
            print >>sys.stderr, MSG
            raise Exception(MSG)
        else:
            namd.communicate()
 
        #Reduce the DCD if requested
        if REDUCE:
            cmd = [CATDCD,
                "-o", NAMDoutputname + "_reduced.dcd",
                "-stride", STRIDE,
                NAMDoutputname + ".dcd"]
            try:
                catdcd = subprocess.Popen(cmd,
                    stderr=open(CATerr, "w"),
                    stdout=open(CATout, "w"))
            except OSError:
                print >>sys.stderr, "Faild to run catdcd"
                raise Exception("Faild to run catdcd")
            else:
                catdcd.communicate()
 
        if check(OLDDIR) == True:
            RESTARTS = 0
            clean(NEWDIR)
            NewRestartTimestep = GetLastTimeStep(NEWDIR)[1]
            mkconf(NewRestartTimestep, NEWDIR)
            OLDDIR = NEWDIR
            run_counter = int(run_counter) + 1
        else:
            NEWDIR = NEWDIR + "_failed_" + str(RESTARTS)
            RESTARTS = RESTARTS + 1
            MSG = "Restarting simulation for %s. time (max %s times)"\
                % (RESTARTS, MAX_RESTARTS)
            print >>sys.stderr, MSG
            clean(NEWDIR)
            if RESTARTS == MAX_RESTARTS:
                MSG = "Reached maximum number of restarts (%s),\
                    terminating job now" % (MAX_RESTARTS)
                print >>sys.stderr, MSG
                raise Exception(MSG)
            else:
                NewRestartTimestep = GetLastTimeStep(OLDDIR)[1]
                mkconf(NewRestartTimestep, OLDDIR)
RunNAMD()