""" Set up of solvated complexes using AmberTools LEaP. @author John D. Chodera """ #=============================================================================== # IMPORTS #=============================================================================== import sys import os #=============================================================================== # TEMPLATES #=============================================================================== leap_templates = dict() leap_templates['getcharge'] = """ # Set up complex for explicit solvent simulation. # Load AMBER ff99sb-ildn forcefield for protein. source leaprc.%(protein_forcefield)s # Load GAFF parameters for ligand. source leaprc.gaff # Load in protein (with all residues modeled in). receptor = loadPdb "%(receptor_pdb_filename)s" %(receptor_extra)s # Load ligand and parameters. loadAmberParams "%(ligand_frcmod_filename)s" ligand = loadMol2 "%(ligand_mol2_filename)s" # Create complex. complex = combine { receptor ligand } # Get complex charge charge complex # Exit quit """ leap_templates['explicit'] = """ # Set up complex for explicit solvent simulation. # Load AMBER ff99sb-ildn forcefield for protein. source leaprc.%(protein_forcefield)s # Load GAFF parameters for ligand. source leaprc.gaff # Load in protein (with all residues modeled in). receptor = loadPdb "%(receptor_pdb_filename)s" %(receptor_extra)s # Add counterions addions receptor Na+ %(ncations)d addions receptor Cl- %(nanions)d # Load ligand and parameters. loadAmberParams "%(ligand_frcmod_filename)s" ligand = loadMol2 "%(ligand_mol2_filename)s" # Create complex. complex = combine { receptor ligand } # Check complex. check complex # Report on net charge. charge ligand charge receptor charge complex # Solvate ligand, receptor, and complex. solvatebox ligand TIP3PBOX %(solvent_clearance)f %(closeness)f iso solvatebox receptor TIP3PBOX %(solvent_clearance)f %(closeness)f iso solvatebox complex TIP3PBOX %(solvent_clearance)f %(closeness)f iso # Write parameters. saveAmberParm ligand "ligand.prmtop" "ligand.inpcrd" saveAmberParm receptor "receptor.prmtop" "receptor.inpcrd" saveAmberParm complex "complex.prmtop" "complex.inpcrd" # Exit quit """ leap_templates['implicit'] = """ # Set up complex for GBSA simulation with OBC model. # Set GB radii to recommended values for OBC. set default PBRadii mbondi2 # Load AMBER ff99sb-ildn forcefield for protein. source leaprc.%(protein_forcefield)s # Load GAFF parameters for ligand. source leaprc.gaff # Load in protein (with all residues modeled in). receptor = loadPdb %(receptor_filename)s # Load ligand and parameters. loadAmberParams %(ligand_frcmod_filename)s ligand = loadMol2 %(ligand_mol2_filename)s # Create complex. complex = combine { receptor ligand } # Check complex. check complex # Report on net charge. charge ligand charge receptor charge complex # Write parameters. saveAmberParm ligand "ligand.prmtop" "ligand.inpcrd" saveAmberParm receptor "receptor.prmtop" "receptor.inpcrd" saveAmberParm complex "complex.prmtop" "complex.inpcrd" # Exit quit """ #=============================================================================== # SUBROUTINES #=============================================================================== def write_file(filename, contents): """Write the specified contents to a file. ARGUMENTS filename (string) - the file to be written contents (string) - the contents of the file to be written """ outfile = open(filename, 'w') if type(contents) == list: for line in contents: outfile.write(line) elif type(contents) == str: outfile.write(contents) else: raise "Type for 'contents' not supported: " + repr(type(contents)) outfile.close() return def read_file(filename): """Read contents of the specified file. ARGUMENTS filename (string) - the name of the file to be read RETURNS lines (list of strings) - the contents of the file, split by line """ infile = open(filename, 'r') lines = infile.readlines() infile.close() return lines def create_system(receptor_pdb_filename, ligand_frcmod_filename, ligand_mol2_filename, prefix, solvent='explicit', protein_forcefield='ff99SBildn', solvent_clearance=10.0, ncations=0, nanions=0, otherpdb=None, autocharge=False, leap_output_filename=None, verbose=False): """ Set up AMBER system using LEaP. ARGUMENTS receptor_pdb_filename (string) - source receptor filename (which must have crystallographic waters included already, if desired) ligand_frcmod_filename (string) - ligand frcmod filename ligand_mol2_filename (string) - ligand GAFF mol2 filename prefix (string) - prefix for writing out AMBER prmtop/inpcrd files (prefix-{receptor,ligand,complex}.{prmtop,inpcrd}) OPTIONAL ARGUMENTS solvent (string) - either 'explicit' or 'implicit' (default: 'explicit') otherpdb (string) - if specified, add these (e.g. crystallographic waters) to receptor (default: None) autocharge (boolean) - if True, will determine number of monovalent counterions to add automatically through an initial LEaP pass (default: False) leap_output_filename (string) - if specified, leap output will be written to this directory verbose (boolean) - if True, will print verbose output (default: False) TODO * Make this more general. * Allow leap input file to be written to output destination directory. """ # Get original filenames. import os.path receptor_pdb_filename = os.path.abspath(receptor_pdb_filename) ligand_frcmod_filename = os.path.abspath(ligand_frcmod_filename) ligand_mol2_filename = os.path.abspath(ligand_mol2_filename) if otherpdb: otherpdb = os.path.abspath(otherpdb) if leap_output_filename: leap_output_filename = os.path.abspath(leap_output_filename) prefix = os.path.abspath(prefix) # Change to a temporary working directory. import os, tempfile original_directory = os.getcwd() temporary_directory = tempfile.mkdtemp() os.chdir(temporary_directory) if verbose: print "Using temporary working directory: %s" % temporary_directory # Copy files into temporary directory. import shutil shutil.copyfile(receptor_pdb_filename, 'receptor.pdb') shutil.copyfile(ligand_frcmod_filename, 'ligand.frcmod') shutil.copyfile(ligand_mol2_filename, 'ligand.mol2') # Assign parameters parameters = dict() parameters['protein_forcefield'] = protein_forcefield parameters['solvent_clearance'] = solvent_clearance parameters['closeness'] = 0.75 # to permit crystallographic waters parameters['receptor_pdb_filename'] = 'receptor.pdb' parameters['ligand_frcmod_filename'] = 'ligand.frcmod' parameters['ligand_mol2_filename'] = 'ligand.mol2' parameters['prefix'] = prefix if autocharge: charge = get_charge('receptor.pdb', 'ligand.frcmod', 'ligand.mol2', prefix, solvent=solvent, protein_forcefield=protein_forcefield, solvent_clearance=solvent_clearance, ncations=ncations, nanions=nanions, otherpdb=otherpdb) # Round charge. charge = int(round(charge)) parameters['ncations'] = max(0, -charge) parameters['nanions'] = max(0, charge) else: # TODO: Automate charge / counterion determination. parameters['ncations'] = ncations # number of cations to add parameters['nanions'] = nanions # number of anions to add # Add other PDB file to receptor if specified. parameters['receptor_extra'] = '' if otherpdb: parameters['receptor_extra'] = "other = loadPdb %s ; receptor = combine { receptor other }" % otherpdb print parameters # Write LEaP input to temporary file. tleap_input_filename = 'tleap.in' contents = leap_templates[solvent] % parameters write_file(tleap_input_filename, contents) # Run LEaP. if verbose: print "Running LEaP..." import commands command = 'rm -f leap.log ; tleap -f %s' % tleap_input_filename output = commands.getoutput(command) if verbose: print output if leap_output_filename: write_file(leap_output_filename, output) # Generate PDB files. if verbose: print "Generating PDB files..." for system in ['receptor', 'ligand', 'complex']: command = 'cat %(system)s.inpcrd | ambpdb -p %(system)s.prmtop > %(system)s.pdb' % vars() print command output = commands.getoutput(command) print output # Copy files out. import shutil for system in ['receptor', 'ligand', 'complex']: for suffix in ['inpcrd', 'prmtop', 'pdb']: filename = '%s.%s' % (system, suffix) shutil.copyfile(filename, prefix + '.' + filename) # Change back to original directory. os.chdir(original_directory) # Clean up temporary directory. # for filename in os.listdir(temporary_directory): # os.unlink(os.path.join(temporary_directory, filename)) # os.removedirs(temporary_directory) if verbose: print "Done." print "" return def get_charge(receptor_pdb_filename, ligand_frcmod_filename, ligand_mol2_filename, prefix, solvent='explicit', protein_forcefield='ff99SBildn', solvent_clearance=10.0, ncations=0, nanions=0, otherpdb=None): """ Get the net charge of an AMBER system to be set up. ARGUMENTS receptor_pdb_filename (string) - source receptor filename (which must have crystallographic waters included already, if desired) ligand_frcmod_filename (string) - ligand frcmod filename ligand_mol2_filename (string) - ligand GAFF mol2 filename prefix (string) - prefix for writing out AMBER prmtop/inpcrd files (prefix-{receptor,ligand,complex}.{prmtop,inpcrd}) OPTIONAL ARGUMENTS solvent (string) - either 'explicit' or 'implicit' (default: 'explicit') otherpdb (string) - if specified, add these (e.g. crystallographic waters) to receptor (default: None) TODO * Make this more general. * Automatically determine number of anions and cations based on system charge. """ print "Constructing parameters..." parameters = dict() parameters['protein_forcefield'] = protein_forcefield parameters['solvent_clearance'] = solvent_clearance parameters['closeness'] = 0.75 # to permit crystallographic waters parameters['receptor_pdb_filename'] = receptor_pdb_filename parameters['ligand_frcmod_filename'] = ligand_frcmod_filename parameters['ligand_mol2_filename'] = ligand_mol2_filename parameters['prefix'] = prefix # TODO: Automate charge / counterion determination. parameters['ncations'] = ncations # number of cations to add parameters['nanions'] = nanions # number of anions to add # Add other PDB file to receptor if specified. parameters['receptor_extra'] = '' if otherpdb: parameters['receptor_extra'] = "other = loadPdb %s ; receptor = combine { receptor other }" % otherpdb print parameters # Write LEaP input to temporary file. print "Writing tleap.in..." tleap_input_filename = 'tleap.in' contents = leap_templates["getcharge"] % parameters write_file(tleap_input_filename, contents) # Run LEaP. print "Running LEaP..." import commands command = 'rm -f leap.log ; tleap -f %s' % tleap_input_filename output = commands.getoutput(command) print output # Get charge from output. import re match = re.search('Total unperturbed charge:\s+(-?\d+\.\d+)', output) charge = float(match.group(1)) return charge #=============================================================================== # MAIN #=============================================================================== if __name__ == '__main__': # Test on some data. receptor_pdb_filename = 'src-reference-model.pdb' crystallographic_waters_pdb_filename = '../modeller/closewater.pdb' ligand_frcmod_filename = '../antechamber/oechem/bosutinib.gaff.frcmod' ligand_mol2_filename = '../antechamber/oechem/bosutinib.gaff.mol2' prefix = 'test' #ncations = 7 #nanions = 0 #create_system(receptor_pdb_filename, ligand_frcmod_filename, ligand_mol2_filename, prefix, ncations=ncations, nanions=nanions, otherpdb=crystallographic_waters_pdb_filename) create_system(receptor_pdb_filename, ligand_frcmod_filename, ligand_mol2_filename, prefix, otherpdb=crystallographic_waters_pdb_filename, autocharge=True)