pKa/pdb2pka_interface.py
def pdb2pka_desolv_backgr(self,residue):
        """Calculate the desolvation and background interaction energy for a residue"""
        protein,routines,forcefield,apbs_setup,lig_titgrps = pdb2pka.pre_init(pdbfilename=self.pdbfile,
                                                                              ff='parse',
                                                                              ligand=None,
                                                                              verbose=1)
        mypkaRoutines = pdb2pka.pKaRoutines(protein, routines, forcefield,apbs_setup)
        #
        # Find our group
        #
        sp=residue.split(':')
        chainid=sp[0]
        #if chainid!='':
        #    raise Exception,'pKD cannot handle PDB files with ChainIDs!'
        resnum=int(sp[1])
        target=sp[2]
        mypkaRoutines.findTitratableGroups()
        this_pKa=None
        for pKa in mypkaRoutines.pKas:
            print pKa.residue
            print pKa.uniqueid
            print pKa.residue.resSeq,resnum, pKa.residue.chainID,'CID',chainid,pKa.residue.name,target,pKa.pKaGroup.name,target
            #if pKa.residue.resSeq==resnum and pKa.residue.chainID==chainid and pKa.residue.name==target and pKa.pKaGroup.name==target:
            if pKa.residue.resSeq==resnum and pKa.residue.chainID==chainid and pKa.residue.name==target and pKa.pKaGroup.name==target:
                this_pKa=pKa
                break
        if not this_pKa:
            raise Exception,'Could not find inserted titratable group'
        mypkaRoutines.calculateBackground(onlypKa=pKa)
        mypkaRoutines.calculateDesolvation(onlypKa=pKa)
        #
        # Get the intrinsic pKa
        #
        pKaGroup=pKa.pKaGroup
        Gtype=pKa.pKaGroup.type
        #
        # We measure intrinsic pKa values against a single reference state
        #
        desolv=[]
        backgr=[]
        import math
        ln10=math.log(10)
        for titration in pKaGroup.DefTitrations:
            #
            # Find an uncharged reference state
            #
            ref_state=mypkaRoutines.neutral_ref_state[pKa][titration]
            all_states=titration.allstates
            all_states.sort()
            for state in all_states:
                if mypkaRoutines.is_charged(pKa,titration,state)==1:
                    dpKa_desolv=(pKa.desolvation[state]-pKa.desolvation[ref_state])/ln10
                    dpKa_backgr=(pKa.background[state]-pKa.background[ref_state])/ln10
                    #
                    # Make acid and base modifications
                    #
                    if Gtype=='base':
                        dpKa_desolv=-dpKa_desolv
                        dpKa_backgr=-dpKa_backgr
                    #
                    # Now calculate intrinsic pKa
                    #
                    backgr.append(dpKa_backgr)
                    desolv.append(dpKa_desolv)
                    intpKa=titration.modelpKa+dpKa_desolv+dpKa_backgr
                    #print 'Energy difference for %s -> %s [reference state] is %5.2f pKa units' %(state,ref_state,intpKa)
                    pKa.intrinsic_pKa[state]=intpKa
                else:
                    #
                    # Neutral states - why do we treat them differently?
                    #
                    dpKa_desolv=(pKa.desolvation[state]-pKa.desolvation[ref_state])/ln10
                    dpKa_backgr=(pKa.background[state]-pKa.background[ref_state])/ln10
                    #
                    # Make acid and base modifications
                    #
                    if Gtype=='base':
                        dpKa_desolv=-dpKa_desolv
                        dpKa_backgr=-dpKa_backgr
                    backgr.append(dpKa_backgr)
                    desolv.append(dpKa_desolv)
                    dpKa=dpKa_desolv+dpKa_backgr
                    #print 'Energy difference for %s -> %s [reference state] is %5.2f kT' %(state,ref_state,dpKa)
                    pKa.intrinsic_pKa[state]=dpKa
        #print '-----------------'
        #
        # One value is zero, so just get the avg of the rest
        #
        des_sum=0.0
        for val in desolv:
            des_sum=des_sum+val
        des_sum=des_sum/float(len(desolv)-1)
        #
        bac_sum=0.0
        for val in backgr:
            bac_sum=bac_sum+val
        bac_sum=bac_sum/float(len(backgr)-1)
        return  {residue:des_sum},{residue:bac_sum}
