Changeset 1613

Timestamp:

Dec 19, 2014 9:04:21 AM (8 years ago)

Author:

vondreele

Message:

missing menu items in PDF controls
trap PDF plotting if no PDF to plot
some more on SS refinement constraints
Begin implementation of SStructureFactor & SStructureFactorDerv

Location:

trunk

Files:

• GSASIIgrid.py (modified) (1 diff)
• GSASIIplot.py (modified) (1 diff)
• GSASIIpwdGUI.py (modified) (1 diff)
• GSASIIspc.py (modified) (2 diffs)
• GSASIIstrMath.py (modified) (3 diffs)

trunk/GSASIIgrid.py

@@ -3553 +3553 @@
         self.PDFEdit.Append(help='Copy PDF controls', id=wxID_PDFCOPYCONTROLS, kind=wx.ITEM_NORMAL,
             text='Copy controls')
-        #        self.PDFEdit.Append(help='Load PDF controls from file',id=wxID_PDFLOADCONTROLS, kind=wx.ITEM_NORMAL,
-        #            text='Load Controls')
-        #        self.PDFEdit.Append(help='Save PDF controls to file', id=wxID_PDFSAVECONTROLS, kind=wx.ITEM_NORMAL,
-        #            text='Save controls')
+        self.PDFEdit.Append(help='Load PDF controls from file',id=wxID_PDFLOADCONTROLS, kind=wx.ITEM_NORMAL,
+            text='Load Controls')
+        self.PDFEdit.Append(help='Save PDF controls to file', id=wxID_PDFSAVECONTROLS, kind=wx.ITEM_NORMAL,
+            text='Save controls')
         self.PDFEdit.Append(help='Compute PDF', id=wxID_PDFCOMPUTE, kind=wx.ITEM_NORMAL,
             text='Compute PDF')

trunk/GSASIIplot.py

@@ -1923 +1923 @@
     lenX = 0
     for Pattern in PlotList:
-        xye = Pattern[1]
+        try:
+            xye = Pattern[1]
+        except IndexError:
+            return
         Ymax = max(Ymax,max(xye[1]))
     offset = G2frame.Offset[0]*Ymax/100.0

trunk/GSASIIpwdGUI.py

@@ -4271 +4271 @@
                 
     def OnSavePDFControls(event):
-        print 'save PDF controls?'
+        print 'save PDF controls? TBD'
         
     def OnLoadPDFControls(event):
-        print 'Load PDF controls?'
+        print 'Load PDF controls? TBD'
         
     def OnAddElement(event):

trunk/GSASIIspc.py

@@ -1426 +1426 @@
         'Sadp':[[0,1,2,3,4,5, 0,1,2,3,4,5],[1.,1.,1.,1.,1.,1., 1.,1.,1.,1.,1.,1.]],
         'Smag':[[0,1,2, 0,1,2],[1.,1.,1., 1.,1.,1.]]}
-    deltx = np.eye((3))*.001
+    deltx = np.ones((3,4))*.01
+    deltx[:3,:3] = np.eye((3))*.001
     deltu = np.eye((6))*.0001
     xyz = np.array(XYZ)%1.
+    xyzt = np.array(XYZ+[0,])%1.
     uij = np.array(UIJ)
     SGOps = SGData['SGOps']
@@ -1524 +1526 @@
 #    elif SGData['SGLaue'] in ['6/m','6/mmm']:
 #        
-    xsin = np.zeros(3)
-    xcos = np.zeros(3)
-    usin = np.zeros(6)
-    ucos = np.zeros(6)
+    xsin = np.zeros(3,dtype='i')
+    xcos = np.zeros(3,dtype='i')
+    usin = np.zeros(6,dtype='i')
+    ucos = np.zeros(6,dtype='i')
+    csi = np.ones((6),dtype='i')*-1
     for i,idelt in enumerate(deltx):
-        nxyz = (np.inner(sop[0],(xyz+idelt))+sop[1])%1.
-        print 'nxyz',i,nxyz
-        xsin[i] = np.equal((xyz-idelt)%1.,nxyz)[i]
-        print 'sin',(xyz-idelt)%1.
-        xcos[i] = np.equal((xyz+idelt)%1.,nxyz)[i]
-        print 'cos',(xyz+idelt)%1.
+        print 'idelt:',idelt
+        nxyzt = np.inner(ssop[0],(xyzt+idelt))+ssop[1]
+        nxyzt[3] -= ssop[1][3]
+        print 'nxyz',nxyzt
+        xsin[i] = np.allclose((xyzt-idelt),nxyzt,1.e-6)
+        print 'sin ',(xyzt-idelt),xsin[i]
+        xcos[i] = np.allclose((xyzt+idelt),nxyzt,1.e-6)
+        print 'cos ',(xyzt+idelt),xcos[i]
+    n = -1
+    for i,isin in enumerate(xsin):
+        if isin:
+            n += 1
+            csi[i] = n
+    for i,icos in enumerate(xcos):
+        if icos:
+            n += 1
+            csi[i+3] = n
+    print csi
     print CSI['Spos'][0]
     print xsin,xcos

