Changeset 2111


Ignore:
Timestamp:
Jan 4, 2016 1:36:49 PM (6 years ago)
Author:
vondreele
Message:

fix to mod vec edit & search
comment out unused StructureFactorDerv? versions

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIindex.py

    r2090 r2111  
    313313    for v,r in zip(ssopt['ModVec'],Vref):
    314314        if r:
    315             ranges += [slice(0.05,1.98,.05),]
     315            ranges += [slice(0.02,0.98,.02),]
    316316            values += [v,]
    317317    dmin = getDmin(peaks)-0.005
  • trunk/GSASIIphsGUI.py

    r2110 r2111  
    16911691                if nextName:
    16921692                    nextneigh = G2mth.FindNeighbors(data,nextName,AtNames,notName=neigh[0])
    1693                     neigh[1][1].append(nextneigh[1][1][0])
     1693                    if nextneigh[0]:
     1694                        neigh[1][1].append(nextneigh[1][1][0])
    16941695                neigh[2] = max(0,nH)  #set expected no. H's needed
    16951696                if len(neigh[1][0]):
  • trunk/GSASIIpwdGUI.py

    r2109 r2111  
    24702470        Id = Indx[ObjId]
    24712471        try:
    2472             value = min(1.0,max(-1.0,float(Obj.GetValue())))
     2472            value = min(0.98,max(-0.98,float(Obj.GetValue())))
    24732473        except ValueError:
    24742474            value = ssopt['ModVec'][Id]
     
    24812481        ObjId = Obj.GetId()
    24822482        Id,valObj = Indx[ObjId]
    2483         move = Obj.GetValue()*0.001
     2483        move = Obj.GetValue()*0.01
    24842484        Obj.SetValue(0)
    2485         value = min(1.0,max(-1.0,float(valObj.GetValue())+move))
     2485        value = min(0.98,max(-0.98,float(valObj.GetValue())+move))
    24862486        valObj.SetValue('%.4f'%(value))
    24872487        ssopt['ModVec'][Id] = value
  • trunk/GSASIIstrMath.py

    r2110 r2111  
    771771#    print ' %d sf time %.4f\r'%(nRef,time.time()-time0)
    772772   
    773 def StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
    774     '''Compute structure factor derivatives on single reflections - keep as it works for twins
    775     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
    789     '''
    790     phfx = pfx.split(':')[0]+hfx
    791     ast = np.sqrt(np.diag(G))
    792     Mast = twopisq*np.multiply.outer(ast,ast)
    793     SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
    794     SGT = np.array([ops[1] for ops in SGData['SGOps']])
    795     FFtables = calcControls['FFtables']
    796     BLtables = calcControls['BLtables']
    797     TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])
    798     TwDict = refDict.get('TwDict',{})           
    799     if 'S' in calcControls[hfx+'histType']:
    800         NTL = calcControls[phfx+'NTL']
    801         NM = calcControls[phfx+'TwinNMN']+1
    802         TwinLaw = calcControls[phfx+'TwinLaw']
    803         TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
    804         TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
    805     nTwin = len(TwinLaw)       
    806     nRef = len(refDict['RefList'])
    807     Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
    808     mSize = len(Mdata)
    809     FF = np.zeros(len(Tdata))
    810     if 'NC' in calcControls[hfx+'histType']:
    811         FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
    812     elif 'X' in calcControls[hfx+'histType']:
    813         FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
    814         FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
    815     Uij = np.array(G2lat.U6toUij(Uijdata))
    816     bij = Mast*Uij.T
    817     dFdvDict = {}
    818     if nTwin > 1:
    819         dFdfr = np.zeros((nRef,nTwin,mSize))
    820         dFdx = np.zeros((nRef,nTwin,mSize,3))
    821         dFdui = np.zeros((nRef,nTwin,mSize))
    822         dFdua = np.zeros((nRef,nTwin,mSize,6))
    823         dFdbab = np.zeros((nRef,nTwin,2))
    824         dFdfl = np.zeros((nRef,nTwin))
    825         dFdtw = np.zeros((nRef,nTwin))
    826     else:
    827         dFdfr = np.zeros((nRef,mSize))
    828         dFdx = np.zeros((nRef,mSize,3))
    829         dFdui = np.zeros((nRef,mSize))
    830         dFdua = np.zeros((nRef,mSize,6))
    831         dFdbab = np.zeros((nRef,2))
    832         dFdfl = np.zeros((nRef))
    833         dFdtw = np.zeros((nRef))
    834     Flack = 1.0
    835     if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
    836         Flack = 1.-2.*parmDict[phfx+'Flack']
    837     time0 = time.time()
    838     nref = len(refDict['RefList'])/100   
    839     for iref,refl in enumerate(refDict['RefList']):
    840         if 'T' in calcControls[hfx+'histType']:
    841             FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12])
    842         H = np.array(refl[:3])
    843         H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(3,nTwins) or (3)
    844         TwMask = np.any(H,axis=-1)
    845         if TwinLaw.shape[0] > 1 and TwDict:
    846             if iref in TwDict:
    847                 for i in TwDict[iref]:
    848                     for n in range(NTL):
    849                         H[i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
    850             TwMask = np.any(H,axis=-1)
    851         SQ = 1./(2.*refl[4])**2             # or (sin(theta)/lambda)**2
    852         SQfactor = 8.0*SQ*np.pi**2
    853         dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
    854         Bab = parmDict[phfx+'BabA']*dBabdA
    855         Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
    856         FF = refDict['FF']['FF'][iref].T[Tindx].T
    857         Uniq = np.inner(H,SGMT)             # array(nSGOp,3) or (nTwin,nSGOp,3)
    858         Phi = np.inner(H,SGT)
    859         phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
    860         sinp = np.sin(phase)
    861         cosp = np.cos(phase)
    862         occ = Mdata*Fdata/len(SGT)
    863         biso = -SQfactor*Uisodata[:,nxs]
    864         Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*nTwin,axis=1)
    865         HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
    866         Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])
    867         Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6)))
    868         Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
    869         Tcorr = (np.reshape(Tiso,Tuij.shape)*Tuij).T*occ
    870         fot = (FF+FP-Bab)*Tcorr
    871         fotp = FPP*Tcorr       
    872         fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
    873         fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
    874 #        GSASIIpath.IPyBreak()
    875         fas = np.sum(np.sum(fa,axis=-1),axis=-1)      #real sum over atoms & unique hkl array(2,nTwins)
    876         fbs = np.sum(np.sum(fb,axis=-1),axis=-1)      #imag sum over atoms & uniq hkl
    877         fax = np.array([-fot*sinp,-fotp*cosp])   #positions array(2,ntwi,nEqv,nAtoms)
    878         fbx = np.array([fot*cosp,-fotp*sinp])
    879         #sum below is over Uniq
    880         dfadfr = np.sum(fa/occ,axis=-2)        #array(2,ntwin,nAtom) Fdata != 0 avoids /0. problem
    881         dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
    882         dfadui = np.sum(-SQfactor*fa,axis=-2)
    883         if nTwin > 1:
    884             dfadx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
    885             dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fa,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
    886             # array(nTwin,2,nAtom,3) & array(nTwin,2,nAtom,6)
    887         else:
    888             dfadx = np.sum(twopi*Uniq*np.swapaxes(fax,-2,-1)[:,:,:,nxs],axis=-2)
    889             dfadua = np.sum(-Hij*np.swapaxes(fa,-2,-1)[:,:,:,nxs],axis=-2)
    890             # array(2,nAtom,3) & array(2,nAtom,6)
    891         if not SGData['SGInv']:
    892             dfbdfr = np.sum(fb/occ,axis=-2)        #Fdata != 0 avoids /0. problem
    893             dfadba /= 2.
    894             dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)/2.
    895             dfbdui = np.sum(-SQfactor*fb,axis=-2)
    896             if len(TwinLaw) > 1:
    897                 dfbdx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)])           
    898                 dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fb,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)])
    899             else:
    900                 dfadfl = np.sum(-FPP*Tcorr*sinp)
    901                 dfbdfl = np.sum(FPP*Tcorr*cosp)
    902                 dfbdx = np.sum(twopi*Uniq*np.swapaxes(fbx,-2,-1)[:,:,:,nxs],axis=2)           
    903                 dfbdua = np.sum(-Hij*np.swapaxes(fb,-2,-1)[:,:,:,nxs],axis=2)
    904         else:
    905             dfbdfr = np.zeros_like(dfadfr)
    906             dfbdx = np.zeros_like(dfadx)
    907             dfbdui = np.zeros_like(dfadui)
    908             dfbdua = np.zeros_like(dfadua)
    909             dfbdba = np.zeros_like(dfadba)
    910             dfadfl = 0.0
    911             dfbdfl = 0.0
    912         #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
    913         SA = fas[0]+fas[1]
    914         SB = fbs[0]+fbs[1]
    915         if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro
    916             dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)+   \
    917                 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
    918             dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])+  \
    919                 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
    920             dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])+   \
    921                 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
    922             dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])+   \
    923                 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
    924         else:
    925             if nTwin > 1:
    926                 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)]
    927                 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)]
    928                 dFdui[iref] = [2.*TwMask[it]*(SA[it]*dfadui[it][0]+SA[it]*dfadui[it][1]+SB[it]*dfbdui[it][0]+SB[it]*dfbdui[it][1]) for it in range(nTwin)]
    929                 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)]
    930                 dFdtw[iref] = np.sum(TwMask*fas,axis=0)**2+np.sum(TwMask*fbs,axis=0)**2
    931             else:   #these are good for no twin single crystals
    932                 dFdfr[iref] = (2.*SA*(dfadfr[0]+dfadfr[1])+2.*SB*(dfbdfr[0]+dfbdfr[1]))*Mdata/len(Uniq)
    933                 dFdx[iref] = 2.*SA*(dfadx[0]+dfadx[1])+2.*SB*(dfbdx[0]+dfbdx[1])
    934                 dFdui[iref] = 2.*SA*(dfadui[0]+dfadui[1])+2.*SB*(dfbdui[0]+dfbdui[1])
    935                 dFdua[iref] = 2.*SA*(dfadua[0]+dfadua[1])+2.*SB*(dfbdua[0]+dfbdua[1])
    936                 dFdfl[iref] = -SA*dfadfl-SB*dfbdfl  #array(nRef,)
    937         dFdbab[iref] = fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
    938             fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
    939 #        GSASIIpath.IPyBreak()
    940         if not iref%100 :
    941             print ' %d derivative time %.4f\r'%(iref,time.time()-time0),
    942     print ' %d derivative time %.4f\r'%(len(refDict['RefList']),time.time()-time0)
    943         #loop over atoms - each dict entry is list of derivatives for all the reflections
    944     if nTwin > 1:
    945         for i in range(len(Mdata)):     #these all OK?
    946             dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,nxs],axis=0)
    947             dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,nxs],axis=0)
    948             dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,nxs],axis=0)
    949             dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,nxs],axis=0)
    950             dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,nxs],axis=0)
    951             dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,nxs],axis=0)
    952             dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,nxs],axis=0)
    953             dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,nxs],axis=0)
    954             dFdvDict[pfx+'AU12:'+str(i)] = 2.*np.sum(dFdua.T[3][i]*TwinFr[:,nxs],axis=0)
    955             dFdvDict[pfx+'AU13:'+str(i)] = 2.*np.sum(dFdua.T[4][i]*TwinFr[:,nxs],axis=0)
    956             dFdvDict[pfx+'AU23:'+str(i)] = 2.*np.sum(dFdua.T[5][i]*TwinFr[:,nxs],axis=0)
    957     else:
    958         for i in range(len(Mdata)):
    959             dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
    960             dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
    961             dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
    962             dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
    963             dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
    964             dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
    965             dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
    966             dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
    967             dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
    968             dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
    969             dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
    970         dFdvDict[phfx+'Flack'] = 4.*dFdfl.T
    971     dFdvDict[phfx+'BabA'] = dFdbab.T[0]
    972     dFdvDict[phfx+'BabU'] = dFdbab.T[1]
    973     if nTwin > 1:
    974         for i in range(nTwin):
    975             dFdvDict[phfx+'TwinFr:'+str(i)] = dFdtw.T[i]
    976     return dFdvDict
     773#def StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
     774#    '''Compute structure factor derivatives on single reflections - keep as it works for twins
     775#    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
     789#    '''
     790#    phfx = pfx.split(':')[0]+hfx
     791#    ast = np.sqrt(np.diag(G))
     792#    Mast = twopisq*np.multiply.outer(ast,ast)
     793#    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
     794#    SGT = np.array([ops[1] for ops in SGData['SGOps']])
     795#    FFtables = calcControls['FFtables']
     796#    BLtables = calcControls['BLtables']
     797#    TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])
     798#    TwDict = refDict.get('TwDict',{})           
     799#    if 'S' in calcControls[hfx+'histType']:
     800#        NTL = calcControls[phfx+'NTL']
     801#        NM = calcControls[phfx+'TwinNMN']+1
     802#        TwinLaw = calcControls[phfx+'TwinLaw']
     803#        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
     804#        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
     805#    nTwin = len(TwinLaw)       
     806#    nRef = len(refDict['RefList'])
     807#    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
     808#    mSize = len(Mdata)
     809#    FF = np.zeros(len(Tdata))
     810#    if 'NC' in calcControls[hfx+'histType']:
     811#        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
     812#    elif 'X' in calcControls[hfx+'histType']:
     813#        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
     814#        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
     815#    Uij = np.array(G2lat.U6toUij(Uijdata))
     816#    bij = Mast*Uij.T
     817#    dFdvDict = {}
     818#    if nTwin > 1:
     819#        dFdfr = np.zeros((nRef,nTwin,mSize))
     820#        dFdx = np.zeros((nRef,nTwin,mSize,3))
     821#        dFdui = np.zeros((nRef,nTwin,mSize))
     822#        dFdua = np.zeros((nRef,nTwin,mSize,6))
     823#        dFdbab = np.zeros((nRef,nTwin,2))
     824#        dFdfl = np.zeros((nRef,nTwin))
     825#        dFdtw = np.zeros((nRef,nTwin))
     826#    else:
     827#        dFdfr = np.zeros((nRef,mSize))
     828#        dFdx = np.zeros((nRef,mSize,3))
     829#        dFdui = np.zeros((nRef,mSize))
     830#        dFdua = np.zeros((nRef,mSize,6))
     831#        dFdbab = np.zeros((nRef,2))
     832#        dFdfl = np.zeros((nRef))
     833#        dFdtw = np.zeros((nRef))
     834#    Flack = 1.0
     835#    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
     836#        Flack = 1.-2.*parmDict[phfx+'Flack']
     837#    time0 = time.time()
     838#    nref = len(refDict['RefList'])/100   
     839#    for iref,refl in enumerate(refDict['RefList']):
     840#        if 'T' in calcControls[hfx+'histType']:
     841#            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12])
     842#        H = np.array(refl[:3])
     843#        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(3,nTwins) or (3)
     844#        TwMask = np.any(H,axis=-1)
     845#        if TwinLaw.shape[0] > 1 and TwDict:
     846#            if iref in TwDict:
     847#                for i in TwDict[iref]:
     848#                    for n in range(NTL):
     849#                        H[i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
     850#            TwMask = np.any(H,axis=-1)
     851#        SQ = 1./(2.*refl[4])**2             # or (sin(theta)/lambda)**2
     852#        SQfactor = 8.0*SQ*np.pi**2
     853#        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
     854#        Bab = parmDict[phfx+'BabA']*dBabdA
     855#        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
     856#        FF = refDict['FF']['FF'][iref].T[Tindx].T
     857#        Uniq = np.inner(H,SGMT)             # array(nSGOp,3) or (nTwin,nSGOp,3)
     858#        Phi = np.inner(H,SGT)
     859#        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
     860#        sinp = np.sin(phase)
     861#        cosp = np.cos(phase)
     862#        occ = Mdata*Fdata/len(SGT)
     863#        biso = -SQfactor*Uisodata[:,nxs]
     864#        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*nTwin,axis=1)
     865#        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
     866#        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])
     867#        Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6)))
     868#        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
     869#        Tcorr = (np.reshape(Tiso,Tuij.shape)*Tuij).T*occ
     870#        fot = (FF+FP-Bab)*Tcorr
     871#        fotp = FPP*Tcorr       
     872#        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
     873#        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
     874##        GSASIIpath.IPyBreak()
     875#        fas = np.sum(np.sum(fa,axis=-1),axis=-1)      #real sum over atoms & unique hkl array(2,nTwins)
     876#        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)      #imag sum over atoms & uniq hkl
     877#        fax = np.array([-fot*sinp,-fotp*cosp])   #positions array(2,ntwi,nEqv,nAtoms)
     878#        fbx = np.array([fot*cosp,-fotp*sinp])
     879#        #sum below is over Uniq
     880#        dfadfr = np.sum(fa/occ,axis=-2)        #array(2,ntwin,nAtom) Fdata != 0 avoids /0. problem
     881#        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
     882#        dfadui = np.sum(-SQfactor*fa,axis=-2)
     883#        if nTwin > 1:
     884#            dfadx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
     885#            dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fa,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
     886#            # array(nTwin,2,nAtom,3) & array(nTwin,2,nAtom,6)
     887#        else:
     888#            dfadx = np.sum(twopi*Uniq*np.swapaxes(fax,-2,-1)[:,:,:,nxs],axis=-2)
     889#            dfadua = np.sum(-Hij*np.swapaxes(fa,-2,-1)[:,:,:,nxs],axis=-2)
     890#            # array(2,nAtom,3) & array(2,nAtom,6)
     891#        if not SGData['SGInv']:
     892#            dfbdfr = np.sum(fb/occ,axis=-2)        #Fdata != 0 avoids /0. problem
     893#            dfadba /= 2.
     894#            dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)/2.
     895#            dfbdui = np.sum(-SQfactor*fb,axis=-2)
     896#            if len(TwinLaw) > 1:
     897#                dfbdx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)])           
     898#                dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fb,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)])
     899#            else:
     900#                dfadfl = np.sum(-FPP*Tcorr*sinp)
     901#                dfbdfl = np.sum(FPP*Tcorr*cosp)
     902#                dfbdx = np.sum(twopi*Uniq*np.swapaxes(fbx,-2,-1)[:,:,:,nxs],axis=2)           
     903#                dfbdua = np.sum(-Hij*np.swapaxes(fb,-2,-1)[:,:,:,nxs],axis=2)
     904#        else:
     905#            dfbdfr = np.zeros_like(dfadfr)
     906#            dfbdx = np.zeros_like(dfadx)
     907#            dfbdui = np.zeros_like(dfadui)
     908#            dfbdua = np.zeros_like(dfadua)
     909#            dfbdba = np.zeros_like(dfadba)
     910#            dfadfl = 0.0
     911#            dfbdfl = 0.0
     912#        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
     913#        SA = fas[0]+fas[1]
     914#        SB = fbs[0]+fbs[1]
     915#        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro
     916#            dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)+   \
     917#                2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
     918#            dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])+  \
     919#                2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
     920#            dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])+   \
     921#                2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
     922#            dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])+   \
     923#                2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
     924#        else:
     925#            if nTwin > 1:
     926#                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)]
     927#                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)]
     928#                dFdui[iref] = [2.*TwMask[it]*(SA[it]*dfadui[it][0]+SA[it]*dfadui[it][1]+SB[it]*dfbdui[it][0]+SB[it]*dfbdui[it][1]) for it in range(nTwin)]
     929#                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)]
     930#                dFdtw[iref] = np.sum(TwMask*fas,axis=0)**2+np.sum(TwMask*fbs,axis=0)**2
     931#            else:   #these are good for no twin single crystals
     932#                dFdfr[iref] = (2.*SA*(dfadfr[0]+dfadfr[1])+2.*SB*(dfbdfr[0]+dfbdfr[1]))*Mdata/len(Uniq)
     933#                dFdx[iref] = 2.*SA*(dfadx[0]+dfadx[1])+2.*SB*(dfbdx[0]+dfbdx[1])
     934#                dFdui[iref] = 2.*SA*(dfadui[0]+dfadui[1])+2.*SB*(dfbdui[0]+dfbdui[1])
     935#                dFdua[iref] = 2.*SA*(dfadua[0]+dfadua[1])+2.*SB*(dfbdua[0]+dfbdua[1])
     936#                dFdfl[iref] = -SA*dfadfl-SB*dfbdfl  #array(nRef,)
     937#        dFdbab[iref] = fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
     938#            fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
     939##        GSASIIpath.IPyBreak()
     940#        if not iref%100 :
     941#            print ' %d derivative time %.4f\r'%(iref,time.time()-time0),
     942#    print ' %d derivative time %.4f\r'%(len(refDict['RefList']),time.time()-time0)
     943#        #loop over atoms - each dict entry is list of derivatives for all the reflections
     944#    if nTwin > 1:
     945#        for i in range(len(Mdata)):     #these all OK?
     946#            dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,nxs],axis=0)
     947#            dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,nxs],axis=0)
     948#            dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,nxs],axis=0)
     949#            dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,nxs],axis=0)
     950#            dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,nxs],axis=0)
     951#            dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,nxs],axis=0)
     952#            dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,nxs],axis=0)
     953#            dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,nxs],axis=0)
     954#            dFdvDict[pfx+'AU12:'+str(i)] = 2.*np.sum(dFdua.T[3][i]*TwinFr[:,nxs],axis=0)
     955#            dFdvDict[pfx+'AU13:'+str(i)] = 2.*np.sum(dFdua.T[4][i]*TwinFr[:,nxs],axis=0)
     956#            dFdvDict[pfx+'AU23:'+str(i)] = 2.*np.sum(dFdua.T[5][i]*TwinFr[:,nxs],axis=0)
     957#    else:
     958#        for i in range(len(Mdata)):
     959#            dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
     960#            dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
     961#            dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
     962#            dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
     963#            dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
     964#            dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
     965#            dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
     966#            dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
     967#            dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
     968#            dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
     969#            dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
     970#        dFdvDict[phfx+'Flack'] = 4.*dFdfl.T
     971#    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
     972#    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
     973#    if nTwin > 1:
     974#        for i in range(nTwin):
     975#            dFdvDict[phfx+'TwinFr:'+str(i)] = dFdtw.T[i]
     976#    return dFdvDict
    977977   
    978978def StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
     
    11341134    return dFdvDict
    11351135   
    1136 def StructureFactorDervTw(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
    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
    1151     '''
    1152     phfx = pfx.split(':')[0]+hfx
    1153     ast = np.sqrt(np.diag(G))
    1154     Mast = twopisq*np.multiply.outer(ast,ast)
    1155     SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
    1156     SGT = np.array([ops[1] for ops in SGData['SGOps']])
    1157     FFtables = calcControls['FFtables']
    1158     BLtables = calcControls['BLtables']
    1159     TwDict = refDict.get('TwDict',{})           
    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))
    1165     nTwin = len(TwinLaw)       
    1166     nRef = len(refDict['RefList'])
    1167     Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
    1168     mSize = len(Mdata)
    1169     FF = np.zeros(len(Tdata))
    1170     if 'NC' in calcControls[hfx+'histType']:
    1171         FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
    1172     elif 'X' in calcControls[hfx+'histType']:
    1173         FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
    1174         FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
    1175     Uij = np.array(G2lat.U6toUij(Uijdata))
    1176     bij = Mast*Uij.T
    1177     dFdvDict = {}
    1178     dFdfr = np.zeros((nRef,nTwin,mSize))
    1179     dFdx = np.zeros((nRef,nTwin,mSize,3))
    1180     dFdui = np.zeros((nRef,nTwin,mSize))
    1181     dFdua = np.zeros((nRef,nTwin,mSize,6))
    1182     dFdbab = np.zeros((nRef,nTwin,2))
    1183     dFdtw = np.zeros((nRef,nTwin))
    1184     time0 = time.time()
    1185     nref = len(refDict['RefList'])/100   
    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])
    1190         H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(3,nTwins) or (3)
    1191         TwMask = np.any(H,axis=-1)
    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
    1198         SQfactor = 8.0*SQ*np.pi**2
    1199         dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
    1200         Bab = parmDict[phfx+'BabA']*dBabdA
    1201         Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
    1202         FF = refDict['FF']['FF'][iref].T[Tindx].T
    1203         Uniq = np.inner(H,SGMT)             # array(nSGOp,3) or (nTwin,nSGOp,3)
    1204         Phi = np.inner(H,SGT)
    1205         phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
    1206         sinp = np.sin(phase)
    1207         cosp = np.cos(phase)
    1208         occ = Mdata*Fdata/len(SGT)
    1209         biso = -SQfactor*Uisodata[:,nxs]
    1210         Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*nTwin,axis=1)
    1211         HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
    1212         Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])
    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
    1217         fotp = FPP*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])
    1220 #        GSASIIpath.IPyBreak()
    1221         fas = np.sum(np.sum(fa,axis=-1),axis=-1)      #real sum over atoms & unique hkl array(2,nTwins)
    1222         fbs = np.sum(np.sum(fb,axis=-1),axis=-1)      #imag sum over atoms & uniq hkl
    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)
    1227         fbx = np.array([fot*cosp,-fotp*sinp])
    1228         #sum below is over Uniq
    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)
    1235         if not SGData['SGInv']:
    1236             dfbdfr = np.sum(fb/occ,axis=-2)        #Fdata != 0 avoids /0. problem
    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)])
    1242         else:
    1243             dfbdfr = np.zeros_like(dfadfr)
    1244             dfbdx = np.zeros_like(dfadx)
    1245             dfbdui = np.zeros_like(dfadui)
    1246             dfbdua = np.zeros_like(dfadua)
    1247             dfbdba = np.zeros_like(dfadba)
    1248         SA = fas[0]+fas[1]
    1249         SB = fbs[0]+fbs[1]
    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
    1260 #        GSASIIpath.IPyBreak()
    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
    1263     for i in range(len(Mdata)):     #these all OK?
    1264         dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,nxs],axis=0)
    1265         dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,nxs],axis=0)
    1266         dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,nxs],axis=0)
    1267         dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,nxs],axis=0)
    1268         dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,nxs],axis=0)
    1269         dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,nxs],axis=0)
    1270         dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,nxs],axis=0)
    1271         dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,nxs],axis=0)
    1272         dFdvDict[pfx+'AU12:'+str(i)] = 2.*np.sum(dFdua.T[3][i]*TwinFr[:,nxs],axis=0)
    1273         dFdvDict[pfx+'AU13:'+str(i)] = 2.*np.sum(dFdua.T[4][i]*TwinFr[:,nxs],axis=0)
    1274         dFdvDict[pfx+'AU23:'+str(i)] = 2.*np.sum(dFdua.T[5][i]*TwinFr[:,nxs],axis=0)
    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
     1136#def StructureFactorDervTw(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
     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
     1151#    '''
     1152#    phfx = pfx.split(':')[0]+hfx
     1153#    ast = np.sqrt(np.diag(G))
     1154#    Mast = twopisq*np.multiply.outer(ast,ast)
     1155#    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
     1156#    SGT = np.array([ops[1] for ops in SGData['SGOps']])
     1157#    FFtables = calcControls['FFtables']
     1158#    BLtables = calcControls['BLtables']
     1159#    TwDict = refDict.get('TwDict',{})           
     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))
     1165#    nTwin = len(TwinLaw)       
     1166#    nRef = len(refDict['RefList'])
     1167#    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
     1168#    mSize = len(Mdata)
     1169#    FF = np.zeros(len(Tdata))
     1170#    if 'NC' in calcControls[hfx+'histType']:
     1171#        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
     1172#    elif 'X' in calcControls[hfx+'histType']:
     1173#        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
     1174#        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
     1175#    Uij = np.array(G2lat.U6toUij(Uijdata))
     1176#    bij = Mast*Uij.T
     1177#    dFdvDict = {}
     1178#    dFdfr = np.zeros((nRef,nTwin,mSize))
     1179#    dFdx = np.zeros((nRef,nTwin,mSize,3))
     1180#    dFdui = np.zeros((nRef,nTwin,mSize))
     1181#    dFdua = np.zeros((nRef,nTwin,mSize,6))
     1182#    dFdbab = np.zeros((nRef,nTwin,2))
     1183#    dFdtw = np.zeros((nRef,nTwin))
     1184#    time0 = time.time()
     1185#    nref = len(refDict['RefList'])/100   
     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])
     1190#        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(3,nTwins) or (3)
     1191#        TwMask = np.any(H,axis=-1)
     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
     1198#        SQfactor = 8.0*SQ*np.pi**2
     1199#        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
     1200#        Bab = parmDict[phfx+'BabA']*dBabdA
     1201#        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
     1202#        FF = refDict['FF']['FF'][iref].T[Tindx].T
     1203#        Uniq = np.inner(H,SGMT)             # array(nSGOp,3) or (nTwin,nSGOp,3)
     1204#        Phi = np.inner(H,SGT)
     1205#        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
     1206#        sinp = np.sin(phase)
     1207#        cosp = np.cos(phase)
     1208#        occ = Mdata*Fdata/len(SGT)
     1209#        biso = -SQfactor*Uisodata[:,nxs]
     1210#        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*nTwin,axis=1)
     1211#        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
     1212#        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])
     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
     1217#        fotp = FPP*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])
     1220##        GSASIIpath.IPyBreak()
     1221#        fas = np.sum(np.sum(fa,axis=-1),axis=-1)      #real sum over atoms & unique hkl array(2,nTwins)
     1222#        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)      #imag sum over atoms & uniq hkl
     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)
     1227#        fbx = np.array([fot*cosp,-fotp*sinp])
     1228#        #sum below is over Uniq
     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)
     1235#        if not SGData['SGInv']:
     1236#            dfbdfr = np.sum(fb/occ,axis=-2)        #Fdata != 0 avoids /0. problem
     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)])
     1242#        else:
     1243#            dfbdfr = np.zeros_like(dfadfr)
     1244#            dfbdx = np.zeros_like(dfadx)
     1245#            dfbdui = np.zeros_like(dfadui)
     1246#            dfbdua = np.zeros_like(dfadua)
     1247#            dfbdba = np.zeros_like(dfadba)
     1248#        SA = fas[0]+fas[1]
     1249#        SB = fbs[0]+fbs[1]
     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
     1260##        GSASIIpath.IPyBreak()
     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
     1263#    for i in range(len(Mdata)):     #these all OK?
     1264#        dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,nxs],axis=0)
     1265#        dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,nxs],axis=0)
     1266#        dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,nxs],axis=0)
     1267#        dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,nxs],axis=0)
     1268#        dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,nxs],axis=0)
     1269#        dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,nxs],axis=0)
     1270#        dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,nxs],axis=0)
     1271#        dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,nxs],axis=0)
     1272#        dFdvDict[pfx+'AU12:'+str(i)] = 2.*np.sum(dFdua.T[3][i]*TwinFr[:,nxs],axis=0)
     1273#        dFdvDict[pfx+'AU13:'+str(i)] = 2.*np.sum(dFdua.T[4][i]*TwinFr[:,nxs],axis=0)
     1274#        dFdvDict[pfx+'AU23:'+str(i)] = 2.*np.sum(dFdua.T[5][i]*TwinFr[:,nxs],axis=0)
     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
    12801280   
    12811281def StructureFactorDervTw2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
Note: See TracChangeset for help on using the changeset viewer.