Calculate the desolvation and background interaction energy for a residue
Calculate the desolvation and background interaction energy for a residue
Calculate the desolvation and background interaction energy for a residue
def pdb2pka_desolv_backgr(self,residue):
        protein,routines,forcefield,apbs_setup,lig_titgrps = pdb2pka.pre_init(pdbfilename=self.pdbfile,
                                                                              ff='parse',
                                                                              ligand=None,
                                                                              verbose=1)
        mypkaRoutines = pdb2pka.pKaRoutines(protein, routines, forcefield,apbs_setup)
         Find our group
        sp=residue.split(':')
        chainid=sp[0]
        if chainid!='':
            raise Exception,'pKD cannot handle PDB files with ChainIDs!'
        resnum=int(sp[1])
        target=sp[2]
        mypkaRoutines.findTitratableGroups()
        this_pKa=None
        for pKa in mypkaRoutines.pKas:
            print pKa.residue
            print pKa.uniqueid
            print pKa.residue.resSeq,resnum, pKa.residue.chainID,'CID',chainid,pKa.residue.name,target,pKa.pKaGroup.name,target
            if pKa.residue.resSeq==resnum and pKa.residue.chainID==chainid and pKa.residue.name==target and pKa.pKaGroup.name==target:
            if pKa.residue.resSeq==resnum and pKa.residue.chainID==chainid and pKa.residue.name==target and pKa.pKaGroup.name==target:
                this_pKa=pKa
                break
        if not this_pKa:
            raise Exception,'Could not find inserted titratable group'
        mypkaRoutines.calculateBackground(onlypKa=pKa)
        mypkaRoutines.calculateDesolvation(onlypKa=pKa)
         Get the intrinsic pKa
        pKaGroup=pKa.pKaGroup
        Gtype=pKa.pKaGroup.type
         We measure intrinsic pKa values against a single reference state
        desolv=[]
        backgr=[]
        import math
        ln10=math.log(10)
        for titration in pKaGroup.DefTitrations:
             Find an uncharged reference state
            ref_state=mypkaRoutines.neutral_ref_state[pKa][titration]
            all_states=titration.allstates
            all_states.sort()
            for state in all_states:
                if mypkaRoutines.is_charged(pKa,titration,state)==1:
                    dpKa_desolv=(pKa.desolvation[state]-pKa.desolvation[ref_state])/ln10
                    dpKa_backgr=(pKa.background[state]-pKa.background[ref_state])/ln10
                     Make acid and base modifications
                    if Gtype=='base':
                        dpKa_desolv=-dpKa_desolv
                        dpKa_backgr=-dpKa_backgr
                     Now calculate intrinsic pKa
                    backgr.append(dpKa_backgr)
                    desolv.append(dpKa_desolv)
                    intpKa=titration.modelpKa+dpKa_desolv+dpKa_backgr
                    print 'Energy difference for %s -> %s [reference state] is %5.2f pKa units' %(state,ref_state,intpKa)
                    pKa.intrinsic_pKa[state]=intpKa
                else:
                     Neutral states - why do we treat them differently?
                    dpKa_desolv=(pKa.desolvation[state]-pKa.desolvation[ref_state])/ln10
                    dpKa_backgr=(pKa.background[state]-pKa.background[ref_state])/ln10
                     Make acid and base modifications
                    if Gtype=='base':
                        dpKa_desolv=-dpKa_desolv
                        dpKa_backgr=-dpKa_backgr
                    backgr.append(dpKa_backgr)
                    desolv.append(dpKa_desolv)
                    dpKa=dpKa_desolv+dpKa_backgr
                    print 'Energy difference for %s -> %s [reference state] is %5.2f kT' %(state,ref_state,dpKa)
                    pKa.intrinsic_pKa[state]=dpKa
        print '-----------------'
         One value is zero, so just get the avg of the rest
        des_sum=0.0
        for val in desolv:
            des_sum=des_sum+val
        des_sum=des_sum/float(len(desolv)-1)
        bac_sum=0.0
        for val in backgr:
            bac_sum=bac_sum+val
        bac_sum=bac_sum/float(len(backgr)-1)
        return  {residue:des_sum},{residue:bac_sum}
