#!/usr/bin/env python

########################################################################
#                                                                      #
# MMPBSA.py is a script for performing (M)olecular (M)echanics         #
# (P)oisson (B)oltzmann (S)urface (A)rea to find free energies of      #
# binding. Refer to the AMBER manual and/or relevant literature for a  #
# more thorough overview of the method. This implementation uses AMBER #
# executables to find energies using either Poisson Boltzmann or       #
# Generalized Born implicit solvent models of a complex of a receptor  #
# with a bound ligand. This script must be included with the python    #
# modules inputparse.py(c), utils.py(c), and alamdcrd.py(c). This      #
# script was written by Dwight McGee, Billy Miller III, and Jason      #
# Swails in Adrian Roitberg's research group at the Quantum Theory     #
# Project at the University of Florida.                                #
#                                                                      #
#       Last updated: 05/21/2010                                       #
#                                                                      #
########################################################################

########################## GPL LICENSE INFO ############################

#  Copyright (C) 2009  Dwight McGee, Billy Miller III, and Jason Swails

#  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.

#  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.


# import modules (1) exist on system and (2) are included with MMPBSA.py
import sys, os, time, math         # (1)
import inputparse, utils, alamdcrd # (2)

# import mpi4py, but check to make sure that it exists on the system
try:
   from mpi4py import MPI
except ImportError:
   print >> sys.stderr, 'Error: You must install the mpi4py package to use MMPBSA.py.MPI!'
   sys.exit()

########################################################################
#
#  BEGIN IMPORTANT NOTES: please add comments if you think there are
#  important comments to be made that will make understanding the code/
#  program flow easier.
#
########################################################################
 
#  utils.remove(debug,0) quits program as long as debug=-2
#  all functions return -1 in the case of an error
#  Temperature is hard-coded into ptraj and nab for entropy calculations, 
#   so you must change temperature in these places (and here)
#   and re-compile if you plan to calculate entropy not at 298.15 K
 
########################################################################
#
#  END Important notes
#
########################################################################

########################################################################
#
# BEGIN Initialize variables needed throughout script -- changing these
# default values may break the script (since the default values for
# some are necessary)
#
########################################################################

# MPI variables
communicator = MPI.COMM_WORLD
mpi_size     = communicator.Get_size()
rank         = communicator.Get_rank()
master       = rank == 0

# Integers (and float)
debug             = -2      # debug value for utils.remove()
x                 = 0       # top level loop counter
y                 = 0       # 2nd level loop counter
ligstart          = -1      # the residue in complex where ligand starts
numres_tot        = -1      # num residues in solvated prmtop
numres_com        = 0       # num residues in complex prmtop
numres_rec        = 0       # num residues in receptor prmtop
numres_lig        = 0       # num residues in ligand prmtop
numres_mut_com    = 0       # num residues in mutant complex prmtop
numres_mut_rec    = 0       # num residues in mutant receptor prmtop
numres_mut_lig    = 0       # num residues in mutant ligand prmtop
numframes         = 0       # num frames processed by ptraj 
numframesnmode    = 0       # num frames that were processed by nmode
numframesnmodechk = 0       # checks the above variable
isinerr           = 0       # variable used to test for errors
testdummy1        = 0       # checks for dummy com inpcrd file existence
testdummy2        = 0       # checks for dummy rec inpcrd file existence
testdummy3        = 0       # checks for dummy lig inpcrd file existence
mutant_residue    = 0       # index of mutant res in alanine scanning
topnumber         = 0       # the lowest frame processed by nmode
success           = 0       # catches return value of functions
entropy_holder    = 0.0     # holder for entropy val added to array
mpi_startframe    = 0       # the first frame this MPI thread works on
mpi_endframe      = 0       # the last frame this MPI thread works on
mpi_interval      = 0       # the how many frames each MPI thread works on
temp              = 298.15  # temperature for entropy unit conversion 
# Arrays
warnings          = []      # store all of the warnings for final output
INPUT             = []      # store all input variables described below
checkradiis       = []      # prmtops to be checked for consistent radii
hasmpiar          = []      # stores words in mpi_cmd
alascan           = []      # deltaGavg, deltaGstdev from ala. scanning
noalascan         = []      # deltaGavg, deltaGstdev not from ala. scan.
mdcrdholder       = []      # holder for mdcrds parsed from command line
suffix            = []      # stores suffices for frames used for nmode
entropyvalues     = []      # stores the values for the entropies
maskholder        = ['','',''] # 3-element array that holds lig & rec masks
# Logical variables
one_trajectory    = True    # Using single trajectory protocol (or TTP)
make_mdins        = False   # default value for make-mdins flag
use_mdins         = False   # default value for use-mdins flag
gbrun             = False   # variable if we are doing GB calc
pbrun             = False   # variable if we are doing PB calc
alarun            = False   # variable if we are doing ala scan
nmoderun          = False   # variable if we are doing nmode calc
decomprun         = False   # variable if we are doing decomp calc
overwrite         = False   # overwrite output file or not?
rewrite_output    = False   # varaible if we want to compile output
# Strings
cmdarg            = ''      # holds a command-line argument
inputfile_name    = 'none'  # default input file name -- there is none
ptraj             = ''      # absolute path for ptraj executable
sander            = ''      # absolute path for sander executable
sanderpb          = ''      # absolute path for sander to be used for pb
nmode             = ''      # absolute path for nmode executable
hasmpi            = ''      # absolute path of first word in mpi_cmd
mutstring         = ''      # mutation description (i.e. D152A)
# Files
centertraj        = ''      # process initial mdcrd 
ligandtraj        = ''      # extract ligand mdcrd (and nmode frames)
receptortraj      = ''      # extract receptor mdcrd (and nmode frames)
complextraj       = ''      # extract nmode complex snapshots
mutligandtraj     = ''      # extract mutant ligand mdcrd nmode frames
mutreceptortraj   = ''      # extract mutant receptor mdcrd nmode frames
mutcomplextraj    = ''      # extract mutant complex mdcrd nmode frames
getcrd            = ''      # extract top frame for dummy inpcrd
ptrajentropy      = ''      # ptraj script for quasi-harmonic NMA
grpfile           = ''      # write groupfile for multisander
avgpdb            = ''      # ptraj script to get avg. pdb 
# Timer variables 
tglobalstart      = time.time() #global timer
tglobalend        = 0.0
tchkmutprmstart   = 0.0     # check mutant prmtop timer
tchkmutprmend     = 0.0
tmdcrdmutstart    = 0.0     # mutate mdcrd for alanine scanning timer
tmdcrdmutend      = 0.0
tptrajnormstart   = 0.0     # normal ptraj system calls timer
tptrajnormend     = 0.0
tgbstart          = 0.0     # GB calculation timer
tgbend            = 0.0
tpbstart          = 0.0     # PB calculation timer
tpbend            = 0.0
tptrajentstart    = 0.0     # ptraj entropy calculation timer
tptrajentend      = 0.0
tnmodestart       = 0.0     # nmode entropy calculation timer
tnmodeend         = 0.0
toutputstart      = 0.0     # output timer for parsing/writing output
toutputend        = 0.0

########################################################################
#
#                END initializing variables 
#
########################################################################

# Open unbuffered stdout and stderr so redirected output gets written to
# stdout immediately

if master:
   sys.stdout = os.fdopen(sys.stdout.fileno(),'w',0)
   sys.stderr = os.fdopen(sys.stderr.fileno(),'w',0)
else: # non-masters will write everything to null device to suppress redundant output
   sys.stdout = open(os.devnull,'w',0)
   sys.stderr = open(os.devnull,'w',0)

#########################################################################
#
#                 BEGIN Explanation of input structure
#
########################################################################

# There is a single INPUT array whose indices are detailed below.
# This method allows all input variables to be passed at once to
# a subroutine to use them as they're needed. It is also the method
# I chose to facilitate a dynamic input file, such that an ordered
# input file is unnecessary.
#
# 0 ............. calc_type (see legend below)
# 1 ............. solvated_prmtop
# 2 ............. complex_prmtop
# 3 ............. receptor_prmtop
# 4 ............. ligand_prmtop
# 5 ............. startframe
# 6 ............. endframe
# 7 ............. interval
# 8 ............. mdcrd
# 9 ............. igb
# 10 ............ gbsa
# 11 ............ saltcon
# 12 ............ extdiel
# 13 ............ intdiel
# 14 ............ surften
# 15 ............ surfoff
# 16 ............ indi
# 17 ............ exdi
# 18 ............ scale
# 19 ............ linit
# 20 ............ prbrad
# 21 ............ istrng
# 22 ............ npopt
# 23 ............ cavity_surften
# 24 ............ cavity_offset
# 25 ............ fillratio
# 26 ............ nproc
# 27 ............ mpi_cmd
# 28 ............ radiopt
# 29 ............ verbose
# 30 ............ strip_mdcrd
# 31 ............ mutant_residue
# 32 ............ mutant_complex_prmtop
# 33 ............ mutant_receptor_prmtop
# 34 ............ mutant_ligand_prmtop
# 35 ............ receptor_mask
# 36 ............ ligand_mask
# 37 ............ sander_apbs
# 38 ............ mutant_only
# 39 ............ entropy
# 40 ............ dielc
# 41 ............ drms
# 42 ............ maxcyc
# 43 ............ nmstartframe (nmode)
# 44 ............ nmendframe (nmode)
# 45 ............ nminterval (nmode)
# 46 ............ receptor_mdcrd
# 47 ............ ligand_mdcrd
# 48 ............ solvated_receptor_prmtop
# 49 ............ solvated_ligand_prmtop
# 50 ............ ala_entropy
# 51 ............ keep_files
# 52 ............ strip_mask
# 53 ............ output_file
# 54 ............ idecomp
# 55 ............ print_res
# 56 ............ dec_verbose
# 57 ............ decompout
# 58 ............ nmode_igb
# 59 ............ nmode_istrng