trunk/GSASIIstrMath.py

@@ -612 +612 @@
         refl[10] = atan2d(fbs[0],fas[0])
     
+def SStructureFactor(refDict,im,G,hfx,pfx,SGData,calcControls,parmDict):
+    ''' 
+    Compute super structure factors for all h,k,l,m for phase
+    puts the result, F^2, in each ref[8+im] in refList
+    input:
+    
+    :param dict refDict: where
+        'RefList' list where each ref = h,k,l,m,d,...
+        'FF' dict of form factors - filed in below
+    :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)
+    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
+    SGT = np.array([ops[1] for ops in SGData['SGOps']])
+    FFtables = calcControls['FFtables']
+    BLtables = calcControls['BLtables']
+    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
+    FF = np.zeros(len(Tdata))
+    if 'NC' in calcControls[hfx+'histType']:
+        FP,FPP = G2el.BlenResCW(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
+    if not len(refDict['FF']):
+        if 'N' in calcControls[hfx+'histType']:
+            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
+        else:
+            dat = G2el.getFFvalues(FFtables,0.)        
+        refDict['FF']['El'] = dat.keys()
+        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))   
+    for iref,refl in enumerate(refDict['RefList']):
+        if 'NT' in calcControls[hfx+'histType']:
+            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl[14+im])
+        fbs = np.array([0,0])
+        H = refl[:4]
+        if refl[3]:
+            continue
+        SQ = 1./(2.*refl[4+im])**2
+        SQfactor = 4.0*SQ*twopisq
+        Bab = parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor)
+        if not np.any(refDict['FF']['FF'][iref]):                #no form factors - 1st time thru StructureFactor
+            if 'N' in calcControls[hfx+'histType']:
+                dat = G2el.getBLvalues(BLtables)
+                refDict['FF']['FF'][iref] = dat.values()
+            else:       #'X'
+                dat = G2el.getFFvalues(FFtables,SQ)
+                refDict['FF']['FF'][iref] = dat.values()
+        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
+        FF = refDict['FF']['FF'][iref][Tindx]
+        Uniq = np.inner(H[:3],SGMT)
+        Phi = np.inner(H[:3],SGT)
+        phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+Phi[:,np.newaxis])
+        sinp = np.sin(phase)
+        cosp = np.cos(phase)
+        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*Mdata*Fdata/len(Uniq)
+        fa = np.array([(FF+FP-Bab)*cosp*Tcorr,-FPP*sinp*Tcorr])
+        fas = np.sum(np.sum(fa,axis=1),axis=1)        #real
+        if not SGData['SGInv']:
+            fb = np.array([(FF+FP-Bab)*sinp*Tcorr,FPP*cosp*Tcorr])
+            fbs = np.sum(np.sum(fb,axis=1),axis=1)
+        fasq = fas**2
+        fbsq = fbs**2        #imaginary
+        refl[9+im] = np.sum(fasq)+np.sum(fbsq)
+        refl[10+im] = atan2d(fbs[0],fas[0])
+    
 def StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
     ''' Compute structure factors for all h,k,l for phase
@@ -735 +815 @@
         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
+        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
+        FF = refDict['FF']['FF'][iref].T[Tindx]
+        Uniq = np.inner(H,SGMT)
+        Phi = np.inner(H,SGT)
+        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)        #Fdata != 0 ever avoids /0. problem
+        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 for al2O3!    
+        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] = 2.*fas[0]*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)        #Fdata != 0 ever avoids /0. problem
+            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] += 2.*fbs[0]*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 SStructureFactorDerv(refDict,im,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)
+    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
+    SGT = np.array([ops[1] for ops in SGData['SGOps']])
+    FFtables = calcControls['FFtables']
+    BLtables = calcControls['BLtables']
+    nRef = len(refDict['RefList'])
+    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
+    mSize = len(Mdata)
+    FF = np.zeros(len(Tdata))
+    if 'NC' in calcControls[hfx+'histType']:
+        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
+    elif 'X' in calcControls[hfx+'histType']:
+        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((nRef,mSize))
+    dFdx = np.zeros((nRef,mSize,3))
+    dFdui = np.zeros((nRef,mSize))
+    dFdua = np.zeros((nRef,mSize,6))
+    dFdbab = np.zeros((nRef,2))
+    for iref,refl in enumerate(refDict['RefList']):
+        if 'T' in calcControls[hfx+'histType']:
+            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12+im])
+        H = np.array(refl[:4])
+        if H[3]:
+            continue
+        SQ = 1./(2.*refl[4+im])**2             # or (sin(theta)/lambda)**2
         SQfactor = 8.0*SQ*np.pi**2
         dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
@@ -2352 +2534 @@
             refDict = Histogram['Data']
             time0 = time.time()
-            StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
+            if im:
+                SStructureFactor(refDict,im,G,hfx,pfx,SGData,calcControls,parmDict)
+            else:
+                StructureFactor(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
+#                StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
 #            print 'sf-calc time: %.3f'%(time.time()-time0)
             df = np.zeros(len(refDict['RefList']))