monomer_classes_Cl.py

# Classes for MPD and TMC monomers, can be used for fragments within polymer chain

import MDAnalysis as mda
import numpy as np
    
class MPD():
    def __init__(self, N):
        '''
        N should be one of the 2 nitrogens on the MPD monomer
        
        Positions on MPD ring
               _1_
              4    2
              |    |
         N2 - 5_  _3 - N1
                6

        '''

        # Save input nitrogen and molecule information
        if N.type == '7': 
            self.N1 = N # N1 is always crosslinked
            self.atoms = [self.N1]
        else:
            self.N2 = N
            self.atoms = [self.N2]
        self.mol = N.residue

        # Define the aromatic carbons by diagram above
        self.ca3 = [atom for atom in N.bonded_atoms if atom.type == '2'][0]
        self.atoms.append(self.ca3)
        
        ca26 = [atom for atom in self.ca3.bonded_atoms if atom.type == '2']
        for ca in ca26:
            ca5 = [atom for atom in ca.bonded_atoms if atom != self.ca3 and ('7' in atom.bonded_atoms.types or '8' in atom.bonded_atoms.types)]
            if len(ca5) > 0:
                self.ca5 = ca5[0]
                self.atoms.append(self.ca5)
                self.ca6 = ca
                self.atoms.append(self.ca6)
            else:
                self.ca1 = [atom for atom in ca.bonded_atoms if atom.type == '2' and atom != self.ca3][0]
                self.atoms.append(self.ca1)
                self.ca2 = ca
                self.atoms.append(self.ca2)

        self.ca4 = [atom for atom in self.ca5.bonded_atoms if atom.type == '2' and atom != self.ca6][0]
        self.atoms.append(self.ca4)

        # Save other nitrogen
        if N.type == '7':
            self.N2 = [atom for atom in self.ca5.bonded_atoms if atom.type == '7' or atom.type == '8'][0]
            self.atoms.append(self.N2)
        else:
            self.N1 = [atom for atom in self.ca5.bonded_atoms if atom.type == '7' or atom.type == '8'][0]
            self.atoms.append(self.N1)

        # Save aromatic hydrogens
        self.ha1 = [atom for atom in self.ca1.bonded_atoms if atom.type == '3'][0]
        self.atoms.append(self.ha1)
        self.ha2 = [atom for atom in self.ca2.bonded_atoms if atom.type == '3'][0]
        self.atoms.append(self.ha2)
        self.ha4 = [atom for atom in self.ca4.bonded_atoms if atom.type == '3'][0]
        self.atoms.append(self.ha4)
        self.ha6 = [atom for atom in self.ca6.bonded_atoms if atom.type == '3'][0]
        self.atoms.append(self.ha6)
    
        # Save hydrogens on nitrogens
        hn1 = [atom for atom in self.N1.bonded_atoms if atom.type == '4']
        self.hn1a = hn1[0]
        self.atoms.append(self.hn1a)
        if len(hn1) == 2:
            self.hn1b = hn1[1]
            self.atoms.append(self.hn1b)

        hn2 = [atom for atom in self.N2.bonded_atoms if atom.type == '4']
        self.hn2a = hn2[0]
        self.atoms.append(self.hn2a)
        if len(hn2) == 2:
            self.hn2b = hn2[1]
            self.atoms.append(self.hn2b)

        # Classify fragment as monomer, linear, or terminated
        n_hn = len(hn1) + len(hn2)
        if n_hn == 4:
            self.frag_type = 'mono'
            self.name = 'MPD'
        elif n_hn == 3:
            self.frag_type = 'term'
            self.name = 'MPD_T'
        elif n_hn == 2:
            self.frag_type = 'lin'
            self.name = 'MPD_L'
        else:
            raise AttributeError(f'cannot have number of hn type = {n_hn}')
        
        n_cl = len([atom for atom in self.mol.atoms if atom.type == '6'])
        if n_cl > 0:
            self.name += '_Cl'

    
    def assign_charges(self, charges):
        '''Assign charges to the fragment from the charges dictionary'''

        for a_type in charges.keys():
            charge = round(charges[a_type]['charge'], 6)
            # charge = charges[a_type]['charge']
            atom = getattr(self, a_type)
            atom.charge = charge        



