changed openbabel legacy functions and made modifactions to indivdual…

…s and openmm minimization routines

.gitignore

@@ -1,2 +1,3 @@
 *.pyc 
 *.bak
+.idea

evaluators.py

@@ -376,7 +376,7 @@ def openmm_energy_minimize(settings, individual):
         dielectric = 80.0
         gpu_openmm = True
         gpudeviceindex = 0
-        simAnneal = True 
+        simAnneal = False
         md = False
 
     Parameters
@@ -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.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)
+        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)
         state = sim.context.getState(getEnergy=True,getPositions=True)
         finalEnergy=state.getPotentialEnergy().value_in_unit(u.kilocalories_per_mole)
 
@@ -585,10 +585,15 @@ def helical_stability(settings, individuals, fitness_index, pop_start=0):
             print("Already computed: " , i, " -> member ", already_done)
             individuals[i].mol = openbabel.OBMol(individuals[already_done].mol)
             individuals[i].fitnesses = individuals[already_done].fitnesses
+            individuals[i].energies = individuals[already_done].energies
             continue     
 
         add = 0.0
-        negate = settings.initial_energy
+        negate = 0.0
+
+        #Entropy corrections
+        entropy_add=0.0
+        entropy_negate=0.0
 
         print("Minimising: ", i, [mi.getResType(individuals[i].mol.GetResidue(j)) for j in settings.composition_residue_indexes])
         print("Rotamers: ", [cnts.selected_rotamers[v] for v in individuals[i].composition])
@@ -606,12 +611,29 @@ def helical_stability(settings, individuals, fitness_index, pop_start=0):
 
             add += constants.energies[str(settings.originalResidues[j])][settings.helical_dielectric]
             negate += constants.energies[res][settings.helical_dielectric]
-            print('add',str(settings.originalResidues[j]), constants.energies[str(settings.originalResidues[j])][settings.helical_dielectric])
-            print('negate',res,'-(', settings.initial_energy,'+', constants.energies[res][settings.helical_dielectric],')')
-
+            if settings.entropy_correction:
+                entropy_add+=constants.entropy_corrections[str(settings.originalResidues[j])][settings.entropy_correction]
+                entropy_negate+=constants.entropy_corrections[res][settings.entropy_correction]
+
+            #print('add',str(settings.originalResidues[j]), constants.energies[str(settings.originalResidues[j])][settings.helical_dielectric])
+            #print('negate',res,'-(', settings.initial_energy,'+', constants.energies[res][settings.helical_dielectric],')')
+
+        if settings.write_energies:
+            individuals[i].setEnergies(settings.initial_energy, add, finalEnergy, negate)
+
+        negate +=settings.initial_energy
+
         print("min_energy:", finalEnergy, add, negate)
         finalEnergy += (add - negate)
 
+        if settings.entropy_correction:
+            print("Entropy Correction: ", (entropy_add - entropy_negate))
+            finalEnergy+=  (entropy_add - entropy_negate)
+
+
+
+
+
         individuals[i].setFitness(fitness_index, finalEnergy)
 
         if (settings.solution_min_fitness is not None):

main.py

@@ -29,7 +29,7 @@
 from src.gaapi import generators
 from src.gaapi import replacers
 from src import MoleculeInfo as mi
-
+import src.outprocesses as op
 from collections import Counter 
 import evaluators as evals
 # import  printResidueInfo as pri
@@ -416,6 +416,9 @@ def printIterationInfo(settings, curr_iteration, pop, ending):
         print(pop[min_indiv].phi_dihedrals, pop[min_indiv].psi_dihedrals)
 
     print("Best ind - Fitness(es):", pop[min_indiv].fitnesses)
+
+    if settings.write_energies:
+        op.writeEnergyLog(path=settings.output_path, energies=pop[min_indiv].energies, iteration=curr_iteration)
 
     obConversion = openbabel.OBConversion()
     obConversion.SetOutFormat("pdb")

src/JobConfig.py