# There are 15 possible types of calculations that can be done involving
# all combinations of gb, pb, alanine scanning, and nmode

# a 1 in the 1's position means gb is on
# a 1 in the 10's position means pb is on
# a 1 in the 100's position means ala is on
# a 1 in the 1000's position means nmode is on

# For example, 100 is alanine scanning only, 1 is gb only, 10 is pb only
# 111 is all three, 1000 is nmode only, etc. Note that a solvent method
# must be chosen for alanine scanning, so if only alanine scanning is 
# specified, INPUT[0] is adjusted to 101.

########################################################################
#
#                 END explanation of input structure
#
########################################################################

########################################################################
#
#                 BEGIN load default input values
#
########################################################################

INPUT.append('1')                       # calc_type
INPUT.append('solvated_prmtop')         # solvated_prmtop
INPUT.append('complex_prmtop')          # complex_prmtop
INPUT.append('receptor_prmtop')         # receptor_prmtop
INPUT.append('ligand_prmtop')           # ligand_prmtop
INPUT.append('1')                       # startframe
INPUT.append('1000000000')              # endframe
INPUT.append('1')                       # interval
INPUT.append([])                        # mdcrd
INPUT.append('5')                       # igb
INPUT.append('2')                       # gbsa
INPUT.append('0.0')                     # saltcon
INPUT.append('78.3')                    # extdiel
INPUT.append('1')                       # intdiel
INPUT.append('0.0072')                  # surften
INPUT.append('0.00')                    # surfoff
INPUT.append('1')                       # indi
INPUT.append('80')                      # exdi
INPUT.append('2.0')                     # scale
INPUT.append('1000')                    # linit
INPUT.append('1.4')                     # prbrad
INPUT.append('0.0')                     # istrng
INPUT.append('1')                       # npopt
INPUT.append('0.00542')                 # cavity_surften
INPUT.append('-1.008')                  # cavity_offset
INPUT.append('4')                       # fillratio
INPUT.append('3')                       # nproc
INPUT.append('none')                    # mpi_cmd
INPUT.append('0')                       # radiopt
INPUT.append('1')                       # verbose
INPUT.append('1')                       # strip_mdcrd
INPUT.append('1')                       # mutant_residue
INPUT.append('mutant_complex_prmtop')   # mutant_complex_prmtop
INPUT.append('none')                    # mutant_receptor_prmtop
INPUT.append('none')                    # mutant_ligand_prmtop
INPUT.append(' ')                       # receptor_mask
INPUT.append(' ')                       # ligand_mask
INPUT.append('0')                       # sander_apbs
INPUT.append('0')                       # mutant_only
INPUT.append('0')                       # entropy
INPUT.append('4')                       # dielc
INPUT.append('0.001')                   # drms
INPUT.append('10000')                   # maxcyc
INPUT.append('1')                       # nmstartframe
INPUT.append('1000000000')              # nmendframe
INPUT.append('1')                       # nminterval
INPUT.append([])                        # receptor mdcrd
INPUT.append([])                        # ligand mdcrd
INPUT.append('none')                    # solvated receptor prmtop
INPUT.append('none')                    # solvated ligand prmtop
INPUT.append('1')                       # ala_entropy
INPUT.append('1')                       # keep_files
INPUT.append(':WAT:Cl-:CIO:Cs+:IB:K+:Li+:MG2:Na+:Rb+') # strip_mask
INPUT.append('FINAL_RESULTS_MMPBSA.dat')# output_file
INPUT.append('0')                       # idecomp
INPUT.append('none')                    # print_res
INPUT.append('0')                       # dec_verbose
INPUT.append('FINAL_DECOMP_MMPBSA.dat') # decompout
INPUT.append('1')                       # nmode_igb
INPUT.append('0')                       # nmode_istrng

########################################################################
#
#                 END loading default input values
#
########################################################################


########################################################################
#
#                 BEGIN Read command line arguments
#
########################################################################

# Any call for help will be met with usage instructions
for x in range(len(sys.argv)):
   cmdarg = sys.argv[x].lower()
   if cmdarg == '-help' or cmdarg == '--help' or cmdarg == '-h':
      utils.printusage()
      sys.exit()

# All files created by this script are removed
for x in range(len(sys.argv)):
   cmdarg = sys.argv[x].lower()
   if '-clear' in cmdarg or '-clean' in cmdarg:
      if master:
         utils.remove(0,mpi_size)
      sys.exit() 

# print details about the parallelism
print >> sys.stdout, 'MMPBSA.py.MPI being run on {0} processors'.format(mpi_size)

# Read all input file names for all prmtops, the mdcrd, and the input 
# file. It also tests for out-of-bounds entries so the script crashes 
# gracefully, informing the user that they made a command line error.  
# If not all files are given on command line, the defaults are used. If
# any future files are necessary to add to the command line, it is easy
# to do here by following the convention that is already present.

try:
   for x in range(len(sys.argv)):
      if sys.argv[x] == '-i':            # input file flag
         inputfile_name = sys.argv[x+1]
      elif sys.argv[x] == '-sp':         # solvated_prmtop flag
         INPUT[1] = sys.argv[x+1]
      elif sys.argv[x] == '-cp':         # complex prmtop flag
         INPUT[2] = sys.argv[x+1]
      elif sys.argv[x] == '-rp':         # receptor_prmtop flag
         INPUT[3] = sys.argv[x+1]
      elif sys.argv[x] == '-lp':         # ligand_prmtop flag
         INPUT[4] = sys.argv[x+1]
      elif sys.argv[x] == '-mc':         # mutant_complex_prmtop flag
         INPUT[32] = sys.argv[x+1]
      elif sys.argv[x] == '-mr':         # mutant_receptor_prmtop flag
         INPUT[33] = sys.argv[x+1]
      elif sys.argv[x] == '-ml':         # mutant_ligand_prmtop flag
         INPUT[34] = sys.argv[x+1]
      elif sys.argv[x] == '-srp':        # solvated_receptor_prmtop flag
         INPUT[48] = sys.argv[x+1]
      elif sys.argv[x] == '-slp':        # solvated_ligand_prmtop flag
         INPUT[49] = sys.argv[x+1]
      elif sys.argv[x] == '-make-mdins': # create mdin files
         make_mdins = True
      elif sys.argv[x] == '-use-mdins':  # use existing mdin files flag
         use_mdins = True
      elif sys.argv[x] == '-rewrite-output': # just compile output
         rewrite_output = True
      elif sys.argv[x] == '-o':          # output_file
         INPUT[53] = sys.argv[x+1]
      elif sys.argv[x] == '-do':         # decomp output file
         INPUT[57] = sys.argv[x+1]       
      elif sys.argv[x] == '-O':          # overwrite flag
         overwrite = True
# Get arbitrary num of mdcrd files comma or white-space delimited
      elif sys.argv[x] == '-y':           # mdcrd flag
         x = x + 1
         while x < len(sys.argv) and not sys.argv[x].startswith("-"):
            mdcrdholder = sys.argv[x].split(",")
            for y in range(len(mdcrdholder)):
               if len(mdcrdholder[y].strip()) != 0:
                  INPUT[8].append(mdcrdholder[y].strip())
            x = x + 1
# while loop may have been triggered by a flag; rewind to read flag
         x = x - 1 
      elif sys.argv[x] == '-yr':          # receptor_mdcrd flag
         x = x + 1
         while x < len(sys.argv) and not sys.argv[x].startswith("-"):
            mdcrdholder = sys.argv[x].split(",")
            for y in range(len(mdcrdholder)):
               if len(mdcrdholder[y].strip()) != 0:
                  INPUT[46].append(mdcrdholder[y].strip())
            x = x + 1
# while loop may have been triggered by a flag; rewind to read flag
         x = x - 1
      elif sys.argv[x] == '-yl':          # ligand_mdcrd flag
         x = x + 1
         while x < len(sys.argv) and not sys.argv[x].startswith("-"):
            mdcrdholder = sys.argv[x].split(",")
            for y in range(len(mdcrdholder)):
               if len(mdcrdholder[y].strip()) != 0:
                  INPUT[47].append(mdcrdholder[y].strip())
            x = x + 1
# while loop may have been triggered by a flag; rewind to read flag
         x = x - 1 
# End getting mdcrd files
      elif sys.argv[x].startswith('-'):
         print >> sys.stderr, 'Error: Command line argument "{0}" not recognized!'.format(sys.argv[x])
         utils.printusage()
         sys.exit()

except IndexError: # IndexError occurs if flag is last item on cmd line
   print >> sys.stderr, 'Error: command line error!\n'
   utils.printusage()
   sys.exit()

# Give default mdcrd if no mdcrd was specified
if len(INPUT[8]) == 0:
   INPUT[8].append('mdcrd')

if make_mdins and use_mdins: # These flags are mutually exclusive
   print >> sys.stderr, 'Error: -use-mdins and -make-mdins cannot be used simultaneously!'
   utils.printusage()
   sys.exit()

if rewrite_output and not master: # output file is not parallel
   sys.exit()

if (make_mdins or use_mdins) and rewrite_output: # these options are mutually exclusive
   print >> sys.stderr, 'Error: -rewrite-output cannot be used with -make/use-mdins!'
   utils.printusage()
   sys.exit()

