Ignore:
Timestamp:
Dec 16, 2015 2:04:20 PM (7 years ago)
Author:
vondreele
Message:

implement reflection block processing for sf derivatives for non twinned data
block size optimized at 32 for sf derivatives & 100 for sf calcs

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIstrMath.py

    r2089 r2094  
    3838twopi = 2.0*np.pi
    3939twopisq = 2.0*np.pi**2
     40nxs = np.newaxis
    4041
    4142################################################################################
    … …  
    686687    Uij = np.array(G2lat.U6toUij(Uijdata))
    687688    bij = Mast*Uij.T
    688     blkSize = 100       #no. of reflections in a block
     689    blkSize = 100       #no. of reflections in a block - size seems optimal
    689690    nRef = refDict['RefList'].shape[0]
    690691    if not len(refDict['FF']):                #no form factors - 1st time thru StructureFactor
    … …  
    703704#reflection processing begins here - big arrays!
    704705    iBeg = 0
     706    time0 = time.time()
    705707    while iBeg < nRef:
    706708        iFin = min(iBeg+blkSize,nRef)
    … …  
    734736        sinp = np.sin(phase)
    735737        cosp = np.cos(phase)
    736         biso = -SQfactor*Uisodata[:,np.newaxis]
     738        biso = -SQfactor*Uisodata[:,nxs]
    737739        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*len(TwinLaw),axis=1).T
    738740        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
    … …  
    755757            if len(TwinLaw) > 1:
    756758                refl.T[9] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2   #FcT from primary twin element
    757                 refl.T[7] = np.sum(TwinFr*np.sum(TwMask[np.newaxis,:,:]*fas,axis=0)**2,axis=-1)+   \
    758                     np.sum(TwinFr*np.sum(TwMask[np.newaxis,:,:]*fbs,axis=0)**2,axis=-1)                        #Fc sum over twins
     759                refl.T[7] = np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fas,axis=0)**2,axis=-1)+   \
     760                    np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fbs,axis=0)**2,axis=-1)                        #Fc sum over twins
    759761                refl.T[10] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f"
    760762            else:
    … …  
    764766#        refl.T[10] = atan2d(np.sum(fbs,axis=0),np.sum(fas,axis=0)) #include f' & f"
    765767        iBeg += blkSize
     768#    print ' %d sf time %.4f\r'%(nRef,time.time()-time0)
    766769   
    767770def StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
    768     'Needs a doc string'
     771    '''Compute structure factor derivatives on blocks of reflections - keep as it works for twins
     772    but is slower for powders/nontwins
     773    '''
    769774    phfx = pfx.split(':')[0]+hfx
    770775    ast = np.sqrt(np.diag(G))
    … …  
    840845        cosp = np.cos(phase)
    841846        occ = Mdata*Fdata/len(SGT)
    842         biso = -SQfactor*Uisodata[:,np.newaxis]
     847        biso = -SQfactor*Uisodata[:,nxs]
    843848        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*nTwin,axis=1)
    844849        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
    … …  
    851856        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
    852857        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
     858#        GSASIIpath.IPyBreak()
    853859        fas = np.sum(np.sum(fa,axis=-1),axis=-1)      #real sum over atoms & unique hkl array(2,nTwins)
    854860        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)      #imag sum over atoms & uniq hkl
    … …  
    857863        #sum below is over Uniq
    858864        dfadfr = np.sum(fa/occ,axis=-2)        #array(2,ntwin,nAtom) Fdata != 0 avoids /0. problem
    859         dfadba = np.sum(-cosp*Tcorr[:,np.newaxis],axis=1)
     865        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
    860866        dfadui = np.sum(-SQfactor*fa,axis=-2)
    861867        if nTwin > 1:
    862             dfadx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fax,-2,-1)[:,it,:,:,np.newaxis],axis=-2) for it in range(nTwin)])
    863             dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fa,-2,-1)[:,it,:,:,np.newaxis],axis=-2) for it in range(nTwin)])
     868            dfadx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
     869            dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fa,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
    864870            # array(nTwin,2,nAtom,3) & array(nTwin,2,nAtom,6)
    865871        else:
    866             dfadx = np.sum(twopi*Uniq*np.swapaxes(fax,-2,-1)[:,:,:,np.newaxis],axis=-2)
    867             dfadua = np.sum(-Hij*np.swapaxes(fa,-2,-1)[:,:,:,np.newaxis],axis=-2)
     872            dfadx = np.sum(twopi*Uniq*np.swapaxes(fax,-2,-1)[:,:,:,nxs],axis=-2)
     873            dfadua = np.sum(-Hij*np.swapaxes(fa,-2,-1)[:,:,:,nxs],axis=-2)
    868874            # array(2,nAtom,3) & array(2,nAtom,6)
    869875        if not SGData['SGInv']:
    870876            dfbdfr = np.sum(fb/occ,axis=-2)        #Fdata != 0 avoids /0. problem
    871             dfbdba = np.sum(-sinp*Tcorr[:,np.newaxis],axis=1)
     877            dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)
    872878            dfbdui = np.sum(-SQfactor*fb,axis=-2)
    873879            if len(TwinLaw) > 1:
    874                 dfbdx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fbx,-2,-1)[:,it,:,:,np.newaxis],axis=2) for it in range(nTwin)])           
    875                 dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fb,-2,-1)[:,it,:,:,np.newaxis],axis=2) for it in range(nTwin)])
     880                dfbdx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)])           
     881                dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fb,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)])
    876882            else:
    877883                dfadfl = np.sum(-FPP*Tcorr*sinp)
    878884                dfbdfl = np.sum(FPP*Tcorr*cosp)
    879                 dfbdx = np.sum(twopi*Uniq*np.swapaxes(fbx,-2,-1)[:,:,:,np.newaxis],axis=2)           
    880                 dfbdua = np.sum(-Hij*np.swapaxes(fb,-2,-1)[:,:,:,np.newaxis],axis=2)
     885                dfbdx = np.sum(twopi*Uniq*np.swapaxes(fbx,-2,-1)[:,:,:,nxs],axis=2)           
     886                dfbdua = np.sum(-Hij*np.swapaxes(fb,-2,-1)[:,:,:,nxs],axis=2)
    881887        else:
    882888            dfbdfr = np.zeros_like(dfadfr)
    … …  
    918924        if not iref%100 :
    919925            print ' %d derivative time %.4f\r'%(iref,time.time()-time0),
     926    print ' %d derivative time %.4f\r'%(len(refDict['RefList']),time.time()-time0)
    920927        #loop over atoms - each dict entry is list of derivatives for all the reflections
    921928    if nTwin > 1:
    922929        for i in range(len(Mdata)):     #these all OK?
    923             dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,np.newaxis],axis=0)
    924             dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,np.newaxis],axis=0)
    925             dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,np.newaxis],axis=0)
    926             dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,np.newaxis],axis=0)
    927             dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,np.newaxis],axis=0)
    928             dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,np.newaxis],axis=0)
    929             dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,np.newaxis],axis=0)
    930             dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,np.newaxis],axis=0)
    931             dFdvDict[pfx+'AU12:'+str(i)] = 2.*np.sum(dFdua.T[3][i]*TwinFr[:,np.newaxis],axis=0)
    932             dFdvDict[pfx+'AU13:'+str(i)] = 2.*np.sum(dFdua.T[4][i]*TwinFr[:,np.newaxis],axis=0)
    933             dFdvDict[pfx+'AU23:'+str(i)] = 2.*np.sum(dFdua.T[5][i]*TwinFr[:,np.newaxis],axis=0)
     930            dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,nxs],axis=0)
     931            dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,nxs],axis=0)
     932            dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,nxs],axis=0)
     933            dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,nxs],axis=0)
     934            dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,nxs],axis=0)
     935            dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,nxs],axis=0)
     936            dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,nxs],axis=0)
     937            dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,nxs],axis=0)
     938            dFdvDict[pfx+'AU12:'+str(i)] = 2.*np.sum(dFdua.T[3][i]*TwinFr[:,nxs],axis=0)
     939            dFdvDict[pfx+'AU13:'+str(i)] = 2.*np.sum(dFdua.T[4][i]*TwinFr[:,nxs],axis=0)
     940            dFdvDict[pfx+'AU23:'+str(i)] = 2.*np.sum(dFdua.T[5][i]*TwinFr[:,nxs],axis=0)
    934941    else:
    935942        for i in range(len(Mdata)):
    … …  
    953960    return dFdvDict
    954961   
     962def StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
     963    '''Compute structure factor derivatives on blocks of reflections- works for powders/nontwins
     964    faster than StructureFactorDerv
     965    '''
     966    phfx = pfx.split(':')[0]+hfx
     967    ast = np.sqrt(np.diag(G))
     968    Mast = twopisq*np.multiply.outer(ast,ast)
     969    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
     970    SGT = np.array([ops[1] for ops in SGData['SGOps']])
     971    FFtables = calcControls['FFtables']
     972    BLtables = calcControls['BLtables']
     973    TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])
     974    TwDict = refDict.get('TwDict',{})           
     975    if 'S' in calcControls[hfx+'histType']:
     976        NTL = calcControls[phfx+'NTL']
     977        NM = calcControls[phfx+'TwinNMN']+1
     978        TwinLaw = calcControls[phfx+'TwinLaw']
     979        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
     980        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
     981    nTwin = len(TwinLaw)       
     982    nRef = len(refDict['RefList'])
     983    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
     984    mSize = len(Mdata)
     985    FF = np.zeros(len(Tdata))
     986    if 'NC' in calcControls[hfx+'histType']:
     987        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
     988    elif 'X' in calcControls[hfx+'histType']:
     989        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
     990        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
     991    Uij = np.array(G2lat.U6toUij(Uijdata))
     992    bij = Mast*Uij.T
     993    dFdvDict = {}
     994    if nTwin > 1:
     995        dFdfr = np.zeros((nRef,nTwin,mSize))
     996        dFdx = np.zeros((nRef,nTwin,mSize,3))
     997        dFdui = np.zeros((nRef,nTwin,mSize))
     998        dFdua = np.zeros((nRef,nTwin,mSize,6))
     999        dFdbab = np.zeros((nRef,nTwin,2))
     1000        dFdfl = np.zeros((nRef,nTwin))
     1001        dFdtw = np.zeros((nRef,nTwin))
     1002    else:
     1003        dFdfr = np.zeros((nRef,mSize))
     1004        dFdx = np.zeros((nRef,mSize,3))
     1005        dFdui = np.zeros((nRef,mSize))
     1006        dFdua = np.zeros((nRef,mSize,6))
     1007        dFdbab = np.zeros((nRef,2))
     1008        dFdfl = np.zeros((nRef))
     1009        dFdtw = np.zeros((nRef))
     1010    Flack = 1.0
     1011    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
     1012        Flack = 1.-2.*parmDict[phfx+'Flack']
     1013    time0 = time.time()
     1014    nref = len(refDict['RefList'])/100   
     1015#reflection processing begins here - big arrays!
     1016    iBeg = 0
     1017    blkSize = 32       #no. of reflections in a block - optimized for speed
     1018    while iBeg < nRef:
     1019        iFin = min(iBeg+blkSize,nRef)
     1020        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
     1021        H = refl.T[:3]
     1022        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(3,nTwins) or (3)
     1023        TwMask = np.any(H,axis=-1)
     1024        if TwinLaw.shape[0] > 1 and TwDict:
     1025            for ir in range(blkSize):
     1026                iref = ir+iBeg
     1027                if iref in TwDict:
     1028                    for i in TwDict[iref]:
     1029                        for n in range(NTL):
     1030                            H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
     1031            TwMask = np.any(H,axis=-1)
     1032        SQ = 1./(2.*refl.T[4])**2             # or (sin(theta)/lambda)**2
     1033        SQfactor = 8.0*SQ*np.pi**2
     1034        if 'T' in calcControls[hfx+'histType']:
     1035            if 'P' in calcControls[hfx+'histType']:
     1036                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
     1037            else:
     1038                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
     1039            FP = np.repeat(FP.T,len(SGT)*len(TwinLaw),axis=0)
     1040            FPP = np.repeat(FPP.T,len(SGT)*len(TwinLaw),axis=0)
     1041        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
     1042        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT)*len(TwinLaw))
     1043        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
     1044        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT)*len(TwinLaw),axis=0)
     1045        Uniq = np.inner(H,SGMT)             # array(nSGOp,3) or (nTwin,nSGOp,3)
     1046        Phi = np.inner(H,SGT)
     1047        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
     1048        sinp = np.sin(phase)
     1049        cosp = np.cos(phase)
     1050        occ = Mdata*Fdata/len(SGT)
     1051        biso = -SQfactor*Uisodata[:,nxs]
     1052        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*len(TwinLaw),axis=1).T
     1053        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
     1054        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
     1055        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/len(SGMT)
     1056        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])
     1057        if nTwin > 1:
     1058            Hij = np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(-1,nTwin,len(SGT),6))
     1059        else:
     1060            Hij = np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(-1,len(SGT),6))
     1061        fot = np.reshape(((FF+FP).T-Bab).T,cosp.shape)*Tcorr
     1062        fotp = FPP*Tcorr       
     1063        if 'T' in calcControls[hfx+'histType']:
     1064            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
     1065            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr])
     1066        else:
     1067            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
     1068            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
     1069#        GSASIIpath.IPyBreak()
     1070        fas = np.sum(np.sum(fa,axis=-1),axis=-1)      #real sum over atoms & unique hkl array(2,refBlk,nTwins)
     1071        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)      #imag sum over atoms & uniq hkl
     1072        fax = np.array([-fot*sinp,-fotp*cosp])   #positions array(2,refBlk,ntwi,nEqv,nAtoms)
     1073        fbx = np.array([fot*cosp,-fotp*sinp])
     1074        #sum below is over Uniq
     1075        dfadfr = np.sum(fa/occ,axis=-2)        #array(2,refBlk,ntwin,nAtom) Fdata != 0 avoids /0. problem
     1076        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=-2)  #array(refBlk,nTwin,nOps,nAtom)
     1077        if nTwin > 1:
     1078            dfadfr = np.swapaxes(dfadfr,0,1)
     1079            dfadx = np.swapaxes(np.array([np.sum(twopi*Uniq[it,nxs,:,nxs,:,:]*np.swapaxes(fax,-2,-1)[it,:,:,:,:,nxs],axis=-2) for it in range(nTwin)]),0,1)
     1080            dfadui = np.swapaxes(np.sum(-SQfactor[nxs,:,nxs,nxs,nxs]*fa,axis=-2),0,1) #array(Ops,refBlk,nTw,nAtoms)
     1081            dfadua = np.swapaxes(np.array([np.sum(-Hij[it,nxs,:,nxs,:,:]*np.swapaxes(fa,-2,-1)[it,:,:,:,:,nxs],axis=-2) for it in range(nTwin)]),0,1)
     1082            # array(nTwin,2,refBlk,nAtom,3) & array(nTwin,2,refBlk,nAtom,6)
     1083        else:
     1084            dfadx = np.sum(twopi*Uniq[nxs,:,nxs,:,:]*np.swapaxes(fax,-2,-1)[:,:,:,:,nxs],axis=-2)
     1085            dfadui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fa,axis=-2) #array(Ops,refBlk,nTw,nAtoms)
     1086            dfadua = np.sum(-Hij[nxs,:,nxs,:,:]*np.swapaxes(fa,-2,-1)[:,:,:,:,nxs],axis=-2)
     1087            # array(2,refBlk,nAtom,3) & array(2,refBlk,nAtom,6)
     1088        if not SGData['SGInv']:
     1089            dfbdfr = np.sum(fb/occ,axis=-2)        #Fdata != 0 avoids /0. problem
     1090            dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=-2)
     1091            if len(TwinLaw) > 1:
     1092                dfbdfr = np.swapaxes(dfbdfr,0,1)
     1093                dfbdx = np.swapaxes(np.array([np.sum(twopi*Uniq[it,nxs,:,nxs,:,:]*np.swapaxes(fbx,-2,-1)[it,:,:,:,:,nxs],axis=-2) for it in range(nTwin)]),0,1)         
     1094                dfbdui = np.swapaxes(np.sum(-SQfactor[nxs,:,nxs,nxs,nxs]*fb,axis=-2),0,1)      #array(Ops,refBlk,nTw,nAtoms)
     1095                dfbdua = np.swapaxes(np.array([np.sum(-Hij[it,nxs,:,nxs,:,:]*np.swapaxes(fb,-2,-1)[it,:,:,:,:,nxs],axis=-2) for it in range(nTwin)]),0,1)
     1096            else:
     1097                dfadfl = np.sum(np.sum(-FPP*Tcorr*sinp,axis=-1),axis=-1)
     1098                dfbdfl = np.sum(np.sum(FPP*Tcorr*cosp,axis=-1),axis=-1)
     1099                dfbdx = np.sum(twopi*Uniq[nxs,:,nxs,:,:]*np.swapaxes(fbx,-2,-1)[:,:,:,:,nxs],axis=-2)           
     1100                dfbdui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fb,axis=-2)
     1101                dfbdua = np.sum(-Hij[nxs,:,nxs,:,:]*np.swapaxes(fb,-2,-1)[:,:,:,:,nxs],axis=-2)
     1102        else:
     1103            dfbdfr = np.zeros_like(dfadfr)
     1104            dfbdx = np.zeros_like(dfadx)
     1105            dfbdui = np.zeros_like(dfadui)
     1106            dfbdua = np.zeros_like(dfadua)
     1107            dfbdba = np.zeros_like(dfadba)
     1108            dfadfl = 0.0
     1109            dfbdfl = 0.0
     1110        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
     1111        SA = fas[0]+fas[1]
     1112        SB = fbs[0]+fbs[1]
     1113        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro
     1114            dFdfr[iBeg:iFin] = 2.*(fas[0,:,nxs]*dfadfr[0]+fas[1,:,nxs]*dfadfr[1])*Mdata/len(SGMT)+   \
     1115                2.*(fbs[0,:,nxs]*dfbdfr[0]-fbs[1,:,nxs]*dfbdfr[1])*Mdata/len(SGMT)
     1116            dFdx[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadx[0]+fas[1,:,nxs,nxs]*dfadx[1])+  \
     1117                2.*(fbs[0,:,nxs,nxs]*dfbdx[0]+fbs[1,:,nxs,nxs]*dfbdx[1])
     1118            dFdui[iBeg:iFin] = 2.*(fas[0,:,nxs]*dfadui[0]+fas[1,:,nxs]*dfadui[1])+   \
     1119                2.*(fbs[0,:,nxs]*dfbdui[0]-fbs[1,:,nxs]*dfbdui[1])
     1120            dFdua[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadua[0]+fas[1,:,nxs,nxs]*dfadua[1])+   \
     1121                2.*(fbs[0,:,nxs,nxs]*dfbdua[0]+fbs[1,:,nxs,nxs]*dfbdua[1])
     1122        else:
     1123            if nTwin > 1:
     1124                dFdfr[iBeg:iFin] = np.swapaxes(np.array([2.*TwMask[:,it,nxs]*(SA[:,it,nxs]*(np.sum(dfadfr[:,:,it],axis=1))+ \
     1125                    SB[:,it,nxs]*(np.sum(dfbdfr[:,:,it],axis=1)))*Mdata/len(SGMT) for it in range(nTwin)]),0,1)
     1126                dFdx[iBeg:iFin] = np.swapaxes(np.array([2.*TwMask[:,it,nxs,nxs]*(SA[:,it,nxs,nxs]*np.sum(dfadx[:,:,it],axis=1)+    \
     1127                    SB[:,it,nxs,nxs]*np.sum(dfbdx[:,:,it],axis=1)) for it in range(nTwin)]),0,1)
     1128                dFdui[iBeg:iFin] = np.swapaxes(np.array([2.*TwMask[:,it,nxs]*(SA[:,it,nxs]*np.sum(dfadui[:,:,it],axis=1)+  \
     1129                    SB[:,it,nxs]*np.sum(dfbdui[:,:,it],axis=1)) for it in range(nTwin)]),0,1)
     1130                dFdua[iBeg:iFin] = np.swapaxes(np.array([2.*TwMask[:,it,nxs,nxs]*(SA[:,it,nxs,nxs]*np.sum(dfadua[:,:,it],axis=1)+  \
     1131                    SB[:,it,nxs,nxs]*np.sum(dfbdua[:,:,it],axis=1)) for it in range(nTwin)]),0,1)
     1132                dFdtw[iBeg:iFin] = np.sum(TwMask*fas,axis=0)**2+np.sum(TwMask*fbs,axis=0)**2
     1133               
     1134            else:   #these are good for no twin single crystals
     1135                dFdfr[iBeg:iFin] = (2.*SA[:,nxs]*(dfadfr[0]+dfadfr[1])+2.*SB[:,nxs]*(dfbdfr[0]+dfbdfr[1]))*Mdata/len(SGMT)
     1136                dFdx[iBeg:iFin] = 2.*SA[:,nxs,nxs]*(dfadx[0]+dfadx[1])+2.*SB[:,nxs,nxs]*(dfbdx[0]+dfbdx[1])
     1137                dFdui[iBeg:iFin] = 2.*SA[:,nxs]*(dfadui[0]+dfadui[1])+2.*SB[:,nxs]*(dfbdui[0]+dfbdui[1])
     1138                dFdua[iBeg:iFin] = 2.*SA[:,nxs,nxs]*(dfadua[0]+dfadua[1])+2.*SB[:,nxs,nxs]*(dfbdua[0]+dfbdua[1])
     1139                dFdfl[iBeg:iFin] = -SA*dfadfl-SB*dfbdfl  #array(nRef,)
     1140#        dFdbab[iBeg:iFin] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
     1141#            2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
     1142#        GSASIIpath.IPyBreak()
     1143        iBeg += blkSize
     1144    print ' %d derivative time %.4f\r'%(nRef,time.time()-time0)
     1145        #loop over atoms - each dict entry is list of derivatives for all the reflections
     1146    if nTwin > 1:
     1147        for i in range(len(Mdata)):     #these all OK?
     1148            dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,nxs],axis=0)
     1149            dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,nxs],axis=0)
     1150            dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,nxs],axis=0)
     1151            dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,nxs],axis=0)
     1152            dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,nxs],axis=0)
     1153            dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,nxs],axis=0)
     1154            dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,nxs],axis=0)
     1155            dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,nxs],axis=0)
     1156            dFdvDict[pfx+'AU12:'+str(i)] = 2.*np.sum(dFdua.T[3][i]*TwinFr[:,nxs],axis=0)
     1157            dFdvDict[pfx+'AU13:'+str(i)] = 2.*np.sum(dFdua.T[4][i]*TwinFr[:,nxs],axis=0)
     1158            dFdvDict[pfx+'AU23:'+str(i)] = 2.*np.sum(dFdua.T[5][i]*TwinFr[:,nxs],axis=0)
     1159    else:
     1160        for i in range(len(Mdata)):
     1161            dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
     1162            dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
     1163            dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
     1164            dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
     1165            dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
     1166            dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
     1167            dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
     1168            dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
     1169            dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
     1170            dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
     1171            dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
     1172        dFdvDict[phfx+'Flack'] = 4.*dFdfl.T
     1173    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
     1174    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
     1175    if nTwin > 1:
     1176        for i in range(nTwin):
     1177            dFdvDict[phfx+'TwinFr:'+str(i)] = dFdtw.T[i]
     1178    return dFdvDict
     1179   
    9551180def SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
    9561181    '''
    … …  
    9691194
    9701195    '''
    971     nxs = np.newaxis       
    9721196    phfx = pfx.split(':')[0]+hfx
    9731197    ast = np.sqrt(np.diag(G))
    … …  
    9781202    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
    9791203    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
    980 #    eps = SSGMT[:,3,3]
    9811204    FFtables = calcControls['FFtables']
    9821205    BLtables = calcControls['BLtables']
    … …  
    10421265        UniqP = np.inner(HP.T,SGMT)
    10431266        Phi = np.inner(H.T,SSGT)
    1044 #        PhiP = np.inner(HP.T,SGT)
    10451267        if SGInv:   #if centro - expand HKL sets
    10461268            Uniq = np.hstack((Uniq,-Uniq))
    10471269            Phi = np.hstack((Phi,-Phi))
    10481270            UniqP = np.hstack((UniqP,-UniqP))
    1049 #            PhiP = np.hstack((PhiP,-PhiP))
    10501271        if 'T' in calcControls[hfx+'histType']:
    10511272            if 'P' in calcControls[hfx+'histType']:
    … …  
    10851306            if len(TwinLaw) > 1:
    10861307                refl.T[10] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2                  #FcT from primary twin element
    1087                 refl.T[8] = np.sum(TwinFr*np.sum(TwMask[np.newaxis,:,:]*fas,axis=0)**2,axis=-1)+   \
    1088                     np.sum(TwinFr*np.sum(TwMask[np.newaxis,:,:]*fbs,axis=0)**2,axis=-1)                 #Fc sum over twins
     1308                refl.T[8] = np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fas,axis=0)**2,axis=-1)+   \
     1309                    np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fbs,axis=0)**2,axis=-1)                 #Fc sum over twins
    10891310                refl.T[11] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f"
    10901311            else:
    … …  
    10971318def SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
    10981319    'Needs a doc string'
    1099     nxs = np.newaxis       
    11001320    phfx = pfx.split(':')[0]+hfx
    11011321    ast = np.sqrt(np.diag(G))
    … …  
    24422662                dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
    24432663            else:
    2444                 dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
     2664                dFdvDict = StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
    24452665#            print 'sf-derv time %.3fs'%(time.time()-time0)
    24462666            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
    … …  
    26892909    A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
    26902910    G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
     2911    TwinLaw = calcControls[phfx+'TwinLaw']
    26912912    refDict = Histogram['Data']
    26922913    if parmDict[phfx+'Scale'] < 0.:
    … …  
    26952916        dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
    26962917    else:
    2697         dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
     2918        if len(TwinLaw) > 1:
     2919            dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
     2920        else:
     2921            dFdvDict = StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
    26982922    ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
    26992923    dMdvh = np.zeros((len(varylist),len(refDict['RefList'])))
    … …  
    28493073            dMdvh = getPowderProfileDerv(parmDict,x[xB:xF],
    28503074                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
    2851             Wt = ma.sqrt(W[xB:xF])[np.newaxis,:]
    2852             Dy = dy[xB:xF][np.newaxis,:]
     3075            Wt = ma.sqrt(W[xB:xF])[nxs,:]
     3076            Dy = dy[xB:xF][nxs,:]
    28533077            dMdvh *= Wt
    28543078            if dlg:
Note: See TracChangeset for help on using the changeset viewer.