def pdb2pka_desolv_backgr ( self , residue ) : protein , routines , forcefield , apbs_setup , lig_titgrps = pdb2pka . pre_init ( pdbfilename = self . pdbfile , ff = 'parse' , ligand = None , verbose = 1 ) mypkaRoutines = pdb2pka . pKaRoutines ( protein , routines , forcefield , apbs_setup ) sp = residue . split ( ':' ) chainid = sp [ 0 ] resnum = int ( sp [ 1 ] ) target = sp [ 2 ] mypkaRoutines . findTitratableGroups ( ) this_pKa = None for pKa in mypkaRoutines . pKas : print pKa . residue print pKa . uniqueid print pKa . residue . resSeq , resnum , pKa . residue . chainID , 'CID' , chainid , pKa . residue . name , target , pKa . pKaGroup . name , target if pKa . residue . resSeq == resnum and pKa . residue . chainID == chainid and pKa . residue . name == target and pKa . pKaGroup . name == target : this_pKa = pKa break if not this_pKa : raise Exception , 'Could not find inserted titratable group' mypkaRoutines . calculateBackground ( onlypKa = pKa ) mypkaRoutines . calculateDesolvation ( onlypKa = pKa ) pKaGroup = pKa . pKaGroup Gtype = pKa . pKaGroup . type desolv = [ ] backgr = [ ] import math ln10 = math . log ( 10 ) for titration in pKaGroup . DefTitrations : ref_state = mypkaRoutines . neutral_ref_state [ pKa ] [ titration ] all_states = titration . allstates all_states . sort ( ) for state in all_states : if mypkaRoutines . is_charged ( pKa , titration , state ) == 1 : dpKa_desolv = ( pKa . desolvation [ state ] - pKa . desolvation [ ref_state ] ) / ln10 dpKa_backgr = ( pKa . background [ state ] - pKa . background [ ref_state ] ) / ln10 if Gtype == 'base' : dpKa_desolv = - dpKa_desolv dpKa_backgr = - dpKa_backgr backgr . append ( dpKa_backgr ) desolv . append ( dpKa_desolv ) intpKa = titration . modelpKa + dpKa_desolv + dpKa_backgr pKa . intrinsic_pKa [ state ] = intpKa else : dpKa_desolv = ( pKa . desolvation [ state ] - pKa . desolvation [ ref_state ] ) / ln10 dpKa_backgr = ( pKa . background [ state ] - pKa . background [ ref_state ] ) / ln10 if Gtype == 'base' : dpKa_desolv = - dpKa_desolv dpKa_backgr = - dpKa_backgr backgr . append ( dpKa_backgr ) desolv . append ( dpKa_desolv ) dpKa = dpKa_desolv + dpKa_backgr pKa . intrinsic_pKa [ state ] = dpKa des_sum = 0.0 for val in desolv : des_sum = des_sum + val des_sum = des_sum / float ( len ( desolv ) - 1 ) bac_sum = 0.0 for val in backgr : bac_sum = bac_sum + val bac_sum = bac_sum / float ( len ( backgr ) - 1 ) return { residue : des_sum } , { residue : bac_sum }
Calculate the desolvation and background interaction energy for a residue
Calculate the desolvation and background interaction energy for a residue
"""Calculate the desolvation and background interaction energy for a residue""" # # Find our group # #if chainid!='': #    raise Exception,'pKD cannot handle PDB files with ChainIDs!' #if pKa.residue.resSeq==resnum and pKa.residue.chainID==chainid and pKa.residue.name==target and pKa.pKaGroup.name==target: # # Get the intrinsic pKa # # # We measure intrinsic pKa values against a single reference state # # # Find an uncharged reference state # # # Make acid and base modifications # # # Now calculate intrinsic pKa # #print 'Energy difference for %s -> %s [reference state] is %5.2f pKa units' %(state,ref_state,intpKa) # # Neutral states - why do we treat them differently? # # # Make acid and base modifications # #print 'Energy difference for %s -> %s [reference state] is %5.2f kT' %(state,ref_state,dpKa) #print '-----------------' # # One value is zero, so just get the avg of the rest # #
self
residue