########################################################################
#
#                 END read command line arguments
#
########################################################################


########################################################################
#
# BEGIN Send input file and INPUT array to parser to be initialized
#
########################################################################

# If inputfile_name == 'none', then nothing was specified on command 
# line so only the default values will be used, which means only gb run

# inputparse.InputParse function returns a value that specifies the type
# of calculation to be run, namely INPUT[0] described above

if inputfile_name != 'none':
   INPUT[0] = inputparse.InputParse(inputfile_name, INPUT, warnings)
else:
# if make-mdins was specified, create all mdins if no input file given
   if make_mdins:
      INPUT[0] = 1111 

# inputparse returned error but already printed message; just exit.
if INPUT[0] == -1: 
   sys.exit()

########################################################################
#
#             END input parsing and INPUT initialization
#
########################################################################


########################################################################
#
#    BEGIN Transmit INPUT array to more descriptive variables to be used
#    throughout the rest of the program. This also allows me to check
#    that all input variables are of the correct type.
#
########################################################################
try:
   solvated_prmtop           =       INPUT[1]
   complex_prmtop            =       INPUT[2]
   receptor_prmtop           =       INPUT[3]
   ligand_prmtop             =       INPUT[4]
   startframe                =   int(INPUT[5])
   endframe                  =   int(INPUT[6])
   interval                  =   int(INPUT[7])
   mdcrd                     =       INPUT[8]
   igb                       =   int(INPUT[9])
   gbsa                      =   int(INPUT[10])
   saltcon                   = float(INPUT[11])
   extdiel                   = float(INPUT[12])
   intdiel                   = float(INPUT[13])
   surften                   = float(INPUT[14])
   surfoff                   = float(INPUT[15])
   indi                      = float(INPUT[16])
   exdi                      = float(INPUT[17])
   scale                     = float(INPUT[18])
   linit                     = float(INPUT[19])
   prbrad                    = float(INPUT[20])
   istrng                    = float(INPUT[21])
   npopt                     =   int(INPUT[22])
   cavity_surften            = float(INPUT[23])
   cavity_offset             = float(INPUT[24])
   fillratio                 = float(INPUT[25])
   nproc                     =   int(INPUT[26])
   mpi_cmd                   =       INPUT[27]
   radiopt                   = float(INPUT[28])
   verbose                   =   int(INPUT[29])
   strip_mdcrd               =   int(INPUT[30])
#  MUTANT_RESIDUE IS NOT AN INPUTTED VARIABLE ANYMORE
   mutant_complex_prmtop     =       INPUT[32]
   mutant_receptor_prmtop    =       INPUT[33]
   mutant_ligand_prmtop      =       INPUT[34]
   receptor_mask             =       INPUT[35]
   ligand_mask               =       INPUT[36]
   sander_apbs               =   int(INPUT[37])
   mutant_only               =   int(INPUT[38])
   entropy                   =   int(INPUT[39])
   dielc                     = float(INPUT[40])
   drms                      = float(INPUT[41])
   maxcyc                    =   int(INPUT[42])
   nmstartframe              =   int(INPUT[43])
   nmendframe                =   int(INPUT[44])
   nminterval                =   int(INPUT[45])
   receptor_mdcrd            =       INPUT[46]
   ligand_mdcrd              =       INPUT[47]
   solvated_receptor_prmtop  =       INPUT[48]
   solvated_ligand_prmtop    =       INPUT[49]
   ala_entropy               =   int(INPUT[50])
   keep_files                =   int(INPUT[51])
   strip_mask                =       INPUT[52]
   output_file               =       INPUT[53]
   idecomp                   =   int(INPUT[54])
   print_res                 =       INPUT[55]
   dec_verbose               =   int(INPUT[56])
   decompout                 =       INPUT[57]
   nmode_igb                 =   int(INPUT[58])
   nmode_istrng              = float(INPUT[59])
   
except ValueError:
   print >> sys.stderr, '\nInput Error: Check to make sure you have' + \
           ' proper input variable types (integer, string, float, etc.)'
   sys.exit()

# Calculate some variables needed from inputs provided

# invert the scale and put it back into INPUT as a string
scale = 1 / scale
INPUT[18] = str(scale)

# NOTE: the three mandatory ala variables are taken care of later in the
#       event that alanine scanning is asked for.

########################################################################
#
#             END variable initialization
#
########################################################################

# If we've made it this far, you're about to start writing files. Erase 
# all outputs and inputs from  previous calculations (assuming nobody in
# their right minds would create files starting with _MMPBSA_ and 
# FINAL_RESULTS_MMPBSA on their own).

if use_mdins and master:
   utils.remove(-3,mpi_size) # remove everything except the mdins.
elif not rewrite_output and master:
   utils.remove(0,mpi_size) # remove all files and back up old final results

# Before we start writing any files, we will check to see if the output
# file already exists. If it exists, print a message and quit in error

if not overwrite and not make_mdins:
   if utils.fileexists_noprint(output_file) != -1:
      print >> sys.stderr, "Error: Output file {0} already exists! Use -O to overwrite".format(output_file)
      sys.exit()
   if utils.fileexists_noprint(decompout) != -1:
      print >> sys.stderr, "Error: Decomp output file {0} already exists! Use -O flag to overwrite".\
               format(decompout)
      sys.exit()

########################################################################
#
# BEGIN Sets whether gb, pb, ala, nmode, and decomp will be run (0 no, 
# 1 yes) based on value of INPUT[0].  See above (near top) for legend.
#
########################################################################

# Detect which calculations are being requested by looking at the digits
# of INPUT[0].

if utils.digit(INPUT[0],4) == 1:
   if idecomp == 0:
      print >> sys.stderr, 'Error: You must specify idecomp in the &decomp namelist!'
      utils.remove(debug,0)
   decomprun = True
if not decomprun:                    # decomp requires different mdins
   if utils.digit(INPUT[0],0) == 1:
      if not use_mdins and master:
         utils.gbmdin(INPUT)
      elif use_mdins:
         if utils.fileexists('_MMPBSA_gb.mdin') == -1:
            utils.remove(debug,0)
      gbrun = True                  # turn gbrun on
   if utils.digit(INPUT[0],1) == 1:
      if not use_mdins and master: 
         utils.pbmdin(INPUT)
   #     utils.pbmdin_old(INPUT) # uncomment this for amber10 or older for amber pbsa only!
      elif use_mdins:
         if utils.fileexists('_MMPBSA_pb.mdin') == -1:
            utils.remove(debug,0)
      pbrun = True                  # turn pbrun on
   if utils.digit(INPUT[0],2) == 1: 
      alarun = True                 # turn alanine scanning on
   
   # For alanine scanning, set non-specified mutant prmtop files to the 
   # normal prmtop files, thereby assuming un-specified prmtop files are
   # unchanged in the mutant
   if alarun and mutant_receptor_prmtop == 'none':
      mutant_receptor_prmtop = receptor_prmtop
      INPUT[33] = mutant_receptor_prmtop
   if alarun and mutant_ligand_prmtop == 'none':
      mutant_ligand_prmtop = ligand_prmtop
      INPUT[34] = mutant_ligand_prmtop
   
   # If they requested -make-mdins, then quit now after all of the mdins 
   # have been created.
   if make_mdins:
      print >> sys.stdout, 'mdin Files Created! Exiting...'
      sys.exit()
else:
   if utils.digit(INPUT[0],0) == 1:
      gbrun = True
   if utils.digit(INPUT[0],1) == 1:
      pbrun = True
   if utils.digit(INPUT[0],2) == 1:
      alarun = True

if utils.digit(INPUT[0],3) == 1: 
   nmoderun = True               # turn nmode on

########################################################################
#
#          END setting gbrun, pbrun, alarun, and nmoderun
#
########################################################################


########################################################################
#
#          BEGIN Find and initialize required executables
#
########################################################################

# look for ptraj
if not rewrite_output:
   ptraj = utils.which('ptraj')  
   if ptraj == 'none': 
      print >> sys.stderr, "\nError: ptraj is needed for MMPBSA"
      utils.remove(debug,0)
   else:
      print >> sys.stdout, "ptraj found! Using " + ptraj
   
   # look for sander and sander.APBS if requested
   sander = utils.which('sander')
   if sander == 'none':
      print >> sys.stderr, "\nError: sander is needed for MMPBSA"
      utils.remove(debug,0)
   else:
      print >> sys.stdout, "sander found! Using " + sander

   if sander_apbs == 1:
      sanderpb = utils.which('sander.APBS')
      if sanderpb == 'none':
         print >> sys.stderr, 'Error: sander.APBS could not be found!' +  \
                  ' Either place sander.APBS in your path or remove'
         print >> sys.stderr, '       sander_apbs=1 from your input file and rerun.'
         utils.remove(debug,0)
      else:
         print >> sys.stdout, 'sander.APBS found! Using {0} for PB calculations'.format(sanderpb)
   else:
      sanderpb = sander 
   # The PB input file has changed from amber10 to amber11, so check
   # if sanderpb executable is in an amber10 directory. If it is, overwrite
   # _MMPBSA_pb.mdin with a copy for the older version of sander.
      if "amber10" in sanderpb.lower() or "amber9" in sanderpb.lower() and master:
         print >> sys.stderr, "Assuming " + sanderpb + " is part of\n" +  \
               "amber9 or amber10. Using old PB input file."
         os.system("rm -f _MMPBSA_pb.mdin")
         utils.pbmdin_old(INPUT)
   
   # look for nmode if requested
   if nmoderun:
      nmode = utils.which('mmpbsa_py_nabnmode')
      if nmode == 'none':
         print >> sys.stderr, 'Error: You specified an nmode calculation, but mmpbsa_py_nabnmode cannot be found!'
         utils.remove(debug,0)
      print >> sys.stdout, 'nmode program found! Using ' + nmode
   
