Changeset 1103

Timestamp:

Oct 13, 2013 7:55:07 AM (10 years ago)

Author:

vondreele

Message:

fixes to protein restraints & rigid body GUI issues

Location:

trunk

Files:

• GSASIIrestrGUI.py (modified) (3 diffs)
• GSASIIstrMath.py (modified) (1 diff)

trunk/GSASIIrestrGUI.py

@@ -592 +592 @@
                                         continue
                         else:
-                            for res,name,id in Chains[chain][residue]:
-                                try:
-                                    ipos = Atms.index(name)
-                                    ids[ipos] = id
-                                except ValueError:
-                                    continue
-                            for res,name,id in Chains[chain][residue+1]:
-                                try:
-                                    ipos = pAtms.index(name)
-                                    ids[ipos] = id
-                                except ValueError:
-                                    continue
+                            try:
+                                for res,name,id in Chains[chain][residue]:
+                                    try:
+                                        ipos = Atms.index(name)
+                                        ids[ipos] = id
+                                    except ValueError:
+                                        continue
+                                for res,name,id in Chains[chain][residue+1]:
+                                    try:
+                                        ipos = pAtms.index(name)
+                                        ids[ipos] = id
+                                    except ValueError:
+                                        continue
+                            except KeyError:
+                                continue
                         if np.all(ids):
                             torsion = [list(ids),['1','1','1','1'],Name,Esd]
@@ -667 +670 @@
                     else:
                         for residue in residues[1:-1]:
-                            for res,name,id in Chains[chain][residue-1]:
-                                try:
-                                    ipos = mAtms.index(name)
-                                    ids[ipos] = id
-                                except ValueError:
-                                    continue
-                            for res,name,id in Chains[chain][residue+1]:
-                                try:
-                                    ipos = pAtms.index(name)
-                                    ids[ipos] = id
-                                except ValueError:
-                                    continue
-                            for res,name,id in Chains[chain][residue]:
-                                if Res == res:
+                            try:
+                                for res,name,id in Chains[chain][residue-1]:
                                     try:
-                                        ipos = Atms.index(name)
+                                        ipos = mAtms.index(name)
                                         ids[ipos] = id
                                     except ValueError:
                                         continue
-                            if np.all(ids):
-                                rama = [list(ids),['1','1','1','1','1'],Name,Esd]
-                                if rama not in ramaRestData['Ramas']:
-                                    ramaRestData['Ramas'].append(rama)
-                                ids = np.array([0,0,0,0,0])
+                                for res,name,id in Chains[chain][residue+1]:
+                                    try:
+                                        ipos = pAtms.index(name)
+                                        ids[ipos] = id
+                                    except ValueError:
+                                        continue
+                                for res,name,id in Chains[chain][residue]:
+                                    if Res == res:
+                                        try:
+                                            ipos = Atms.index(name)
+                                            ids[ipos] = id
+                                        except ValueError:
+                                            continue
+                                if np.all(ids):
+                                    rama = [list(ids),['1','1','1','1','1'],Name,Esd]
+                                    if rama not in ramaRestData['Ramas']:
+                                        ramaRestData['Ramas'].append(rama)
+                                    ids = np.array([0,0,0,0,0])
+                            except KeyError:
+                                continue
             macStr = macro.readline()
         macro.close()
@@ -1813 +1819 @@
             ramaRestData = restrData['Rama']
             UpdateRamaRestr(ramaRestData)
-            G2plt.PlotRama(G2frame,phaseName,rama,ramaName)
+            wx.CallAfter(G2plt.PlotRama,G2frame,phaseName,rama,ramaName)
         elif text == 'Chem. comp. restraints':
             G2gd.SetDataMenuBar(G2frame,G2frame.dataFrame.RestraintMenu)

trunk/GSASIIstrMath.py

@@ -591 +591 @@
         refl[10] = atan2d(fbs[0],fas[0])
     
+def StructureFactor2(refList,G,hfx,pfx,SGData,calcControls,parmDict):
+    ''' Compute structure factors for all h,k,l for phase
+    puts the result, F^2, in each ref[8] in refList - want to do this for blocks of reflections
+    & not one by one.
+    input:
+    
+    :param list refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv] 
+    :param np.array G:      reciprocal metric tensor
+    :param str pfx:    phase id string
+    :param dict SGData: space group info. dictionary output from SpcGroup
+    :param dict calcControls:
+    :param dict ParmDict:
+
+    '''        
+    twopi = 2.0*np.pi
+    twopisq = 2.0*np.pi**2
+    phfx = pfx.split(':')[0]+hfx
+    ast = np.sqrt(np.diag(G))
+    Mast = twopisq*np.multiply.outer(ast,ast)
+    FFtables = calcControls['FFtables']
+    BLtables = calcControls['BLtables']
+    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
+    FF = np.zeros(len(Tdata))
+    if 'N' in calcControls[hfx+'histType']:
+        FP,FPP = G2el.BlenRes(Tdata,BLtables,parmDict[hfx+'Lam'])
+    else:
+        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
+        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
+    Uij = np.array(G2lat.U6toUij(Uijdata))
+    bij = Mast*Uij.T
+    for refl in refList:
+        fbs = np.array([0,0])
+        H = refl[:3]
+        SQ = 1./(2.*refl[4])**2
+        SQfactor = 4.0*SQ*twopisq
+        Bab = parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor)
+        if not len(refl[-1]):                #no form factors
+            if 'N' in calcControls[hfx+'histType']:
+                refl[-1] = G2el.getBLvalues(BLtables)
+            else:       #'X'
+                refl[-1] = G2el.getFFvalues(FFtables,SQ)
+        for i,El in enumerate(Tdata):
+            FF[i] = refl[-1][El]           
+        Uniq = refl[11]
+        phi = refl[12]
+        phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+phi[:,np.newaxis])
+        sinp = np.sin(phase)
+        cosp = np.cos(phase)
+        occ = Mdata*Fdata/len(Uniq)
+        biso = -SQfactor*Uisodata
+        Tiso = np.where(biso<1.,np.exp(biso),1.0)
+        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
+        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
+        Tcorr = Tiso*Tuij
+        fa = np.array([(FF+FP-Bab)*occ*cosp*Tcorr,-FPP*occ*sinp*Tcorr])
+        fas = np.sum(np.sum(fa,axis=1),axis=1)        #real
+        if not SGData['SGInv']:
+            fb = np.array([(FF+FP-Bab)*occ*sinp*Tcorr,FPP*occ*cosp*Tcorr])
+            fbs = np.sum(np.sum(fb,axis=1),axis=1)
+        fasq = fas**2
+        fbsq = fbs**2        #imaginary
+        refl[9] = np.sum(fasq)+np.sum(fbsq)
+        refl[10] = atan2d(fbs[0],fas[0])
+    
 def StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict):
     'Needs a doc string'
