import getopt, sys, os
import subprocess
def writeSoluteConstrain(f, wt=500,ntr=False):
r"""
Write the lines to apply constraints to solute
Parameters
----------
f: file, Required, default: None
File of Amber input
ntr: bool, Optional, default: False
Use the ntr restraint or ibelly constrain
wt: float, Optional, default: 500
Weight of restraint if using the Amber ntr restraint
Returns
-------
None
"""
if ntr:
f.write("ntr=1,\n")
f.write("restraintmask=\':1\',\n")
f.write("restraint_wt={:f},\n".format(wt))
else:
f.write("ibelly=1,\n")
f.write("bellymask = \':2-\',\n")
[docs]def runMM(filename='water_solvated', stepsmmheat=10000, stepsmmnve=0, stepsmmnpt=300000, srun_use=False, pmemduse=False, dryrun=False, freeze_solute=False):
r"""
Equilibrate with MM
Parameters
----------
filename : str, Optional, default: 'water_solvated'
Filename prefix for .prmtop and .inpcrd files
stepsmmheat : int, Optional, default: 10000
MM heating steps
stepsmmnve : int, Optional, default: 0
MM NVE steps
stepsmmnpt : int, Optional, default: 300000
MM NPT steps
srun_use : bool, Optional, default: False
Run all commands with a srun prefix.
pmemduse : bool, Optional, default: False
Use pmemd.CUDA instead of sander
dryrun : bool, Optional, default: False
Dry run mode: only generate the commands to run MD programs and save them into a file without executing the commands
Returns
-------
None
Results stored as .netcdf files and log files
"""
frun = None
if dryrun:
frun = open('runMM.sh','w')
print('MM Energy minimization')
cmd=' -O -i mmmin.in -o mmmin.out -p '+filename+'.prmtop -c '+filename+'.inpcrd -r mm.ncrst -inf mmmin.info'
if freeze_solute:
cmd = cmd + ' -ref '+filename+'.inpcrd '
if pmemduse:
cmd= 'pmemd.cuda' +cmd
else:
cmd= 'sander'+ cmd
if srun_use:
cmd='srun -n 1 '+cmd
if dryrun:
frun.write(cmd+'\n')
else:
subprocess.call(cmd, shell=True)
if stepsmmheat>0:
print('MM Heating')
cmd=' -O -i mmheat.in -o mmheat.out -p '+filename+'.prmtop -c mm.ncrst -r mm.ncrst -x '+filename+'-heat.netcdf -inf mmheat.info'
if freeze_solute:
cmd = cmd + ' -ref '+filename+'.inpcrd '
if pmemduse:
cmd= 'pmemd.cuda' +cmd
else:
cmd= 'sander'+ cmd
if srun_use:
cmd='srun -n 1 '+cmd
if dryrun:
frun.write(cmd+'\n')
else:
subprocess.call(cmd, shell=True)
if stepsmmnve>0:
print('MM NVE equilibration')
cmd=' -O -i mmnve.in -o mmnve.out -p '+filename+'.prmtop -c mm.ncrst -r mm.ncrst -x '+filename+'-mmnve.netcdf -inf mmnve.info'
if freeze_solute:
cmd = cmd + ' -ref '+filename+'.inpcrd '
if pmemduse:
cmd= 'pmemd.cuda' +cmd
else:
cmd= 'sander'+ cmd
if srun_use:
cmd='srun -n 1 '+cmd
if dryrun:
frun.write(cmd+'\n')
else:
subprocess.call(cmd, shell=True)
if stepsmmnpt>0:
print('MM NPT equilibration')
cmd=' -O -i mmnpt.in -o mmnpt.out -p '+filename+'.prmtop -c mm.ncrst -r mm.ncrst -x '+filename+'-mmnpt.netcdf -inf mmnpt.info'
if freeze_solute:
cmd = cmd + ' -ref '+filename+'.inpcrd '
if pmemduse:
cmd= 'pmemd.cuda' +cmd
else:
cmd= 'sander'+ cmd
if srun_use:
cmd='srun -n 1 '+cmd
if dryrun:
frun.write(cmd+'\n')
else:
subprocess.call(cmd, shell=True)
if dryrun:
frun.close()
[docs]def writeQMMMTemplate(spinmult=1,charge=0,functional="b3lyp"):
r"""
Write Terachem file-based interface to Amber
Parameters
----------
spinmult : float, default: 1
spin multiplicity of solvated system
charge : int, Optional, default: 0
charge of solvated system
Returns
-------
None
Result stored as tc_job.tpl
"""
f = open("tc_job.tpl","w")
f.write("basis lacvps_ecp\n")
if spinmult==1:
f.write("method "+functional+"\n")
else:
f.write("method u"+functional+"\n")
f.write("dispersion yes\n")
f.write("scf diis+a\n")
f.write("threall 1e-13\n")
f.write("convthre 1e-6\n")
f.write("xtol 1e-4\n")
f.write("maxit 200\n")
f.write("dftgrid 1\n")
f.write("charge "+str(charge)+"\n")
f.write("spinmult "+str(spinmult)+"\n")
f.write("scrdir ./scr\n")
f.write("keep_scr yes\n")
f.write("end\n")
f.close()
[docs]def runQMMM(filename='water_solvated', spinmult=1, srun_use=False, stepsqmmmmin=250, stepsqmmmheat=1000, stepsqmmmnve=0, stepsqmmmnvt=10000, dryrun=False):
r"""
Run QMMM minimization, heating and NVT trajectory run
Parameters
----------
filename : str, Optional, default: 'water_solvated'
Filename prefix for .prmtop input file and .netcdf output file
spinmult : int, Required
Spin multiplicity of system
srun_use : bool, Optional, default: False
Run all commands with a srun prefix.
stepsqmmmmin : int, Optional, default: 250
Number of QMMM minimization steps
stepsqmmmheat : int, Optional, default: 1000
Number of QMMM heating steps
stepsqmmmnve : int, Optional, default: 0
Number of QMMM NVE steps
stepsqmmmnvt : int, Optional, default: 10000
Number of QMMM NVT trajectory steps
dryrun : bool, Optional, default: False
Dry run mode: only generate the commands to run MD programs and save them into a file without executing the commands
Returns
-------
None
Results stored in .netcdf files and log files
"""
frun = None
if dryrun:
frun = open('runQMMMM.sh','w')
if stepsqmmmmin>0:
print('QMMM Energy minimization')
cmd='sander -O -i qmmmmin.in -o qmmmmin.out -p '+filename+'.prmtop -c mm.ncrst -r qmmm.ncrst -inf qmmmmin.info -x '+filename+'-qmmmmin.netcdf'
if srun_use:
cmd='srun -n 1 '+cmd
if dryrun:
frun.write(cmd+'\n')
else:
subprocess.call(cmd, shell=True)
if spinmult>1:
print('Adjusting terachem input file for higher Spin multiplicity')
cmd="sed -i '3 a guess ./scr/ca0 ./scr/cb0' tc_job.tpl"
if srun_use:
cmd='srun -n 1 '+cmd
if dryrun:
frun.write(cmd+'\n')
else:
subprocess.call(cmd, shell=True)
cmd='sander -O -i qmmmmin.in -o qmmmmin.out -p '+filename+'.prmtop -c qmmm.ncrst -r qmmm.ncrst -inf qmmmmin.info -x '+filename+'-qmmmmin.netcdf'
if srun_use:
cmd='srun -n 1 '+cmd
if dryrun:
frun.write(cmd+'\n')
else:
subprocess.call(cmd, shell=True)
if stepsqmmmheat>0:
print('QMMM Heating')
cmd='sander -O -i qmmmheat.in -o qmmmheat.out -p '+filename+'.prmtop -c qmmm.ncrst -r qmmm.ncrst -inf qmmmheat.info -x '+filename+'-qmmmheat.netcdf'
if srun_use:
cmd='srun -n 1 '+cmd
if dryrun:
frun.write(cmd+'\n')
else:
subprocess.call(cmd, shell=True)
if stepsqmmmnve>0:
print('QMMM NVE Run')
cmd='sander -O -i qmmmnve.in -o qmmmnve.out -p '+filename+'.prmtop -c qmmm.ncrst -r qmmm.ncrst -inf qmmmnve.info -x '+filename+'-qmmmnve.netcdf'
if srun_use:
cmd='srun -n 1 '+cmd
if dryrun:
frun.write(cmd+'\n')
else:
subprocess.call(cmd, shell=True)
if stepsqmmmnvt>0:
print('QMMM NVT Run')
cmd='sander -O -i qmmmnvt.in -o qmmmnvt.out -p '+filename+'.prmtop -c qmmm.ncrst -r qmmm.ncrst -inf qmmmnvt.info -x '+filename+'-qmmmnvt.netcdf'
if srun_use:
cmd='srun -n 1 '+cmd
if dryrun:
frun.write(cmd+'\n')
else:
subprocess.call(cmd, shell=True)
if dryrun:
frun.close()
[docs]def startmd(argumentList):
r"""
Wrap function that parses command line options for autosolvate clustergen,
generates inputfiles for Amber and TeraChem,
runs MM and QMMM stages.
.. warning::
Currently, some simulation parameters cannot be set by the user, for example:
* simulation time step
* integrator type
* nonbonded cutoff
* thermostat type
* Langevin collision frequency
* barostat type
* pressure relaxation time
* frequency of trajectory writing
Parameters
----------
argumentList: list
The list contains the command line options to specify MM and QMMM stage options.
Command line options:
-f, --filename prefix of .prmtop and .inpcrd files
-t, --temp temperature in Kelvin to equilibrate
-p, --pressure pressure in bar to equilibrate during MM NPT step
-i, --stepsmmmin Number of MM minimization steps
-m, --stepsmmheat Number of MM heating steps, setting to 0 skips the MM heating step
-b, --stepsmmnve Number of MM NVE steps, setting to 0 skips the MM NVE step
-n, --stepsmmnpt Number of MM NPT steps, setting to 0 skips the MM NPT step
-l, --stepsqmmmmin Number of QMMM minimization steps, setting to 0 skips the QMMM minimization step
-o, --stepsqmmmheat Number of QMMM heating steps, setting to 0 skips the QMMM heating step
-v, --stepsqmmmnve Number of QMMM NVE steps, setting to 0 skips the QMMM NVE step
-s, --stepsqmmmnvt Number of QMMM NVT steps, setting to 0 skips the QMMM NVT step
-q, --charge Total charge of system
-u, --spinmultiplicity Spin multiplicity of whole system
-k, --functional DFT functional to use for the QM part in QM/MM
-r, --srunuse option to run inside a slurm job
-x, --pmemduse Speed up MM with pmemd.CUDA instead of sander
-d, --dryrun Dry run mode: only generate the commands to run MD programs and save them into a file without executing the command
-z, --freezesolute Freeze the solute structure throughout all MD steps
-h, --help short usage description
Returns
-------
None
Generate MD simulation input files and execute MD programs, or same the MD program execution commands in runMM.sh and runQMMM.sh
"""
#print(argumentList)
options = "hf:t:p:i:m:b:n:l:o:v:s:q:u:k:rxdz"
long_options = ["help", "filename", "temp", "pressure", "stepsmmmin", "stepsmmheat", "stepsmmnve", "stepsmmnpt", "stepsqmmmmin", "stepsqmmmheat", "stepsqmmmnve", "stepsqmmmnvt", "charge", "spinmultiplicity","functional", "srunuse", "pmemduse","dryrun","freezesolute"]
arguments, values = getopt.getopt(argumentList, options, long_options)
srun_use=False
temperature=300
pressure=1
stepsmmmin=2000
stepsmmheat=10000
stepsmmnve=0
stepsmmnpt=300000
stepsqmmmmin=250
stepsqmmmheat=1000
stepsqmmmnve=0
stepsqmmmnvt=10000
charge=0
spinmult=1
functional="b3lyp"
srun_use=False
pmemduse=False
dryrun=False
freeze_solute=False
for currentArgument, currentValue in arguments:
if currentArgument in ("-h", "-help"):
print('Usage: autosolvate mdrun [OPTIONS]')
print(' -f, --filename prefix of .prmtop and .inpcrd files')
print(' -t, --temp temperature in Kelvin to equilibrate')
print(' -p, --pressure pressure in bar to equilibrate during MM NPT step')
print(' -i, --stepsmmmin Number of MM minimization steps')
print(' -m, --stepsmmheat Number of MM heating steps, setting to 0 skips the MM heating step')
print(' -b, --stepsmmnve Number of MM NVE steps, setting to 0 skips the MM NVE step')
print(' -n, --stepsmmnpt Number of MM NPT steps, setting to 0 skips the MM NPT step')
print(' -l, --stepsqmmmmin Number of QMMM minimization steps, setting to 0 skips the QMMM minimization step')
print(' -o, --stepsqmmmheat Number of QMMM heating steps, setting to 0 skips the QMMM heating step')
print(' -v, --stepsqmmmnve Number of QMMM NVE steps, setting to 0 skips the QMMM NVE step')
print(' -s, --stepsqmmmnvt Number of QMMM NVT steps, setting to 0 skips the QMMM NVT step')
print(' -q, --charge Total charge of system')
print(' -u, --spinmultiplicity Spin multiplicity of whole system')
print(' -k, --functional DFT functional to use for the QM part in QM/MM')
print(' -r, --srunuse option to run inside a slurm job')
print(' -x, --pmemduse Speed up MM with pmemd.CUDA instead of sander')
print(' -d, --dryrun Dry run mode')
print(' -z, --freezesolute Freeze the solute')
print(' -h, --help short usage description')
exit()
elif currentArgument in ("-f", "-filename"):
print ("Filename:", currentValue)
filename=str(currentValue)
elif currentArgument in ("-t", "-temp"):
print ("Temperature in K:", currentValue)
temperature=float(currentValue)
elif currentArgument in ("-p", "-pressure"):
print ("Pressure in bar:", currentValue)
pressure=float(currentValue)
elif currentArgument in ("-i","-stepsmmmin"):
print ("Steps MM min:", currentValue)
stepsmmmin=int(currentValue)
elif currentArgument in ("-m","-stepsmmheat"):
print ("Steps MM heat:", currentValue)
stepsmmheat=int(currentValue)
elif currentArgument in ("-b","-stepsmmnve"):
print ("Steps MM NVE:", currentValue)
stepsmmnve=int(currentValue)
elif currentArgument in ("-n", "-stepsmmnpt"):
print ("Steps MM NPT:", currentValue)
stepsmmnpt=int(currentValue)
elif currentArgument in ("-l","-stepsqmmmmin"):
print ("Steps QMMM min:", currentValue)
stepsqmmmmin=int(currentValue)
elif currentArgument in ("-o","-stepsqmmmheat"):
print ("Steps QMMM heat:", currentValue)
stepsqmmmheat=int(currentValue)
elif currentArgument in ("-v","-stepsqmmmnve"):
print ("Steps QMMM NVE:", currentValue)
stepsqmmmnve=int(currentValue)
elif currentArgument in ("-s", "-stepsqmmmnvt"):
print ("Steps QMMM NPT:", currentValue)
stepsqmmmnvt=int(currentValue)
elif currentArgument in ("-q", "-charge"):
print("Charge:", currentValue)
charge=int(currentValue)
elif currentArgument in ("-u", "-spinmultiplicity"):
print ("Spinmultiplicity:", currentValue)
spinmult=int(currentValue)
elif currentArgument in ("-k", "-functional"):
print ("DFT functional:", currentValue)
functional=currentValue
elif currentArgument in ("-r", "-srunuse"):
print("using srun")
srun_use=True
elif currentArgument in ("-x", "-pmemduse"):
print("using pmemd.cuda instead of sander")
pmemduse=True
elif currentArgument in ("-d", "-dryrun"):
print("Dry run mode: only generate the commands to run MD programs and save them into a file without executing the commands")
dryrun=True
elif currentArgument in ("-z", "-freezesolute"):
print("Freeze the solute while running MD")
freeze_solute=True
print("Ignoring all QM/MM options. QM/MM will not run")
stepsqmmmmin = 0
stepsqmmmheat = 0
stepsqmmmnve = 0
stepsqmmmnvt = 0
writeMMminInput(stepsmmmin=stepsmmmin,freeze_solute=freeze_solute)
writeMMheatInput(temperature=temperature, stepsmmheat=stepsmmheat, freeze_solute=freeze_solute)
writeMMNVEInput(stepsmmnve=stepsmmnve, freeze_solute=freeze_solute)
writeMMNPTInput(temperature=temperature, pressure=pressure, stepsmmnpt=stepsmmnpt, freeze_solute=freeze_solute)
writeQMMMMinInput(stepsqmmmmin=stepsqmmmmin)
writeQMMMTemplate(spinmult=spinmult, charge=charge, functional=functional)
writeQMMMInput(temperature=temperature, stepsqmmm=stepsqmmmheat, charge=charge, infilename='qmmmheat.in')
writeQMMMInput(stepsqmmm=stepsqmmmnve, charge=charge, infilename='qmmmnve.in', nve=True)
writeQMMMInput(temperature=temperature, stepsqmmm=stepsqmmmnvt, charge=charge, infilename='qmmmnvt.in')
runMM(filename=filename, stepsmmheat=stepsmmheat, stepsmmnpt=stepsmmnpt, stepsmmnve=stepsmmnve, srun_use=srun_use, pmemduse=pmemduse, dryrun=dryrun, freeze_solute=freeze_solute)
runQMMM(filename=filename, spinmult=spinmult, srun_use=srun_use, stepsqmmmmin=stepsqmmmmin, stepsqmmmheat=stepsqmmmheat, stepsqmmmnve=stepsqmmmnve, stepsqmmmnvt=stepsqmmmnvt, dryrun=dryrun)
if __name__ == '__main__':
argumentList = sys.argv[1:]
startmd(argumentList)