Implementation of Centroid degree of freedom. Centeroid of the molecule is determined by some point in space.

Issues: creation of the molecule is sensitive to the order of applying of degrees of freedom.
Need to add more comments for creation of new DOF object for future generetions.

Author: maksimovdmitrii

examples/ga.py

@@ -5,13 +5,15 @@
 from fafoom import MoleculeDescription, Structure, selection, print_output,\
     remover_dir, set_default, file2dict
 import fafoom.run_utilities as run_util
-from utilities import atoms_positions, centreofthebox
 from visual import draw_picture
 
 
 from rdkit import Chem
 from rdkit.Chem import AllChem
 
+from measure import centroid_measure
+from utilities import sdf2xyz
+
 # Decide for restart or a simple run.
 opt = run_util.simple_or_restart()
 p_file = sys.argv[1]
@@ -31,7 +33,6 @@
 min_energy = []
 
 
-
 #***********************************************************************
 """
 Creation of the folders for valid and invalid structures. 
@@ -72,6 +73,7 @@
         print_output("New trial")
         str3d = Structure(mol)
         str3d.generate_structure()
+        print '\n{}'.format(sdf2xyz(str3d.sdf_string))
         aims_object = AimsObject(os.path.join(os.getcwd(),'adds')) #Need for creation of the input file. Does not affect the algoritm.
         if not str3d.is_geometry_valid():
             print_output("The geometry of "+str(str3d)+" is invalid. Copied to /invalid")
@@ -85,16 +87,6 @@
             os.mkdir(os.path.join(os.getcwd(),'valid',str(cnt)+'_geometry')) # creates the folder for particular structure inside th "valid" folder
             shutil.copy('geometry.in',os.path.join(os.getcwd(), 'valid', str(cnt)+'_geometry','geometry.in')) # copy input to self-titled folder
 ############            draw_picture(os.path.join(os.getcwd(), 'valid', str(cnt)+'_geometry','geometry.in'), image_write = 'yes') # Part of post-processing module. Under construction. Produce nice image with PyMol
-
-            print str3d.sdf_string
-            #~ print str3d
-#            print atoms_positions(str3d.sdf_string)
-#            print centreofthebox(atoms_positions(str3d.sdf_string))
-#            print atoms_shift(atoms_positions(str3d.sdf_string), centreofthebox(atoms_positions(str3d.sdf_string)))
-#            print centreofthebox(atoms_shift(atoms_positions(str3d.sdf_string), centreofthebox(atoms_positions(str3d.sdf_string))))
-            new = apply_coord_translation(str3d.sdf_string, atoms_shift(atoms_positions(str3d.sdf_string), centreofthebox(atoms_positions(str3d.sdf_string))))
-            str3d.sdf_string = new
-            print str3d.sdf_string
             name = "initial_%d" % (len(population))
             # Perform the local optimization
             run_util.optimize(str3d, energy_function, params, name)

examples/parameters_aims.txt

@@ -1,15 +1,15 @@
 [Molecule]
 
 smiles="CC(=O)N[C@H](C(=O)NC)C"
-optimize_torsion=True
 optimize_cistrans=True
+optimize_torsion=True
 smarts_torsion= "[*]~[!$(*#*)&!D1]-&!@[!$(*#*)&!D1]~[*]"
 smarts_cistrans= "C~[$(C=O)]-[$(NC)]~[C]"
 filter_smarts_torsion= "C~[$(C=O)]-[$(NC)]~[C]"
 rmsd_type="cartesian"
 rmsd_cutoff_uniq=0.25
 chiral=True
-
+optimize_centroid=True
 [GA settings]
 
 energy_var=0.001

fafoom/deg_of_freedom.py

@@ -26,11 +26,19 @@
     dihedral_measure,
     dihedral_set,
     pyranosering_measure,
-    pyranosering_set
+    pyranosering_set,
+    centroid_measure,
+    centroid_set
 )
 
 from genetic_operations import mutation
 
+from rdkit import Chem
+from rdkit.Chem import AllChem
+
+from operator import itemgetter
+from rdkit.Chem import rdMolTransforms
+from utilities import *
 
 class DOF:
 
@@ -40,6 +48,77 @@ def __init__(self, name):
     def common_function():
         pass
 
+class Centroid(DOF):
+    '''Find and handle centre of the molecule. '''
+    values_options = range(-10, 10, 1)
+    @staticmethod
+    def find(smiles, positions=None):
+        if positions is None:
+            mol = Chem.MolFromSmiles(smiles)
+            if mol is None:
+                raise ValueError("The smiles is invalid")
+            pattern_cent = Chem.MolFromSmarts(smiles)
+            cent = list(mol.GetSubstructMatches(pattern_cent))
+            positions = cleaner(cent)
+        return positions
+
+    def __init__(self, positions):
+        """Initialaize the Centroid object from the positions."""
+        self.name = 'Centroid'        
+        self.type = "centroid"
+        self.positions = positions
+
+    def apply_on_string(self, string, values_to_set=None):
+        if values_to_set is not None:
+            self.values = values_to_set
+        string = centroid_set(string, self.values)
+        return string
+
+    def get_random_values(self):
+        """Generate a random value for position of the Centroid object"""
+        self.values = [choice(Centroid.values_options) for i in range(3)]
+
+    def get_weighted_values(self, weights):
+        if len(weights) == len(Centroid.values_options):
+            self.values = [Centroid.values_options[find_one_in_list(sum(
+                           weights), weights)]
+                           for i in range(len(self.positions))]
+        else:
+            self.values = [choice(Centroid.values_options)
+                           for i in range(len(self.positions))]        
+
+    def mutate_values(self, max_mutations=None, weights=None):
+
+        if max_mutations is None:
+            max_mutations = max(1, int(math.ceil(len(self.values)/2.0)))
+
+        self.values = mutation(self.values, max_mutations,
+                               Centroid.values_options, weights, periodic=False)
+
+    def update_values(self, string):
+        self.values = centroid_measure(string)
+
+    def is_equal(self, other, threshold, chiral=False):
+        values = []
+        tmp = []
+        for i in get_vec(self.values, other.values):
+            if i == 0:
+                tmp.append(0)
+            else:
+                tmp.append(1)
+        values.append(sum(tmp)/len(tmp))
+        if hasattr(other, "initial_values"):
+            tmp = []
+            for i in get_vec(self.values, other.initial_values):
+                if i == 0:
+                    tmp.append(0)
+                else:
+                    tmp.append(1)
+            values.append(sum(tmp)/len(tmp))
+        if min(values) > threshold:
+            return False
+        else:
+            return True
 
 class Torsion(DOF):
     """ Find, create and handle rotatable bonds"""
@@ -93,6 +172,7 @@ def find(smiles, smarts_torsion="[*]~[!$(*#*)&!D1]-&!@[!$(*#*)&!D1]~[*]",
 
     def __init__(self, positions):
         """Initialaize the Torsion object from the positions."""
+        self.name = 'Torsion'
         self.type = "torsion"
         self.positions = positions
 
@@ -283,6 +363,7 @@ def find(smiles, pyranosering_pattern="C1(CCCCO1)O", positions=None):
         return positions
 
     def __init__(self, positions):
+        self.name = 'PyranoseRing'
         self.type = "pyranosering"
         self.positions = positions
 
@@ -361,9 +442,10 @@ def find(smiles, smarts_cistrans=None, positions=None):
             pattern_cistrans = Chem.MolFromSmarts(smarts_cistrans)
             cistrans = list(mol.GetSubstructMatches(pattern_cistrans))
             positions = cleaner(cistrans)
-        return positions
+        return positions      
 
     def __init__(self, positions):
+        self.name = 'CisTrans'
         self.type = "cistrans"
         self.positions = positions
 
@@ -421,3 +503,33 @@ def is_equal(self, other, threshold, chiral=True):
             return False
         else:
             return True
+#=======================================================================
+#~ '''
+#~ For test of the module only
+#~ '''
+
+#~ smiles = 'CC(=O)N[C@H](C(=O)NC)C'
+#~ obj = Centroid(smiles)
+#~ mol = Chem.MolFromSmiles(obj.positions)
+#~ AllChem.EmbedMolecule(mol)
+#~ string = Chem.MolToMolBlock(mol)
+
+#~ print 'Initial coordinates:'
+#~ print sdf2xyz(string)
+#~ print 'Initial centroid:'
+#~ print centroid_measure(string)
+#~ obj.get_random_values()
+#~ print 'Centroid will be set to:'
+#~ print obj.values
+#~ string = obj.apply_on_string(string, obj.values)
+#~ print 'Final coordinates:'
+#~ print sdf2xyz(string)
+#~ print 'Final centroid:'
+#~ print centroid_measure(string)
+#***********************************************************************
+
+
+
+
+
+

fafoom/get_parameters.py

@@ -18,7 +18,7 @@
 from __future__ import division
 from rdkit import Chem
 from rdkit.Chem import AllChem
-from deg_of_freedom import Torsion, CisTrans, PyranoseRing
+from deg_of_freedom import Torsion, CisTrans, PyranoseRing, Centroid
 from utilities import check_geo_sdf
 
 
@@ -84,6 +84,13 @@ def get_positions(type_of_deg, smiles, **kwargs):
         else:
             return PyranoseRing.find(smiles)
 
+    if type_of_deg == "centroid":
+        if 'list_of_centroid' in kwargs:
+            return Centroid.find(smiles,
+                                     positions=kwargs['list_of_centroid'])
+        else:
+            return Centroid.find(smiles)
+
 
 def create_dof_object(type_of_deg, positions):
     """Initialize the degree of freedom from the positions