class TMC():
    def __init__(self, C):
        '''
        C should be one of the 3 carbonyl carbons on TMC
        
        Positions on TMC ring
            C1
            |
           _1_
          4    2
          |    |
        C3-5_  _3-C2
            6
                
        '''

        # Save input carbon and molecule information
        self.C1 = C
        self.mol = C.residue
        self.atoms = [self.C1]

        # Define the aromatic carbons by diagram above
        self.ca1 = [atom for atom in C.bonded_atoms if atom.type == '2'][0]
        self.atoms.append(self.ca1)
        self.ca2 = [atom for atom in self.ca1.bonded_atoms if atom.type == '2'][0]
        self.atoms.append(self.ca2)
        self.ca3 = [atom for atom in self.ca2.bonded_atoms if atom.type == '2' and atom != self.ca1][0]
        self.atoms.append(self.ca3)
        self.ca6 = [atom for atom in self.ca3.bonded_atoms if atom.type == '2' and atom != self.ca2][0]
        self.atoms.append(self.ca6)
        self.ca5 = [atom for atom in self.ca6.bonded_atoms if atom.type == '2' and atom != self.ca3][0]
        self.atoms.append(self.ca5)
        self.ca4 = [atom for atom in self.ca5.bonded_atoms if atom.type == '2' and atom != self.ca6][0]
        self.atoms.append(self.ca4)

        # Save other carbonyl carbons
        self.C2 = [atom for atom in self.ca3.bonded_atoms if atom.type == '1'][0]
        self.atoms.append(self.C2)
        self.C3 = [atom for atom in self.ca5.bonded_atoms if atom.type == '1'][0]
        self.atoms.append(self.C3)

        # Save aromatic hydrogens
        self.ha2 = [atom for atom in self.ca2.bonded_atoms if atom.type == '3'][0]
        self.atoms.append(self.ha2)
        self.ha4 = [atom for atom in self.ca4.bonded_atoms if atom.type == '3'][0]
        self.atoms.append(self.ha4)
        self.ha6 = [atom for atom in self.ca6.bonded_atoms if atom.type == '3'][0]
        self.atoms.append(self.ha6)

        # Save carbonyl oxygens
        self.O1 = [atom for atom in self.C1.bonded_atoms if atom.type == '9'][0]
        self.atoms.append(self.O1)
        self.O2 = [atom for atom in self.C2.bonded_atoms if atom.type == '9'][0]
        self.atoms.append(self.O2)
        self.O3 = [atom for atom in self.C3.bonded_atoms if atom.type == '9'][0]
        self.atoms.append(self.O3)

        # Classify fragment as monomer, terminated, linear, or crosslinked and save hydroxide groups
        n_oh = 0
        n_cl = 0
        c_types = []
        for i,c in enumerate([self.C1, self.C2, self.C3]):
            oh = [atom for atom in c.bonded_atoms if atom.type == '10']
            cl = [atom for atom in c.bonded_atoms if atom.type == '6']
            if len(oh) > 0:
                c_types.append('term')
                n_oh += 1
                setattr(self, f'oh{i+1}', oh[0])
                self.atoms.append(getattr(self, f'oh{i+1}'))
                ho = [atom for atom in oh[0].bonded_atoms if atom.type == '5']
                setattr(self, f'ho{i+1}', ho[0])
                self.atoms.append(getattr(self, f'ho{i+1}'))
            elif len(cl) > 0:
                c_types.append('unreact')
                n_cl += 1
                setattr(self, f'Cl{i+1}', cl[0])
                self.atoms.append(getattr(self, f'Cl{i+1}'))
            else:
                c_types.append('xlink')
    
        if n_oh == 0 and n_cl == 0:
            self.frag_type = 'xlink'
            self.name = 'TMC_C'
        elif n_oh == 0 and n_cl == 1:
            self.frag_type = 'lin'
            if c_types[0] == 'unreact':
                self.name = 'TMC_L_Cl_1'
            elif c_types[1] == 'unreact':
                self.name = 'TMC_L_Cl_2'
            elif c_types[2] == 'unreact':
                self.name = 'TMC_L_Cl_3'
        elif n_oh == 0 and n_cl == 2:
            self.frag_type = 'term'
            if c_types[0] == 'xlink':
                self.name = 'TMC_T_Cl_1'
            elif c_types[1] == 'xlink':
                self.name = 'TMC_T_Cl_2'
            elif c_types[2] == 'xlink':
                self.name = 'TMC_T_Cl_3'
        elif n_oh == 0 and n_cl == 3:
            self.frag_type = 'mono'
            self.name = 'TMC_Cl'
        elif n_oh == 1:
            self.frag_type = 'lin'
            if c_types[0] == 'term':
                self.name = 'TMC_L_1'
            elif c_types[1] == 'term':
                self.name = 'TMC_L_2'
            elif c_types[2] == 'term':
                self.name = 'TMC_L_3'
        elif n_oh == 2:
            self.frag_type = 'term'
            if c_types[0]== 'xlink':
                self.name = 'TMC_T_1'
            elif c_types[1] == 'xlink':
                self.name = 'TMC_T_2'
            elif c_types[2] == 'xlink':
                self.name = 'TMC_T_3'
        elif n_oh == 3:
            self.frag_type = 'mono'
            self.name = 'TMC'
        else:
            raise AttributeError(f'cannot have number of oh type = {n_oh}')


    def assign_charges(self, charges):
        '''Assign charges to the fragment from the charges dictionary'''

        for a_type in charges.keys():
            charge = round(charges[a_type]['charge'], 6)
            # charge = charges[a_type]['charge']
            atom = getattr(self, a_type)
            atom.charge = charge    
            


