Changeset 1103


Ignore:
Timestamp:
Oct 13, 2013 7:55:07 AM (8 years ago)
Author:
vondreele
Message:

fixes to protein restraints & rigid body GUI issues

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIrestrGUI.py

    r1099 r1103  
    592592                                        continue
    593593                        else:
    594                             for res,name,id in Chains[chain][residue]:
    595                                 try:
    596                                     ipos = Atms.index(name)
    597                                     ids[ipos] = id
    598                                 except ValueError:
    599                                     continue
    600                             for res,name,id in Chains[chain][residue+1]:
    601                                 try:
    602                                     ipos = pAtms.index(name)
    603                                     ids[ipos] = id
    604                                 except ValueError:
    605                                     continue
     594                            try:
     595                                for res,name,id in Chains[chain][residue]:
     596                                    try:
     597                                        ipos = Atms.index(name)
     598                                        ids[ipos] = id
     599                                    except ValueError:
     600                                        continue
     601                                for res,name,id in Chains[chain][residue+1]:
     602                                    try:
     603                                        ipos = pAtms.index(name)
     604                                        ids[ipos] = id
     605                                    except ValueError:
     606                                        continue
     607                            except KeyError:
     608                                continue
    606609                        if np.all(ids):
    607610                            torsion = [list(ids),['1','1','1','1'],Name,Esd]
     
    667670                    else:
    668671                        for residue in residues[1:-1]:
    669                             for res,name,id in Chains[chain][residue-1]:
    670                                 try:
    671                                     ipos = mAtms.index(name)
    672                                     ids[ipos] = id
    673                                 except ValueError:
    674                                     continue
    675                             for res,name,id in Chains[chain][residue+1]:
    676                                 try:
    677                                     ipos = pAtms.index(name)
    678                                     ids[ipos] = id
    679                                 except ValueError:
    680                                     continue
    681                             for res,name,id in Chains[chain][residue]:
    682                                 if Res == res:
     672                            try:
     673                                for res,name,id in Chains[chain][residue-1]:
    683674                                    try:
    684                                         ipos = Atms.index(name)
     675                                        ipos = mAtms.index(name)
    685676                                        ids[ipos] = id
    686677                                    except ValueError:
    687678                                        continue
    688                             if np.all(ids):
    689                                 rama = [list(ids),['1','1','1','1','1'],Name,Esd]
    690                                 if rama not in ramaRestData['Ramas']:
    691                                     ramaRestData['Ramas'].append(rama)
    692                                 ids = np.array([0,0,0,0,0])
     679                                for res,name,id in Chains[chain][residue+1]:
     680                                    try:
     681                                        ipos = pAtms.index(name)
     682                                        ids[ipos] = id
     683                                    except ValueError:
     684                                        continue
     685                                for res,name,id in Chains[chain][residue]:
     686                                    if Res == res:
     687                                        try:
     688                                            ipos = Atms.index(name)
     689                                            ids[ipos] = id
     690                                        except ValueError:
     691                                            continue
     692                                if np.all(ids):
     693                                    rama = [list(ids),['1','1','1','1','1'],Name,Esd]
     694                                    if rama not in ramaRestData['Ramas']:
     695                                        ramaRestData['Ramas'].append(rama)
     696                                    ids = np.array([0,0,0,0,0])
     697                            except KeyError:
     698                                continue
    693699            macStr = macro.readline()
    694700        macro.close()
     
    18131819            ramaRestData = restrData['Rama']
    18141820            UpdateRamaRestr(ramaRestData)
    1815             G2plt.PlotRama(G2frame,phaseName,rama,ramaName)
     1821            wx.CallAfter(G2plt.PlotRama,G2frame,phaseName,rama,ramaName)
    18161822        elif text == 'Chem. comp. restraints':
    18171823            G2gd.SetDataMenuBar(G2frame,G2frame.dataFrame.RestraintMenu)
  • trunk/GSASIIstrMath.py

    r1101 r1103  
    591591        refl[10] = atan2d(fbs[0],fas[0])
    592592   
     593def StructureFactor2(refList,G,hfx,pfx,SGData,calcControls,parmDict):
     594    ''' Compute structure factors for all h,k,l for phase
     595    puts the result, F^2, in each ref[8] in refList - want to do this for blocks of reflections
     596    & not one by one.
     597    input:
     598   
     599    :param list refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv]
     600    :param np.array G:      reciprocal metric tensor
     601    :param str pfx:    phase id string
     602    :param dict SGData: space group info. dictionary output from SpcGroup
     603    :param dict calcControls:
     604    :param dict ParmDict:
     605
     606    '''       
     607    twopi = 2.0*np.pi
     608    twopisq = 2.0*np.pi**2
     609    phfx = pfx.split(':')[0]+hfx
     610    ast = np.sqrt(np.diag(G))
     611    Mast = twopisq*np.multiply.outer(ast,ast)
     612    FFtables = calcControls['FFtables']
     613    BLtables = calcControls['BLtables']
     614    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
     615    FF = np.zeros(len(Tdata))
     616    if 'N' in calcControls[hfx+'histType']:
     617        FP,FPP = G2el.BlenRes(Tdata,BLtables,parmDict[hfx+'Lam'])
     618    else:
     619        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
     620        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
     621    Uij = np.array(G2lat.U6toUij(Uijdata))
     622    bij = Mast*Uij.T
     623    for refl in refList:
     624        fbs = np.array([0,0])
     625        H = refl[:3]
     626        SQ = 1./(2.*refl[4])**2
     627        SQfactor = 4.0*SQ*twopisq
     628        Bab = parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor)
     629        if not len(refl[-1]):                #no form factors
     630            if 'N' in calcControls[hfx+'histType']:
     631                refl[-1] = G2el.getBLvalues(BLtables)
     632            else:       #'X'
     633                refl[-1] = G2el.getFFvalues(FFtables,SQ)
     634        for i,El in enumerate(Tdata):
     635            FF[i] = refl[-1][El]           
     636        Uniq = refl[11]
     637        phi = refl[12]
     638        phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+phi[:,np.newaxis])
     639        sinp = np.sin(phase)
     640        cosp = np.cos(phase)
     641        occ = Mdata*Fdata/len(Uniq)
     642        biso = -SQfactor*Uisodata
     643        Tiso = np.where(biso<1.,np.exp(biso),1.0)
     644        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
     645        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
     646        Tcorr = Tiso*Tuij
     647        fa = np.array([(FF+FP-Bab)*occ*cosp*Tcorr,-FPP*occ*sinp*Tcorr])
     648        fas = np.sum(np.sum(fa,axis=1),axis=1)        #real
     649        if not SGData['SGInv']:
     650            fb = np.array([(FF+FP-Bab)*occ*sinp*Tcorr,FPP*occ*cosp*Tcorr])
     651            fbs = np.sum(np.sum(fb,axis=1),axis=1)
     652        fasq = fas**2
     653        fbsq = fbs**2        #imaginary
     654        refl[9] = np.sum(fasq)+np.sum(fbsq)
     655        refl[10] = atan2d(fbs[0],fas[0])
     656   
    593657def StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict):
    594658    'Needs a doc string'
     659    twopi = 2.0*np.pi
     660    twopisq = 2.0*np.pi**2
     661    phfx = pfx.split(':')[0]+hfx
     662    ast = np.sqrt(np.diag(G))
     663    Mast = twopisq*np.multiply.outer(ast,ast)
     664    FFtables = calcControls['FFtables']
     665    BLtables = calcControls['BLtables']
     666    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
     667    FF = np.zeros(len(Tdata))
     668    if 'N' in calcControls[hfx+'histType']:
     669        FP = 0.
     670        FPP = 0.
     671    else:
     672        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
     673        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
     674    Uij = np.array(G2lat.U6toUij(Uijdata))
     675    bij = Mast*Uij.T
     676    dFdvDict = {}
     677    dFdfr = np.zeros((len(refList),len(Mdata)))
     678    dFdx = np.zeros((len(refList),len(Mdata),3))
     679    dFdui = np.zeros((len(refList),len(Mdata)))
     680    dFdua = np.zeros((len(refList),len(Mdata),6))
     681    dFdbab = np.zeros((len(refList),2))
     682    for iref,refl in enumerate(refList):
     683        H = np.array(refl[:3])
     684        SQ = 1./(2.*refl[4])**2             # or (sin(theta)/lambda)**2
     685        SQfactor = 8.0*SQ*np.pi**2
     686        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
     687        Bab = parmDict[phfx+'BabA']*dBabdA
     688        for i,El in enumerate(Tdata):           
     689            FF[i] = refl[-1][El]           
     690        Uniq = refl[11]
     691        phi = refl[12]
     692        phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+phi[np.newaxis,:])
     693        sinp = np.sin(phase)
     694        cosp = np.cos(phase)
     695        occ = Mdata*Fdata/len(Uniq)
     696        biso = -SQfactor*Uisodata
     697        Tiso = np.where(biso<1.,np.exp(biso),1.0)
     698        HbH = -np.inner(H,np.inner(bij,H))
     699        Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
     700        Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
     701        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
     702        Tcorr = Tiso*Tuij
     703        fot = (FF+FP-Bab)*occ*Tcorr
     704        fotp = FPP*occ*Tcorr
     705        fa = np.array([fot[:,np.newaxis]*cosp,fotp[:,np.newaxis]*cosp])       #non positions
     706        fb = np.array([fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])
     707       
     708        fas = np.sum(np.sum(fa,axis=1),axis=1)
     709        fbs = np.sum(np.sum(fb,axis=1),axis=1)
     710        fax = np.array([-fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])   #positions
     711        fbx = np.array([fot[:,np.newaxis]*cosp,-fot[:,np.newaxis]*cosp])
     712        #sum below is over Uniq
     713        dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)
     714        dfadx = np.sum(twopi*Uniq*fax[:,:,:,np.newaxis],axis=2)
     715        dfadui = np.sum(-SQfactor*fa,axis=2)
     716        dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
     717        dfadba = np.sum(-cosp*(occ*Tcorr)[:,np.newaxis],axis=1)
     718        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for     
     719        dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)
     720        dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])
     721        dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])
     722        dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])
     723        dFdbab[iref] = np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
     724        if not SGData['SGInv']:
     725            dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2)        #problem here if occ=0 for some atom
     726            dfbdx = np.sum(twopi*Uniq*fbx[:,:,:,np.newaxis],axis=2)         
     727            dfbdui = np.sum(-SQfactor*fb,axis=2)
     728            dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
     729            dfbdba = np.sum(-sinp*(occ*Tcorr)[:,np.newaxis],axis=1)
     730            dFdfr[iref] += 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
     731            dFdx[iref] += 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
     732            dFdui[iref] += 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
     733            dFdua[iref] += 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
     734            dFdbab[iref] += np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
     735        #loop over atoms - each dict entry is list of derivatives for all the reflections
     736    for i in range(len(Mdata)):     
     737        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
     738        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
     739        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
     740        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
     741        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
     742        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
     743        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
     744        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
     745        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
     746        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
     747        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
     748        dFdvDict[pfx+'BabA'] = dFdbab.T[0]
     749        dFdvDict[pfx+'BabU'] = dFdbab.T[1]
     750    return dFdvDict
     751   
     752def StructureFactorDerv2(refList,G,hfx,pfx,SGData,calcControls,parmDict):
     753    '''Needs a doc string - want to do this for blocks of reflections
     754    & not one by one.'''
    595755    twopi = 2.0*np.pi
    596756    twopisq = 2.0*np.pi**2
Note: See TracChangeset for help on using the changeset viewer.