@@ -100,6 +107,8 @@ def create_dof_object(type_of_deg, positions):
         return CisTrans(positions)
     if type_of_deg == "pyranosering":
         return PyranoseRing(positions)
+    if type_of_deg == "centroid":
+        return Centroid(positions)
 
 #~ , distance_cutoff_1, distance_cutoff_2
 def template_sdf(smiles):

fafoom/measure.py

@@ -22,13 +22,33 @@
 from rdkit import Chem
 from rdkit.Chem import rdMolTransforms
 
-from utilities import get_vec, tor_rmsd, xyz2sdf
+from utilities import get_vec, tor_rmsd, xyz2sdf, sdf2xyz
 
 
 def ig(x):
     return itemgetter(x)
 
+def centroid_measure(sdf_string):
+    mol = Chem.MolFromMolBlock(sdf_string, removeHs=False)
+    pos = mol.GetConformer()
+    centroid = rdMolTransforms.ComputeCentroid(pos, ignoreHs=True) 
+    return [centroid.x, centroid.y, centroid.z]
 
+def centroid_set(sdf_string, values_to_set):
+    atoms_list = []
+    mol = Chem.MolFromMolBlock(sdf_string, removeHs=False)
+    for i in range(0, mol.GetNumAtoms()):
+        pos = mol.GetConformer().GetAtomPosition(i)
+        atoms_list.append([pos.x, pos.y, pos.z])
+    new_coordinates = []
+    shift = [values_to_set[0] - centroid_measure(sdf_string)[0], values_to_set[1] - centroid_measure(sdf_string)[1], values_to_set[2] - centroid_measure(sdf_string)[2]]
+    for i in range(len(atoms_list)):
+        new_coordinates.append([atoms_list[i][0] + shift[0], atoms_list[i][1] + shift[1], atoms_list[i][2] + shift[2]])
+    for i in range(0, mol.GetNumAtoms()):
+        mol.GetConformer().SetAtomPosition(i, new_coordinates[i])
+    sdf_string = Chem.MolToMolBlock(mol)
+    return sdf_string
+    
 def dihedral_measure(sdf_string, position):
     """ Measure the dihedral angle.
 

fafoom/structure.py

@@ -58,7 +58,8 @@ def __init__(self, parameter_file=None, **kwargs):
             new_names_dict = {'smile': 'smiles',
                               'smart_torsion': 'smarts_torsion',
                               'filter_smart_torsion': 'filter_smarts_torsion',
-                              'smart_cistrans': 'smarts_cistrans'}
+                              'smart_cistrans': 'smarts_cistrans',
+                              'centroid': 'centroids'}
             for key in new_names_dict:
                 try:
                     params[new_names_dict[key]] = params.pop(key)
@@ -269,11 +270,12 @@ def generate_structure(self, values={}):
                     dof.get_weighted_values(weights)
                 else:
                     dof.get_random_values()
+                    print 'Initial values for {} are {}'.format(dof.name ,dof.values) 
                 new_string = dof.apply_on_string(new_string)
         self.sdf_string = new_string
         for dof in self.dof:
             dof.update_values(self.sdf_string)
-
+            print 'Updated values for {} are {}'.format(dof.name, dof.values)
     def is_geometry_valid(self):
         """Return True if the geometry is valid."""
         #~ , self.mol_info.distance_cutoff_1, self.mol_info.distance_cutoff_2