+    twopi = 2.0*np.pi
+    twopisq = 2.0*np.pi**2
+    phfx = pfx.split(':')[0]+hfx
+    ast = np.sqrt(np.diag(G))
+    Mast = twopisq*np.multiply.outer(ast,ast)
+    FFtables = calcControls['FFtables']
+    BLtables = calcControls['BLtables']
+    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
+    FF = np.zeros(len(Tdata))
+    if 'N' in calcControls[hfx+'histType']:
+        FP = 0.
+        FPP = 0.
+    else:
+        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
+        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
+    Uij = np.array(G2lat.U6toUij(Uijdata))
+    bij = Mast*Uij.T
+    dFdvDict = {}
+    dFdfr = np.zeros((len(refList),len(Mdata)))
+    dFdx = np.zeros((len(refList),len(Mdata),3))
+    dFdui = np.zeros((len(refList),len(Mdata)))
+    dFdua = np.zeros((len(refList),len(Mdata),6))
+    dFdbab = np.zeros((len(refList),2))
+    for iref,refl in enumerate(refList):
+        H = np.array(refl[:3])
+        SQ = 1./(2.*refl[4])**2             # or (sin(theta)/lambda)**2
+        SQfactor = 8.0*SQ*np.pi**2
+        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
+        Bab = parmDict[phfx+'BabA']*dBabdA
+        for i,El in enumerate(Tdata):            
+            FF[i] = refl[-1][El]           
+        Uniq = refl[11]
+        phi = refl[12]
+        phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+phi[np.newaxis,:])
+        sinp = np.sin(phase)
+        cosp = np.cos(phase)
+        occ = Mdata*Fdata/len(Uniq)
+        biso = -SQfactor*Uisodata
+        Tiso = np.where(biso<1.,np.exp(biso),1.0)
+        HbH = -np.inner(H,np.inner(bij,H))
+        Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
+        Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
+        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
+        Tcorr = Tiso*Tuij
+        fot = (FF+FP-Bab)*occ*Tcorr
+        fotp = FPP*occ*Tcorr
+        fa = np.array([fot[:,np.newaxis]*cosp,fotp[:,np.newaxis]*cosp])       #non positions
+        fb = np.array([fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])
+        
+        fas = np.sum(np.sum(fa,axis=1),axis=1)
+        fbs = np.sum(np.sum(fb,axis=1),axis=1)
+        fax = np.array([-fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])   #positions
+        fbx = np.array([fot[:,np.newaxis]*cosp,-fot[:,np.newaxis]*cosp])
+        #sum below is over Uniq
+        dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)
+        dfadx = np.sum(twopi*Uniq*fax[:,:,:,np.newaxis],axis=2)
+        dfadui = np.sum(-SQfactor*fa,axis=2)
+        dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
+        dfadba = np.sum(-cosp*(occ*Tcorr)[:,np.newaxis],axis=1)
+        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for      
+        dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)
+        dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])
+        dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])
+        dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])
+        dFdbab[iref] = np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
+        if not SGData['SGInv']:
+            dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2)        #problem here if occ=0 for some atom
+            dfbdx = np.sum(twopi*Uniq*fbx[:,:,:,np.newaxis],axis=2)          
+            dfbdui = np.sum(-SQfactor*fb,axis=2)
+            dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
+            dfbdba = np.sum(-sinp*(occ*Tcorr)[:,np.newaxis],axis=1)
+            dFdfr[iref] += 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
+            dFdx[iref] += 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
+            dFdui[iref] += 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
+            dFdua[iref] += 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
+            dFdbab[iref] += np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
+        #loop over atoms - each dict entry is list of derivatives for all the reflections
+    for i in range(len(Mdata)):     
+        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
+        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
+        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
+        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
+        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
+        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
+        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
+        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
+        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
+        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
+        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
+        dFdvDict[pfx+'BabA'] = dFdbab.T[0]
+        dFdvDict[pfx+'BabU'] = dFdbab.T[1]
+    return dFdvDict
+    
+def StructureFactorDerv2(refList,G,hfx,pfx,SGData,calcControls,parmDict):
+    '''Needs a doc string - want to do this for blocks of reflections
+    & not one by one.'''
     twopi = 2.0*np.pi
     twopisq = 2.0*np.pi**2