########################################################################
#
#          END find and initialize required executables
#
########################################################################

   # Check to make sure the mdcrd file(s) exist
   for x in range(len(mdcrd)):
      if utils.fileexists(mdcrd[x]) == -1:
         utils.printusage()
         utils.remove(debug,0)

########################################################################
#
#   BEGIN Get numbers of residues from the topology files. Also checks
#   for the existence of all topology files and outputs errors for them.
#
########################################################################

if strip_mdcrd == 1:
   numres_tot = utils.resnum(solvated_prmtop)
numres_com = utils.resnum(complex_prmtop)
checkradiis.append(complex_prmtop)
numres_rec = utils.resnum(receptor_prmtop)
checkradiis.append(receptor_prmtop)
numres_lig = utils.resnum(ligand_prmtop)
checkradiis.append(ligand_prmtop)

if alarun:
   numres_mut_com = utils.resnum(mutant_complex_prmtop)
   checkradiis.append(mutant_complex_prmtop)
   numres_mut_rec = utils.resnum(mutant_receptor_prmtop)
   checkradiis.append(mutant_receptor_prmtop)
   numres_mut_lig = utils.resnum(mutant_ligand_prmtop)
   checkradiis.append(mutant_ligand_prmtop)

# test for any incompatibilities in the inputs
if utils.CheckIncomp(INPUT, numres_com) == -1:
   utils.remove(debug,0)

# set the default receptor and ligand masks
if ligand_mask != ' ' and receptor_mask == ' ':
   print >> sys.stderr, 'Warning: You only defined ligand_mask! This will be overridden by the default ligand_mask.'
   warnings.append('You only defined ligand_mask! This will be overridden by the default ligand_mask.')
   maskholder = utils.GetMasks(complex_prmtop, receptor_prmtop, ligand_prmtop, debug)

   ligand_mask   = maskholder[0]
   receptor_mask = maskholder[1]
   ligstart      = maskholder[2]

elif ligand_mask == ' ' and receptor_mask != ' ':
   print >> sys.stderr, 'Warning: You only defined receptor_mask! This will be overridden by the default receptor_mask.'
   warnings.append('You only defined receptor_mask! This will be overridden by the default ligand_mask.')
   maskholder = utils.GetMasks(complex_prmtop, receptor_prmtop, ligand_prmtop, debug)

   ligand_mask   = maskholder[0]
   receptor_mask = maskholder[1]
   ligstart      = maskholder[2]

elif ligand_mask == ' ' and receptor_mask == ' ':
   maskholder = utils.GetMasks(complex_prmtop, receptor_prmtop, ligand_prmtop, debug)

   ligand_mask   = maskholder[0]
   receptor_mask = maskholder[1]
   ligstart      = maskholder[2]

# Put the masks back into INPUT
INPUT[35] = receptor_mask
INPUT[36] = ligand_mask
# end set default receptor and ligand masks

# test for valid topology files (mutant prmtop check done later)
isinerr = 0
if numres_tot == -1 and strip_mdcrd == 1:
   isinerr = 1
if numres_com == -1 or numres_rec == -1 or numres_lig == -1:
   isinerr = 1

if isinerr == 1:
   utils.remove(debug,0)

# Check radii across different prmtops
if gbrun or pbrun:
   if utils.CheckRadiis(igb, checkradiis, warnings) == -1:
      utils.remove(debug,0)

if alarun and decomprun:
   print >> sys.stderr, "Error: Alanine scanning with decomposition is not supported!"
   utils.remove(debug,0)

# now that we have already checked the prmtops and figured out the mask
if decomprun:
   if use_mdins:
      if gbrun:
         if utils.fileexists("_MMPBSA_gb_decomp_com.mdin") + utils.fileexists("_MMPBSA_gb_decomp_rec.mdin") + \
               utils.fileexists("_MMPBSA_gb_decomp_lig.mdin") <= -1:
            utils.remove(debug,0)
      if pbrun:
         if utils.fileexists("_MMPBSA_pb_decomp_com.mdin") + utils.fileexists("_MMPBSA_pb_decomp_rec.mdin") + \
               utils.fileexists("_MMPBSA_pb_decomp_lig.mdin") <= -1:
            utils.remove(debug,0)
   else:
      if ligstart == -1:
         print >> sys.stderr, "WARNING: Running decomp with your own masks requires more care!"
         print >> sys.stderr, "         Consider removing your mask definitions and using defaults -or-\n"
         print >> sys.stderr, "Fix RES cards in decomp inputs for your system and re-run with the -use-mdins flag.\n"
         make_mdins = True
      if gbrun and master:
         utils.decomp_gbmdin(INPUT, ligstart)
      if pbrun and master:
         utils.decomp_pbmdin(INPUT, ligstart)

   if make_mdins:
      print >> sys.stdout, 'mdin Files Created! Exiting...'
      sys.exit()
      
# Now test the mutant topology files in more detail
if alarun:
   tchkmutprmstart = time.time()
   print >> sys.stdout, 'Checking mutant topology files...\n'
   if numres_mut_com == -1 or numres_mut_rec == -1 or numres_mut_lig == -1:
      print >> sys.stderr, 'Error: Invalid mutant topology file(s)!'
      utils.remove(debug,0)

# Check the validity of the mutant topology files
   isinerr = utils.CheckRecLigParms(receptor_prmtop,                   \
         mutant_receptor_prmtop, ligand_prmtop, mutant_ligand_prmtop)

   if isinerr == -1:
      print >> sys.stderr, 'Error: Invalid alanine scanning topology file(s)!'
      utils.remove(debug,0)
   tchkmutprmend = time.time()


########################################################################
#
#          END residue number acquisition and prmtop checks
#
########################################################################
 
########################################################################
#
# BEGIN Write prtaj input files to process the original mdcrd file.
# Complex, receptor, and ligand trajectory files are extracted via these
# ptraj calls. Also generate dummy inpcrd files for imin=5 sander runs.
#
########################################################################

if not rewrite_output:
   print >> sys.stdout, "\nPreparing trajectories with ptraj... "
   
   tptrajnormstart = time.time()
   
   if master:
      centertraj = open('_MMPBSA_cenptraj.in', 'w')
      for x in range(len(mdcrd)):
         centertraj.write("trajin {0} {1} {2} {3}\n".format(mdcrd[x],startframe,endframe,interval))
      if strip_mdcrd == 1:
         centertraj.write('strip ' + strip_mask + '\n')
         centertraj.write('center ' + receptor_mask + ' mass origin\n')
         centertraj.write('image origin center\n')
         centertraj.write('center :1-{0} mass origin\n'.format(numres_com))
         centertraj.write('image origin center\n')
      centertraj.write('rms first mass ' + receptor_mask + '\n')
      centertraj.write("average _MMPBSA_avgcomplex.pdb pdb\n")
      centertraj.write('trajout _MMPBSA_complex.mdcrd nobox\n')
      centertraj.close()

      ligandtraj = open('_MMPBSA_ligtraj.in', 'w')
      if len(ligand_mdcrd) != 0:
         one_trajectory = False
         for x in range(len(ligand_mdcrd)):
            if utils.fileexists(ligand_mdcrd[x]) == -1:
               utils.remove(debug,0)
            ligandtraj.write('trajin {0} {1} {2} {3}\n'.format(ligand_mdcrd[x], startframe, endframe, interval))
         if strip_mdcrd == 1:
            ligandtraj.write('strip {0}\n'.format(strip_mask))
      else:
         ligandtraj.write('trajin _MMPBSA_complex.mdcrd\n')
         ligandtraj.write('strip {0}\n'.format(receptor_mask))
      ligandtraj.write('trajout _MMPBSA_ligand.mdcrd nobox\n')
      ligandtraj.close()
   
   
      receptortraj = open('_MMPBSA_rectraj.in', 'w')
      if len(receptor_mdcrd) != 0:
         one_trajectory = False
         for x in range(len(receptor_mdcrd)):
            if utils.fileexists(receptor_mdcrd[x]) == -1:
               utils.remove(debug,0)
            receptortraj.write('trajin {0} {1} {2} {3}\n'.format(receptor_mdcrd[x],startframe,endframe,interval))
         if strip_mdcrd == 1:
            receptortraj.write('strip ' + strip_mask + '\n')
      else:
         receptortraj.write('trajin _MMPBSA_complex.mdcrd\n')
         receptortraj.write('strip '+ ligand_mask + '\n')
      receptortraj.write('trajout _MMPBSA_receptor.mdcrd nobox\n')
      receptortraj.close()

   
   # Extract only those frames specified in &nmode to do NMA on
   if nmoderun and master:
      complextraj = open('_MMPBSA_complex_nm.in', 'w')
      complextraj.write("trajin _MMPBSA_complex.mdcrd {0} {1} {2}\n".format(nmstartframe,nmendframe,nminterval))
      complextraj.write('trajout _MMPBSA_complex_nm.pdb pdb\n')
      complextraj.close()
      
      ligandtraj = open('_MMPBSA_ligtraj_nm.in', 'w')
      ligandtraj.write("trajin _MMPBSA_ligand.mdcrd {0} {1} {2}\n".format(nmstartframe,nmendframe,nminterval))
      ligandtraj.write('trajout _MMPBSA_ligand_nm.pdb pdb\n')
      ligandtraj.close()
      
      
      receptortraj = open('_MMPBSA_rectraj_nm.in', 'w')
      receptortraj.write("trajin _MMPBSA_receptor.mdcrd {0} {1} {2}\n".format(nmstartframe,nmendframe,nminterval))
      receptortraj.write('trajout _MMPBSA_receptor_nm.pdb pdb\n')
      receptortraj.close()
   
   if master:
      getcrd = open('_MMPBSA_complexinpcrd.in','w')
      getcrd.write("trajin _MMPBSA_complex.mdcrd 1 1 1\n")
      getcrd.write("trajout _MMPBSA_dummycomplex.inpcrd restart\n")
      getcrd.close()
   
   
      getcrd = open('_MMPBSA_receptorinpcrd.in','w')
      getcrd.write("""trajin _MMPBSA_receptor.mdcrd 1 1 1
trajout _MMPBSA_dummyreceptor.inpcrd restart
   """)
      getcrd.close()
   
   
      getcrd = open('_MMPBSA_ligandinpcrd.in','w')
      getcrd.write("""trajin _MMPBSA_ligand.mdcrd 1 1 1
trajout _MMPBSA_dummyligand.inpcrd restart
   """)
      getcrd.close()
      
      if entropy == 1 and mutant_only != 1:
         ptrajentropy = open('_MMPBSA_ptrajentropy.in','w')
         ptrajentropy.write('trajin _MMPBSA_complex.mdcrd\nreference _MMPBSA_avgcomplex.pdb\n\n')
         ptrajentropy.write('rms mass reference :1-{0}\n'.format(numres_com))
         ptrajentropy.write('matrix mwcovar name comp.matrix :1-{0}\n'.format(numres_com))
         ptrajentropy.write('analyze matrix comp.matrix out _MMPBSA_complex_entropy.out thermo reduce\n\n\n')
         ptrajentropy.write('rms mass reference ' + receptor_mask + '\n')
         ptrajentropy.write('matrix mwcovar name rec.matrix {0}\n'.format(receptor_mask))
         ptrajentropy.write('analyze matrix rec.matrix out _MMPBSA_receptor_entropy.out thermo reduce\n\n\n')
         ptrajentropy.write('rms mass reference {0}\n'.format(ligand_mask))
         ptrajentropy.write('matrix mwcovar name lig.matrix {0}\n'.format(ligand_mask))
         ptrajentropy.write('analyze matrix lig.matrix out _MMPBSA_ligand_entropy.out thermo reduce\n\n\n')
         ptrajentropy.close()
      
