diff --git a/evaluators.py b/evaluators.py index b2ba470..26dcc14 100755 --- a/evaluators.py +++ b/evaluators.py @@ -312,7 +312,7 @@ def amber_energy_minimize(settings, individual): obconv.SetOutFormat("pdb") obconv.WriteFile(individual.mol, directory + "/mol.pdb") - op.runtleap(work_dir=directory + "/", mol_file='mol.pdb', tleaptemp=settings.tleap_template, tleapin="leap.in") + op.runtleap(work_dir=directory + "/", mol_file='mol2.pdb', tleaptemp=settings.tleap_template, tleapin="leap.in") op.runPMEMD(work_dir=directory + "/", np=settings.mpi_procs, amberin=settings.amber_params) @@ -415,7 +415,7 @@ def openmm_energy_minimize(settings, individual): obconv.SetOutFormat("pdb") obconv.WriteFile(individual.mol, directory + "/mol.pdb") - op.runtleap(work_dir=directory + "/", mol_file='mol.pdb', tleaptemp=settings.tleap_template, tleapin="leap.in") + op.runtleap(work_dir=directory + "/", mol_file='mol_amber.pdb', tleaptemp=settings.tleap_template, tleapin="leap.in") tleapfiles = load_file(os.getcwd()+"/amber_run/mol.prmtop",os.getcwd()+"/amber_run/mol.inpcrd") system = tleapfiles.createSystem(nonbondedMethod=app.NoCutoff,implicitSolvent=app.OBC2,) @@ -452,10 +452,10 @@ def openmm_energy_minimize(settings, individual): integrator = openmmtools.integrators.GradientDescentMinimizationIntegrator(initial_step_size=0.04 * u.angstroms) sim = app.Simulation(tleapfiles.topology, system, integrator, platform, prop) sim.context.setPositions(tleapfiles.positions) - sim.minimizeEnergy(maxIterations=1000) - sim.reporters.append(StateDataReporter(os.devnull, 1000, step=True, potentialEnergy=True,kineticEnergy=True, temperature=True)) - sim.reporters.append(app.PDBReporter(directory +'/min.pdb', 1000)) - sim.step(1000) + #sim.minimizeEnergy(maxIterations=1000) + sim.reporters.append(StateDataReporter(os.devnull, 15000, step=True, potentialEnergy=True,kineticEnergy=True, temperature=True)) + sim.reporters.append(app.PDBReporter(directory +'/min.pdb', 15000)) + sim.step(15000) state = sim.context.getState(getEnergy=True,getPositions=True) finalEnergy=state.getPotentialEnergy().value_in_unit(u.kilocalories_per_mole) @@ -608,7 +608,6 @@ def helical_stability(settings, individuals, fitness_index, pop_start=0): #add += constants.energies['ALA'][settings.helical_dielectric] # TODO this won't work for a non poly ala protein!! #negate += constants.energies[res][settings.helical_dielectric] res = mi.getResType(individuals[i].mol.GetResidue(settings.composition_residue_indexes[j])) - add += constants.energies[str(settings.originalResidues[j])][settings.helical_dielectric] negate += constants.energies[res][settings.helical_dielectric] if settings.entropy_correction: @@ -844,13 +843,10 @@ def mmpbsa_multi(settings, individuals, fitness_index, pop_start=0): op.runtleapMultiInd(work_dir=work_dir, leap_inputfile=settings.tleap_template, n_frames=settings.no_frames) # run pmemd using multipmemd version. Need to use pmemd.MPI because of force truncation in pmemd.cuda! minReturnCode = op.runPMEMDMultiInd(work_dir=work_dir, np=settings.mpi_procs, amberin=settings.amber_params, - n_frames=settings.no_frames, compute_cluster=settings.use_compute_cluster, - nodes=settings.compute_cluster_nodes, - ntasks=settings.compute_cluster_ntasks, - queuename=settings.compute_cluster_queuename) + n_frames=settings.no_frames) logfile = open(work_dir + "/evaluator_mmpbsa_multi.log", "a") if minReturnCode == 0: - print("Setting indivdual fitness to 0 due to error in the minimization (Seg fault)") + print("Setting individual fitness to 0 due to error in the minimization (Seg fault)") finalEnergy = 0 pass else: @@ -860,10 +856,6 @@ def mmpbsa_multi(settings, individuals, fitness_index, pop_start=0): # Do mmpbsa calculation after frames have been minimized and converted to a solvent free pdb file. MMPBSA calculation can run on as many cores as there are frames. finalEnergy = op.runMMPBSAMultiInd(work_dir=work_dir, mmpbsa_inputfile=settings.mmpbsa_params, n_frames=settings.no_frames, - compute_cluster=settings.use_compute_cluster, - nodes=settings.compute_cluster_nodes, - ntasks=settings.compute_cluster_ntasks, - queuename=settings.compute_cluster_queuename, energy_evaluator=settings.energy_calculator) logfile = open(work_dir + "/evaluator_mmpbsa_multi.log", "a") logfile.write("finished MMPBSA min\n") diff --git a/main.py b/main.py index 4bb1282..1636931 100755 --- a/main.py +++ b/main.py @@ -66,6 +66,7 @@ def mainLoop(settings): parent_population = [] ga = initialise_ga(settings) + # if ("helical_stability" in settings.evaluators): # print('Fitness of {}:'.format(settings.initial_molecule_path)) @@ -418,6 +419,7 @@ def printIterationInfo(settings, curr_iteration, pop, ending): print("Best ind - Fitness(es):", pop[min_indiv].fitnesses) if settings.write_energies: + print("writing energies to file") op.writeEnergyLog(path=settings.output_path, energies=pop[min_indiv].energies, iteration=curr_iteration) obConversion = openbabel.OBConversion() diff --git a/src/JobConfig.py b/src/JobConfig.py index e31490a..70e39ab 100755 --- a/src/JobConfig.py +++ b/src/JobConfig.py @@ -272,6 +272,7 @@ def __init__(self, configFilePath): self.update_files = True if "stability_multi" in self.evaluators or 'helical_stability' in self.evaluators: self.initial_fitness_computation = True + self.write_energies = self.config.getboolean('EVALUATOR', 'write_energies', fallback=False) if "helical_stability" in self.evaluators: self.tleap_template = self.config.get('EVALUATOR', 'tleap_template') self.amber_params = self.config.get('EVALUATOR', 'amber_params') @@ -283,7 +284,7 @@ def __init__(self, configFilePath): self.gpudeviceindex = self.config.get('EVALUATOR', 'gpudeviceindex', fallback=0) self.simAnneal = self.config.getboolean('EVALUATOR', 'simAnneal', fallback=False) - self.write_energies = self.config.getboolean('EVALUATOR', 'writeEnergies', fallback=False) + self.write_energies = self.config.getboolean('EVALUATOR', 'write_energies', fallback=False) self.entropy_correction = self.config.get('EVALUATOR', 'entropy_correction', fallback=None) diff --git a/src/MoleculeCreator.py b/src/MoleculeCreator.py index c53f304..3b6290f 100755 --- a/src/MoleculeCreator.py +++ b/src/MoleculeCreator.py @@ -173,18 +173,25 @@ def swapsidechain(settings, mol, res_index, aa_mol): mol.BeginModify() aa_mol.BeginModify() - curr = mol.GetResidue(res_index) - new = aa_mol.GetResidue(0) + #mol.SetChainsPerceived(False) + + curr = mol.GetResidue(res_index) + + new = aa_mol.GetResidue(0) mol_CA = mi.getAlphaCarbon(curr) mol_CB = mi.getBetaAtom(curr) mol_N = mi.getBBNitrogen(curr) + mol_NplusOne = mi.getPlus1BBNitrogen(curr) + mol_C = mi.getBBCarboxyl(curr) aa_CA = mi.getAlphaCarbon(new) aa_CB = mi.getBetaAtom(new) aa_bb_nitrogen = mi.getBBNitrogen(new) aa_gamma_atom = mi.getChi1DihedralAtom(new) + aa_bb_oxygen = mi.getBBCarboxyl(new) + frag_res = new frag_name = frag_res.GetName() @@ -193,12 +200,22 @@ def swapsidechain(settings, mol, res_index, aa_mol): res_name =mi.getResType(curr) aa_tor = 0 + aa_tor2 = 0 if (aa_gamma_atom is not None): aa_tor = aa_mol.GetTorsion(aa_gamma_atom, aa_CA, aa_CB, aa_bb_nitrogen) + - mol_cb_vec = np.asarray([mol_CB.GetX(), mol_CB.GetY(), mol_CB.GetZ()]) - mol_ca_vec = np.asarray([mol_CA.GetX(), mol_CA.GetY(), mol_CA.GetZ()]) + psi_tor = mol.GetTorsion(mol_N, mol_CA, mol_C, mol_NplusOne) + + + try: + mol_cb_vec = np.asarray([mol_CB.GetX(), mol_CB.GetY(), mol_CB.GetZ()]) + mol_ca_vec = np.asarray([mol_CA.GetX(), mol_CA.GetY(), mol_CA.GetZ()]) + except: + print('could not get coordinates from the molecule. EVOLVE needs a protonated structure') + exit() + frag_cb_vec = np.asarray([aa_CB.GetX(), aa_CB.GetY(), aa_CB.GetZ()]) frag_ca_vec = np.asarray([aa_CA.GetX(), aa_CA.GetY(), aa_CA.GetZ()]) @@ -350,17 +367,16 @@ def swapsidechain(settings, mol, res_index, aa_mol): if curr.GetAtomID(obatom) in ["HE3", "HE4", "HE5", "HE6", "HE7", "HE8"]: curr.SetAtomID(obatom, "HE1") - + mol.SetAromaticPerceived(False) # Fix to get around openbabel screwing up the structure when calling EndModify() # This is necessary because we need to set the Chi1 dihedral correctly obConversion = openbabel.OBConversion() obConversion.SetInAndOutFormats("pdb", "pdb") obConversion.WriteFile(mol, "temp.pdb") - obConversion = openbabel.OBConversion() obConversion.SetInAndOutFormats("pdb", "pdb") obConversion.ReadFile(mol, "temp.pdb") - os.remove('temp.pdb') + # os.remove('temp.pdb') @@ -371,10 +387,23 @@ def swapsidechain(settings, mol, res_index, aa_mol): bb_nitrogen = mi.getBBNitrogen(curr) beta_atom = mi.getBetaAtom(curr) chi_atom = mi.getChi1DihedralAtom(curr) - + bb_carbonyl = mi.getBBCarboxyl(curr) + bb_nitrogenPlusOne = mi.getPlus1BBNitrogen(curr) + + psi_tor_before = mol.GetTorsion(bb_nitrogen, alpha_carbon, bb_carbonyl, bb_nitrogenPlusOne) + # print('BB oxygen', aa_mol.GetTorsion(aa_gamma_atom, aa_CA, aa_CB, aa_bb_oxygen)) + + ## Here the correct dihedral needs to be set! - if (chi_atom is not None and aa_gamma_atom is not None): - mol.SetTorsion(bb_nitrogen, alpha_carbon, beta_atom, chi_atom, aa_tor * (old_div(np.pi, 180.0))) + if (chi_atom is not None and aa_gamma_atom is not None ): + # mol.SetTorsion(bb_carbonyl, alpha_carbon, beta_atom, chi_atom, aa_tor2 * np.pi/180.0) + mol.SetTorsion(bb_nitrogen, alpha_carbon, beta_atom, chi_atom, aa_tor * np.pi/180.0) + # mol.SetTorsion(bb_nitrogen, alpha_carbon, bb_carbonyl, bb_nitrogenPlusOne, psi_tor*-180* np.pi/180.0) + psi_tor_after = mol.GetTorsion(bb_nitrogen, alpha_carbon, bb_carbonyl, bb_nitrogenPlusOne) + + if psi_tor_before != psi_tor_after: + print('something went wrong for this individual') + diff --git a/src/MoleculeInfo.py b/src/MoleculeInfo.py index 3d095b6..9b13339 100755 --- a/src/MoleculeInfo.py +++ b/src/MoleculeInfo.py @@ -122,7 +122,6 @@ def getAlphaCarbon(obres): if (carboxl_carbon is not None): # now also find the nitrogen amide_nitro = getConnectedAmideNitrogen(obatom) - if (amide_nitro is not None): return obatom return None @@ -191,7 +190,7 @@ def getBetaAtom(obres): bbcarboxyl = getBBCarboxyl(obres) bbnitrogen = getBBNitrogen(obres) - + if alpha_carbon is None: return None @@ -206,7 +205,6 @@ def getBetaAtom(obres): carboxyl_vec = np.asarray([bbcarboxyl.GetX(), bbcarboxyl.GetY(), bbcarboxyl.GetZ()]) res = getResType(obres) - for obatom in openbabel.OBAtomAtomIter(alpha_carbon): if (res == "GLY"): if (obatom.GetType() == 'H'): diff --git a/src/constants.py b/src/constants.py index 64a951d..b418b30 100755 --- a/src/constants.py +++ b/src/constants.py @@ -59,7 +59,7 @@ 'M01', 'M02', 'M03', 'M04', 'M05', 'M06', 'M07', 'M08', 'M09', 'M10', \ 'M11', 'M12', 'M13', 'P01', 'P02', 'P03', 'P04', 'S01', 'S02', 'S03', \ 'T01', 'T02', 'T03', 'W01', 'W02', 'W03', 'W04', 'W05', 'W06', 'W07', \ - 'Y01', 'Y02', 'Y03', 'Y04', 'V01', 'V02', 'V03' + 'Y01', 'Y02', 'Y03', 'Y04', 'V01', 'V02', 'V03', 'HIN' ] selected_rotamers = rotamers @@ -82,11 +82,11 @@ 'MET', 'MET', 'MET', 'MET', 'MET', 'MET', 'MET', 'MET', 'MET', 'MET', \ 'MET', 'MET', 'MET', 'PHE', 'PHE', 'PHE', 'PHE', 'SER', 'SER', 'SER', \ 'THR', 'THR', 'THR', 'TRP', 'TRP', 'TRP', 'TRP', 'TRP', 'TRP', 'TRP', \ - 'TYR', 'TYR', 'TYR', 'TYR', 'VAL', 'VAL', 'VAL' + 'TYR', 'TYR', 'TYR', 'TYR', 'VAL', 'VAL', 'VAL', 'HIN' ] selected_rotamer_types = rotamer_types -allowed_residue_types = ['GLY', 'ALA', 'ARG', 'ASP', 'ASH', 'ASN', 'CYS', 'GLU', 'GLH', 'GLN', 'HID', 'HIE', 'HIP', 'ILE', 'LEU', 'LYS', 'LYN', 'MET', 'PHE', 'SER', 'THR', 'TRP', 'TYR', 'VAL'] +allowed_residue_types = ['GLY', 'ALA', 'ARG', 'ASP', 'ASH', 'ASN', 'CYS', 'GLU', 'GLH', 'GLN', 'HID', 'HIE', 'HIP', 'ILE', 'LEU', 'LYS', 'LYN', 'MET', 'PHE', 'SER', 'THR', 'TRP', 'TYR', 'VAL', 'HIN'] # Vacic et al, BMC Bioinformatics 211(8), 2007, https://doi.org/10.1186/1471-2105-8-211 swissprot_probabilities = {'ALA': 0.0789, 'ARG': 0.054, 'ASH': 0.02675, 'ASP': 0.02675, 'ASN': 0.0413, \ @@ -99,7 +99,7 @@ # 10, 50, 80 energies = {'ALA': (13.03, 12.14 , 12.06), 'ARG': (-170.77, -176.22 , -176.73), 'ASH': (-43.74, -45.64, -45.82), 'ASP': (-52.27, -58.92, -59.55), 'ASN': (-72.72, -74.29, -74.44), \ 'CYS': (13.55, 12.58, 12.49), 'GLH': (-36.86, -38.93, -39.13), 'GLU': (-47.00, -53.61, -54.30), 'GLN': (-56.93, -58.79, -58.97), 'GLY':(7.11, 6.15, 6.06), \ - 'HID': (12.43, 10.83, 10.68), 'HIE':(8.74, 7.22, 7.08), 'HIP': (22.66, 17.17, 16.62), 'ILE': (10.74, 9.93, 9.85), 'LEU': (-10.04, -10.88, -10.96), \ + 'HIN': (12.43, 10.83, 10.68), 'HID': (12.43, 10.83, 10.68), 'HIE':(8.74, 7.22, 7.08), 'HIP': (22.66, 17.17, 16.62), 'ILE': (10.74, 9.93, 9.85), 'LEU': (-10.04, -10.88, -10.96), \ 'LYN': (-4.36, -5.73, -5.85), 'LYS': (-9.60, -15.72, -16.29), 'MET': (8.36, 7.47, 7.29), 'PHE': (11.05, 10.04, 9.95), 'SER': (-0.71, -2.32, -2.47), \ 'THR': (-19.37, -20.45, -20.59), 'TRP':(14.43, 13.07, 12.94), 'TYR': (-13.28, -14.89, -15.05), 'VAL':(-7.58, -8.41, -8.49) } @@ -108,7 +108,7 @@ 'GLN':(-2.11,-2.02,-1.29), 'GLU':(-1.81,-1.65,-1.31), 'GLY':(0,0,0), 'HIS':(-0.96,-0.99,-0.92), 'ILE':(-0.89,-0.75,-0.94),\ 'LEU':(-0.78,-0.75,-0.94), 'LYS':(-1.94,-2.21,-1.63), 'MET':(-1.61,-1.53,-1.24), 'PHE':(-0.58,-0.58,-0.65), 'PRO':(0,0,-0.3),\ 'SER':(-1.71,-1.19,-0.43), 'THR':(-1.63,-1.12,-0.57), 'TRP':(-0.97,-0.97,-1.14), 'TYR':(-0.98,-0.99,-1.07),\ - 'VAL':(-0.51,-0.5,-0.62),'HID':(-0.96,-0.99,-0.92),'HIP':(-0.96,-0.99,-0.92),'HIE':(-0.96,-0.99,-0.92),\ + 'VAL':(-0.51,-0.5,-0.62), 'HIN': (-0.96,-0.99,-0.92),'HID':(-0.96,-0.99,-0.92),'HIP':(-0.96,-0.99,-0.92),'HIE':(-0.96,-0.99,-0.92),\ 'ASH':(-1.25,-0.61,-0.65),'LYN':(-1.94,-2.21,-1.63),'GLH':(-1.81,-1.65,-1.31),} def subselect_rotamers(type_list): indices = [i for i, x in enumerate(rotamer_types) if (x in type_list)] diff --git a/src/gaapi/Individual.py b/src/gaapi/Individual.py index 35d7a85..55ee559 100755 --- a/src/gaapi/Individual.py +++ b/src/gaapi/Individual.py @@ -307,6 +307,8 @@ def applyComposition(self, settings): return print(self.composition_residue_indexes) cmp = [] + if (settings.verbose): + print("Created Individual:") for i in range(0, len(self.composition_residue_indexes)): if (settings.multi_individual): for x in range(0, settings.no_frames): @@ -318,7 +320,10 @@ def applyComposition(self, settings): rotamer_type = constants.rotamers[self.composition[i]] obConversion.ReadFile(frag, settings.composition_library + "/" + rotamer_type + ".mol2") - mcr.swapsidechain(settings, self.mol[x], self.composition_residue_indexes[i], frag) + try: + mcr.swapsidechain(settings, self.mol[x], self.composition_residue_indexes[i], frag) + except AttributeError as e: + print(f"{i} went wrong, ignore it") pass cmp.append(rotamer_type) # append only after loop because all individuals are identical else: @@ -328,11 +333,12 @@ def applyComposition(self, settings): rotamer_type = constants.selected_rotamers[self.composition[i]] cmp.append(rotamer_type) + if (settings.verbose): + print(cmp) obConversion.ReadFile(frag, settings.composition_library + "/" + rotamer_type + ".mol2") mcr.swapsidechain(settings, self.mol, self.composition_residue_indexes[i], frag) - if (settings.verbose): - print("Created Individual:", cmp) + diff --git a/src/outprocesses.py b/src/outprocesses.py index cc9427f..4b46820 100644 --- a/src/outprocesses.py +++ b/src/outprocesses.py @@ -21,7 +21,7 @@ import time -def runtleap(work_dir="", mol_file="mol.pdb", tleaptemp="tleap_template.in", tleapin="leap.in", inpcrd_out="mol.inpcrd", prmtop_out="mol.prmtop"): +def runtleap(work_dir="", mol_file="mol_amber.pdb", tleaptemp="tleap_template.in", tleapin="leap.in", inpcrd_out="mol.inpcrd", prmtop_out="mol.prmtop"): ''' runs tleap to process a pdb file into a valid amber prmtop file. @@ -49,7 +49,24 @@ def runtleap(work_dir="", mol_file="mol.pdb", tleaptemp="tleap_template.in", tle os.remove(work_dir + inpcrd_out) if (os.path.exists(work_dir + prmtop_out)): os.remove(work_dir + prmtop_out) - + + try: + f = open(work_dir+mol_file, "w") + FNULL = open(os.devnull, 'w') + pPDB4Amber = subprocess.Popen(["pdb4amber", work_dir+'mol.pdb', ], universal_newlines=True, stdout=f, + stderr=FNULL) + pPDB4Amber.wait() + print('running pdb4amber') + FNULL.close() + f.close() + except IOError as e: + sys.exit("I/O error on '%s': %s" % (e.filename, e.strerror)) + except subprocess.CalledProcessError as e: + sys.exit("pdb4amber failed processing mutated.pdb, returned code %d " % (e.returncode)) + except OSError as e: + sys.exit("failed to execute pdb4amber: %s" % (str(e))) + pass + f = open(tleaptemp, "r") g = f.readlines() filename = (work_dir + mol_file).split('.')[0]