#!/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 your@mail.tld #$ -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 specifies 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 specified 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 specifies 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 affected by Langevin dynamics will be adjusted. langevinTemp $$temp # If langevinDamping is set then setting langevinHydrogen to off will turn off # Langevin dynamics for hydrogen atoms. This parameter has no effect if Langevin # coupling coefficients are read from a PDB file. langevinHydrogen off # # Constant Pressure Control (variable volume) # # Is required in conjunction with rigidBonds (SHAKE). useGroupPressure yes # Anisotropic cell fluctuations useFlexibleCell yes # When enabled, NAMD keeps the dimension of the unit cell in the x-y plane # constant, while allowing fluctuations 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 fluctuations along all axes. The useFlexibleCell # option is required for this option. useConstantRatio no # Use Langevin piston pressure control? LangevinPiston on # Specifies 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 # fluctuations in the cell volume. A reasonable choice for the piston period # would be 200 fs. LangevinPistonPeriod 200 # Specifies barostat damping time scale for Langevin piston method # (<= LangevinPistonPeriod) LangevinPistonDecay 100 # Values: positive decimal Specifies barostat noise temperature for Langevin # piston method (should be ~temperature) LangevinPistonTemp $$temp # Specifies surface tension target. Must be used with useFlexibleCell and # periodic boundary conditions (dyn/cm). The pressure specified 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()