########################################################################
#
#          END ptraj input file creation for normal trajectory
#
########################################################################

# Generate the dry complex trajectory
      if strip_mdcrd == 1:
         os.system("{0} {1} _MMPBSA_cenptraj.in > _MMPBSA_ptraj1.out 2>&1".format(ptraj, solvated_prmtop))
         if len(ligand_mdcrd) > 0:
            os.system('{0} {1} _MMPBSA_ligtraj.in > _MMPBSA_ptraj2.out 2>&1'.format(ptraj, solvated_ligand_prmtop))
         else:
            os.system('{0} {1} _MMPBSA_ligtraj.in > _MMPBSA_ptraj2.out 2>&1'.format(ptraj, complex_prmtop))
         if len(receptor_mdcrd) > 0:
            os.system('{0} {1} _MMPBSA_rectraj.in > _MMPBSA_ptraj3.out 2>&1'.format(ptraj, solvated_receptor_prmtop))
         else:
            os.system('{0} {1} _MMPBSA_rectraj.in > _MMPBSA_ptraj3.out 2>&1'.format(ptraj, complex_prmtop))
      else:
         os.system('{0} {1} _MMPBSA_cenptraj.in > _MMPBSA_ptraj1.out 2>&1'.format(ptraj,complex_prmtop))
         if len(ligand_mdcrd) > 0:
            os.system('{0} {1} _MMPBSA_ligtraj.in > _MMPBSA_ptraj2.out 2>&1'.format(ptraj, ligand_prmtop))
         else:
            os.system('{0} {1} _MMPBSA_ligtraj.in > _MMPBSA_ptraj2.out 2>&1'.format(ptraj, complex_prmtop))
         if len(receptor_mdcrd) > 0:
            os.system('{0} {1} _MMPBSA_rectraj.in > _MMPBSA_ptraj3.out 2>&1'.format(ptraj, receptor_prmtop))
         else:
            os.system('{0} {1} _MMPBSA_rectraj.in > _MMPBSA_ptraj3.out 2>&1'.format(ptraj, complex_prmtop))
   # Generate snapshots for nmode if requested
      if nmoderun:
         os.system('{0} {1} _MMPBSA_complex_nm.in > _MMPBSA_ptraj12.out 2>&1'.format(ptraj,complex_prmtop))
         os.system('{0} {1} _MMPBSA_ligtraj_nm.in > _MMPBSA_ptraj13.out 2>&1'.format(ptraj,ligand_prmtop))
         os.system('{0} {1} _MMPBSA_rectraj_nm.in > _MMPBSA_ptraj14.out 2>&1'.format(ptraj,receptor_prmtop))

         if alarun: # now we have to get the mutant snapshots
            mutant_residue = alamdcrd.MutateMdcrd("_MMPBSA_complex.mdcrd","_MMPBSA_mutant_complex.mdcrd", 
                                                  complex_prmtop, mutant_complex_prmtop)

            if receptor_prmtop != mutant_receptor_prmtop:
               mutant_residue = alamdcrd.MutateMdcrd("_MMPBSA_receptor.mdcrd","_MMPBSA_mutant_receptor.mdcrd", 
                                                     receptor_prmtop, mutant_receptor_prmtop)
            else:
               os.system('cp _MMPBSA_receptor.mdcrd _MMPBSA_mutant_receptor.mdcrd')

            if ligand_prmtop != mutant_ligand_prmtop:
               mutant_residue = alamdcrd.MutateMdcrd("_MMPBSA_ligand.mdcrd","_MMPBSA_mutant_ligand.mdcrd", 
                                                     ligand_prmtop, mutant_ligand_prmtop)
            else:
               os.system('cp _MMPBSA_ligand.mdcrd _MMPBSA_mutant_ligand.mdcrd')

            mutcomplextraj = open('_MMPBSA_mutant_comtraj_nm.in','w')
            mutcomplextraj.write("trajin _MMPBSA_mutant_complex.mdcrd {0} {1} {2}\n".format(
                              nmstartframe, nmendframe, nminterval))
            mutcomplextraj.write('trajout _MMPBSA_mutant_complex_nm.pdb pdb\n')
            mutcomplextraj.close()
   
            mutligandtraj = open('_MMPBSA_mutant_ligtraj_nm.in', 'w')
            mutligandtraj.write("trajin _MMPBSA_mutant_ligand.mdcrd {0} {1} {2}\n".format(
                       nmstartframe, nmendframe, nminterval))
            mutligandtraj.write('trajout _MMPBSA_mutant_ligand_nm.pdb pdb\n')
            mutligandtraj.close()
   
            mutreceptortraj = open('_MMPBSA_mutant_rectraj_nm.in', 'w')
            mutreceptortraj.write("trajin _MMPBSA_mutant_receptor.mdcrd {0} {1} {2}\n".format(
                         nmstartframe, nmendframe, nminterval))
            mutreceptortraj.write('trajout _MMPBSA_mutant_receptor_nm.pdb pdb\n')
            mutreceptortraj.close()

            os.system('{0} {1} _MMPBSA_mutant_comtraj_nm.in > _MMPBSA_ptraj15.out 2>&1'.format(ptraj, mutant_complex_prmtop))
            os.system('{0} {1} _MMPBSA_mutant_ligtraj_nm.in > _MMPBSA_ptraj16.out 2>&1'.format(ptraj, mutant_ligand_prmtop))
            os.system('{0} {1} _MMPBSA_mutant_rectraj_nm.in > _MMPBSA_ptraj17.out 2>&1'.format(ptraj, mutant_receptor_prmtop))


########################################################################
#
#       BEGIN find the number of frames processed by ptraj
#
########################################################################
communicator.Barrier() # let all processors catch up
# now we quit if creation of decomp input files failed
if decomprun and gbrun and (utils.fileexists_noprint('_MMPBSA_gb_decomp_com.mdin') + \
     utils.fileexists_noprint('_MMPBSA_gb_decomp_rec.mdin') + utils.fileexists_noprint('_MMPBSA_gb_decomp_rec.mdin') < 0):
   utils.remove(debug,0)
if decomprun and pbrun and (utils.fileexists_noprint('_MMPBSA_pb_decomp_com.mdin') + \
     utils.fileexists_noprint('_MMPBSA_pb_decomp_rec.mdin') + utils.fileexists_noprint('_MMPBSA_pb_decomp_rec.mdin') < 0):
   utils.remove(debug,0)

# check to make sure we're not trying to do decomp with amber10 or older
if ("amber10" in sanderpb.lower() or "amber9" in sanderpb.lower() or "amber10" in sander.lower() \
     or "amber9" in sander.lower()) and decomprun:
   print >> sys.stderr, 'Error: Energy decomposition only works with Amber11 or later for MMPBSA.py!'
   utils.remove(debug,0)

numframes = int(utils.PtrajNumFrames('_MMPBSA_ptraj1.out'))
if numframes <= 0:
   print >> sys.stderr, 'Error: Either no frames were read in from original mdcrd or there was a ptraj error!'
   utils.remove(debug,0)
