Changeset 2110


Ignore:
Timestamp:
Jan 3, 2016 12:40:04 PM (7 years ago)
Author:
vondreele
Message:

add menu item for global setting of wave vary parameters - TBD
split StructureFactorDeriv? over twins & incommensurate
do block refl for the nontwin/powder one

Location:
trunk
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIgrid.py

    r2109 r2110  
    6868    wxID_RELOADDRAWATOMS,wxID_ATOMSDISAGL,wxID_ATOMMOVE,wxID_MAKEMOLECULE,
    6969    wxID_ASSIGNATMS2RB,wxID_ATOMSPDISAGL, wxID_ISODISP,wxID_ADDHATOM,wxID_UPDATEHATOM,
    70 ] = [wx.NewId() for item in range(17)]
     70    wxID_WAVEVARY,
     71] = [wx.NewId() for item in range(18)]
    7172
    7273[ wxID_DRAWATOMSTYLE, wxID_DRAWATOMLABEL, wxID_DRAWATOMCOLOR, wxID_DRAWATOMRESETCOLOR,
    … …  
    12921293        self.PostfillDataMenu()
    12931294       
    1294         # Phase / Imcommensurate "waves" tab
     1295        # Phase / Imcommensurate "waves" tab 
    12951296        self.WavesData = wx.MenuBar()
    12961297        self.PrefillDataMenu(self.WavesData,helpType='Wave Data', helpLbl='Imcommensurate wave data')
    12971298        self.WavesData.Append(menu=wx.Menu(title=''),title='Select tab')
    1298         self.WavesDataCompute = wx.Menu(title='')
     1299        self.WavesDataEdit = wx.Menu(title='')
     1300        self.WavesData.Append(menu=self.WavesDataEdit, title='Edit')
     1301        self.WavesDataEdit.Append(id=wxID_WAVEVARY, kind=wx.ITEM_NORMAL,text='Global wave vary',
     1302            help='Global setting of wave vary flags')
    12991303        self.PostfillDataMenu()
    13001304                 

  • trunk/GSASIImath.py

    r2089 r2110  
    10821082    '''
    10831083    H: array ops X hklt proj to hkl
    1084     HP: array nRefBlk x ops X hklt proj to hkl
     1084    HP: array ops X hklt proj to hkl
    10851085    Hij: array 2pi^2[a*^2h^2 b*^2k^2 c*^2l^2 a*b*hk a*c*hl b*c*kl] of projected hklm to hkl space
    10861086    '''
    … …  
    11331133    dGdMzSt = np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(dHdXZt*np.cos(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
    11341134    dGdMzSx = np.sum((Fmod*HbH)[:,:,:,nxs]*(dHdXZx*np.cos(HdotXD)[:,:,:,nxs])*glWt[nxs,nxs,:,nxs],axis=-2)
     1135    dGdMzS = np.concatenate((np.sum(dGdMzSt,axis=-1),dGdMzSx),axis=-1)
     1136#    GSASIIpath.IPyBreak()
     1137    return [dGdMfC,dGdMfS],[dGdMxC,dGdMxS],[dGdMuC,dGdMuS],[dGdMzC,dGdMzS]
     1138   
     1139def ModulationDerv2(H,HP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt):
     1140    '''
     1141    H: array refBlk x ops X hklt proj to hkl
     1142    HP: array refBlk x ops X hklt proj to hkl
     1143    Hij: array 2pi^2[a*^2h^2 b*^2k^2 c*^2l^2 a*b*hk a*c*hl b*c*kl] of projected hklm to hkl space
     1144    '''
     1145   
     1146    Mf = [H.shape[0],]+list(waveShapes[0])    #=[ops,atoms,waves,2] (sin+cos frac mods)
     1147    dGdMfC = np.zeros(Mf)
     1148    dGdMfS = np.zeros(Mf)
     1149    Mx = [H.shape[0],]+list(waveShapes[1])   #=[ops,atoms,waves,6] (sin+cos pos mods)
     1150    dGdMxC = np.zeros(Mx)
     1151    dGdMxS = np.zeros(Mx)
     1152    Mu = [H.shape[0],]+list(waveShapes[2])    #=[ops,atoms,waves,12] (sin+cos Uij mods)
     1153    dGdMuC = np.zeros(Mu)
     1154    dGdMuS = np.zeros(Mu)
     1155   
     1156    D = twopi*H[:,:,3,nxs]*glTau[nxs,nxs,:]              #m*e*tau; refBlk x ops X ngl
     1157    HdotX = twopi*np.inner(HP,Xmod)        #refBlk x ops x atoms X ngl
     1158    HdotXD = HdotX+D[:,:,nxs,:]
     1159    if nWaves[2]:
     1160        Umod = np.swapaxes((UmodAB),2,4)      #atoms x waves x ngl x 3x3 (symmetric so I can do this!)
     1161        HuH = np.sum(HP[:,:,nxs,nxs,nxs]*np.inner(HP,Umod),axis=-1)    #refBlk x ops x atoms x waves x ngl
     1162        HuH = np.sum(HP[:,:,nxs,nxs,nxs]*np.inner(HP,Umod),axis=-1)    #refBlk x ops x atoms x waves x ngl
     1163        HbH = np.exp(-np.sum(HuH,axis=-2)) #refBlk x ops x atoms x ngl; sum waves - OK vs Modulation version
     1164        part1 = -np.exp(-HuH)*Fmod    #refBlk x ops x atoms x waves x ngl
     1165        dUdAu = Hij[:,:,nxs,nxs,nxs,:]*np.rollaxis(G2lat.UijtoU6(SCtauU[0]),0,4)[nxs,nxs,:,:,:,:]    #ops x atoms x waves x ngl x 6sinUij
     1166        dUdBu = Hij[:,:,nxs,nxs,nxs,:]*np.rollaxis(G2lat.UijtoU6(SCtauU[1]),0,4)[nxs,nxs,:,:,:,:]    #ops x atoms x waves x ngl x 6cosUij
     1167        dGdMuCa = np.sum(part1[:,:,:,:,:,nxs]*dUdAu*np.cos(HdotXD)[:,:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)   #ops x atoms x waves x 6uij; G-L sum
     1168        dGdMuCb = np.sum(part1[:,:,:,:,:,nxs]*dUdBu*np.cos(HdotXD)[:,:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)
     1169        dGdMuC = np.concatenate((dGdMuCa,dGdMuCb),axis=-1)   #ops x atoms x waves x 12uij
     1170        dGdMuSa = np.sum(part1[:,:,:,:,:,nxs]*dUdAu*np.sin(HdotXD)[:,:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)   #ops x atoms x waves x 6uij; G-L sum
     1171        dGdMuSb = np.sum(part1[:,:,:,:,:,nxs]*dUdBu*np.sin(HdotXD)[:,:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)
     1172        dGdMuS = np.concatenate((dGdMuSa,dGdMuSb),axis=-1)   #ops x atoms x waves x 12uij
     1173    else:
     1174        HbH = np.ones_like(HdotX)
     1175    dHdXA = twopi*HP[:,:,nxs,nxs,nxs,:]*np.swapaxes(SCtauX[0],-1,-2)[nxs,nxs,:,:,:,:]    #ops x atoms x sine waves x ngl x xyz
     1176    dHdXB = twopi*HP[:,:,nxs,nxs,nxs,:]*np.swapaxes(SCtauX[1],-1,-2)[nxs,nxs,:,:,:,:]    #ditto - cos waves
     1177# ops x atoms x waves x 2xyz - real part - good
     1178    dGdMxCa = -np.sum((Fmod*HbH)[:,:,:,nxs,:,nxs]*(dHdXA*np.sin(HdotXD)[:,:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)
     1179    dGdMxCb = -np.sum((Fmod*HbH)[:,:,:,nxs,:,nxs]*(dHdXB*np.sin(HdotXD)[:,:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)
     1180    dGdMxC = np.concatenate((dGdMxCa,dGdMxCb),axis=-1)
     1181# ops x atoms x waves x 2xyz - imag part - good
     1182    dGdMxSa = np.sum((Fmod*HbH)[:,:,:,nxs,:,nxs]*(dHdXA*np.cos(HdotXD)[:,:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)   
     1183    dGdMxSb = np.sum((Fmod*HbH)[:,:,:,nxs,:,nxs]*(dHdXB*np.cos(HdotXD)[:,:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)   
     1184    dGdMxS = np.concatenate((dGdMxSa,dGdMxSb),axis=-1)
     1185# ZigZag/Block waves - problems here?
     1186    dHdXZt = -twopi*HP[:,:,nxs,nxs,nxs,:]*np.swapaxes(SCtauX[2],-1,-2)[nxs,nxs,:,:,:,:]          #ops x atoms x ngl x 2(ZigZag/Block Tminmax)
     1187    dHdXZx = twopi*HP[:,:,nxs,nxs,:]*np.swapaxes(SCtauX[3],-1,-2)[nxs,nxs,:,:,:]          #ops x atoms x ngl x 3(ZigZag/Block XYZmax)
     1188    dGdMzCt = -np.sum((Fmod*HbH)[:,:,:,nxs,:,nxs]*(dHdXZt*np.sin(HdotXD)[:,:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)
     1189    dGdMzCx = -np.sum((Fmod*HbH)[:,:,:,:,nxs]*(dHdXZx*np.sin(HdotXD)[:,:,:,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
     1190    dGdMzC = np.concatenate((np.sum(dGdMzCt,axis=-1),dGdMzCx),axis=-1)
     1191    dGdMzSt = np.sum((Fmod*HbH)[:,:,:,nxs,:,nxs]*(dHdXZt*np.cos(HdotXD)[:,:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)
     1192    dGdMzSx = np.sum((Fmod*HbH)[:,:,:,:,nxs]*(dHdXZx*np.cos(HdotXD)[:,:,:,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
    11351193    dGdMzS = np.concatenate((np.sum(dGdMzSt,axis=-1),dGdMzSx),axis=-1)
    11361194#    GSASIIpath.IPyBreak()

  • trunk/GSASIIphsGUI.py

    r2109 r2110  
    25472547        G2frame.bottomSizer = ShowAtomInfo()
    25482548        mainSizer.Add(G2frame.bottomSizer)
    2549        
     2549        #wxID_WAVEVARY
    25502550        SetPhaseWindow(G2frame.dataFrame,G2frame.waveData,mainSizer,Scroll)
     2551   
     2552    def OnWaveVary(event):
     2553        print 'set vary flags for all waves - TBD'
    25512554
    25522555################################################################################
    … …  
    63976400        if data['General']['Type'] in ['modulated','magnetic']:
    63986401            FillSelectPageMenu(TabSelectionIdDict, G2frame.dataFrame.WavesData)
     6402            G2frame.dataFrame.Bind(wx.EVT_MENU, OnWaveVary, id=G2gd.wxID_WAVEVARY)           
    63996403        # Draw Options
    64006404        FillSelectPageMenu(TabSelectionIdDict, G2frame.dataFrame.DataDrawOptions)

  • trunk/GSASIIstrIO.py

    r2070 r2110  
    24072407            for i,name in enumerate(Snames):
    24082408                ptlbls += '%12s' % (name)
    2409                 ptstr += '%12.6f' % (hapData[4][i])
     2409                ptstr += '%12.1f' % (hapData[4][i])
    24102410                if mustrainSig[1][i]:
    24112411                    refine = True
    2412                     sigstr += '%12.6f' % (mustrainSig[1][i])
     2412                    sigstr += '%12.1f' % (mustrainSig[1][i])
    24132413                else:
    24142414                    sigstr += 12*' '

  • trunk/GSASIIstrMath.py

    r2097 r2110  
    632632    return np.array(waveTypes),FSSdata,XSSdata,USSdata,MSSdata
    633633   
    634 def GetSSTauM(SGOps,SSOps,pfx,calcControls,XData):
    635    
    636     Natoms = calcControls['Natoms'][pfx]
    637     maxSSwave = calcControls['maxSSwave'][pfx]
    638     Smult = np.zeros((Natoms,len(SGOps)))
    639     TauT = np.zeros((Natoms,len(SGOps)))
    640     for ix,xyz in enumerate(XData.T):
    641         for isym,(sop,ssop) in enumerate(zip(SGOps,SSOps)):
    642             sdet,ssdet,dtau,dT,tauT = G2spc.getTauT(0,sop,ssop,xyz)
    643             Smult[ix][isym] = sdet
    644             TauT[ix][isym] = tauT
    645     return Smult,TauT
    646    
     634#def GetSSTauM(SGOps,SSOps,pfx,calcControls,XData):
     635#   
     636#    Natoms = calcControls['Natoms'][pfx]
     637#    maxSSwave = calcControls['maxSSwave'][pfx]
     638#    Smult = np.zeros((Natoms,len(SGOps)))
     639#    TauT = np.zeros((Natoms,len(SGOps)))
     640#    for ix,xyz in enumerate(XData.T):
     641#        for isym,(sop,ssop) in enumerate(zip(SGOps,SSOps)):
     642#            sdet,ssdet,dtau,dT,tauT = G2spc.getTauT(0,sop,ssop,xyz)
     643#            Smult[ix][isym] = sdet
     644#            TauT[ix][isym] = tauT
     645#    return Smult,TauT
     646#   
    647647def StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
    648648    ''' Compute structure factors for all h,k,l for phase
    … …  
    742742        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
    743743        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/len(SGMT)
    744         if 'T' in calcControls[hfx+'histType']:
     744        if 'T' in calcControls[hfx+'histType']: #fa,fb are 2 X blkSize X nTwin X nOps x nAtoms
    745745            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
    746746            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr])
    … …  
    748748            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
    749749            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
    750         fas = np.sum(np.sum(fa,axis=-1),axis=-1)       #real sum over atoms & unique hkl
    751         fbs = np.sum(np.sum(fb,axis=-1),axis=-1)  #imag sum over atoms & uniq hkl
     750        fas = np.sum(np.sum(fa,axis=-1),axis=-1)  #real 2 x blkSize x nTwin; sum over atoms & uniq hkl
     751        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)  #imag
    752752        if SGData['SGInv']: #centrosymmetric; B=0
    753753            fbs[0] *= 0.
    754         if 'P' in calcControls[hfx+'histType']:
    755             refl.T[9] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0)
     754            fas[1] *= 0.
     755        if 'PWDR' in calcControls[hfx+'histType']:     #PWDR: F^2 = A[0]^2 + A[1]^2 + B[0]^2 + B[1]^2
     756            refl.T[9] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0)     
    756757            refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
    757         else:
     758        else:                                       #HKLF: F^2 = (A[0]+A[1])^2 + (B[0]+B[1])^2
    758759            if len(TwinLaw) > 1:
    759760                refl.T[9] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2   #FcT from primary twin element
    … …  
    761762                    np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fbs,axis=0)**2,axis=-1)                        #Fc sum over twins
    762763                refl.T[10] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f"
    763             else:
     764#                GSASIIpath.IPyBreak()
     765            else:   # checked correct!!
    764766                refl.T[9] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2
    765767                refl.T[7] = np.copy(refl.T[9])               
    766768                refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
    767 #        refl.T[10] = atan2d(np.sum(fbs,axis=0),np.sum(fas,axis=0)) #include f' & f"
     769#                refl.T[10] = atan2d(np.sum(fbs,axis=0),np.sum(fas,axis=0)) #include f' & f"
    768770        iBeg += blkSize
    769771#    print ' %d sf time %.4f\r'%(nRef,time.time()-time0)
    … …  
    772774    '''Compute structure factor derivatives on single reflections - keep as it works for twins
    773775    but is slower for powders/nontwins
     776    input:
     777   
     778    :param dict refDict: where
     779        'RefList' list where each ref = h,k,l,it,d,...
     780        'FF' dict of form factors - filled in below
     781    :param np.array G:      reciprocal metric tensor
     782    :param str hfx:    histogram id string
     783    :param str pfx:    phase id string
     784    :param dict SGData: space group info. dictionary output from SpcGroup
     785    :param dict calcControls:
     786    :param dict ParmDict:
     787   
     788    :returns: dict dFdvDict: dictionary of derivatives
    774789    '''
    775790    phfx = pfx.split(':')[0]+hfx
    … …  
    914929                dFdua[iref] = [2.*TwMask[it]*(SA[it]*dfadua[it][0]+SA[it]*dfadua[it][1]+SB[it]*dfbdua[it][0]+SB[it]*dfbdua[it][1]) for it in range(nTwin)]
    915930                dFdtw[iref] = np.sum(TwMask*fas,axis=0)**2+np.sum(TwMask*fbs,axis=0)**2
    916                
    917931            else:   #these are good for no twin single crystals
    918932                dFdfr[iref] = (2.*SA*(dfadfr[0]+dfadfr[1])+2.*SB*(dfbdfr[0]+dfbdfr[1]))*Mdata/len(Uniq)
    … …  
    964978def StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
    965979    '''Compute structure factor derivatives on blocks of reflections - for powders/nontwins only
    966     faster than StructureFactorDerv
     980    faster than StructureFactorDerv - correct for powders/nontwins!!
     981    input:
     982   
     983    :param dict refDict: where
     984        'RefList' list where each ref = h,k,l,it,d,...
     985        'FF' dict of form factors - filled in below
     986    :param np.array G:      reciprocal metric tensor
     987    :param str hfx:    histogram id string
     988    :param str pfx:    phase id string
     989    :param dict SGData: space group info. dictionary output from SpcGroup
     990    :param dict calcControls:
     991    :param dict parmDict:
     992   
     993    :returns: dict dFdvDict: dictionary of derivatives
    967994    '''
    968995    phfx = pfx.split(':')[0]+hfx
    … …  
    10331060        fotp = FPP*Tcorr       
    10341061        if 'T' in calcControls[hfx+'histType']:
    1035             fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
    1036             fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr])
     1062            fa = np.array([fot*cosp,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
     1063            fb = np.array([fot*sinp,np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr])
    10371064        else:
    1038             fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
    1039             fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
     1065            fa = np.array([fot*cosp,-Flack*FPP*sinp*Tcorr])
     1066            fb = np.array([fot*sinp,Flack*FPP*cosp*Tcorr])
    10401067#        GSASIIpath.IPyBreak()
    10411068        fas = np.sum(np.sum(fa,axis=-1),axis=-1)      #real sum over atoms & unique hkl array(2,refBlk,nTwins)
    … …  
    11081135   
    11091136def StructureFactorDervTw(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
    1110     '''Compute structure factor derivatives on blocks of reflections - for twins only
    1111     faster than StructureFactorDerv - not working yet
     1137    '''Compute structure factor derivatives on single reflections - for twins only
     1138    input:
     1139   
     1140    :param dict refDict: where
     1141        'RefList' list where each ref = h,k,l,it,d,...
     1142        'FF' dict of form factors - filled in below
     1143    :param np.array G:      reciprocal metric tensor
     1144    :param str hfx:    histogram id string
     1145    :param str pfx:    phase id string
     1146    :param dict SGData: space group info. dictionary output from SpcGroup
     1147    :param dict calcControls:
     1148    :param dict parmDict:
     1149   
     1150    :returns: dict dFdvDict: dictionary of derivatives
    11121151    '''
    11131152    phfx = pfx.split(':')[0]+hfx
    … …  
    11181157    FFtables = calcControls['FFtables']
    11191158    BLtables = calcControls['BLtables']
    1120     TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])
    11211159    TwDict = refDict.get('TwDict',{})           
    1122     if 'S' in calcControls[hfx+'histType']:
    1123         NTL = calcControls[phfx+'NTL']
    1124         NM = calcControls[phfx+'TwinNMN']+1
    1125         TwinLaw = calcControls[phfx+'TwinLaw']
    1126         TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
    1127         TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
     1160    NTL = calcControls[phfx+'NTL']
     1161    NM = calcControls[phfx+'TwinNMN']+1
     1162    TwinLaw = calcControls[phfx+'TwinLaw']
     1163    TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
     1164    TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
    11281165    nTwin = len(TwinLaw)       
    11291166    nRef = len(refDict['RefList'])
    … …  
    11441181    dFdua = np.zeros((nRef,nTwin,mSize,6))
    11451182    dFdbab = np.zeros((nRef,nTwin,2))
    1146     dFdfl = np.zeros((nRef,nTwin))
    11471183    dFdtw = np.zeros((nRef,nTwin))
    1148     Flack = 1.0
    1149     if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
    1150         Flack = 1.-2.*parmDict[phfx+'Flack']
    11511184    time0 = time.time()
    11521185    nref = len(refDict['RefList'])/100   
    1153 #reflection processing begins here - big arrays!
    1154     iBeg = 0
    1155     blkSize = 32       #no. of reflections in a block - optimized for speed
    1156     while iBeg < nRef:
    1157         iFin = min(iBeg+blkSize,nRef)
    1158         refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
    1159         H = refl.T[:3]
     1186    for iref,refl in enumerate(refDict['RefList']):
     1187        if 'T' in calcControls[hfx+'histType']:
     1188            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12])
     1189        H = np.array(refl[:3])
    11601190        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(3,nTwins) or (3)
    11611191        TwMask = np.any(H,axis=-1)
    1162         for ir in range(blkSize):
    1163             iref = ir+iBeg
    1164             if iref in TwDict:
    1165                 for i in TwDict[iref]:
    1166                     for n in range(NTL):
    1167                         H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
    1168         TwMask = np.any(H,axis=-1)
    1169         SQ = 1./(2.*refl.T[4])**2             # or (sin(theta)/lambda)**2
     1192        if iref in TwDict:
     1193            for i in TwDict[iref]:
     1194                for n in range(NTL):
     1195                    H[i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
     1196            TwMask = np.any(H,axis=-1)
     1197        SQ = 1./(2.*refl[4])**2             # or (sin(theta)/lambda)**2
    11701198        SQfactor = 8.0*SQ*np.pi**2
    1171         if 'T' in calcControls[hfx+'histType']:
    1172             if 'P' in calcControls[hfx+'histType']:
    1173                 FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
    1174             else:
    1175                 FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
    1176             FP = np.repeat(FP.T,len(SGT)*len(TwinLaw),axis=0)
    1177             FPP = np.repeat(FPP.T,len(SGT)*len(TwinLaw),axis=0)
    11781199        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
    1179         Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT)*len(TwinLaw))
     1200        Bab = parmDict[phfx+'BabA']*dBabdA
    11801201        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
    1181         FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT)*len(TwinLaw),axis=0)
     1202        FF = refDict['FF']['FF'][iref].T[Tindx].T
    11821203        Uniq = np.inner(H,SGMT)             # array(nSGOp,3) or (nTwin,nSGOp,3)
    11831204        Phi = np.inner(H,SGT)
    11841205        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
    1185         sinp = np.sin(phase)        #refBlk x nOps x nAtoms
     1206        sinp = np.sin(phase)
    11861207        cosp = np.cos(phase)
    11871208        occ = Mdata*Fdata/len(SGT)
    11881209        biso = -SQfactor*Uisodata[:,nxs]
    1189         Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*len(TwinLaw),axis=1).T
     1210        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*nTwin,axis=1)
    11901211        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
    1191         Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
    1192         Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/len(SGMT)
    11931212        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])
    1194         Hij = np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(-1,nTwin,len(SGT),6))
    1195         fot = np.reshape(((FF+FP).T-Bab).T,cosp.shape)*Tcorr
     1213        Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6)))
     1214        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
     1215        Tcorr = (np.reshape(Tiso,Tuij.shape)*Tuij).T*occ
     1216        fot = (FF+FP-Bab)*Tcorr
    11961217        fotp = FPP*Tcorr       
    1197         if 'T' in calcControls[hfx+'histType']:
    1198             fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
    1199             fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr])
    1200         else:
    1201             fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
    1202             fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
     1218        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-FPP*sinp*Tcorr])
     1219        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,FPP*cosp*Tcorr])
    12031220#        GSASIIpath.IPyBreak()
    1204         fas = np.sum(np.sum(fa,axis=-1),axis=-1)      #real sum over atoms & unique hkl array(2,refBlk,nTwins)
     1221        fas = np.sum(np.sum(fa,axis=-1),axis=-1)      #real sum over atoms & unique hkl array(2,nTwins)
    12051222        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)      #imag sum over atoms & uniq hkl
    1206         fax = np.array([-fot*sinp,-fotp*cosp])   #positions array(2,refBlk,ntwi,nEqv,nAtoms)
     1223        if SGData['SGInv']: #centrosymmetric; B=0
     1224            fbs[0] *= 0.
     1225            fas[1] *= 0.
     1226        fax = np.array([-fot*sinp,-fotp*cosp])   #positions array(2,ntwi,nEqv,nAtoms)
    12071227        fbx = np.array([fot*cosp,-fotp*sinp])
    12081228        #sum below is over Uniq
    1209         dfadfr = np.sum(fa/occ,axis=-2)        #array(2,refBlk,ntwin,nAtom) Fdata != 0 avoids /0. problem
    1210         dfadba = np.sum(-cosp*Tcorr,axis=-2)  #array(refBlk,nAtom)
    1211         dfadfr = np.swapaxes(dfadfr,0,1)
    1212         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)
    1213         dfadui = np.swapaxes(np.sum(-SQfactor[nxs,:,nxs,nxs,nxs]*fa,axis=-2),0,1) #array(Ops,refBlk,nTw,nAtoms)
    1214         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)
    1215         # array(nTwin,2,refBlk,nAtom,3) & array(nTwin,2,refBlk,nAtom,6)
     1229        dfadfr = np.sum(fa/occ,axis=-2)        #array(2,ntwin,nAtom) Fdata != 0 avoids /0. problem
     1230        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
     1231        dfadui = np.sum(-SQfactor*fa,axis=-2)
     1232        dfadx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
     1233        dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fa,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
     1234        # array(nTwin,2,nAtom,3) & array(nTwin,2,nAtom,6)
    12161235        if not SGData['SGInv']:
    12171236            dfbdfr = np.sum(fb/occ,axis=-2)        #Fdata != 0 avoids /0. problem
    1218             dfbdba = np.sum(-sinp*Tcorr,axis=-2)
    1219             dfbdfr = np.swapaxes(dfbdfr,0,1)
    1220             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)         
    1221             dfbdui = np.swapaxes(np.sum(-SQfactor[nxs,:,nxs,nxs,nxs]*fb,axis=-2),0,1)      #array(Ops,refBlk,nTw,nAtoms)
    1222             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)
     1237            dfadba /= 2.
     1238            dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)/2.
     1239            dfbdui = np.sum(-SQfactor*fb,axis=-2)
     1240            dfbdx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)])           
     1241            dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fb,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)])
    12231242        else:
    12241243            dfbdfr = np.zeros_like(dfadfr)
    … …  
    12271246            dfbdua = np.zeros_like(dfadua)
    12281247            dfbdba = np.zeros_like(dfadba)
    1229             dfadfl = 0.0
    1230             dfbdfl = 0.0
    1231         #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
    12321248        SA = fas[0]+fas[1]
    12331249        SB = fbs[0]+fbs[1]
    1234         dFdfr[iBeg:iFin] = np.swapaxes(np.array([2.*TwMask[:,it,nxs]*(SA[:,it,nxs]*(np.sum(dfadfr[:,:,it],axis=1))+ \
    1235             SB[:,it,nxs]*(np.sum(dfbdfr[:,:,it],axis=1)))*Mdata/len(SGMT) for it in range(nTwin)]),0,1)
    1236         dFdx[iBeg:iFin] = np.swapaxes(np.array([2.*TwMask[:,it,nxs,nxs]*(SA[:,it,nxs,nxs]*np.sum(dfadx[:,:,it],axis=1)+    \
    1237             SB[:,it,nxs,nxs]*np.sum(dfbdx[:,:,it],axis=1)) for it in range(nTwin)]),0,1)
    1238         dFdui[iBeg:iFin] = np.swapaxes(np.array([2.*TwMask[:,it,nxs]*(SA[:,it,nxs]*np.sum(dfadui[:,:,it],axis=1)+  \
    1239             SB[:,it,nxs]*np.sum(dfbdui[:,:,it],axis=1)) for it in range(nTwin)]),0,1)
    1240         dFdua[iBeg:iFin] = np.swapaxes(np.array([2.*TwMask[:,it,nxs,nxs]*(SA[:,it,nxs,nxs]*np.sum(dfadua[:,:,it],axis=1)+  \
    1241             SB[:,it,nxs,nxs]*np.sum(dfbdua[:,:,it],axis=1)) for it in range(nTwin)]),0,1)
    1242         dFdtw[iBeg:iFin] = np.sum(TwMask*fas,axis=0)**2+np.sum(TwMask*fbs,axis=0)**2
    1243                
    1244 #        dFdbab[iBeg:iFin] = 2.*(fas[0,nxs]*np.array([np.sum(dfadba.T*dBabdA,axis=0),np.sum(-dfadba.T*parmDict[phfx+'BabA']*SQfactor*dBabdA,axis=0)])+ \
    1245 #                            fbs[0,nxs]*np.array([np.sum(dfbdba.T*dBabdA,axis=0),np.sum(-dfbdba.T*parmDict[phfx+'BabA']*SQfactor*dBabdA,axis=0)])).T
     1250        dFdfr[iref] = [2.*TwMask[it]*(SA[it]*dfadfr[0,it]+SA[it]*dfadfr[1,it]+SB[it]*dfbdfr[0,it]+SB[it]*dfbdfr[1,it])*Mdata/len(Uniq[it]) for it in range(nTwin)]
     1251        dFdx[iref] = [2.*TwMask[it]*(SA[it]*dfadx[it,0]+SA[it]*dfadx[it,1]+SB[it]*dfbdx[it,0]+SB[it]*dfbdx[it,1]) for it in range(nTwin)]
     1252        dFdui[iref] = [2.*TwMask[it]*(SA[it]*dfadui[0,it]+SA[it]*dfadui[1,it]+SB[it]*dfbdui[0,it]+SB[it]*dfbdui[1,it]) for it in range(nTwin)]
     1253        dFdua[iref] = [2.*TwMask[it]*(SA[it]*dfadua[it,0]+SA[it]*dfadua[it,1]+SB[it]*dfbdua[it,0]+SB[it]*dfbdua[it,1]) for it in range(nTwin)]
     1254        if SGData['SGInv']: #centrosymmetric; B=0
     1255            dFdtw[iref] = np.sum(TwMask[nxs,:]*fas,axis=0)**2
     1256        else:               
     1257            dFdtw[iref] = np.sum(TwMask[nxs,:]*fas,axis=0)**2+np.sum(TwMask[nxs,:]*fbs,axis=0)**2
     1258        dFdbab[iref] = fas[0,:,nxs]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
     1259            fbs[0,:,nxs]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
    12461260#        GSASIIpath.IPyBreak()
    1247         iBeg += blkSize
    1248     print ' %d derivative time %.4f\r'%(nRef,time.time()-time0)
    1249         #loop over atoms - each dict entry is list of derivatives for all the reflections
     1261    print ' %d derivative time %.4f\r'%(len(refDict['RefList']),time.time()-time0)
     1262    #loop over atoms - each dict entry is list of derivatives for all the reflections
    12501263    for i in range(len(Mdata)):     #these all OK?
    12511264        dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,nxs],axis=0)
    … …  
    12601273        dFdvDict[pfx+'AU13:'+str(i)] = 2.*np.sum(dFdua.T[4][i]*TwinFr[:,nxs],axis=0)
    12611274        dFdvDict[pfx+'AU23:'+str(i)] = 2.*np.sum(dFdua.T[5][i]*TwinFr[:,nxs],axis=0)
    1262     dFdvDict[phfx+'Flack'] = 4.*dFdfl.T
     1275    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
     1276    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
     1277    for i in range(nTwin):
     1278        dFdvDict[phfx+'TwinFr:'+str(i)] = dFdtw.T[i]
     1279    return dFdvDict
     1280   
     1281def StructureFactorDervTw2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
     1282    '''Compute structure factor derivatives on blocks of reflections - for twins only
     1283    faster than StructureFactorDervTw
     1284    input:
     1285   
     1286    :param dict refDict: where
     1287        'RefList' list where each ref = h,k,l,it,d,...
     1288        'FF' dict of form factors - filled in below
     1289    :param np.array G:      reciprocal metric tensor
     1290    :param str hfx:    histogram id string
     1291    :param str pfx:    phase id string
     1292    :param dict SGData: space group info. dictionary output from SpcGroup
     1293    :param dict calcControls:
     1294    :param dict parmDict:
     1295   
     1296    :returns: dict dFdvDict: dictionary of derivatives
     1297    '''
     1298    phfx = pfx.split(':')[0]+hfx
     1299    ast = np.sqrt(np.diag(G))
     1300    Mast = twopisq*np.multiply.outer(ast,ast)
     1301    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
     1302    SGT = np.array([ops[1] for ops in SGData['SGOps']])
     1303    FFtables = calcControls['FFtables']
     1304    BLtables = calcControls['BLtables']
     1305    TwDict = refDict.get('TwDict',{})           
     1306    NTL = calcControls[phfx+'NTL']
     1307    NM = calcControls[phfx+'TwinNMN']+1
     1308    TwinLaw = calcControls[phfx+'TwinLaw']
     1309    TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
     1310    TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
     1311    nTwin = len(TwinLaw)       
     1312    nRef = len(refDict['RefList'])
     1313    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
     1314    mSize = len(Mdata)
     1315    FF = np.zeros(len(Tdata))
     1316    if 'NC' in calcControls[hfx+'histType']:
     1317        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
     1318    elif 'X' in calcControls[hfx+'histType']:
     1319        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
     1320        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
     1321    Uij = np.array(G2lat.U6toUij(Uijdata))
     1322    bij = Mast*Uij.T
     1323    dFdvDict = {}
     1324    dFdfr = np.zeros((nRef,nTwin,mSize))
     1325    dFdx = np.zeros((nRef,nTwin,mSize,3))
     1326    dFdui = np.zeros((nRef,nTwin,mSize))
     1327    dFdua = np.zeros((nRef,nTwin,mSize,6))
     1328    dFdbab = np.zeros((nRef,nTwin,2))
     1329    dFdtw = np.zeros((nRef,nTwin))
     1330    time0 = time.time()
     1331    nref = len(refDict['RefList'])/100   
     1332    for iref,refl in enumerate(refDict['RefList']):
     1333        if 'T' in calcControls[hfx+'histType']:
     1334            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12])
     1335        H = np.array(refl[:3])
     1336        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(3,nTwins) or (3)
     1337        TwMask = np.any(H,axis=-1)
     1338        if iref in TwDict:
     1339            for i in TwDict[iref]:
     1340                for n in range(NTL):
     1341                    H[i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
     1342            TwMask = np.any(H,axis=-1)
     1343        SQ = 1./(2.*refl[4])**2             # or (sin(theta)/lambda)**2
     1344        SQfactor = 8.0*SQ*np.pi**2
     1345        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
     1346        Bab = parmDict[phfx+'BabA']*dBabdA
     1347        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
     1348        FF = refDict['FF']['FF'][iref].T[Tindx].T
     1349        Uniq = np.inner(H,SGMT)             # array(nSGOp,3) or (nTwin,nSGOp,3)
     1350        Phi = np.inner(H,SGT)
     1351        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
     1352        sinp = np.sin(phase)
     1353        cosp = np.cos(phase)
     1354        occ = Mdata*Fdata/len(SGT)
     1355        biso = -SQfactor*Uisodata[:,nxs]
     1356        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*nTwin,axis=1)
     1357        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
     1358        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])
     1359        Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6)))
     1360        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
     1361        Tcorr = (np.reshape(Tiso,Tuij.shape)*Tuij).T*occ
     1362        fot = (FF+FP-Bab)*Tcorr
     1363        fotp = FPP*Tcorr       
     1364        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-FPP*sinp*Tcorr])
     1365        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,FPP*cosp*Tcorr])
     1366#        GSASIIpath.IPyBreak()
     1367        fas = np.sum(np.sum(fa,axis=-1),axis=-1)      #real sum over atoms & unique hkl array(2,nTwins)
     1368        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)      #imag sum over atoms & uniq hkl
     1369        if SGData['SGInv']: #centrosymmetric; B=0
     1370            fbs[0] *= 0.
     1371            fas[1] *= 0.
     1372        fax = np.array([-fot*sinp,-fotp*cosp])   #positions array(2,ntwi,nEqv,nAtoms)
     1373        fbx = np.array([fot*cosp,-fotp*sinp])
     1374        #sum below is over Uniq
     1375        dfadfr = np.sum(fa/occ,axis=-2)        #array(2,ntwin,nAtom) Fdata != 0 avoids /0. problem
     1376        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
     1377        dfadui = np.sum(-SQfactor*fa,axis=-2)
     1378        dfadx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
     1379        dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fa,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
     1380        # array(nTwin,2,nAtom,3) & array(nTwin,2,nAtom,6)
     1381        if not SGData['SGInv']:
     1382            dfbdfr = np.sum(fb/occ,axis=-2)        #Fdata != 0 avoids /0. problem
     1383            dfadba /= 2.
     1384            dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)/2.
     1385            dfbdui = np.sum(-SQfactor*fb,axis=-2)
     1386            dfbdx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)])           
     1387            dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fb,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)])
     1388        else:
     1389            dfbdfr = np.zeros_like(dfadfr)
     1390            dfbdx = np.zeros_like(dfadx)
     1391            dfbdui = np.zeros_like(dfadui)
     1392            dfbdua = np.zeros_like(dfadua)
     1393            dfbdba = np.zeros_like(dfadba)
     1394        SA = fas[0]+fas[1]
     1395        SB = fbs[0]+fbs[1]
     1396        dFdfr[iref] = [2.*TwMask[it]*(SA[it]*dfadfr[0,it]+SA[it]*dfadfr[1,it]+SB[it]*dfbdfr[0,it]+SB[it]*dfbdfr[1,it])*Mdata/len(Uniq[it]) for it in range(nTwin)]
     1397        dFdx[iref] = [2.*TwMask[it]*(SA[it]*dfadx[it,0]+SA[it]*dfadx[it,1]+SB[it]*dfbdx[it,0]+SB[it]*dfbdx[it,1]) for it in range(nTwin)]
     1398        dFdui[iref] = [2.*TwMask[it]*(SA[it]*dfadui[0,it]+SA[it]*dfadui[1,it]+SB[it]*dfbdui[0,it]+SB[it]*dfbdui[1,it]) for it in range(nTwin)]
     1399        dFdua[iref] = [2.*TwMask[it]*(SA[it]*dfadua[it,0]+SA[it]*dfadua[it,1]+SB[it]*dfbdua[it,0]+SB[it]*dfbdua[it,1]) for it in range(nTwin)]
     1400        if SGData['SGInv']: #centrosymmetric; B=0
     1401            dFdtw[iref] = np.sum(TwMask[nxs,:]*fas,axis=0)**2
     1402        else:               
     1403            dFdtw[iref] = np.sum(TwMask[nxs,:]*fas,axis=0)**2+np.sum(TwMask[nxs,:]*fbs,axis=0)**2
     1404        dFdbab[iref] = fas[0,:,nxs]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
     1405            fbs[0,:,nxs]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
     1406#        GSASIIpath.IPyBreak()
     1407    print ' %d derivative time %.4f\r'%(len(refDict['RefList']),time.time()-time0)
     1408    #loop over atoms - each dict entry is list of derivatives for all the reflections
     1409    for i in range(len(Mdata)):     #these all OK?
     1410        dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,nxs],axis=0)
     1411        dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,nxs],axis=0)
     1412        dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,nxs],axis=0)
     1413        dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,nxs],axis=0)
     1414        dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,nxs],axis=0)
     1415        dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,nxs],axis=0)
     1416        dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,nxs],axis=0)
     1417        dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,nxs],axis=0)
     1418        dFdvDict[pfx+'AU12:'+str(i)] = 2.*np.sum(dFdua.T[3][i]*TwinFr[:,nxs],axis=0)
     1419        dFdvDict[pfx+'AU13:'+str(i)] = 2.*np.sum(dFdua.T[4][i]*TwinFr[:,nxs],axis=0)
     1420        dFdvDict[pfx+'AU23:'+str(i)] = 2.*np.sum(dFdua.T[5][i]*TwinFr[:,nxs],axis=0)
    12631421    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
    12641422    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
    … …  
    13171475    Uij = np.array(G2lat.U6toUij(Uijdata)).T
    13181476    bij = Mast*Uij
    1319     blkSize = 100       #no. of reflections in a block
     1477    blkSize = 32       #no. of reflections in a block
    13201478    nRef = refDict['RefList'].shape[0]
    13211479    if not len(refDict['FF']):
    … …  
    13401498        H = refl.T[:4]                          #array(blkSize,4)
    13411499        HP = H[:3]+modQ[:,nxs]*H[3:]            #projected hklm to hkl
    1342 #        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(blkSize,nTwins,4) or (blkSize,4)
    1343 #        TwMask = np.any(H,axis=-1)
     1500#        HP = np.squeeze(np.inner(HP.T,TwinLaw))   #maybe array(blkSize,nTwins,4) or (blkSize,4)
     1501#        TwMask = np.any(HP,axis=-1)
    13441502#        if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]
    13451503#            for ir in range(blkSize):
    … …  
    13481506#                    for i in TwDict[iref]:
    13491507#                        for n in range(NTL):
    1350 #                            H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
    1351 #            TwMask = np.any(H,axis=-1)
     1508#                            HP[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
     1509#            TwMask = np.any(HP,axis=-1)
    13521510        SQ = 1./(2.*refl.T[5])**2               #array(blkSize)
    13531511        SQfactor = 4.0*SQ*twopisq               #ditto prev.
    … …  
    14071565
    14081566def SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
    1409     'Needs a doc string - no twins'
     1567    '''
     1568    Compute super structure factor derivatives for all h,k,l,m for phase - no twins
     1569    input:
     1570   
     1571    :param dict refDict: where
     1572        'RefList' list where each ref = h,k,l,m,it,d,...
     1573        'FF' dict of form factors - filled in below
     1574    :param int im: = 1 (could be eliminated)
     1575    :param np.array G:      reciprocal metric tensor
     1576    :param str hfx:    histogram id string
     1577    :param str pfx:    phase id string
     1578    :param dict SGData: space group info. dictionary output from SpcGroup
     1579    :param dict SSGData: super space group info.
     1580    :param dict calcControls:
     1581    :param dict ParmDict:
     1582   
     1583    :returns: dict dFdvDict: dictionary of derivatives
     1584    '''
    14101585    phfx = pfx.split(':')[0]+hfx
    14111586    ast = np.sqrt(np.diag(G))
    … …  
    15681743        if not iref%100 :
    15691744            print ' %d derivative time %.4f\r'%(iref,time.time()-time0),
     1745    for i in range(len(Mdata)):     #loop over atoms
     1746        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
     1747        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
     1748        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
     1749        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
     1750        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
     1751        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
     1752        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
     1753        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
     1754        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
     1755        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
     1756        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
     1757        for j in range(FSSdata.shape[1]):        #loop over waves Fzero & Fwid?
     1758            dFdvDict[pfx+'Fsin:'+str(i)+':'+str(j)] = dFdGf.T[0][j][i]
     1759            dFdvDict[pfx+'Fcos:'+str(i)+':'+str(j)] = dFdGf.T[1][j][i]
     1760        nx = 0
     1761        if waveTypes[i] in ['Block','ZigZag']:
     1762            nx = 1
     1763            dFdvDict[pfx+'Tmin:'+str(i)+':0'] = dFdGz.T[0][i]   #ZigZag/Block waves (if any)
     1764            dFdvDict[pfx+'Tmax:'+str(i)+':0'] = dFdGz.T[1][i]
     1765            dFdvDict[pfx+'Xmax:'+str(i)+':0'] = dFdGz.T[2][i]
     1766            dFdvDict[pfx+'Ymax:'+str(i)+':0'] = dFdGz.T[3][i]
     1767            dFdvDict[pfx+'Zmax:'+str(i)+':0'] = dFdGz.T[4][i]
     1768        for j in range(XSSdata.shape[1]-nx):       #loop over waves
     1769            dFdvDict[pfx+'Xsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[0][j][i]
     1770            dFdvDict[pfx+'Ysin:'+str(i)+':'+str(j+nx)] = dFdGx.T[1][j][i]
     1771            dFdvDict[pfx+'Zsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[2][j][i]
     1772            dFdvDict[pfx+'Xcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[3][j][i]
     1773            dFdvDict[pfx+'Ycos:'+str(i)+':'+str(j+nx)] = dFdGx.T[4][j][i]
     1774            dFdvDict[pfx+'Zcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[5][j][i]
     1775        for j in range(USSdata.shape[1]):       #loop over waves
     1776            dFdvDict[pfx+'U11sin:'+str(i)+':'+str(j)] = dFdGu.T[0][j][i]
     1777            dFdvDict[pfx+'U22sin:'+str(i)+':'+str(j)] = dFdGu.T[1][j][i]
     1778            dFdvDict[pfx+'U33sin:'+str(i)+':'+str(j)] = dFdGu.T[2][j][i]
     1779            dFdvDict[pfx+'U12sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[3][j][i]
     1780            dFdvDict[pfx+'U13sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[4][j][i]
     1781            dFdvDict[pfx+'U23sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[5][j][i]
     1782            dFdvDict[pfx+'U11cos:'+str(i)+':'+str(j)] = dFdGu.T[6][j][i]
     1783            dFdvDict[pfx+'U22cos:'+str(i)+':'+str(j)] = dFdGu.T[7][j][i]
     1784            dFdvDict[pfx+'U33cos:'+str(i)+':'+str(j)] = dFdGu.T[8][j][i]
     1785            dFdvDict[pfx+'U12cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[9][j][i]
     1786            dFdvDict[pfx+'U13cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[10][j][i]
     1787            dFdvDict[pfx+'U23cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[11][j][i]
     1788           
     1789#        GSASIIpath.IPyBreak()
     1790    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
     1791    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
     1792    return dFdvDict
     1793
     1794def SStructureFactorDerv2(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
     1795    'Needs a doc string - no twins'
     1796    phfx = pfx.split(':')[0]+hfx
     1797    ast = np.sqrt(np.diag(G))
     1798    Mast = twopisq*np.multiply.outer(ast,ast)
     1799    SGInv = SGData['SGInv']
     1800    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
     1801    SGT = np.array([ops[1] for ops in SGData['SGOps']])
     1802    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
     1803    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
     1804    FFtables = calcControls['FFtables']
     1805    BLtables = calcControls['BLtables']
     1806    nRef = len(refDict['RefList'])
     1807    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
     1808    mSize = len(Mdata)  #no. atoms
     1809    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
     1810    ngl,nWaves,Fmod,Xmod,Umod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,Mast)
     1811    waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,Mast)
     1812    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
     1813    FF = np.zeros(len(Tdata))
     1814    if 'NC' in calcControls[hfx+'histType']:
     1815        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
     1816    elif 'X' in calcControls[hfx+'histType']:
     1817        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
     1818        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
     1819    Uij = np.array(G2lat.U6toUij(Uijdata)).T
     1820    bij = Mast*Uij
     1821    if not len(refDict['FF']):
     1822        if 'N' in calcControls[hfx+'histType']:
     1823            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
     1824        else:
     1825            dat = G2el.getFFvalues(FFtables,0.)       
     1826        refDict['FF']['El'] = dat.keys()
     1827        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))
     1828    dFdvDict = {}
     1829    dFdfr = np.zeros((nRef,mSize))
     1830    dFdx = np.zeros((nRef,mSize,3))
     1831    dFdui = np.zeros((nRef,mSize))
     1832    dFdua = np.zeros((nRef,mSize,6))
     1833    dFdbab = np.zeros((nRef,2))
     1834    dFdfl = np.zeros((nRef))
     1835    dFdtw = np.zeros((nRef))
     1836    dFdGf = np.zeros((nRef,mSize,FSSdata.shape[1],2))
     1837    dFdGx = np.zeros((nRef,mSize,XSSdata.shape[1],6))
     1838    dFdGz = np.zeros((nRef,mSize,5))
     1839    dFdGu = np.zeros((nRef,mSize,USSdata.shape[1],12))
     1840    Flack = 1.0
     1841    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
     1842        Flack = 1.-2.*parmDict[phfx+'Flack']
     1843    time0 = time.time()
     1844    iBeg = 0
     1845    blkSize = 4       #no. of reflections in a block - optimized for speed
     1846    while iBeg < nRef:
     1847        iFin = min(iBeg+blkSize,nRef)
     1848        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
     1849        H = refl.T[:4]
     1850        HP = H[:3].T+modQ*H.T[:,3:]            #projected hklm to hkl
     1851        SQ = 1./(2.*refl.T[4+im])**2             # or (sin(theta)/lambda)**2
     1852        SQfactor = 8.0*SQ*np.pi**2
     1853        if 'T' in calcControls[hfx+'histType']:
     1854            if 'P' in calcControls[hfx+'histType']:
     1855                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[15])
     1856            else:
     1857                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[13])
     1858            FP = np.repeat(FP.T,len(SGT),axis=0)
     1859            FPP = np.repeat(FPP.T,len(SGT),axis=0)
     1860        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
     1861        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT))
     1862        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
     1863        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT),axis=0)
     1864        Uniq = np.inner(H.T,SSGMT)
     1865        Phi = np.inner(H.T,SSGT)
     1866        UniqP = np.inner(HP,SGMT)
     1867        if SGInv:   #if centro - expand HKL sets
     1868            Uniq = np.hstack((Uniq,-Uniq))
     1869            Phi = np.hstack((Phi,-Phi))
     1870            UniqP = np.hstack((UniqP,-UniqP))
     1871            FF = np.vstack((FF,FF))
     1872            Bab = np.concatenate((Bab,Bab))
     1873        phase = twopi*(np.inner(Uniq[:,:,:3],(dXdata+Xdata).T)+Phi[:,:,nxs])
     1874        sinp = np.sin(phase)
     1875        cosp = np.cos(phase)
     1876        occ = Mdata*Fdata/Uniq.shape[1]
     1877        biso = -SQfactor*Uisodata[:,nxs]
     1878        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[1],axis=1).T    #ops x atoms
     1879        HbH = -np.sum(UniqP[:,:,nxs,:3]*np.inner(UniqP[:,:,:3],bij),axis=-1)  #ops x atoms
     1880        Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in np.reshape(UniqP,(-1,3))]) #atoms x 3x3
     1881        Hij = np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(iFin-iBeg,-1,6))                     #atoms x 6
     1882        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)     #ops x atoms
     1883#        GSASIIpath.IPyBreak()
     1884        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[0]  #ops x atoms
     1885        fot = np.reshape(FF+FP[nxs,:]-Bab[:,nxs],cosp.shape)*Tcorr     #ops x atoms
     1886        fotp = FPP*Tcorr            #ops x atoms
     1887        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms
     1888        dGdf,dGdx,dGdu,dGdz = G2mth.ModulationDerv2(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt)
     1889        # GfpuA is 2 x ops x atoms
     1890        # derivs are: ops x atoms x waves x 2,6,12, or 5 parms as [real,imag] parts
     1891        fa = np.array([fot*cosp,-Flack*FPP*sinp*Tcorr]) # array(2,nEqv,nAtoms)
     1892        fb = np.array([fot*sinp,Flack*FPP*cosp*Tcorr])  #or array(2,nEqv,nAtoms)
     1893        fag = fa*GfpuA[0]-fb*GfpuA[1]
     1894        fbg = fb*GfpuA[0]+fa*GfpuA[1]
     1895       
     1896        fas = np.sum(np.sum(fag,axis=-1),axis=-1)     # 2 x refBlk
     1897        fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
     1898        fax = np.array([-fot*sinp,-fotp*cosp])   #positions; 2 x ops x atoms
     1899        fbx = np.array([fot*cosp,-fotp*sinp])
     1900        fax = fax*GfpuA[0]-fbx*GfpuA[1]
     1901        fbx = fbx*GfpuA[0]+fax*GfpuA[1]
     1902        #sum below is over Uniq
     1903        dfadfr = np.sum(fag/occ,axis=-2)        #Fdata != 0 ever avoids /0. problem
     1904        dfbdfr = np.sum(fbg/occ,axis=-2)        #Fdata != 0 avoids /0. problem
     1905        dfadba = np.sum(-cosp*Tcorr,axis=-2)
     1906        dfbdba = np.sum(-sinp*Tcorr,axis=-2)
     1907        dfadui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fag,axis=-2)
     1908        dfbdui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fbg,axis=-2)
     1909        dfadx = np.sum(twopi*Uniq[nxs,:,:,nxs,:3]*fax[:,:,:,:,nxs],axis=-3)  #2 refBlk x x nAtom x 3xyz; sum nOps
     1910        dfbdx = np.sum(twopi*Uniq[nxs,:,:,nxs,:3]*fbx[:,:,:,:,nxs],axis=-3)  #2 refBlk x x nAtom x 3xyz; sum nOps
     1911        dfadua = np.sum(-Hij[nxs,:,:,nxs,:]*fag[:,:,:,:,nxs],axis=-3)             #2 x nAtom x 6Uij; sum nOps
     1912        dfbdua = np.sum(-Hij[nxs,:,:,nxs,:]*fbg[:,:,:,:,nxs],axis=-3)             #2 x nAtom x 6Uij; sum nOps
     1913        # array(2,nAtom,nWave,2) & array(2,nAtom,nWave,6) & array(2,nAtom,nWave,12); sum on nOps
     1914        dfadGf = np.sum(fa[:,:,:,:,nxs,nxs]*dGdf[0][nxs,:,nxs,:,:,:]-fb[:,:,:,:,nxs,nxs]*dGdf[1][nxs,:,nxs,:,:,:],axis=2)
     1915        dfbdGf = np.sum(fb[:,:,:,:,nxs,nxs]*dGdf[0][nxs,:,nxs,:,:,:]+fa[:,:,:,:,nxs,nxs]*dGdf[1][nxs,:,nxs,:,:,:],axis=2)
     1916        dfadGx = np.sum(fa[:,:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:,:]-fb[:,:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:,:],axis=2)
     1917        dfbdGx = np.sum(fb[:,:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:,:]+fa[:,:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:,:],axis=2)
     1918        dfadGz = np.sum(fa[:,:,:,:,nxs]*dGdz[0][nxs,:,:,:,:]-fb[:,:,:,:,nxs]*dGdz[1][nxs,:,:,:,:],axis=2)
     1919        dfbdGz = np.sum(fb[:,:,:,:,nxs]*dGdz[0][nxs,:,:,:,:]+fa[:,:,:,:,nxs]*dGdz[1][nxs,:,:,:,:],axis=2)
     1920        dfadGu = np.sum(fa[:,:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]-fb[:,:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=2)
     1921        dfbdGu = np.sum(fb[:,:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]+fa[:,:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=2)   
     1922        if not SGData['SGInv']:   #Flack derivative
     1923            dfadfl = np.sum(np.sum(-FPP*Tcorr*sinp,axis=-1),axis=-1)
     1924            dfbdfl = np.sum(np.sum(FPP*Tcorr*cosp,axis=-1),axis=-1)
     1925        else:
     1926            dfadfl = 1.0
     1927            dfbdfl = 1.0
     1928        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
     1929        SA = fas[0]+fas[1]      #float = A+A' (might be array[nTwin])
     1930        SB = fbs[0]+fbs[1]      #float = B+B' (might be array[nTwin])
     1931        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro?
     1932            dFdfl[iBeg:iFin] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
     1933            dFdfr[iBeg:iFin] = 2.*(fas[0,:,nxs]*dfadfr[0]+fas[1,:,nxs]*dfadfr[1])*Mdata/len(Uniq)+   \
     1934                2.*(fbs[0,:,nxs]*dfbdfr[0]-fbs[1,:,nxs]*dfbdfr[1])*Mdata/len(Uniq)
     1935            dFdx[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadx[0]+fas[1,:,nxs,nxs]*dfadx[1])+  \
     1936                2.*(fbs[0,:,nxs,nxs]*dfbdx[0]+fbs[1,:,nxs,nxs]*dfbdx[1])
     1937            dFdui[iBeg:iFin] = 2.*(fas[0,:,nxs]*dfadui[0]+fas[1,:,nxs]*dfadui[1])+   \
     1938                2.*(fbs[0,:,nxs]*dfbdui[0]-fbs[1,:,nxs]*dfbdui[1])
     1939            dFdua[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadua[0]+fas[1,:,nxs,nxs]*dfadua[1])+   \
     1940                2.*(fbs[0,:,nxs,nxs]*dfbdua[0]+fbs[1,:,nxs,nxs]*dfbdua[1])
     1941            dFdGf[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs,nxs]*dfadGf[0]+fas[1,:,nxs,nxs,nxs]*dfadGf[1])+  \
     1942                2.*(fbs[0,:,nxs,nxs,nxs]*dfbdGf[0]+fbs[1,:,nxs,nxs,nxs]*dfbdGf[1])
     1943            dFdGx[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs,nxs]*dfadGx[0]+fas[1,:,nxs,nxs,nxs]*dfadGx[1])+  \
     1944                2.*(fbs[0,:,nxs,nxs,nxs]*dfbdGx[0]-fbs[1,:,nxs,nxs,nxs]*dfbdGx[1])
     1945            dFdGz[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadGz[0]+fas[1,:,nxs,nxs]*dfadGz[1])+  \
     1946                2.*(fbs[0,:,nxs,nxs]*dfbdGz[0]+fbs[1,:,nxs,nxs]*dfbdGz[1])
     1947            dFdGu[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs,nxs]*dfadGu[0]+fas[1,:,nxs,nxs,nxs]*dfadGu[1])+  \
     1948                2.*(fbs[0,:,nxs,nxs,nxs]*dfbdGu[0]+fbs[1,:,nxs,nxs,nxs]*dfbdGu[1])
     1949        else:                       #OK, I think
     1950            dFdfr[iBeg:iFin] = 2.*(SA[:,nxs]*(dfadfr[0]+dfadfr[1])+SB[:,nxs]*(dfbdfr[0]+dfbdfr[1]))*Mdata/len(Uniq) #array(nRef,nAtom)
     1951            dFdx[iBeg:iFin] = 2.*(SA[:,nxs,nxs]*(dfadx[0]+dfadx[1])+SB[:,nxs,nxs]*(dfbdx[0]+dfbdx[1]))    #array(nRef,nAtom,3)
     1952            dFdui[iBeg:iFin] = 2.*(SA[:,nxs]*(dfadui[0]+dfadui[1])+SB[:,nxs]*(dfbdui[0]+dfbdui[1]))   #array(nRef,nAtom)
     1953            dFdua[iBeg:iFin] = 2.*(SA[:,nxs,nxs]*(dfadua[0]+dfadua[1])+SB[:,nxs,nxs]*(dfbdua[0]+dfbdua[1]))    #array(nRef,nAtom,6)
     1954            dFdfl[iBeg:iFin] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
     1955                           
     1956            dFdGf[iBeg:iFin] = 2.*(SA[:,nxs,nxs,nxs]*dfadGf[0]+SB[:,nxs,nxs,nxs]*dfbdGf[1])      #array(nRef,natom,nwave,2)
     1957            dFdGx[iBeg:iFin] = 2.*(SA[:,nxs,nxs,nxs]*dfadGx[0]+SB[:,nxs,nxs,nxs]*dfbdGx[1])      #array(nRef,natom,nwave,6)
     1958            dFdGz[iBeg:iFin] = 2.*(SA[:,nxs,nxs]*dfadGz[0]+SB[:,nxs,nxs]*dfbdGz[1])      #array(nRef,natom,5)
     1959            dFdGu[iBeg:iFin] = 2.*(SA[:,nxs,nxs,nxs]*dfadGu[0]+SB[:,nxs,nxs,nxs]*dfbdGu[1])      #array(nRef,natom,nwave,12)
     1960#            GSASIIpath.IPyBreak()
     1961#        dFdbab[iBeg:iFin] = 2.*fas[0,:,nxs]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
     1962#            2.*fbs[0,:,nxs]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
     1963        #loop over atoms - each dict entry is list of derivatives for all the reflections
     1964        print ' %d derivative time %.4f\r'%(iBeg,time.time()-time0),
     1965        iBeg += blkSize
    15701966    for i in range(len(Mdata)):     #loop over atoms
    15711967        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
    … …  
    31623558    else:
    31633559        if len(TwinLaw) > 1:
    3164             dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
    3165         else:
     3560            dFdvDict = StructureFactorDervTw2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
     3561        else:   #correct!!
    31663562            dFdvDict = StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
    31673563    ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
    … …  
    33623758def errRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg=None):       
    33633759    '''Computes the point-by-point discrepancies between every data point in every histogram
    3364     and the observed value
     3760    and the observed value. Used in the Jacobian, Hessian & numeric least-squares to compute function
    33653761   
    33663762    :returns: an np array of differences between observed and computed diffraction values.
Note: See TracChangeset for help on using the changeset viewer.