Changeset 3099


Ignore:
Timestamp:
Sep 27, 2017 10:56:21 AM (4 years ago)
Author:
vondreele
Message:

fixes to protein validation code

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r3097 r3099  
    26202620    for atom in Atoms:
    26212621        if atom[1] in resNames:
     2622            cartAtoms.append(atom[:cx+3])
    26222623            if atom[4].strip() in ['S','Se']:
    26232624                if not old:
    26242625                    continue        #S,Se skipped for erratv2?
    2625                 atom[3] = 'Os'
    2626                 atom[4] = 'O'
    2627             cartAtoms.append(atom[:cx+3])
     2626                cartAtoms[-1][3] = 'Os'
     2627                cartAtoms[-1][4] = 'O'
    26282628            cartAtoms[-1][cx:cx+3] = np.inner(Amat,cartAtoms[-1][cx:cx+3])
    26292629    XYZ = np.array([atom[cx:cx+3] for atom in cartAtoms])
     
    26362636        Boxes[box[0],box[1],box[2],0] += 1
    26372637        Boxes[box[0],box[1],box[2],Boxes[box[0],box[1],box[2],0]] = ib
    2638     #Box content checks with errat.f ibox1 array
     2638    #Box content checks with errat.f $ erratv2.cpp ibox1 arrays
    26392639    indices = (-1,0,1)
    26402640    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
     
    26512651    intact = {'CC':0,'CN':0,'CO':0,'NN':0,'NO':0,'OO':0,'NC':0,'OC':0,'ON':0}
    26522652    for ia,atom in enumerate(cartAtoms):
     2653        jntact = {'CC':0,'CN':0,'CO':0,'NN':0,'NO':0,'OO':0,'NC':0,'OC':0,'ON':0}
    26532654        if atom[2] not in chains:   #get chain id & save residue sequence from last chain
    26542655            chains.append(atom[2])
     
    26872688                tgts = [tgt for tgt in tgts if not (cartAtoms[tgt][3].strip() == 'N' and int(cartAtoms[tgt][0]) in [ires-1,ires+1])]
    26882689        else:
     2690            tgts = [tgt for tgt in tgts if not int(cartAtoms[tgt][0]) in [ires+1,ires+2,ires+3,ires+4,ires+5,ires+6,ires+7,ires+8]]
    26892691            if atom[3].strip() == 'C':
    2690                 tgts = [tgt for tgt in tgts if not (cartAtoms[tgt][3].strip() == 'N' and int(cartAtoms[tgt][0])  in [ires-1,ires+1])]
     2692                tgts = [tgt for tgt in tgts if not (cartAtoms[tgt][3].strip() == 'N' and int(cartAtoms[tgt][0]) == ires+1)]
    26912693            elif atom[3].strip() == 'N':
    2692                 tgts = [tgt for tgt in tgts if not (cartAtoms[tgt][3].strip() == 'C' and int(cartAtoms[tgt][0])  in [ires-1,ires+1])]
     2694                tgts = [tgt for tgt in tgts if not (cartAtoms[tgt][3].strip() == 'C' and int(cartAtoms[tgt][0]) == ires-1)]
    26932695        for tgt in tgts:
    26942696            dsqt = np.sqrt(np.sum((XYZ[ia]-XYZ[tgt])**2))
     
    26972699                mult = 2.*(3.75-dsqt)
    26982700            intype = atom[4].strip()+cartAtoms[tgt][4].strip()
    2699             intact[intype] += mult
     2701            if 'S' not in intype:
     2702                intact[intype] += mult
     2703                jntact[intype] += mult
     2704#        print ia,atom[0]+atom[1]+atom[3],tgts,jntact['CC'],jntact['CN']+jntact['NC'],jntact['CO']+jntact['OC'],jntact['NN'],jntact['NO']+jntact['ON']
    27002705    resNames += resname
    27012706    resIntAct.append(sumintact(intact))
  • trunk/GSASIIphsGUI.py

    r3097 r3099  
    33013301        print 'Ref: Colovos, C. & Yeates, T.O. Protein Science 2, 1511-1519 (1991).'
    33023302        print 'Residue error scores >6 for 5% & >8 for 1% likelihood of being correct'
    3303         print 'Plot 2 is Protein validation based on erratv2.cpp; by D. Obukhov & T. Yeates'
     3303        print 'NB: this calc. matches errat.f result'
     3304        print 'Plot 2 is Protein validation based on erratv2.cpp; by D. Obukhov & T. Yeates (2002)'
    33043305        print 'Ref: Colovos, C. & Yates, T.O. Protein Science 2, 1511-1519 (1991).'
    33053306        print 'Residue error scores >11.5 for 5% & >17.2 for 1% likelihood of being correct'
    3306         print 'NB: this calc. still in error'
     3307        print 'NB: this calc. gives a close approximate to original erratv2 result'
    33073308        G2plt.PlotAAProb(G2frame,resNames,Probs1,Probs2,Title='Error score for %s'%(data['General']['Name']),
    33083309            thresh=[[8.0,6.0],[17.191,11.527]])
Note: See TracChangeset for help on using the changeset viewer.