@@ -282,6 +282,18 @@ def __init__(self, configFilePath):
                 self.replicas = self.config.get('EVALUATOR', 'replicas', fallback=1)
                 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.entropy_correction = self.config.get('EVALUATOR', 'entropy_correction', fallback=None)
+
+                if (self.entropy_correction == 'PicketSternberg'):
+                    self.entropy_correction = 0
+                elif (self.entropy_correction == "AbagyanTotrov"):
+                    self.entropy_correction = 1
+                elif (self.entropy_correction == "KoehlDelarue"):
+                    self.entropy_correction = 2
+
                 if (self.helical_dielectric == 80):
                     self.helical_dielectric = 2
                 elif (self.helical_dielectric == 50):

src/MoleculeCreator.py

@@ -16,6 +16,11 @@
 from openbabel import OBAtom
 from openbabel import OBMol
 from openbabel import OBAtomTyper
+# from openbabel import openbabel
+# from openbabel.openbabel import OBResidue
+# from openbabel.openbabel import OBAtom
+#from openbabel.openbabel import OBMol
+#from openbabel.openbabel import OBAtomTyper
 import numpy as np
 from . import MoleculeInfo as mi
 
@@ -115,9 +120,9 @@ def add_fragment(mol, fragment, id_start):
         new_atom.SetType(frag_atom.GetType())
         new_atom.SetSpinMultiplicity(frag_atom.GetSpinMultiplicity())
         new_atom.SetFormalCharge(frag_atom.GetFormalCharge())
-        new_atom.SetImplicitValence(frag_atom.GetImplicitValence())
+        # new_atom.SetImplicitValence(frag_atom.GetImplicitValence())
         new_atom.SetPartialCharge(frag_atom.GetPartialCharge())
-
+        
         for data in frag_atom.GetData():
             new_atom.SetData(data)
 
@@ -126,15 +131,15 @@ def add_fragment(mol, fragment, id_start):
     for obbond in openbabel.OBMolBondIter(fragment):
         begin_atom = obbond.GetBeginAtom()
         end_atom = obbond.GetEndAtom()
-        mol.AddBond(obbond.GetBeginAtomIdx() + prevAtoms, obbond.GetEndAtomIdx() + prevAtoms, obbond.GetBO(), obbond.GetFlags())
+        mol.AddBond(obbond.GetBeginAtomIdx() + prevAtoms, obbond.GetEndAtomIdx() + prevAtoms, obbond.GetBondOrder(), obbond.GetFlags())
 
 
 def debugAtom(obatom):
     return "AN:", obatom.GetAtomicNum(), "T:", obatom.GetType(), "idx:", obatom.GetIdx(), "id:", obatom.GetId(), "BO(1):", obatom.CountBondsOfOrder(1)
 
 
 def debugBond(obbond):
-    return "begin_idx:", obbond.GetBeginAtomIdx(), "end_idx:", obbond.GetEndAtomIdx(), "norm:", obbond.GetLength(), "BO:", obbond.GetBO()
+    return "begin_idx:", obbond.GetBeginAtomIdx(), "end_idx:", obbond.GetEndAtomIdx(), "norm:", obbond.GetLength(), "BO:", obbond.GetBondOrder()
 
 
 def getAtomByID(mol, atom_id):
@@ -165,13 +170,14 @@ def swapsidechain(settings, mol, res_index, aa_mol):
     mol.BeginModify()
 
     curr = mol.GetResidue(res_index)
-
-    new = aa_mol.GetResidue(0)
 
+    new = aa_mol.GetResidue(0)
+    print(mol, mol.NumResidues())
+    print(mol.GetResidue(4))
     mol_CA = mi.getAlphaCarbon(curr)
     mol_CB = mi.getBetaAtom(curr)
     mol_N = mi.getBBNitrogen(curr)
-   
+
     aa_CA = mi.getAlphaCarbon(new)
     aa_CB = mi.getBetaAtom(new)
     aa_bb_nitrogen = mi.getBBNitrogen(new)
