from openbabel import pybel
import getopt, sys, os
from openbabel import openbabel as ob
import subprocess
import pkg_resources
from autosolvate.globs import keywords_avail, available_qm_programs, available_charge_methods
from autosolvate.resp_classes.resp_factory import resp_factory
from autosolvate.pubchem_api import PubChemAPI
from autosolvate.solute_info import Solute
amber_solv_dict = {'water': [' ','TIP3PBOX '],
'methanol': ['loadOff solvents.lib\n loadamberparams frcmod.meoh\n', 'MEOHBOX '],
'chloroform': ['loadOff solvents.lib\n loadamberparams frcmod.chcl3\n', 'CHCL3BOX '],
'nma': ['loadOff solvents.lib\n loadamberparams frcmod.nma\n', 'NMABOX ']}
custom_solv_dict = {'acetonitrile':'ch3cn'}
[docs]class solventBoxBuilder():
r"""
Solvated molecule in specified solvent.
Parameters
----------
solvent : str, Optional, default: 'water'
Currently implemented solvents are: 'water', 'methanol', 'chloroform', 'nma', 'acetonitrile'
slu_netcharge: int, Optional, default 0
Charge of solute, the solvent box will be neutralized with Cl- and Na+ ions
cube_size: float, Optional, default: 54
Size of MM solvent box
charge_method: str, Optional, default: "resp"
Use 'resp' (quantum mechanical calculation needed) or 'bcc' to estimate partial charges
slu_spinmult: int, Optional, default: 1
Spinmultiplicity of solute
outputFile: str, Optional, default='water_solvated'
Filename-prefix for outputfiles
srun_use: bool, Optional, default='False
Run all commands with a srun prefix
Returns
-------
None
To run solvation, call build function.
"""
[docs] def __init__(self, xyzfile, solvent='water', slu_netcharge=0, cube_size=54,
charge_method="resp", slu_spinmult=1, outputFile="",
srun_use=False, qm_program='gaussian', qm_exe=None, qm_dir=None,
amberhome=None, closeness=0.8,
solvent_off="", solvent_frcmod="",rundir=""):
self.xyz = xyzfile
self.solute = pybel.readfile('xyz', xyzfile).__next__()
self.slu_netcharge = slu_netcharge
self.slu_spinmult = slu_spinmult
# currently hard coded. Can be changed later to be determined automatically based on the density of the solute
self.solvent = solvent
if closeness=='automated':
if self.solvent=='acetonitrile':
self.closeness=1.88
elif self.solvent=='water':
self.closeness=0.50
elif self.solvent=='methanol':
self.closeness=0.60
elif self.solvent=='nma':
self.closeness=0.58
elif self.solvent=='chloroform':
self.closeness=0.58
else:
print("Warning: solvent not supported for automated closeness determination")
else:
self.closeness = closeness
self.waterbox_size = 8.0
# following are for custom organic solvent
self.cube_size = cube_size # in angstrom
self.slu_pos = self.cube_size/2.0
self.slv_count = 210*8 # 210
self.pbcbox_size = self.cube_size+2
self.outputFile=outputFile
if self.outputFile == "":
self.outputFile = self.solvent + "_solvated"
self.srun_use=srun_use
self.charge_method=charge_method
self.qm_program = qm_program
self.qm_exe=qm_exe
self.qm_dir = qm_dir
self.amberhome = amberhome
self.rundir = rundir
self.solvent_off=solvent_off
self.solvent_frcmod=solvent_frcmod
self.is_custom_solvent = False
self.inputCheck()
def inputCheck(self):
if self.charge_method not in available_charge_methods:
print("Error: the requested charge method: {}, is not supported yet".format(self.charge_method))
exit(1)
if self.slu_spinmult > 1:
if self.charge_method != "resp":
print("Error: solute spin multiplicity: ", self.slu_spinmult, " charge method: ", self.charge_method)
print("Error: atomic charge fitting for open-shell system only works for resp charge method")
print("Error: exiting...")
exit(1)
if self.charge_method == "resp":
if self.qm_program not in available_qm_programs:
print("Error! Selectd quantum chemistry package for RESP charge fittig is not supported:", qm_program)
exit()
elif self.qm_program == "gaussian":
if self.qm_exe == None:
print("WARNING: Gaussian executable name is not specified for RESP charge fitting!")
print("WARNING: Using g16 by default. If failed later,\
please rerun with the option -e specified!")
self.qm_exe = 'g16'
if self.qm_dir == None:
print("WARNING: Gaussian executable directory is not specified for RESP charge fitting!")
print("WARNING: Setting to default path: /opt/packages/gaussian/g16RevC.01/g16/")
print("WARNING: If failed later, please rerun with the option -d specified!")
self.qm_dir = '/opt/packages/gaussian/g16RevC.01/g16/'
gspath = os.path.join(self.qm_dir, self.qm_exe)
if not os.path.exists(gspath):
print("Error! Gaussian executable path",gspath)
print("Does not exist! Exiting...")
exit()
elif self.qm_program == "gamess":
if self.qm_exe == None:
print("WARNING: Gamess executable name is not specified for RESP charge fitting!")
print("WARNING: Using rungms by default. If failed later,\
please rerun with the option -e specified!")
self.qm_exe = 'rungms'
if self.qm_dir == None:
print("WARNING: GAMESS executable directory is not specified for RESP charge fitting!")
print("WARNING: Setting to default path: ")
print("WARNING: If failed later, please rerun with the option -d specified!")
self.qm_dir = '/opt/gamess/'
gamess_path = os.path.join(self.qm_dir, self.qm_exe)
if not os.path.exists(gamess_path):
print("Error! Gamess executable path",gamess_path)
print("Does not exist! Exiting...")
exit()
if self.amberhome == None:
print("WARNING: Amber home directory is not specified in input options")
print("WARNING: Checking AMBERHOME environment variable...")
cmd = ["echo", "$AMBERHOME"]
print(cmd)
proc=subprocess.Popen(cmd,
universal_newlines=True,
stdout=subprocess.PIPE)
output=proc.communicate()[0]
if len(output)<2:
print("ERROR: AMBERHOME not defined")
print("ERROR: Please set up AMBERHOME environment or specify through command line option -a")
print("ERROR: Exiting...")
else:
print("WARNING: AMBERHOME detected: ", output)
self.amberhome = output
else:
print("AMBERHOME path provided from input: ", self.amberhome)
print("Validating path...")
if not os.path.isdir(self.amberhome):
print("ERROR: provided AMBERHOME path does not exist! Exiting...")
print("Exporting AMBERHOME environment variable:")
cmd = "export AMBERHOME=" + self.amberhome
print(cmd)
subprocess.call(cmd, shell=True)
print("AMBERHOME environment variable export finished.")
# check whether requested custom solvent is available
if self.solvent not in amber_solv_dict.keys() and self.solvent not in custom_solv_dict.keys():
print("Requested solvent name not contained in AutoSolvate.")
print("Checking available of custom solvent .frcmod and .off files")
if len(self.solvent_frcmod) == 0:
print("ERROR! Custom solvent .frcmod file is not provided!")
exit()
elif len(self.solvent_off) == 0:
print("ERROR! Custom solvent .off library file is not provided!")
exit()
elif not os.path.exists(self.solvent_frcmod):
print("ERROR! Custom solvent .frcmod file ", self.solvent_frcmod,
"ERROR! does not exist!")
exit()
elif not os.path.exists(self.solvent_off):
print("ERROR! Custom solvent .off library file ", self.solvent_frcmod,
"ERROR! does not exist!")
exit()
else:
self.is_custom_solvent = True
[docs] def getSolutePDB(self):
r"""
Convert xyz to pdb
Parameters
----------
None
Returns
-------
None
"""
print("Converting xyz to pdb")
# cmd = "babel " + self.xyz + " solute.pdb"
# subprocess.call(cmd, shell=True)
obConversion = ob.OBConversion()
obConversion.SetInAndOutFormats("xyz", "pdb")
obmol = self.solute.OBMol
obConversion.WriteFile(obmol,os.path.join(self.rundir,'solute.xyz.pdb'))
# Change the residue name from the default UNL to SLU
pdb1 = open('solute.xyz.pdb').readlines()
pdb2 = open('solute.xyz.pdb','w')
for line in pdb1:
newline = line.replace('UNL','SLU')
pdb2.write(newline)
pdb2.close()
def removeConectFromPDB(self):
print("cleaning up solute.xyz.pdb")
pdb1 = open('solute.xyz.pdb').readlines()
pdb2 = open('solute.xyz.pdb','w')
for line in pdb1:
if 'CONECT' not in line:
pdb2.write(line)
pdb2.close()
[docs] def getFrcmod(self):
r"""
Get partial charges and create frcmod
Parameters
----------
None
Returns
-------
None
"""
print("Generate frcmod file for the solute")
print("Remeove CONNECT information from pdb file")
self.removeConectFromPDB()
if self.charge_method == "resp":
myresp = resp_factory(pdbfile="solute.xyz.pdb", charge=self.slu_netcharge,
spinmult=self.slu_spinmult, qm_program=self.qm_program,
qm_exe=self.qm_exe, qm_dir=self.qm_dir,srun_use=self.srun_use,rundir=self.rundir)
myresp.run()
if self.charge_method == "bcc":
print("AnteChamber: Generate mol2 with bcc charge.")
cmd3="$AMBERHOME/bin/antechamber -i solute.xyz.pdb -fi pdb -o solute.mol2 -fo mol2 -c bcc -eq 2 -rn SLU"
cmd3 += " -nc " + str(self.slu_netcharge) + " -m " + str(self.slu_spinmult)
print(cmd3)
if self.srun_use:
cmd3='srun -n 1 '+cmd3
subprocess.call(cmd3, shell=True)
print("Finally generate frcmod with parmchk2")
cmd4 ="$AMBERHOME/bin/parmchk2 -i solute.mol2 -f mol2 -o solute.frcmod"
if self.srun_use:
cmd4='srun -n 1 '+cmd4
subprocess.call(cmd4, shell=True)
[docs] def getHeadTail(self):
r"""
Detect start and end of coordinates
Parameters
----------
None
Returns
-------
None
"""
pdb = open('solute.mol2').readlines()
start =0
end = 0
for i in range(0,len(pdb)):
if "@<TRIPOS>ATOM" in pdb[i]:
start = i+1
if "@<TRIPOS>BOND" in pdb[i]:
end = i
for i in range(start, end):
atom = pdb[i].split()[1]
if "H" not in atom:
self.head = atom
break
for i in reversed(range(start,end)):
atom = pdb[i].split()[1]
if "H" not in atom:
self.tail = atom
break
[docs] def writeTleapcmd1(self):
r"""
Create tleap input file
Parameters
----------
None
Returns
-------
None
"""
self.getHeadTail()
f = open("leap.cmd","w")
f.write("source leaprc.protein.ff14SB\n")
f.write("source leaprc.gaff\n")
f.write("loadamberparams solute.frcmod\n")
f.write("SLU = loadmol2 solute.mol2\n")
f.write("check SLU\n")
f.write("set SLU head SLU.1."+self.head + "\n")
f.write("set SLU tail SLU.1."+self.tail + "\n")
f.write("saveoff SLU solute.lib\n")
f.write("savepdb SLU solute.pdb\n")
f.write("quit\n")
f.close()
[docs] def createLib(self):
r"""
Run tleap
Parameters
----------
None
Returns
-------
None
"""
print("Now create the solute library file")
self.writeTleapcmd1()
cmd ="tleap -s -f leap.cmd > leap_savelib.log"
if self.srun_use:
cmd='srun -n 1 '+cmd
subprocess.call(cmd, shell=True)
if not os.path.isfile('solute.pdb'):
print("tleap failed to generate solute.pdb")
sys.stdout.flush()
sys.exit()
[docs] def writeTleapcmd_add_solvent(self):
r"""
Write tleap input file to add solvent
Parameters
----------
None
Returns
-------
None
"""
if self.solvent in amber_solv_dict:
print("Now add pre-equlibrated solvent box to the solute")
self.getHeadTail()
f = open("leap_add_solventbox.cmd","w")
f.write("source leaprc.protein.ff14SB\n")
f.write("source leaprc.gaff\n")
f.write("source leaprc.water.tip3p\n")
f.write(str(amber_solv_dict[str(self.solvent)][0]))
f.write("loadamberparams solute.frcmod\n")
f.write("mol=loadmol2 solute.mol2\n")
f.write("check mol\n")
f.write("solvatebox mol " + (str(amber_solv_dict[str(self.solvent)][1])) + str(self.slu_pos) + " iso "+str(self.closeness)+" #Solvate the complex with a cubic solvent box\n")
# Notice that we want to add the ion after solvation because we don't want the counter ion to be too close to solute
if self.slu_netcharge != 0:
if self.slu_netcharge > 0:
ion = 'Cl-'
else:
ion = 'Na+'
f.write("addIons2 mol " + ion + " 0\n")
f.write("check mol\n")
f.write("check mol\n")
f.write("savepdb mol " + str(self.outputFile) + ".pdb\n")
f.write("saveamberparm mol " + str(self.outputFile) + ".prmtop " + str(self.outputFile) + ".inpcrd\n")
f.write("quit\n")
f.close()
def writeTleapcmd_add_solvent_custom(self):
r"""
Write tleap input file to add solvent from user provided .off and .frcmod files
Parameters
----------
None
Returns
-------
None
"""
print("Now add custom pre-equlibrated solvent box to the solute")
self.getHeadTail()
f = open("leap_add_solventbox.cmd","w")
f.write("source leaprc.protein.ff14SB\n")
f.write("source leaprc.gaff\n")
f.write("source leaprc.water.tip3p\n")
f.write("loadoff " + self.solvent_off + "\n")
f.write("loadamberparams " + self.solvent_frcmod + "\n")
f.write("loadamberparams solute.frcmod\n")
f.write("mol=loadmol2 solute.mol2\n")
f.write("check mol\n")
f.write("solvatebox mol " + self.solvent + " "
+ str(self.slu_pos) + " iso 0.8 #Solvate the complex with a cubic solvent box\n")
# Notice that we want to add the ion after solvation because we don't want the counter ion to be too close to solute
if self.slu_netcharge != 0:
if self.slu_netcharge > 0:
ion = 'Cl-'
else:
ion = 'Na+'
f.write("addIons2 mol " + ion + " 0\n")
f.write("check mol\n")
f.write("check mol\n")
f.write("savepdb mol " + str(self.outputFile) + ".pdb\n")
f.write("saveamberparm mol " + str(self.outputFile) + ".prmtop " + str(self.outputFile) + ".inpcrd\n")
f.write("quit\n")
f.close()
[docs] def processPackmolPDB(self):
r"""
Convert file for custom solvents like CH3CN
Parameters
----------
None
Returns
-------
None
"""
solvPrefix = custom_solv_dict[self.solvent]
data=open(solvPrefix + '_solvated.packmol.pdb').readlines()
output=open(solvPrefix + '_solvated.processed.pdb','w')
this_resid = 1
last_resid = 1
for line in data:
if 'ATOM' in line:
last_resid = int(this_resid)
this_resid = int(line[22:26])
if last_resid != this_resid:
output.write("TER\n")
output.write(line)
output.close()
[docs] def packSLUSLV(self):
r"""
Write packmol input file for custom solvents
Parameters
----------
None
Returns
-------
None
"""
print("Now use packmol to pack the solute and solvent into a box")
if self.solvent not in custom_solv_dict.keys():
print("The type of solvent is not implemented yet")
return
else:
# Solvent pdb file stored in our package data folder
# Path is too long for Packmol to recognize. Copy to current rundir
solvPrefix = custom_solv_dict[self.solvent]
solvent_pdb = solvPrefix+'.pdb'
solvent_pdb_origin = pkg_resources.resource_filename('autosolvate',
os.path.join('data/', solvPrefix, solvent_pdb))
subprocess.call(['cp',solvent_pdb_origin,solvent_pdb])
output_pdb = solvPrefix + "_solvated.packmol.pdb"
packmol_inp = open('packmol.inp','w')
packmol_inp.write("# All atoms from diferent molecules will be at least %s Angstroms apart\n" % self.closeness)
packmol_inp.write("tolerance %s\n" % self.closeness)
packmol_inp.write("\n")
packmol_inp.write("filetype pdb\n")
packmol_inp.write("\n")
packmol_inp.write("output " + output_pdb + "\n")
packmol_inp.write("\n")
packmol_inp.write("# add the solute\n")
packmol_inp.write("structure solute.pdb\n")
packmol_inp.write(" number 1\n")
packmol_inp.write(" fixed " + " " + str(self.slu_pos) + " "+ str(self.slu_pos) + " " + str(self.slu_pos) + " 0. 0. 0.\n")
packmol_inp.write(" centerofmass\n")
packmol_inp.write(" resnumbers 2 \n")
packmol_inp.write("end structure\n")
packmol_inp.write("\n")
packmol_inp.write("# add first type of solvent molecules\n")
packmol_inp.write("structure "+ solvent_pdb + "\n")
packmol_inp.write(" number " + str(self.slv_count) + " \n")
packmol_inp.write(" inside cube 0. 0. 0. " + str(self.cube_size) + " \n")
packmol_inp.write(" resnumbers 2 \n")
packmol_inp.write("end structure\n")
packmol_inp.close()
cmd ="packmol < packmol.inp > packmol.log"
subprocess.call(cmd, shell=True)
if self.srun_use:
cmd='srun -n 1 '+cmd
self.processPackmolPDB()
[docs] def writeTleapcmd_custom_solvated(self):
r"""
Write tleap file for custom solvents like CH3CN
Parameters
----------
None
Returns
-------
None
"""
solvPrefix = custom_solv_dict[self.solvent]
solvent_frcmod = solvPrefix+'.frcmod'
solvent_frcmod_path = pkg_resources.resource_filename('autosolvate',
os.path.join('data',solvPrefix,solvent_frcmod))
solvent_prep = solvPrefix+'.prep'
solvent_prep_path = pkg_resources.resource_filename('autosolvate',
os.path.join('data',solvPrefix,solvent_prep))
f = open("leap_packmol_solvated.cmd","w")
f.write("source leaprc.protein.ff14SB\n")
f.write("source leaprc.gaff\n")
f.write("source leaprc.water.tip3p\n") # This will load the Ions. Neccessary
f.write("loadamberparams " + solvent_frcmod_path + "\n")
f.write("loadAmberPrep " + solvent_prep_path + "\n")
f.write("loadamberparams solute.frcmod\n")
f.write("loadoff solute.lib\n")
f.write("\n")
f.write("SYS = loadpdb " + solvPrefix + "_solvated.processed.pdb\n")
f.write("check SYS\n")
f.write("\n")
if self.slu_netcharge > 0:
ion = 'Cl-'
else:
ion = 'Na+'
if self.slu_netcharge != 0:
if self.slu_netcharge > 0:
ion = 'Cl-'
else:
ion = 'Na+'
f.write("addIons2 SYS " + ion + " 0\n")
f.write("check SYS\n")
f.write("# set the dimension of the periodic box\n")
f.write("set SYS box {" + str(self.pbcbox_size) +", " + str(self.pbcbox_size) + ", " + str(self.pbcbox_size) + "}\n")
f.write("\n")
f.write("saveamberparm SYS "+str(self.outputFile)+".prmtop "+str(self.outputFile)+".inpcrd #Save AMBER topology and coordinate files\n")
f.write("savepdb SYS "+str(self.outputFile)+".pdb\n")
f.write("quit\n")
f.close()
[docs] def createAmberParm(self):
r"""
Generate Amber parameters with tleap
Parameters
---------
None
Returns
---------
None
Creates files in current directory.
"""
print("Generate Amber parameters for the solvated system")
if self.solvent in amber_solv_dict:
self.writeTleapcmd_add_solvent()
cmd ="tleap -s -f leap_add_solventbox.cmd > leap_add_solventbox.log"
if self.srun_use:
cmd='srun -n 1 '+cmd
subprocess.call(cmd, shell=True)
elif self.solvent in custom_solv_dict.keys():
self.packSLUSLV()
self.writeTleapcmd_custom_solvated()
cmd ="tleap -s -f leap_packmol_solvated.cmd > leap_packmol_solvated.log"
if self.srun_use:
cmd='srun -n 1 '+cmd
subprocess.call(cmd, shell=True)
elif self.is_custom_solvent:
self.writeTleapcmd_add_solvent_custom()
cmd ="tleap -s -f leap_add_solventbox.cmd > leap_add_solventbox.log"
if self.srun_use:
cmd='srun -n 1 '+cmd
subprocess.call(cmd, shell=True)
else:
print('solvent not supported')
[docs] def build(self):
r"""
Run solventBoxBuilder
Parameters
---------
None
Returns
---------
None
Creates files in current directory.
"""
self.getSolutePDB()
self.getFrcmod()
self.createLib()
self.createAmberParm()
print("The script has finished successfully")
[docs]def startboxgen(argumentList):
r"""
Wrap function that parses command line options for autosolvate boxgen,
adds solvent box to a given solute,
and generates related force field parameters.
Parameters
----------
argumentList: list
The list contains the command line options to specify solute, solvent, and other options
related to structure and force field parameter generation.
Command line option definitions
-n, --solutename initialize suggestion for box parameter for a given solute's name
-m, --main solute xyz file
-s, --solvent name of solvent (water, methanol, chloroform, nma)
-o, --output prefix of the output file names
-c, --charge formal charge of solute
-u, --spinmultiplicity spin multiplicity of solute
-g, --chargemethod name of charge fitting method (bcc, resp)
-b, --cubesize size of solvent cube in angstroms
-r, --srunuse option to run inside a slurm job
-q, --qmprogram name of the quantum chemistry program (gaussian, gamess, terachem)
-e, --qmexe name of the quantum chemistry package executable used to generate electrostatic potential needed for RESP charge fitting
-d, --qmdir path to the quantum chemistry package
-a, --amberhome path to the AMBER molecular dynamics package root directory. Definition of the environment variable $AMBERHOME
-t, --closeness Solute-solvent closeness setting, for acetonitrile tolerance parameter in packmol in Å, for water, methanol, nma, chloroform the scaling factor in tleap, setting to 'automated' will automatically set this parameter based on solvent.
-l, --solventoff path to the custom solvent .off library file. Required if the user want to use some custom solvent other than the 5 solvents contained in AutoSolvate (TIP3P water, methanol, NMA, chloroform, MeCN)
-p, --solventfrcmod path to the custom solvent .frcmod file. Required if the user wants to use some custom solvent other than the 5 solvents contained in AutoSolvate.
-v, --validation option to run validation step for given solute's name
-h, --help short usage description
Returns
-------
None
Generates the structure files and save as ```.pdb```. Generates the MD parameter-topology and coordinates files and saves as ```.prmtop``` and ```.inpcrd```
"""
#print(argumentList)
options = "hm:n:s:o:c:b:g:u:rq:e:d:a:t:l:p:vD:"
long_options = ["help", "main","solutename", "solvent", "output", "charge", "cubesize", "chargemethod", "spinmultiplicity", "srunuse","qmprogram","qmexe", "qmdir", "amberhome", "closeness","solventoff","solventfrcmod", "validation", "runningdirectory"]
arguments, values = getopt.getopt(argumentList, options, long_options)
solutename = ""
solutexyz=""
solvent='water'
slu_netcharge=0
cube_size=54
charge_method="bcc"
slu_spinmult=1
mult_given = 0
mult_suggest = 0
outputFile=""
srun_use=False
amberhome=None
qmprogram="gaussian"
qmexe=None
qmdir=None
closeness=0.8
solvent_off=""
solvent_frcmod=""
rundir = ""
for currentArgument, currentValue in arguments:
if currentArgument in ("-h", "--help"):
print('Usage: autosolvate boxgen [OPTIONS]')
print(' -m, --main solute xyz file')
print(' -s, --solvent name of solvent')
print(' -o, --output prefix of the output file names')
print(' -c, --charge formal charge of solute')
print(' -u, --spinmultiplicity spin multiplicity of solute')
print(' -g, --chargemethod name of charge fitting method (bcc, resp)')
print(' -b, --cubesize size of solvent cube in angstroms')
print(' -r, --srunuse option to run inside a slurm job')
print(' -q, --qmprogram name of the quantum chemistry program')
print(' -e, --qmexe name of the quantum chemistry package executable')
print(' -d, --qmdir path to the quantum chemistry package')
print(' -a, --amberhome path to the AMBER molecular dynamics package root directory')
print(' -t, --closeness Solute-solvent closeness setting')
print(' -l, --solventoff path to the custom solvent .off library file')
print(' -D, --rundir running directory where temporary files are stored')
print(' -p, --solventfrcmod path to the custom solvent .frcmod file')
print(' -n, --solutename initialize suggested parameter for given solute')
print(' -v, --validation option to run validation step for input parameters')
print(' -h, --help short usage description')
exit()
elif currentArgument in ("-m", "--main"):
print ("Main/solutexyz", currentValue)
solutexyz=str(currentValue)
elif currentArgument in ("-n", "--solutename"):
if solutexyz == "":
solutename=str(currentValue)
sol=PubChemAPI(solutename)
info=sol.get_info()
solutexyz=str(info[3])
slu_netcharge = info[2]
solS=Solute(info[0], info[1], info[2], info[3])
cube_size = solS.get_box_length()
slu_spinmult = solS.get_spin_multiplicity()
charge_method = solS.get_methods()[0]
else:
solS = Solute("", "", slu_netcharge, solutexyz)
mol = next(pybel.readfile("xyz", solutexyz))
total_electrons = sum(atom.atomicnum for atom in mol.atoms)
slu_spinmult = (total_electrons - slu_netcharge) % 2 + 1
if slu_spinmult > 1: charge_method = 'resp'
cube_size = solS.get_box_length()
elif currentArgument in ("-s", "--solvent"):
print ("Solvent:", currentValue)
solvent=str(currentValue)
elif currentArgument in ("-o", "--output"):
print ("Output:", currentValue)
outputFile=str(currentValue)
elif currentArgument in ("-c", "--charge"):
print ("Charge:", currentValue)
slu_netcharge=int(currentValue)
elif currentArgument in ("-b", "--cubesize"):
print ("Cubesize:", currentValue)
cube_size=float(currentValue)
elif currentArgument in ("-g", "--chargemethod"):
print ("Chargemethod:", currentValue)
charge_method=str(currentValue)
elif currentArgument in ("-u", "--spinmultiplicity"):
print ("Spinmultiplicity:", currentValue)
slu_spinmult=int(currentValue)
elif currentArgument in ("-r", "--srunuse"):
print("usign srun")
srun_use=True
elif currentArgument in ("-q","--qmprogram"):
print("Quantum chemistry program name:", currentValue)
qmprogram = currentValue
elif currentArgument in ("-e","--qmexe"):
print("Quantum chemistry package executable name:", currentValue)
qmexe = currentValue
elif currentArgument in ("-d","--qmdir"):
print("Quantum chemistry package directory:", currentValue)
qmdir = currentValue
elif currentArgument in ("-a","--amberhome"):
print("Amber home directory:", currentValue)
amberhome = currentValue
elif currentArgument in ("-t", "--closeness"):
print("Solute-Solvente closeness parameter", currentValue)
closeness = currentValue
elif currentArgument in ("-l", "--solventoff"):
print("Custom solvent .off library path:", currentValue)
solvent_off = currentValue
elif currentArgument in ("-p", "--solventfrcmod"):
print("Custom solvent .frcmmod file path:", currentValue)
solvent_frcmod = currentValue
elif currentArgument in ("-v", "--validation"):
print('Validating...')
solS = Solute("", "", slu_netcharge, solutexyz)
mol = next(pybel.readfile("xyz", solutexyz))
total_electrons = sum(atom.atomicnum for atom in mol.atoms)
mult_suggest = ((total_electrons - slu_netcharge) % 2 + 1) % 2
mult_given = slu_spinmult % 2
if slu_spinmult == 0 or mult_suggest != mult_given:
raise Exception("Incorrect solute spin multiplicity given, please double check your value or use suggestion function enabled by -n or --solutename")
if slu_spinmult > 1 and charge_method != 'resp':
raise Exception("Incorrect charge method given, please double check your value or use suggestion function enabled by -n or --solutename")
if cube_size < solS.get_box_length():
raise Exception("Solvent box is too small, please increase box length.")
print('Passed')
elif currentArgument in ('-D, --rundir '):
print("Directory for Temporary files:",currentValue)
rundir = currentValue
if solutexyz == "":
print("Error! Solute xyzfile must be provided!\nExiting...")
exit()
elif not os.path.exists(solutexyz):
print("Error! Solute xyzfile path ",solutexyz, " does not exist!\nExiting...")
exit()
try:
pybel.readfile('xyz', solutexyz).__next__()
except:
print("Error! Solute xyzfile format issue!")
print(solutexyz," cannot be opened with openbabel.\n Exiting...")
exit()
builder = solventBoxBuilder(solutexyz, solvent=solvent, slu_netcharge=slu_netcharge,
cube_size=cube_size, charge_method=charge_method,
slu_spinmult=slu_spinmult, outputFile=outputFile, srun_use=srun_use,
qm_program=qmprogram, qm_exe=qmexe, qm_dir=qmdir,
amberhome=amberhome, closeness=closeness,
solvent_off=solvent_off, solvent_frcmod=solvent_frcmod,rundir=rundir)
builder.build()
if __name__ == '__main__':
argumentList = sys.argv[1:]
startboxgen(argumentList)