Python rdkit.Chem.GetAdjacencyMatrix() Examples

def CalculateSchiultz(mol):
    """
    #################################################################
    Calculation of Schiultz number

    ---->Tsch(log value)

    Usage:

        result=CalculateSchiultz(mol)

        Input: mol is a molecule object

        Output: result is a numeric value
    #################################################################
    """
    Distance = numpy.array(Chem.GetDistanceMatrix(mol), "d")
    Adjacent = numpy.array(Chem.GetAdjacencyMatrix(mol), "d")
    VertexDegree = sum(Adjacent)

    return sum(scipy.dot((Distance + Adjacent), VertexDegree)) 

def valences_not_too_large(rdkit_mol):
    valence_dict = {5: 3, 6: 4, 7: 3, 8: 2, 9: 1, 14: 4, 15: 5, 16: 6, 17: 1, 34: 2, 35: 1, 53: 1}
    atomicNumList = [a.GetAtomicNum() for a in rdkit_mol.GetAtoms()]
    valences = [valence_dict[atomic_num] for atomic_num in atomicNumList]
    BO = Chem.GetAdjacencyMatrix(rdkit_mol, useBO=True)
    number_of_bonds_list = BO.sum(axis=1)
    for valence, number_of_bonds in zip(valences, number_of_bonds_list):
        if number_of_bonds > valence:
            return False

    return True


# code modified from https://github.com/haroldsultan/MCTS/blob/master/mcts.py 

def test_smiles_from_adjacent_matrix(smiles):

    charged_fragments = True
    quick = True

    # Cut apart the smiles
    mol = get_mol(smiles)
    atoms = get_atoms(mol)
    charge = Chem.GetFormalCharge(mol)
    adjacent_matrix = Chem.GetAdjacencyMatrix(mol)

    #
    mol = Chem.RemoveHs(mol)
    canonical_smiles = Chem.MolToSmiles(mol)

    # Define new molecule template from atoms
    new_mol = x2m.get_proto_mol(atoms)

    # reconstruct the molecule from adjacent matrix, atoms and total charge
    new_mol = x2m.AC2mol(new_mol, adjacent_matrix, atoms, charge, charged_fragments, quick)
    new_mol = Chem.RemoveHs(new_mol)
    new_mol_smiles = Chem.MolToSmiles(new_mol)

    assert new_mol_smiles == canonical_smiles

    return 

def CalculateBalaban(mol):

    """
    #################################################################
    Calculation of Balaban index in a molecule

    ---->J

    Usage:

        result=CalculateBalaban(mol)

        Input: mol is a molecule object

        Output: result is a numeric value
    #################################################################
    """
    adjMat = Chem.GetAdjacencyMatrix(mol)
    Distance = Chem.GetDistanceMatrix(mol)
    Nbond = mol.GetNumBonds()
    Natom = mol.GetNumAtoms()
    S = numpy.sum(Distance, axis=1)
    mu = Nbond - Natom + 1
    sumk = 0.0
    for i in range(len(Distance)):
        si = S[i]
        for j in range(i, len(Distance)):
            if adjMat[i, j] == 1:
                sumk += 1.0 / numpy.sqrt(si * S[j])
    if mu + 1 != 0:
        J = float(Nbond) / float(mu + 1) * sumk
    else:
        J = 0
    return J 

def create_datasets(task, dataset, radius, device):

    dir_dataset = '../dataset/' + task + '/' + dataset + '/'

    """Initialize x_dict, in which each key is a symbol type
    (e.g., atom and chemical bond) and each value is its index.
    """
    atom_dict = defaultdict(lambda: len(atom_dict))
    bond_dict = defaultdict(lambda: len(bond_dict))
    fingerprint_dict = defaultdict(lambda: len(fingerprint_dict))
    edge_dict = defaultdict(lambda: len(edge_dict))

    def create_dataset(filename):

        print(filename)

        """Load a dataset."""
        with open(dir_dataset + filename, 'r') as f:
            smiles_property = f.readline().strip().split()
            data_original = f.read().strip().split('\n')

        """Exclude the data contains '.' in its smiles."""
        data_original = [data for data in data_original
                         if '.' not in data.split()[0]]

        dataset = []

        for data in data_original:

            smiles, property = data.strip().split()

            """Create each data with the above defined functions."""
            mol = Chem.AddHs(Chem.MolFromSmiles(smiles))
            atoms = create_atoms(mol, atom_dict)
            molecular_size = len(atoms)
            i_jbond_dict = create_ijbonddict(mol, bond_dict)
            fingerprints = extract_fingerprints(radius, atoms, i_jbond_dict,
                                                fingerprint_dict, edge_dict)
            adjacency = Chem.GetAdjacencyMatrix(mol)

            """Transform the above each data of numpy
            to pytorch tensor on a device (i.e., CPU or GPU).
            """
            fingerprints = torch.LongTensor(fingerprints).to(device)
            adjacency = torch.FloatTensor(adjacency).to(device)
            if task == 'classification':
                property = torch.LongTensor([int(property)]).to(device)
            if task == 'regression':
                property = torch.FloatTensor([[float(property)]]).to(device)

            dataset.append((fingerprints, adjacency, molecular_size, property))

        return dataset

    dataset_train = create_dataset('data_train.txt')
    dataset_train, dataset_dev = split_dataset(dataset_train, 0.9)
    dataset_test = create_dataset('data_test.txt')

    N_fingerprints = len(fingerprint_dict)

    return dataset_train, dataset_dev, dataset_test, N_fingerprints 