print >> sys.stdout, '{0} frames were read in and processed by ptraj for use in calculation.'.format(numframes)
if numframes == 1 and entropy == 1: # quasi-harmonic approximation cannot be done with only 1 frame
   print >> sys.stderr, 'Error: Quasi-harmonic approximation cannot be done with only 1 frame!'
   utils.remove(debug,0)

if nmoderun:
   if numframes < nmendframe:
      topnumber = numframes
   else:
      topnumber = nmendframe
   numframesnmode = int(math.floor(float(topnumber - nmstartframe) / float(nminterval))) + 1
   numframesnmodechk = int(utils.PtrajNumFrames('_MMPBSA_ptraj12.out'))

   if numframesnmode < 0 and numframesnmodechk < 0:
      print >> sys.stderr, 'Error: No frames were read in for nmode analysis!'
      utils.remove(debug)
   if numframesnmode != numframesnmodechk:
      print >> sys.stderr, 'Warning: Number of frames processed for' + \
               ' nmode differs from what is expected! {0} vs {1}'.format(numframesnmode, numframesnmodechk)
      warnings.append('Number of frames processed for nmode differs from what is expected!')

   print >> sys.stdout, "Processing {0} frames with normal mode analysis.".format(numframesnmode)

   for x in range(numframesnmode):
      suffix.append(str(nmstartframe + x * nminterval))

   if len(suffix) < mpi_size:
      print >> sys.stderr, "Error: You cannot have more MPI threads than " + \
            "nmode frames!"
      utils.remove(debug,0)

   mpi_interval = numframesnmode/mpi_size # both are integers, so this is "floor"
   mpi_extras = numframesnmode - mpi_interval * mpi_size
   if rank < mpi_extras:
      mpi_startframe = rank * mpi_interval + rank + 1
      mpi_endframe = mpi_interval * (rank + 1) + rank + 1
   else:
      mpi_startframe = rank * mpi_interval + mpi_extras + 1
      mpi_endframe = mpi_interval * (rank + 1) + mpi_extras

   suffix = suffix[mpi_startframe-1:mpi_endframe]

########################################################################
#
#       END find the number of frames processed by ptraj
#
########################################################################

########################################################################
#
#    BEGIN determine which threads handles which frames
#
########################################################################

if numframes < mpi_size: # don't allow more threads than frames
   print >> sys.stderr, "Error: You cannot have more MPI threads than " +\
         "frames!"
   utils.remove(debug,0)

mpi_interval = numframes/mpi_size # both are integers, so this is "floor"
mpi_extras = numframes - mpi_interval * mpi_size
if rank < mpi_extras:
   mpi_startframe = rank * mpi_interval + rank + 1
   mpi_endframe = mpi_interval * (rank + 1) + rank + 1
else:
   mpi_startframe = rank * mpi_interval + mpi_extras + 1
   mpi_endframe = mpi_interval * (rank + 1) + mpi_extras

########################################################################
#
#  BEGIN ptraj input file creation for each MPI thread. First calculate
#  which frames each processor will be working on
#
########################################################################

if not rewrite_output:
   centertraj = open('_MMPBSA_comtraj.in.{0}'.format(rank), 'w')
   centertraj.write("trajin _MMPBSA_complex.mdcrd {0} {1} 1\n".format(
        mpi_startframe,mpi_endframe))
   centertraj.write("trajout _MMPBSA_complex.mdcrd.{0}".format(rank))
   centertraj.close()
   
   ligandtraj = open('_MMPBSA_ligtraj.in.{0}'.format(rank), 'w')
   ligandtraj.write("trajin _MMPBSA_ligand.mdcrd {0} {1} 1\n".format(
        mpi_startframe,mpi_endframe))
   ligandtraj.write("trajout _MMPBSA_ligand.mdcrd.{0}".format(rank))
   ligandtraj.close()

   receptortraj = open('_MMPBSA_rectraj.in.{0}'.format(rank), 'w')
   receptortraj.write("trajin _MMPBSA_receptor.mdcrd {0} {1} 1\n".format(
        mpi_startframe,mpi_endframe))
   receptortraj.write("trajout _MMPBSA_receptor.mdcrd.{0}".format(rank))
   receptortraj.close()

########################################################################
#
# BEGIN ptraj system calls for normal trajectories. Also calls ptraj
# for the quasiharmonic approximation. (That takes awhile).
#
########################################################################

   # Generate the dry ligand trajectory
   os.system('{0} {1} _MMPBSA_comtraj.in.{2} > _MMPBSA_ptraj1.out.{2} 2>&1'.\
            format(ptraj, complex_prmtop, rank))
   os.system(ptraj + ' ' + ligand_prmtop + 
        ' _MMPBSA_ligtraj.in.{0} > _MMPBSA_ptraj2.out.{0} 2>&1'.format(rank))
   # Generate the dry receptor trajectory
   os.system(ptraj + ' ' + receptor_prmtop +
        ' _MMPBSA_rectraj.in.{0} > _MMPBSA_ptraj3.out.{0} 2>&1'.format(rank))
   
   # Generate the dummy input coordinate files for sander
   if master:
      os.system(ptraj + ' ' + complex_prmtop + ' _MMPBSA_complexinpcrd.in >' \
                + ' _MMPBSA_ptraj4.out 2>&1')
      os.system(ptraj + ' ' + receptor_prmtop + ' _MMPBSA_receptorinpcrd.in' \
                + ' > _MMPBSA_ptraj5.out 2>&1')
      os.system(ptraj + ' ' + ligand_prmtop + ' _MMPBSA_ligandinpcrd.in > ' +\
                '_MMPBSA_ptraj6.out 2>&1')
   
   # QUASI-HARMONIC APPROXIMATION PTRAJ CALL.
   if entropy == 1 and mutant_only != 1 and master:
      print >> sys.stdout, 'Beginning quasi-harmonic entropy calculation' \
               + ' with ptraj...'
      tptrajentstart = time.time()
      os.system(ptraj + ' ' + complex_prmtop + ' _MMPBSA_ptrajentropy.in' \
                + ' > _MMPBSA_ptraj_entropy.out')
      tptrajentend = time.time()
   
   # check to make sure that ptraj ran successfully. Easiest way of doing 
   # this is to see if it made the dummy inpcrd files

   communicator.Barrier() # let all threads catch up
   
   testdummy1 = utils.fileexists_noprint('_MMPBSA_dummycomplex.inpcrd.1')
   testdummy2 = utils.fileexists_noprint('_MMPBSA_dummyreceptor.inpcrd.1')
   testdummy3 = utils.fileexists_noprint('_MMPBSA_dummyligand.inpcrd.1')
   
   if testdummy1 == -1 or testdummy2 == -1 or testdummy3 == -1:
      print >> sys.stderr, 'Error! Ptraj failed. Check coordinate and topology files for the complex.'
      utils.remove(debug,0)

   communicator.Barrier() # prevent master from moving files before other threads test their existence

   # remove the .1 from the end of the dummy inpcrd files
   if master:
      os.system('mv _MMPBSA_dummycomplex.inpcrd.1 _MMPBSA_dummycomplex.inpcrd') 
      os.system('mv _MMPBSA_dummyreceptor.inpcrd.1 _MMPBSA_dummyreceptor.inpcrd')
      os.system('mv _MMPBSA_dummyligand.inpcrd.1 _MMPBSA_dummyligand.inpcrd')
      tptrajnormend = time.time()
   
########################################################################
#
#          END ptraj system calls for normal trajectories
#
########################################################################

# BEGIN check for the existence of files for rewrite-outputs

if rewrite_output:
   if utils.fileexists("_MMPBSA_ptraj1.out") == -1:
      utils.remove(debug,0)
   if nmoderun:
      if not mutant_only and (utils.fileexists("_MMPBSA_ptraj12.out") + \
           utils.fileexists("_MMPBSA_complex_nm.out") + utils.fileexists("_MMPBSA_receptor_nm.out") + \
           utils.fileexists("_MMPBSA_ligand_nm.out")) <= -1:
         utils.remove(debug,0)
      if alarun and (utils.fileexists("_MMPBSA_mutant_complex_nm.out") + \
          utils.fileexists("_MMPBSA_mutant_receptor_nm.out") + utils.fileexists("_MMPBSA_mutant_ligand_nm.out")) <= -1:
         utils.remove(debug,0)
   if gbrun:
      if not mutant_only and (utils.fileexists("_MMPBSA_complex_gb.mdout") + \
           utils.fileexists("_MMPBSA_receptor_gb.mdout") + utils.fileexists("_MMPBSA_ligand_gb.mdout")) <= -1:
         utils.remove(debug,0)
      if alarun and (utils.fileexists("_MMPBSA_mutant_complex_gb.mdout") + \
           utils.fileexists("_MMPBSA_mutant_receptor_gb.mdout") + \
           utils.fileexists("_MMPBSA_mutant_ligand_gb.mdout")) <= -1:
         utils.remove(debug,0)
   if pbrun:
      if not mutant_only and (utils.fileexists("_MMPBSA_complex_pb.mdout") + \
           utils.fileexists("_MMPBSA_receptor_pb.mdout") + utils.fileexists("_MMPBSA_ligand_pb.mdout")) <= -1:
         utils.remove(debug,0)
      if alarun and (utils.fileexists("_MMPBSA_mutant_complex_pb.mdout") + \
           utils.fileexists("_MMPBSA_mutant_receptor_pb.mdout") + \
           utils.fileexists("_MMPBSA_mutant_ligand_pb.mdout") ) <= -1:
         utils.remove(debug,0)
   if entropy == 1:
      if not mutant_only and utils.fileexists("_MMPBSA_ptraj_entropy.out") == -1:
         utils.remove(debug,0)
      if alarun and utils.fileexists("_MMPBSA_mutant_ptraj_entropy.out") == -1:
         utils.remove(debug,0)
   if len(receptor_mdcrd) + len(ligand_mdcrd) > 0:
      one_trajectory = False
   