@@ -296,12 +302,13 @@ def swapsidechain(settings, mol, res_index, aa_mol):
     # if old resname = glycine we need to rename the remaining H to HA, the other has been replaced by the sidechain
     if res_name=="GLY":
         for obatom in openbabel.OBResidueAtomIter(curr):
-            if curr.GetAtomID(obatom)=="HA2" or curr.GetAtomID(obatom)=="HA3":
+            if curr.GetAtomID(obatom) in ["HA2", "HA3", " HA2", " HA2 "]:
                 curr.SetAtomID(obatom, "HA")
 
     if frag_res_type=="GLY" and settings.use_res_type:
         for obatom in openbabel.OBResidueAtomIter(curr):
-            if curr.GetAtomID(obatom)=="HA":
+            #print("."+curr.GetAtomID(obatom)+".")
+            if curr.GetAtomID(obatom) in ["HA", " HA", " HA "] :
                 curr.SetAtomID(obatom, "HA2")
 
     if frag_res_type=="HID" and settings.use_res_type:

src/MoleculeInfo.py

@@ -15,6 +15,9 @@
 import openbabel
 from openbabel import OBResidue
 from openbabel import OBAtom
+# from openbabel import openbabel
+# from openbabel.openbabel import OBResidue
+# from openbabel.openbabel import OBAtom
 import numpy as np
 from . import constants
 
