Changeset 1613


Ignore:
Timestamp:
Dec 19, 2014 9:04:21 AM (7 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:
5 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIgrid.py

    r1601 r1613  
    35533553        self.PDFEdit.Append(help='Copy PDF controls', id=wxID_PDFCOPYCONTROLS, kind=wx.ITEM_NORMAL,
    35543554            text='Copy controls')
    3555         #        self.PDFEdit.Append(help='Load PDF controls from file',id=wxID_PDFLOADCONTROLS, kind=wx.ITEM_NORMAL,
    3556         #            text='Load Controls')
    3557         #        self.PDFEdit.Append(help='Save PDF controls to file', id=wxID_PDFSAVECONTROLS, kind=wx.ITEM_NORMAL,
    3558         #            text='Save controls')
     3555        self.PDFEdit.Append(help='Load PDF controls from file',id=wxID_PDFLOADCONTROLS, kind=wx.ITEM_NORMAL,
     3556            text='Load Controls')
     3557        self.PDFEdit.Append(help='Save PDF controls to file', id=wxID_PDFSAVECONTROLS, kind=wx.ITEM_NORMAL,
     3558            text='Save controls')
    35593559        self.PDFEdit.Append(help='Compute PDF', id=wxID_PDFCOMPUTE, kind=wx.ITEM_NORMAL,
    35603560            text='Compute PDF')
  • trunk/GSASIIplot.py

    r1609 r1613  
    19231923    lenX = 0
    19241924    for Pattern in PlotList:
    1925         xye = Pattern[1]
     1925        try:
     1926            xye = Pattern[1]
     1927        except IndexError:
     1928            return
    19261929        Ymax = max(Ymax,max(xye[1]))
    19271930    offset = G2frame.Offset[0]*Ymax/100.0
  • trunk/GSASIIpwdGUI.py

    r1595 r1613  
    42714271               
    42724272    def OnSavePDFControls(event):
    4273         print 'save PDF controls?'
     4273        print 'save PDF controls? TBD'
    42744274       
    42754275    def OnLoadPDFControls(event):
    4276         print 'Load PDF controls?'
     4276        print 'Load PDF controls? TBD'
    42774277       
    42784278    def OnAddElement(event):
  • trunk/GSASIIspc.py

    r1612 r1613  
    14261426        '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.]],
    14271427        'Smag':[[0,1,2, 0,1,2],[1.,1.,1., 1.,1.,1.]]}
    1428     deltx = np.eye((3))*.001
     1428    deltx = np.ones((3,4))*.01
     1429    deltx[:3,:3] = np.eye((3))*.001
    14291430    deltu = np.eye((6))*.0001
    14301431    xyz = np.array(XYZ)%1.
     1432    xyzt = np.array(XYZ+[0,])%1.
    14311433    uij = np.array(UIJ)
    14321434    SGOps = SGData['SGOps']
     
    15241526#    elif SGData['SGLaue'] in ['6/m','6/mmm']:
    15251527#       
    1526     xsin = np.zeros(3)
    1527     xcos = np.zeros(3)
    1528     usin = np.zeros(6)
    1529     ucos = np.zeros(6)
     1528    xsin = np.zeros(3,dtype='i')
     1529    xcos = np.zeros(3,dtype='i')
     1530    usin = np.zeros(6,dtype='i')
     1531    ucos = np.zeros(6,dtype='i')
     1532    csi = np.ones((6),dtype='i')*-1
    15301533    for i,idelt in enumerate(deltx):
    1531         nxyz = (np.inner(sop[0],(xyz+idelt))+sop[1])%1.
    1532         print 'nxyz',i,nxyz
    1533         xsin[i] = np.equal((xyz-idelt)%1.,nxyz)[i]
    1534         print 'sin',(xyz-idelt)%1.
    1535         xcos[i] = np.equal((xyz+idelt)%1.,nxyz)[i]
    1536         print 'cos',(xyz+idelt)%1.
     1534        print 'idelt:',idelt
     1535        nxyzt = np.inner(ssop[0],(xyzt+idelt))+ssop[1]
     1536        nxyzt[3] -= ssop[1][3]
     1537        print 'nxyz',nxyzt
     1538        xsin[i] = np.allclose((xyzt-idelt),nxyzt,1.e-6)
     1539        print 'sin ',(xyzt-idelt),xsin[i]
     1540        xcos[i] = np.allclose((xyzt+idelt),nxyzt,1.e-6)
     1541        print 'cos ',(xyzt+idelt),xcos[i]
     1542    n = -1
     1543    for i,isin in enumerate(xsin):
     1544        if isin:
     1545            n += 1
     1546            csi[i] = n
     1547    for i,icos in enumerate(xcos):
     1548        if icos:
     1549            n += 1
     1550            csi[i+3] = n
     1551    print csi
    15371552    print CSI['Spos'][0]
    15381553    print xsin,xcos
  • trunk/GSASIIstrMath.py

    r1606 r1613  
    612612        refl[10] = atan2d(fbs[0],fas[0])
    613613   
     614def SStructureFactor(refDict,im,G,hfx,pfx,SGData,calcControls,parmDict):
     615    '''
     616    Compute super structure factors for all h,k,l,m for phase
     617    puts the result, F^2, in each ref[8+im] in refList
     618    input:
     619   
     620    :param dict refDict: where
     621        'RefList' list where each ref = h,k,l,m,d,...
     622        'FF' dict of form factors - filed in below
     623    :param np.array G:      reciprocal metric tensor
     624    :param str pfx:    phase id string
     625    :param dict SGData: space group info. dictionary output from SpcGroup
     626    :param dict calcControls:
     627    :param dict ParmDict:
     628
     629    '''       
     630    twopi = 2.0*np.pi
     631    twopisq = 2.0*np.pi**2
     632    phfx = pfx.split(':')[0]+hfx
     633    ast = np.sqrt(np.diag(G))
     634    Mast = twopisq*np.multiply.outer(ast,ast)
     635    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
     636    SGT = np.array([ops[1] for ops in SGData['SGOps']])
     637    FFtables = calcControls['FFtables']
     638    BLtables = calcControls['BLtables']
     639    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
     640    FF = np.zeros(len(Tdata))
     641    if 'NC' in calcControls[hfx+'histType']:
     642        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
     643    else:
     644        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
     645        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
     646    Uij = np.array(G2lat.U6toUij(Uijdata))
     647    bij = Mast*Uij.T
     648    if not len(refDict['FF']):
     649        if 'N' in calcControls[hfx+'histType']:
     650            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
     651        else:
     652            dat = G2el.getFFvalues(FFtables,0.)       
     653        refDict['FF']['El'] = dat.keys()
     654        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))   
     655    for iref,refl in enumerate(refDict['RefList']):
     656        if 'NT' in calcControls[hfx+'histType']:
     657            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl[14+im])
     658        fbs = np.array([0,0])
     659        H = refl[:4]
     660        if refl[3]:
     661            continue
     662        SQ = 1./(2.*refl[4+im])**2
     663        SQfactor = 4.0*SQ*twopisq
     664        Bab = parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor)
     665        if not np.any(refDict['FF']['FF'][iref]):                #no form factors - 1st time thru StructureFactor
     666            if 'N' in calcControls[hfx+'histType']:
     667                dat = G2el.getBLvalues(BLtables)
     668                refDict['FF']['FF'][iref] = dat.values()
     669            else:       #'X'
     670                dat = G2el.getFFvalues(FFtables,SQ)
     671                refDict['FF']['FF'][iref] = dat.values()
     672        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
     673        FF = refDict['FF']['FF'][iref][Tindx]
     674        Uniq = np.inner(H[:3],SGMT)
     675        Phi = np.inner(H[:3],SGT)
     676        phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+Phi[:,np.newaxis])
     677        sinp = np.sin(phase)
     678        cosp = np.cos(phase)
     679        biso = -SQfactor*Uisodata
     680        Tiso = np.where(biso<1.,np.exp(biso),1.0)
     681        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
     682        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
     683        Tcorr = Tiso*Tuij*Mdata*Fdata/len(Uniq)
     684        fa = np.array([(FF+FP-Bab)*cosp*Tcorr,-FPP*sinp*Tcorr])
     685        fas = np.sum(np.sum(fa,axis=1),axis=1)        #real
     686        if not SGData['SGInv']:
     687            fb = np.array([(FF+FP-Bab)*sinp*Tcorr,FPP*cosp*Tcorr])
     688            fbs = np.sum(np.sum(fb,axis=1),axis=1)
     689        fasq = fas**2
     690        fbsq = fbs**2        #imaginary
     691        refl[9+im] = np.sum(fasq)+np.sum(fbsq)
     692        refl[10+im] = atan2d(fbs[0],fas[0])
     693   
    614694def StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
    615695    ''' Compute structure factors for all h,k,l for phase
     
    735815        H = np.array(refl[:3])
    736816        SQ = 1./(2.*refl[4])**2             # or (sin(theta)/lambda)**2
     817        SQfactor = 8.0*SQ*np.pi**2
     818        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
     819        Bab = parmDict[phfx+'BabA']*dBabdA
     820        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
     821        FF = refDict['FF']['FF'][iref].T[Tindx]
     822        Uniq = np.inner(H,SGMT)
     823        Phi = np.inner(H,SGT)
     824        phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+Phi[np.newaxis,:])
     825        sinp = np.sin(phase)
     826        cosp = np.cos(phase)
     827        occ = Mdata*Fdata/len(Uniq)
     828        biso = -SQfactor*Uisodata
     829        Tiso = np.where(biso<1.,np.exp(biso),1.0)
     830        HbH = -np.inner(H,np.inner(bij,H))
     831        Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
     832        Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
     833        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
     834        Tcorr = Tiso*Tuij
     835        fot = (FF+FP-Bab)*occ*Tcorr
     836        fotp = FPP*occ*Tcorr
     837        fa = np.array([fot[:,np.newaxis]*cosp,fotp[:,np.newaxis]*cosp])       #non positions
     838        fb = np.array([fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])
     839       
     840        fas = np.sum(np.sum(fa,axis=1),axis=1)
     841        fbs = np.sum(np.sum(fb,axis=1),axis=1)
     842        fax = np.array([-fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])   #positions
     843        fbx = np.array([fot[:,np.newaxis]*cosp,-fot[:,np.newaxis]*cosp])
     844        #sum below is over Uniq
     845        dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)        #Fdata != 0 ever avoids /0. problem
     846        dfadx = np.sum(twopi*Uniq*fax[:,:,:,np.newaxis],axis=2)
     847        dfadui = np.sum(-SQfactor*fa,axis=2)
     848        dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
     849        dfadba = np.sum(-cosp*(occ*Tcorr)[:,np.newaxis],axis=1)
     850        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for al2O3!   
     851        dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)
     852        dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])
     853        dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])
     854        dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])
     855        dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
     856        if not SGData['SGInv']:
     857            dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2)        #Fdata != 0 ever avoids /0. problem
     858            dfbdx = np.sum(twopi*Uniq*fbx[:,:,:,np.newaxis],axis=2)           
     859            dfbdui = np.sum(-SQfactor*fb,axis=2)
     860            dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
     861            dfbdba = np.sum(-sinp*(occ*Tcorr)[:,np.newaxis],axis=1)
     862            dFdfr[iref] += 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
     863            dFdx[iref] += 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
     864            dFdui[iref] += 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
     865            dFdua[iref] += 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
     866            dFdbab[iref] += 2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
     867        #loop over atoms - each dict entry is list of derivatives for all the reflections
     868    for i in range(len(Mdata)):     
     869        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
     870        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
     871        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
     872        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
     873        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
     874        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
     875        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
     876        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
     877        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
     878        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
     879        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
     880    dFdvDict[pfx+'BabA'] = dFdbab.T[0]
     881    dFdvDict[pfx+'BabU'] = dFdbab.T[1]
     882    return dFdvDict
     883   
     884def SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,calcControls,parmDict):
     885    'Needs a doc string'
     886    twopi = 2.0*np.pi
     887    twopisq = 2.0*np.pi**2
     888    phfx = pfx.split(':')[0]+hfx
     889    ast = np.sqrt(np.diag(G))
     890    Mast = twopisq*np.multiply.outer(ast,ast)
     891    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
     892    SGT = np.array([ops[1] for ops in SGData['SGOps']])
     893    FFtables = calcControls['FFtables']
     894    BLtables = calcControls['BLtables']
     895    nRef = len(refDict['RefList'])
     896    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
     897    mSize = len(Mdata)
     898    FF = np.zeros(len(Tdata))
     899    if 'NC' in calcControls[hfx+'histType']:
     900        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
     901    elif 'X' in calcControls[hfx+'histType']:
     902        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
     903        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
     904    Uij = np.array(G2lat.U6toUij(Uijdata))
     905    bij = Mast*Uij.T
     906    dFdvDict = {}
     907    dFdfr = np.zeros((nRef,mSize))
     908    dFdx = np.zeros((nRef,mSize,3))
     909    dFdui = np.zeros((nRef,mSize))
     910    dFdua = np.zeros((nRef,mSize,6))
     911    dFdbab = np.zeros((nRef,2))
     912    for iref,refl in enumerate(refDict['RefList']):
     913        if 'T' in calcControls[hfx+'histType']:
     914            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12+im])
     915        H = np.array(refl[:4])
     916        if H[3]:
     917            continue
     918        SQ = 1./(2.*refl[4+im])**2             # or (sin(theta)/lambda)**2
    737919        SQfactor = 8.0*SQ*np.pi**2
    738920        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
     
    23522534            refDict = Histogram['Data']
    23532535            time0 = time.time()
    2354             StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
     2536            if im:
     2537                SStructureFactor(refDict,im,G,hfx,pfx,SGData,calcControls,parmDict)
     2538            else:
     2539                StructureFactor(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
     2540#                StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
    23552541#            print 'sf-calc time: %.3f'%(time.time()-time0)
    23562542            df = np.zeros(len(refDict['RefList']))
Note: See TracChangeset for help on using the changeset viewer.