def _GetBurdenMatrix(mol, propertylabel="m"):
    """
    #################################################################
    *Internal used only**

    Calculate Burden matrix and their eigenvalues.
    #################################################################
    """
    mol = Chem.AddHs(mol)
    Natom = mol.GetNumAtoms()

    AdMatrix = Chem.GetAdjacencyMatrix(mol)
    bondindex = numpy.argwhere(AdMatrix)
    AdMatrix1 = numpy.array(AdMatrix, dtype=numpy.float32)

    # The diagonal elements of B, Bii, are either given by
    # the carbon normalized atomic mass,
    # van der Waals volume, Sanderson electronegativity,
    # and polarizability of atom i.

    for i in range(Natom):
        atom = mol.GetAtomWithIdx(i)
        temp = GetRelativeAtomicProperty(
            element=atom.GetSymbol(), propertyname=propertylabel
        )
        AdMatrix1[i, i] = round(temp, 3)

    # The element of B connecting atoms i and j, Bij,
    # is equal to the square root of the bond
    # order between atoms i and j.

    for i in bondindex:
        bond = mol.GetBondBetweenAtoms(int(i[0]), int(i[1]))
        if bond.GetBondType().name == "SINGLE":
            AdMatrix1[i[0], i[1]] = round(numpy.sqrt(1), 3)
        if bond.GetBondType().name == "DOUBLE":
            AdMatrix1[i[0], i[1]] = round(numpy.sqrt(2), 3)
        if bond.GetBondType().name == "TRIPLE":
            AdMatrix1[i[0], i[1]] = round(numpy.sqrt(3), 3)
        if bond.GetBondType().name == "AROMATIC":
            AdMatrix1[i[0], i[1]] = round(numpy.sqrt(1.5), 3)

    ##All other elements of B (corresponding non bonded
    # atom pairs) are set to 0.001
    bondnonindex = numpy.argwhere(AdMatrix == 0)

    for i in bondnonindex:
        if i[0] != i[1]:

            AdMatrix1[i[0], i[1]] = 0.001

    return numpy.real(numpy.linalg.eigvals(AdMatrix1)) 

def read_graph(source_path,MAX_size):
  Vertex = []
  Adj = [] # Normalized adjacency matrix
  mycount=1
  PAD=0
  mydict={}
  max_size=0
  with tf.gfile.GFile(source_path, mode="r") as source_file:
      source = source_file.readline().strip()
      counter = 0
      while source:
        mol = Chem.MolFromSmiles(source)
        atom_list = []
        for a in mol.GetAtoms():
            m = a.GetSymbol()
            if m not in mydict:
              mydict[m]=mycount
              mycount = mycount +1
            
            atom_list.append(mydict[m])

        if len(atom_list) > max_size:
           max_size = len(atom_list)


        if len(atom_list) < MAX_size:
           pad = [PAD] * (MAX_size - len(atom_list))
           atom_list = atom_list+pad

        
        vertex = np.array(atom_list, np.int32)
        Vertex.append(vertex)

        adja_mat = Chem.GetAdjacencyMatrix(mol)
        adj_temp = []
        for adja in adja_mat:
            if len(adja) < MAX_size:
               pad = [PAD]*(MAX_size - len(adja))
               adja = np.array(list(adja)+pad,np.int32)
            adj_temp.append(adja)
      
        cur_len = len(adj_temp)
        for i in range(MAX_size - cur_len):
            adja =np.array( [PAD]*MAX_size,np.int32)
            adj_temp.append(adja)

        adj_temp = adj_temp + np.eye(MAX_size) # A_hat = A + I
        Adj.append(adj_temp) 
        source = source_file.readline().strip()
  return Vertex,Adj,max_size



################ Reading initial states and weigths 

def read_graph(source_path,MAX_size):
  Vertex = []
  Adj = [] # Normalized adjacency matrix
  mycount=1
  PAD=0
  mydict={}
  max_size=0
  with tf.gfile.GFile(source_path, mode="r") as source_file:
      source = source_file.readline().strip()
      counter = 0
      while source:
        mol = Chem.MolFromSmiles(source)
        atom_list = []
        for a in mol.GetAtoms():
            m = a.GetSymbol()
            if m not in mydict:
              mydict[m]=mycount
              mycount = mycount +1
            
            atom_list.append(mydict[m])

        if len(atom_list) > max_size:
           max_size = len(atom_list)


        if len(atom_list) < MAX_size:
           pad = [PAD] * (MAX_size - len(atom_list))
           atom_list = atom_list+pad

        
        vertex = np.array(atom_list, np.int32)
        Vertex.append(vertex)

        adja_mat = Chem.GetAdjacencyMatrix(mol)
        adj_temp = []
        for adja in adja_mat:
            if len(adja) < MAX_size:
               pad = [PAD]*(MAX_size - len(adja))
               adja = np.array(list(adja)+pad,np.int32)
            adj_temp.append(adja)
      
        cur_len = len(adj_temp)
        for i in range(MAX_size - cur_len):
            adja =np.array( [PAD]*MAX_size,np.int32)
            adj_temp.append(adja)

        adj_temp = adj_temp + np.eye(MAX_size) # A_hat = A + I
        Adj.append(adj_temp) 
        source = source_file.readline().strip()
  return Vertex,Adj,max_size



################ Reading initial states and weigths 

def create_adjacency(mol):
    adjacency = Chem.GetAdjacencyMatrix(mol)
    return np.array(adjacency)