@@ -111,12 +114,10 @@ def getAlphaCarbon(obres):
         Query residue
     """
     res = obres.GetName()
-
     for obatom in openbabel.OBResidueAtomIter(obres):
-        
-        if (obatom.IsCarbon()):
+
+        if (obatom.GetAtomicNum()==6):
             carboxl_carbon = getConnectedCarboxylCarbon(obatom)
-
             # identify carboxyl carbon in peptide
             if (carboxl_carbon is not None):
                 # now also find the nitrogen
@@ -142,14 +143,14 @@ def getConnectedAmideNitrogen(calpha_atom):
     #    --> assume that atom is the alpha_carbon
     if (obres.GetName() == "PRO"):
         for obatom in openbabel.OBAtomAtomIter(calpha_atom):
-            if (obatom.IsNitrogen()):
+            if (obatom.GetAtomicNum()==7):
                 return obatom
 
     else:
         for obatom in openbabel.OBAtomAtomIter(calpha_atom):
-            if (obatom.IsNitrogen()):
+            if (obatom.GetAtomicNum()==7):
                 for obatom2 in openbabel.OBAtomAtomIter(obatom):
-                    if (obatom2.IsHydrogen()):
+                    if (obatom2.GetAtomicNum()==1):
                         return obatom
     return None
 
@@ -166,7 +167,7 @@ def getConnectedCarboxylCarbon(atom):
     obres = atom.GetResidue()
 
     for obatom in openbabel.OBAtomAtomIter(atom):
-        if (obatom.IsCarbon()):
+        if (obatom.GetAtomicNum()==6):
 
             for obbond in openbabel.OBAtomBondIter(obatom):
                 if (countBonds(obbond.GetNbrAtom(obatom)) == 1 and obbond.GetNbrAtom(obatom).GetAtomicNum() == 8):
@@ -221,7 +222,7 @@ def getBetaAtom(obres):
                     return obatom
 
         else:
-            if (bbcarboxyl is not None and obatom.IsCarbon() and obatom.GetIdx() != bbcarboxyl.GetIdx()):
+            if (bbcarboxyl is not None and obatom.GetAtomicNum()==6 and obatom.GetIdx() != bbcarboxyl.GetIdx()):
                 return obatom
     return None
 
@@ -242,7 +243,7 @@ def getBBNitrogen(obres):
         return None
 
     for obatom in openbabel.OBAtomAtomIter(alpha_carbon):
-        if (obatom.IsNitrogen()):
+        if (obatom.GetAtomicNum()==7):
             return obatom
     return None
 
@@ -267,7 +268,7 @@ def getNeg1BBCarboxyl(obres):
         if (obatom == ca):
             continue
 
-        if (obatom.IsHydrogen()):
+        if (obatom.GetAtomicNum()==1):
             continue
 
         for obbond in openbabel.OBAtomBondIter(obatom):
@@ -293,7 +294,7 @@ def getBBCarboxyl(obres):
 
     for obatom in openbabel.OBAtomAtomIter(ca_atom):
 
-        if (not obatom.IsCarbon()):
+        if (not obatom.GetAtomicNum()==6):
             continue
 
         for obbond in openbabel.OBAtomBondIter(obatom):
@@ -319,7 +320,7 @@ def getPlus1BBNitrogen(obres):
 
     for obatom in openbabel.OBAtomAtomIter(bb_n):
 
-        if (obatom.IsNitrogen()):
+        if (obatom.GetAtomicNum()==7):
             return obatom
     return None
 
@@ -613,13 +614,13 @@ def getChi1DihedralAtom(obres):
         if (obatom.GetIdx() == alpha_carbon.GetIdx()):
             continue
 
-        if ((res == "SER" or res == "THR") and obatom.IsOxygen()):
+        if ((res == "SER" or res == "THR") and obatom.GetAtomicNum()==8):
             return obatom
-        elif (res == "CYS" and obatom.IsSulfur()):
+        elif (res == "CYS" and obatom.GetAtomicNum()==16):
             return obatom
-        elif (res == "ALA" and obatom.IsHydrogen()):
+        elif (res == "ALA" and obatom.GetAtomicNum()==1):
             return obatom
-        elif((res != "SER" and res != "THR" and res != "CYS" and res != "ALA") and obatom.IsCarbon()):
+        elif((res != "SER" and res != "THR" and res != "CYS" and res != "ALA") and obatom.GetAtomicNum()==6):
             return obatom
 
     return None

src/gaapi/Individual.py

@@ -62,6 +62,8 @@ def initialise_constructor(self, settings):
         self.chi_angles = None
         self.chi_atoms = None
 
+        self.energies = {"wt":None, "wt_corr":None, "mut":None, "mut_corr":None}
+
 
         if (settings.backbone_dihedral_optimization):
 
@@ -121,6 +123,7 @@ def copy_constructor(self, settings, orig):
 
 
         self.fitnesses = copy.deepcopy(orig.fitnesses)
+        self.energies = orig.energies
 
     def saveMol(self, file_path):
         '''
@@ -193,8 +196,24 @@ def setStability(self,value):
         -------
 
         '''
+
         self.stability=value
 
+    def setEnergies(self,wt, wt_corr, mut, mut_corr):
+        '''
+
+        Set energies 
+
+        Parameters
+        ----------
+        value
+
+        Returns
+        -------
+
+        '''
+        self.energies.update({"wt":wt, "wt_corr":wt_corr, "mut":mut, "mut_corr":mut_corr})
+
     def updatePhiPsiDihedrals(self, settings):
         '''
 

src/outprocesses.py

@@ -433,6 +433,39 @@ def writeFrontLog(path="", front='', iteration=0):
     pass
 
 
+def writeEnergyLog(path="", energies=None, iteration=0):
+    '''
+
+    write log of energies for all min_mols as comma seperated values 
+
+    The format is energy wt, correction wt, energy_mut, correction_mut
+
+    Parameters
+    ----------
+    path: str
+        path to the current working directory, not work_dir/
+    energies: list
+        (energy wt, correction wt, energy_mut, correction_mut)
+    iteration: int
+        iteration counter of the GA
+    '''
+    filename = path + '/energy.log'
+
+    print("Writing energies", filename)
+
+    if iteration > 0:
+        append_write = 'a'  # append if already exists
+    else:
+        append_write = 'w'  # make a new file if not
+    string = ''
+    for item in energies.keys():
+        string += str(energies[item]) + ','
+    frontsfile = open(filename, append_write)
+    frontsfile.write(string + '\n')
+    frontsfile.close()
+    pass
+
+
 def writeFrontFile(path="", file='', fitness='', iteration=0, frontid=0):
     '''
 

tests/nsga2/.gitignore

@@ -0,0 +1 @@
+work_dir