# Compile mdouts at this point
if rewrite_output:
   print >> sys.stdout, "Writing output file..."
   utils.PrintFinalResults(output_file, inputfile_name, strip_mdcrd, solvated_prmtop, \
        complex_prmtop, receptor_mdcrd, solvated_receptor_prmtop, ligand_mdcrd,           \
        solvated_ligand_prmtop, receptor_prmtop, ligand_prmtop, mdcrd, mutant_complex_prmtop, \
        mutant_receptor_prmtop, mutant_ligand_prmtop, maskholder, numres_lig, numframes,  \
        nmoderun, numframesnmode, pbrun, sanderpb, entropy, mutant_only, mutstring, warnings, \
        keep_files, ala_entropy, alarun, gbrun, debug, one_trajectory, verbose, sander_apbs,  \
        decompout, idecomp, dec_verbose, ligstart, decomprun, surften, cavity_surften, temp)
   print >> sys.stdout, "\n\nMMPBSA Finished. Thank you for using. " +    \
            "Please send any bugs/suggestions/comments to " +             \
            "amber@ambermd.org"
   utils.remove(keep_files,mpi_size)
   sys.exit()

########################################################################
#
#  BEGIN alanine scanning setup: mutate mdcrd, check topology files
#  for errors and consistency, and set up/complete ptraj calls
#
########################################################################

if alarun:
   tmdcrdmutstart = time.time()
   print >> sys.stdout, "\nMutating complex mdcrd for alanine " +      \
            "scanning...\n"
   mutant_residue = alamdcrd.MutateMdcrd("_MMPBSA_complex.mdcrd.{0}".\
           format(rank),"_MMPBSA_mutant_complex.mdcrd.{0}".format(rank), 
           complex_prmtop,mutant_complex_prmtop)
# mutant_residue has index 31 in INPUT array, but this is where it is 
# determined -- now add it to INPUT[31] as a string
   INPUT[31] = str(mutant_residue)

   if mutant_residue == -1:
      utils.remove(debug,0)

   if ligand_prmtop == mutant_ligand_prmtop:
      os.system("cp _MMPBSA_ligand.mdcrd.{0} _MMPBSA_mutant_ligand.mdcrd.{0}".\
            format(rank))
   else:
      alamdcrd.MutateMdcrd("_MMPBSA_ligand.mdcrd.{0}".format(rank),
       "_MMPBSA_mutant_ligand.mdcrd.{0}".format(rank),ligand_prmtop,
       mutant_ligand_prmtop)

   if receptor_prmtop == mutant_receptor_prmtop:
      os.system("cp _MMPBSA_receptor.mdcrd.{0} _MMPBSA_mutant_receptor.mdcrd.{0}".\
           format(rank))
   else:
      alamdcrd.MutateMdcrd("_MMPBSA_receptor.mdcrd.{0}".format(rank),
                    "_MMPBSA_mutant_receptor.mdcrd.{0}".format(rank),
                    receptor_prmtop,mutant_receptor_prmtop)

   if master:
      getcrd = open('_MMPBSA_mutant_complexinpcrd.in','w')
      getcrd.write('trajin _MMPBSA_mutant_complex.mdcrd.0 1 1 1\n')
      getcrd.write('trajout _MMPBSA_mutant_dummycomplex.inpcrd restart\n')
      getcrd.close()
   
      getcrd = open('_MMPBSA_mutant_receptorinpcrd.in','w')
      getcrd.write('trajin _MMPBSA_mutant_receptor.mdcrd.0 1 1 1\n')
      getcrd.write('trajout _MMPBSA_mutant_dummyreceptor.inpcrd restart\n')
      getcrd.close()
   
      getcrd = open('_MMPBSA_mutant_ligandinpcrd.in','w')
      getcrd.write('trajin _MMPBSA_mutant_ligand.mdcrd.0 1 1 1\n')
      getcrd.write('trajout _MMPBSA_mutant_dummyligand.inpcrd restart\n')
      getcrd.close()

      if ala_entropy == 1 and entropy == 1:
         avgpdb = open('_MMPBSA_mutant_complex_avg.in','w')
         avgpdb.write('trajin _MMPBSA_mutant_complex.mdcrd\n')
         avgpdb.write('average _MMPBSA_mutant_avgcomplex.pdb pdb\n')
         avgpdb.close()

         ptrajentropy = open('_MMPBSA_mutant_ptrajentropy.in','w')
         ptrajentropy.write('trajin _MMPBSA_mutant_complex.mdcrd\n')
         ptrajentropy.write('reference _MMPBSA_mutant_avgcomplex.pdb\n\n')
         ptrajentropy.write('rms mass reference :1-{0}\n'.format(numres_com))
         ptrajentropy.write('matrix mwcovar name comp.matrix :1-{0}\n'.format(numres_com))
         ptrajentropy.write('analyze matrix comp.matrix out _MMPBSA_mutant_complex_entropy.out thermo reduce\n\n\n')
         ptrajentropy.write('rms mass reference {0}\n'.format(receptor_mask))
         ptrajentropy.write('matrix mwcovar name rec.matrix {0}\n'.format(receptor_mask))
         ptrajentropy.write('analyze matrix rec.matrix out _MMPBSA_mutant_receptor_entropy.out thermo reduce\n\n\n')
         ptrajentropy.write('rms mass reference {0}\n'.format(ligand_mask))
         ptrajentropy.write('matrix mwcovar name lig.matrix {0}\n'.format(ligand_mask))
         ptrajentropy.write('analyze matrix lig.matrix out _MMPBSA_mu' + 
                            'tant_ligand_entropy.out thermo reduce\n\n\n')
         ptrajentropy.close()
   # Do the quasi-harmonic ptraj entropy calculation on the mutant now
         os.system(ptraj + ' ' + mutant_complex_prmtop + ' _MMPBSA_mutant_complex_avg.in > _MMPBSA_ptraj12.out 2>&1')
         print >> sys.stdout, 'Beginning quasi-harmonic entropy calculation on mutant system with ptraj...'
         os.system('{0} {1} _MMPBSA_mutant_ptrajentropy.in > _MMPBSA_mutant_ptraj_entropy.out 2>&1'.format(ptraj,mutant_complex_prmtop))
      
      os.system(ptraj + ' ' + mutant_complex_prmtop + 
               ' _MMPBSA_mutant_complexinpcrd.in > _MMPBSA_ptraj9.out 2>&1' )
      os.system(ptraj + ' ' + mutant_receptor_prmtop + 
                ' _MMPBSA_mutant_receptorinpcrd.in > _MMPBSA_ptraj10.out 2>&1' )
      os.system(ptraj + ' ' + mutant_ligand_prmtop +  
                ' _MMPBSA_mutant_ligandinpcrd.in > _MMPBSA_ptraj11.out 2>&1' )
   
   communicator.Barrier() # wait for all threads to catch up

   testdummy1 = utils.fileexists_noprint('_MMPBSA_mutant_dummycomplex.inpcrd.1')
   testdummy2 = utils.fileexists_noprint('_MMPBSA_mutant_dummyreceptor.inpcrd.1')
   testdummy3 = utils.fileexists_noprint('_MMPBSA_mutant_dummyligand.inpcrd.1')
   
   if testdummy1 == -1 or testdummy2 == -1 or testdummy3 == -1:
      print >> sys.stderr, 'Error! Ptraj failed. Check coordinate ' +  \
               'and topology files for the complex.'
      utils.remove(debug,0)
   
   communicator.Barrier() # prevent master from moving files before other threads can test for them

   if master:
      os.system('mv _MMPBSA_mutant_dummycomplex.inpcrd.1 ' +
                '_MMPBSA_mutant_dummycomplex.inpcrd') 
      os.system('mv _MMPBSA_mutant_dummyreceptor.inpcrd.1 ' +
                '_MMPBSA_mutant_dummyreceptor.inpcrd')
      os.system('mv _MMPBSA_mutant_dummyligand.inpcrd.1 ' + 
                '_MMPBSA_mutant_dummyligand.inpcrd')
      tmdcrdmutend = time.time()

########################################################################
#
#                  END alanine scanning setup.
#
########################################################################

########################################################################
#
#   BEGIN actual calculations, including GB, PB, and mutant sander
#   calls in addition to nmode calls.
#
########################################################################

print >> sys.stdout, "\nStarting calculations\n"

communicator.Barrier() # let everyone catch up

#++++++++++++++++++++++++++++++++
# GB sander calculations
#++++++++++++++++++++++++++++++++

if gbrun:

   tgbstart = time.time()
   print >> sys.stdout, "Starting gb calculation...\n"

# do the non-mutant calculation
   if not alarun or mutant_only != 1:
      utils.gbcalc(sander,'.'+str(rank),'_MMPBSA_',complex_prmtop,
            receptor_prmtop,ligand_prmtop,idecomp)
# do the mutant calculation
   if alarun:
      utils.gbcalc(sander,'.'+str(rank),'_MMPBSA_mutant_',
            mutant_complex_prmtop,mutant_receptor_prmtop,
            mutant_ligand_prmtop,idecomp)
   tgbend = time.time()

#++++++++++++++++++++++++++++++++
# Begin PB sander calculations
#++++++++++++++++++++++++++++++++

if pbrun:

   tpbstart = time.time()

   print >> sys.stdout, "Starting pb calculation...\n"

