some updates

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")

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()

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)
 

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')
+
 
 
 

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'):

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)]

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)
+
 
 
 

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]