def check_N(N):

    # input the atom N to check the other N
    #   if the other N is type NH, then MPD-T
    #   if the other N is type LN, then MPD-L

    # Positions on MPD ring
    #       _1_
    #     4    2
    #     |    |
    # N - 5_  _3 - N
    #       6

    if N.type == '8':
        n_NH = 1
    elif N.type == '7':
        n_NH = 0
    else:
        raise ValueError(f'{N} is not a N')
    
    CA3 = [atom for atom in N.bonded_atoms if atom.type == '2'][0]
    CA26 = [atom for atom in CA3.bonded_atoms if atom.type == '2']
    for CA in CA26:
        n_NH += len([atom for atom in CA.bonded_atoms if atom != CA3 and atom.type == '2' and '8' in atom.bonded_atoms.types])

    return n_NH


def check_C(C):

    # input the atom C to check the other Cs
    #   if n_CT = 0, then TMC-C
    #   if n_CT = 1, then TMC-L
    #   if n_CT = 2, then TMC-T

    # Positions on TMC ring
    #        C(14)
    #        |
    #       _1_
    #     4    2
    #     |    |
    # C - 5_  _3 - C
    #       6

    n_CT = len([atom for atom in C.bonded_atoms if atom.type == '10'])
    CA1 = [atom for atom in C.bonded_atoms if atom.type == '2'][0]
    CA2 = [atom for atom in CA1.bonded_atoms if atom.type == '2'][0]
    CA3 = [atom for atom in CA2.bonded_atoms if atom.type == '2' and atom != CA1][0]
    n_CT += len([atom for atom in CA3.bonded_atoms if atom.type == '1' and '10' in atom.bonded_atoms.types])
    CA6 = [atom for atom in CA3.bonded_atoms if atom.type == '2' and atom != CA2][0]
    CA5 = [atom for atom in CA6.bonded_atoms if atom.type == '2' and atom != CA3][0]
    n_CT += len([atom for atom in CA5.bonded_atoms if atom.type == '1' and '10' in atom.bonded_atoms.types])

    return n_CT

        
if __name__ == '__main__':

    u = mda.Universe('testing/2MPD-2TMC.mol2')

    atoms = u.atoms

    # for mol2, convert to numeric atom types
    atom_mapping = {
        'C.2' : '1',
        'C.ar' : '2',
        'H' : {
            'C.ar' : '3',
            'N.am' : '4',
            'N.pl3' : '4',
            'O.3' : '5'
        },
        'Cl' : '6',
        'N.am' : '7',
        'N.pl3' : '8',
        'O.2' : '9',
        'O.3' : '10'   
    }

    num_types = []
    for atom in atoms:
        if atom.type != 'H':
            num_types.append(atom_mapping[atom.type])
        else:
            bonded_type = atom.bonded_atoms[0].type
            num_types.append(atom_mapping[atom.type][bonded_type])

    atoms.types = num_types
    
    # Create MPD and TMC fragments
    mpds = []
    Ns = []
    for N in u.select_atoms('type 7'):
        mpd = MPD(N)
        if mpd.N2 not in Ns:
            mpds.append(mpd)
        
        Ns.append(mpd.N1)
        Ns.append(mpd.N2)

    tmcs = []
    Cs = []
    for C in u.select_atoms('type 1'):
        tmc = TMC(C)
        if tmc.C2 not in Cs and tmc.C3 not in Cs:
            tmcs.append(tmc)

        Cs.append(tmc.C1)
        Cs.append(tmc.C2)
        Cs.append(tmc.C3)


    # Assign charges to MPD and TMC fragments
    import yaml
    
    with open('charges.yaml', 'r') as file: # read in charges from yaml file
        charges_dict = yaml.safe_load(file)

    total_charge = 0
    for mono in mpds + tmcs:
        print(mono.name)
        mono.assign_charges(charges_dict[mono.name])
        frag_charge = np.array([a.charge for a in mono.atoms]).sum()
        total_charge += frag_charge


    print('New total charge in system: {:.4f}'.format(total_charge))