# do the non-mutant calculation
   if not alarun or mutant_only != 1:
      utils.pbcalc(sanderpb,'.'+str(rank),'_MMPBSA_',complex_prmtop,
            receptor_prmtop,ligand_prmtop,idecomp)
# do the mutant calculation
   if alarun:
      utils.pbcalc(sanderpb,'.'+str(rank),'_MMPBSA_mutant_',
            mutant_complex_prmtop,mutant_receptor_prmtop,
            mutant_ligand_prmtop,idecomp)
   tpbend = time.time()

#++++++++++++++++++++++++++++++++
# nmode calculations
#++++++++++++++++++++++++++++++++

if nmoderun:
   print >> sys.stdout, '\nStarting nmode calculations...'
   tnmodestart = time.time()

   utils.nmodecalc("_MMPBSA_",nmode_istrng,maxcyc,drms,nmode_igb,dielc, \
                   complex_prmtop,receptor_prmtop,ligand_prmtop,suffix,\
                   nmode,'.'+str(rank), temp)

   if alarun and ala_entropy == 1:
      utils.nmodecalc("_MMPBSA_mutant_",nmode_istrng,maxcyc,drms,nmode_igb,dielc, \
                   mutant_complex_prmtop,mutant_receptor_prmtop,mutant_ligand_prmtop,suffix,\
                   nmode,'.'+str(rank), temp)

   tnmodeend = time.time()

# Now it's time to let all threads catch up and kill all but master
# since only master outputs. Then we have to combine all output files
# so they can all be parsed
communicator.Barrier()
print >> sys.stdout, '\nCalculations complete. Writing output file(s)...'
if not master:
   sys.exit()

toutputstart = time.time()

for x in range(mpi_size):
   if gbrun and not mutant_only:
      os.system('cat _MMPBSA_complex_gb.mdout.{0}'.format(x) + 
         ' >> _MMPBSA_complex_gb.mdout; rm ' +
         ' _MMPBSA_complex_gb.mdout.{0}'.format(x))
      os.system('cat _MMPBSA_ligand_gb.mdout.{0}'.format(x) + 
         ' >> _MMPBSA_ligand_gb.mdout; rm ' +
         ' _MMPBSA_ligand_gb.mdout.{0}'.format(x))
      os.system('cat _MMPBSA_receptor_gb.mdout.{0}'.format(x) + 
         ' >> _MMPBSA_receptor_gb.mdout; rm ' +
         ' _MMPBSA_receptor_gb.mdout.{0}'.format(x))
   if gbrun and alarun:
      os.system('cat _MMPBSA_mutant_complex_gb.mdout.{0}'.format(x) + 
         ' >> _MMPBSA_mutant_complex_gb.mdout; rm ' +
         ' _MMPBSA_mutant_complex_gb.mdout.{0}'.format(x))
      os.system('cat _MMPBSA_mutant_ligand_gb.mdout.{0}'.format(x) + 
         ' >> _MMPBSA_mutant_ligand_gb.mdout; rm ' +
         ' _MMPBSA_mutant_ligand_gb.mdout.{0}'.format(x))
      os.system('cat _MMPBSA_mutant_receptor_gb.mdout.{0}'.format(x) + 
         ' >> _MMPBSA_mutant_receptor_gb.mdout; rm ' +
         ' _MMPBSA_mutant_receptor_gb.mdout.{0}'.format(x))
   if pbrun and not mutant_only:
      os.system('cat _MMPBSA_complex_pb.mdout.{0}'.format(x) + 
         ' >> _MMPBSA_complex_pb.mdout; rm ' +
         ' _MMPBSA_complex_pb.mdout.{0}'.format(x))
      os.system('cat _MMPBSA_ligand_pb.mdout.{0}'.format(x) + 
         ' >> _MMPBSA_ligand_pb.mdout; rm ' +
         ' _MMPBSA_ligand_pb.mdout.{0}'.format(x))
      os.system('cat _MMPBSA_receptor_pb.mdout.{0}'.format(x) + 
         ' >> _MMPBSA_receptor_pb.mdout; rm ' +
         ' _MMPBSA_receptor_pb.mdout.{0}'.format(x))
   if pbrun and alarun:
      os.system('cat _MMPBSA_mutant_complex_pb.mdout.{0}'.format(x) + 
         ' >> _MMPBSA_mutant_complex_pb.mdout; rm ' +
         ' _MMPBSA_mutant_complex_pb.mdout.{0}'.format(x))
      os.system('cat _MMPBSA_mutant_ligand_pb.mdout.{0}'.format(x) + 
         ' >> _MMPBSA_mutant_ligand_pb.mdout; rm ' +
         ' _MMPBSA_mutant_ligand_pb.mdout.{0}'.format(x))
      os.system('cat _MMPBSA_mutant_receptor_pb.mdout.{0}'.format(x) + 
         ' >> _MMPBSA_mutant_receptor_pb.mdout; rm ' +
         ' _MMPBSA_mutant_receptor_pb.mdout.{0}'.format(x))
   if nmoderun and not mutant_only:
      os.system(('cat _MMPBSA_complex_nm.out.{0} >> _MMPBSA_complex_nm.out' +
         '; rm _MMPBSA_complex_nm.out.{0}').format(x))
      os.system(('cat _MMPBSA_receptor_nm.out.{0} >> _MMPBSA_receptor_nm.out' +
         '; rm _MMPBSA_receptor_nm.out.{0}').format(x))
      os.system(('cat _MMPBSA_ligand_nm.out.{0} >> _MMPBSA_ligand_nm.out' +
         '; rm _MMPBSA_ligand_nm.out.{0}').format(x))
   if nmoderun and alarun:
      os.system(('cat _MMPBSA_mutant_complex_nm.out.{0} >> _MMPBSA_mutant_complex_nm.out' +
         '; rm _MMPBSA_mutant_complex_nm.out.{0}').format(x))
      os.system(('cat _MMPBSA_mutant_receptor_nm.out.{0} >> _MMPBSA_mutant_receptor_nm.out' +
         '; rm _MMPBSA_mutant_receptor_nm.out.{0}').format(x))
      os.system(('cat _MMPBSA_mutant_ligand_nm.out.{0} >> _MMPBSA_mutant_ligand_nm.out' +
         '; rm _MMPBSA_mutant_ligand_nm.out.{0}').format(x))

if alarun:
   mutstring= utils.getressymbol(utils.getresinfo(mutant_residue,      \
              complex_prmtop,'RESIDUE_LABEL')) + str(mutant_residue) + \
               utils.getressymbol('ALA') + ' MUTANT:'

########################################################################
#
#     Write some details about the calculation to the final output file
#
########################################################################

keep_files = utils.PrintFinalResults(output_file, inputfile_name, strip_mdcrd, solvated_prmtop, \
                 complex_prmtop, receptor_mdcrd, solvated_receptor_prmtop, ligand_mdcrd,           \
                 solvated_ligand_prmtop, receptor_prmtop, ligand_prmtop, mdcrd, mutant_complex_prmtop, \
                 mutant_receptor_prmtop, mutant_ligand_prmtop, maskholder, numres_lig, numframes,  \
                 nmoderun, numframesnmode, pbrun, sanderpb, entropy, mutant_only, mutstring, warnings, \
                 keep_files, ala_entropy, alarun, gbrun, debug, one_trajectory, verbose, sander_apbs,  \
                 decompout, idecomp, dec_verbose, ligstart, decomprun, surften, cavity_surften, temp)

toutputend = time.time()
tglobalend = time.time()

########################################################################
#
#        END print out results to final output file
#
########################################################################

########################################################################
#
#            BEGIN Print out the final timing for the script
#
########################################################################

print >> sys.stdout, "Timing:"
if alarun:
   print >> sys.stdout, "Checking Mutant Topology Files:        " +    \
            utils.Format((tchkmutprmend - tchkmutprmstart)/60,4,3) +   \
            ' min.'
   print >> sys.stdout, "Mutating Complex For Alanine Scanning: " +    \
            utils.Format((tmdcrdmutend - tmdcrdmutstart)/60,4,3) +     \
            ' min.'

print >> sys.stdout, "Processing Trajectories With Ptraj:    " +       \
         utils.Format((tptrajnormend - tptrajnormstart)/60,4,3) +      \
         ' min.'

if gbrun:
   print >> sys.stdout, "Total GB Calculation Time (sander):    " +    \
            utils.Format((tgbend - tgbstart)/60,4,3) + ' min.'

if pbrun:
   print >> sys.stdout, "Total PB Calculation Time (sander):    " +    \
            utils.Format((tpbend - tpbstart)/60,4,3) + ' min.'

if entropy == 1:
   print >> sys.stdout, "Total Quasi-harmonic Calculation Time: " +    \
            utils.Format((tptrajentend - tptrajentstart)/60,4,3) +     \
            ' min.'

if nmoderun:
   print >> sys.stdout, "Total Harmonic nmode Calculation Time: " +    \
            utils.Format((tnmodeend - tnmodestart)/60,4,3) + ' min.'

print >> sys.stdout, "Output File Writing Time:              " +       \
         utils.Format((toutputend - toutputstart)/60,4,3) + ' min.'
print >> sys.stdout, "\nTotal Time Taken:                      " +     \
         utils.Format((tglobalend - tglobalstart)/60,4,3) + ' min.'

########################################################################
#
#              END print out the final timing for the script
#
########################################################################

utils.remove(keep_files, mpi_size)
print >> sys.stdout, "\n\nMMPBSA Finished. Thank you for using. " +    \
         "Please send any bugs/suggestions/comments to " +             \
         "amber@ambermd.org"