Changeset 2097


Ignore:
Timestamp:
Dec 18, 2015 10:58:57 AM (6 years ago)
Author:
vondreele
Message:

name fix in AtomTypeSelect? - remove ( & )
split StructureFactorDeriv? (SS & non SS versions) into twin & nontwin/PWDR versions
non twin versions all test OK for refl blocks, twin versions need work (use single refl version for now). SS versions do single refl for now.

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIphsGUI.py

    r2077 r2097  
    13921392                        name = atomData[r][c]
    13931393                        if len(name) in [2,4]:
    1394                             atomData[r][c-1] = name[:2]+'(%d)'%(r+1)
     1394                            atomData[r][c-1] = name[:2]+'%d'%(r+1)
    13951395                        else:
    1396                             atomData[r][c-1] = name[:1]+'(%d)'%(r+1)
     1396                            atomData[r][c-1] = name[:1]+'%d'%(r+1)
    13971397                PE.Destroy()
    13981398                SetupGeneral()
    … …  
    20532053            else:   
    20542054                data['Atoms'] = result
     2055            Atoms.ClearSelection()
     2056            data['Drawing']['Atoms'] = []
    20552057            OnReloadDrawAtoms(event)           
    20562058            FillAtomsGrid(Atoms)
    … …  
    20642066    def OnDistAnglePrt(event):
    20652067        'save distances and angles to a file'   
    2066         fp = file(os.path.abspath(os.path.splitext(G2frame.GSASprojectfile
    2067                                                    )[0]+'.disagl'),'w')
     2068        fp = file(os.path.abspath(os.path.splitext(G2frame.GSASprojectfile)[0]+'.disagl'),'w')
    20682069        OnDistAngle(event,fp=fp)
    20692070        fp.close()

  • trunk/GSASIIstrMath.py

    r2096 r2097  
    648648    ''' Compute structure factors for all h,k,l for phase
    649649    puts the result, F^2, in each ref[8] in refList
     650    operates on blocks of 100 reflections for speed
    650651    input:
    651652   
    652653    :param dict refDict: where
    653         'RefList' list where each ref = h,k,l,m,d,...
     654        'RefList' list where each ref = h,k,l,it,d,...
    654655        'FF' dict of form factors - filed in below
    655656    :param np.array G:      reciprocal metric tensor
    … …  
    769770   
    770771def StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
    771     '''Compute structure factor derivatives on blocks of reflections - keep as it works for twins
     772    '''Compute structure factor derivatives on single reflections - keep as it works for twins
    772773    but is slower for powders/nontwins
    773774    '''
    … …  
    962963   
    963964def StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
    964     '''Compute structure factor derivatives on blocks of reflections- works for powders/nontwins
     965    '''Compute structure factor derivatives on blocks of reflections - for powders/nontwins only
    965966    faster than StructureFactorDerv
     967    '''
     968    phfx = pfx.split(':')[0]+hfx
     969    ast = np.sqrt(np.diag(G))
     970    Mast = twopisq*np.multiply.outer(ast,ast)
     971    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
     972    SGT = np.array([ops[1] for ops in SGData['SGOps']])
     973    FFtables = calcControls['FFtables']
     974    BLtables = calcControls['BLtables']
     975    nRef = len(refDict['RefList'])
     976    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
     977    mSize = len(Mdata)
     978    FF = np.zeros(len(Tdata))
     979    if 'NC' in calcControls[hfx+'histType']:
     980        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
     981    elif 'X' in calcControls[hfx+'histType']:
     982        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
     983        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
     984    Uij = np.array(G2lat.U6toUij(Uijdata))
     985    bij = Mast*Uij.T
     986    dFdvDict = {}
     987    dFdfr = np.zeros((nRef,mSize))
     988    dFdx = np.zeros((nRef,mSize,3))
     989    dFdui = np.zeros((nRef,mSize))
     990    dFdua = np.zeros((nRef,mSize,6))
     991    dFdbab = np.zeros((nRef,2))
     992    dFdfl = np.zeros((nRef))
     993    dFdtw = np.zeros((nRef))
     994    Flack = 1.0
     995    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
     996        Flack = 1.-2.*parmDict[phfx+'Flack']
     997    time0 = time.time()
     998    nref = len(refDict['RefList'])/100   
     999#reflection processing begins here - big arrays!
     1000    iBeg = 0
     1001    blkSize = 32       #no. of reflections in a block - optimized for speed
     1002    while iBeg < nRef:
     1003        iFin = min(iBeg+blkSize,nRef)
     1004        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
     1005        H = refl.T[:3]
     1006        SQ = 1./(2.*refl.T[4])**2             # or (sin(theta)/lambda)**2
     1007        SQfactor = 8.0*SQ*np.pi**2
     1008        if 'T' in calcControls[hfx+'histType']:
     1009            if 'P' in calcControls[hfx+'histType']:
     1010                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
     1011            else:
     1012                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
     1013            FP = np.repeat(FP.T,len(SGT),axis=0)
     1014            FPP = np.repeat(FPP.T,len(SGT),axis=0)
     1015        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
     1016        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT))
     1017        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
     1018        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT),axis=0)
     1019        Uniq = np.inner(H.T,SGMT)             # array(nSGOp,3) or (nTwin,nSGOp,3)
     1020        Phi = np.inner(H.T,SGT)
     1021        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
     1022        sinp = np.sin(phase)        #refBlk x nOps x nAtoms
     1023        cosp = np.cos(phase)
     1024        occ = Mdata*Fdata/len(SGT)
     1025        biso = -SQfactor*Uisodata[:,nxs]
     1026        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT),axis=1).T
     1027        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
     1028        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
     1029        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/len(SGMT)
     1030        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])
     1031        Hij = np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(-1,len(SGT),6))
     1032        fot = np.reshape(((FF+FP).T-Bab).T,cosp.shape)*Tcorr
     1033        fotp = FPP*Tcorr       
     1034        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])
     1037        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])
     1040#        GSASIIpath.IPyBreak()
     1041        fas = np.sum(np.sum(fa,axis=-1),axis=-1)      #real sum over atoms & unique hkl array(2,refBlk,nTwins)
     1042        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)      #imag sum over atoms & uniq hkl
     1043        fax = np.array([-fot*sinp,-fotp*cosp])   #positions array(2,refBlk,ntwi,nEqv,nAtoms)
     1044        fbx = np.array([fot*cosp,-fotp*sinp])
     1045        #sum below is over Uniq
     1046        dfadfr = np.sum(fa/occ,axis=-2)        #array(2,refBlk,ntwin,nAtom) Fdata != 0 avoids /0. problem
     1047        dfadba = np.sum(-cosp*Tcorr,axis=-2)  #array(refBlk,nAtom)
     1048        dfadx = np.sum(twopi*Uniq[nxs,:,nxs,:,:]*np.swapaxes(fax,-2,-1)[:,:,:,:,nxs],axis=-2)
     1049        dfadui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fa,axis=-2) #array(Ops,refBlk,nTw,nAtoms)
     1050        dfadua = np.sum(-Hij[nxs,:,nxs,:,:]*np.swapaxes(fa,-2,-1)[:,:,:,:,nxs],axis=-2)
     1051        # array(2,refBlk,nAtom,3) & array(2,refBlk,nAtom,6)
     1052        if not SGData['SGInv']:
     1053            dfbdfr = np.sum(fb/occ,axis=-2)        #Fdata != 0 avoids /0. problem
     1054            dfbdba = np.sum(-sinp*Tcorr,axis=-2)
     1055            dfadfl = np.sum(np.sum(-FPP*Tcorr*sinp,axis=-1),axis=-1)
     1056            dfbdfl = np.sum(np.sum(FPP*Tcorr*cosp,axis=-1),axis=-1)
     1057            dfbdx = np.sum(twopi*Uniq[nxs,:,nxs,:,:]*np.swapaxes(fbx,-2,-1)[:,:,:,:,nxs],axis=-2)           
     1058            dfbdui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fb,axis=-2)
     1059            dfbdua = np.sum(-Hij[nxs,:,nxs,:,:]*np.swapaxes(fb,-2,-1)[:,:,:,:,nxs],axis=-2)
     1060        else:
     1061            dfbdfr = np.zeros_like(dfadfr)
     1062            dfbdx = np.zeros_like(dfadx)
     1063            dfbdui = np.zeros_like(dfadui)
     1064            dfbdua = np.zeros_like(dfadua)
     1065            dfbdba = np.zeros_like(dfadba)
     1066            dfadfl = 0.0
     1067            dfbdfl = 0.0
     1068        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
     1069        SA = fas[0]+fas[1]
     1070        SB = fbs[0]+fbs[1]
     1071        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro
     1072            dFdfr[iBeg:iFin] = 2.*(fas[0,:,nxs]*dfadfr[0]+fas[1,:,nxs]*dfadfr[1])*Mdata/len(SGMT)+   \
     1073                2.*(fbs[0,:,nxs]*dfbdfr[0]-fbs[1,:,nxs]*dfbdfr[1])*Mdata/len(SGMT)
     1074            dFdx[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadx[0]+fas[1,:,nxs,nxs]*dfadx[1])+  \
     1075                2.*(fbs[0,:,nxs,nxs]*dfbdx[0]+fbs[1,:,nxs,nxs]*dfbdx[1])
     1076            dFdui[iBeg:iFin] = 2.*(fas[0,:,nxs]*dfadui[0]+fas[1,:,nxs]*dfadui[1])+   \
     1077                2.*(fbs[0,:,nxs]*dfbdui[0]-fbs[1,:,nxs]*dfbdui[1])
     1078            dFdua[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadua[0]+fas[1,:,nxs,nxs]*dfadua[1])+   \
     1079                2.*(fbs[0,:,nxs,nxs]*dfbdua[0]+fbs[1,:,nxs,nxs]*dfbdua[1])
     1080        else:
     1081            dFdfr[iBeg:iFin] = (2.*SA[:,nxs]*(dfadfr[0]+dfadfr[1])+2.*SB[:,nxs]*(dfbdfr[0]+dfbdfr[1]))*Mdata/len(SGMT)
     1082            dFdx[iBeg:iFin] = 2.*SA[:,nxs,nxs]*(dfadx[0]+dfadx[1])+2.*SB[:,nxs,nxs]*(dfbdx[0]+dfbdx[1])
     1083            dFdui[iBeg:iFin] = 2.*SA[:,nxs]*(dfadui[0]+dfadui[1])+2.*SB[:,nxs]*(dfbdui[0]+dfbdui[1])
     1084            dFdua[iBeg:iFin] = 2.*SA[:,nxs,nxs]*(dfadua[0]+dfadua[1])+2.*SB[:,nxs,nxs]*(dfbdua[0]+dfbdua[1])
     1085            dFdfl[iBeg:iFin] = -SA*dfadfl-SB*dfbdfl  #array(nRef,)
     1086        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)])+ \
     1087                            fbs[0,nxs]*np.array([np.sum(dfbdba.T*dBabdA,axis=0),np.sum(-dfbdba.T*parmDict[phfx+'BabA']*SQfactor*dBabdA,axis=0)])).T
     1088#        GSASIIpath.IPyBreak()
     1089        iBeg += blkSize
     1090    print ' %d derivative time %.4f\r'%(nRef,time.time()-time0)
     1091        #loop over atoms - each dict entry is list of derivatives for all the reflections
     1092    for i in range(len(Mdata)):
     1093        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
     1094        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
     1095        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
     1096        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
     1097        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
     1098        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
     1099        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
     1100        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
     1101        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
     1102        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
     1103        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
     1104    dFdvDict[phfx+'Flack'] = 4.*dFdfl.T
     1105    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
     1106    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
     1107    return dFdvDict
     1108   
     1109def 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
    9661112    '''
    9671113    phfx = pfx.split(':')[0]+hfx
    … …  
    9931139    bij = Mast*Uij.T
    9941140    dFdvDict = {}
    995     if nTwin > 1:
    996         dFdfr = np.zeros((nRef,nTwin,mSize))
    997         dFdx = np.zeros((nRef,nTwin,mSize,3))
    998         dFdui = np.zeros((nRef,nTwin,mSize))
    999         dFdua = np.zeros((nRef,nTwin,mSize,6))
    1000         dFdbab = np.zeros((nRef,nTwin,2))
    1001         dFdfl = np.zeros((nRef,nTwin))
    1002         dFdtw = np.zeros((nRef,nTwin))
    1003     else:
    1004         dFdfr = np.zeros((nRef,mSize))
    1005         dFdx = np.zeros((nRef,mSize,3))
    1006         dFdui = np.zeros((nRef,mSize))
    1007         dFdua = np.zeros((nRef,mSize,6))
    1008         dFdbab = np.zeros((nRef,2))
    1009         dFdfl = np.zeros((nRef))
    1010         dFdtw = np.zeros((nRef))
     1141    dFdfr = np.zeros((nRef,nTwin,mSize))
     1142    dFdx = np.zeros((nRef,nTwin,mSize,3))
     1143    dFdui = np.zeros((nRef,nTwin,mSize))
     1144    dFdua = np.zeros((nRef,nTwin,mSize,6))
     1145    dFdbab = np.zeros((nRef,nTwin,2))
     1146    dFdfl = np.zeros((nRef,nTwin))
     1147    dFdtw = np.zeros((nRef,nTwin))
    10111148    Flack = 1.0
    10121149    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
    … …  
    10231160        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(3,nTwins) or (3)
    10241161        TwMask = np.any(H,axis=-1)
    1025         if TwinLaw.shape[0] > 1 and TwDict:
    1026             for ir in range(blkSize):
    1027                 iref = ir+iBeg
    1028                 if iref in TwDict:
    1029                     for i in TwDict[iref]:
    1030                         for n in range(NTL):
    1031                             H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
    1032             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)
    10331169        SQ = 1./(2.*refl.T[4])**2             # or (sin(theta)/lambda)**2
    10341170        SQfactor = 8.0*SQ*np.pi**2
    … …  
    10561192        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/len(SGMT)
    10571193        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])
    1058         if nTwin > 1:
    1059             Hij = np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(-1,nTwin,len(SGT),6))
    1060         else:
    1061             Hij = np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(-1,len(SGT),6))
     1194        Hij = np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(-1,nTwin,len(SGT),6))
    10621195        fot = np.reshape(((FF+FP).T-Bab).T,cosp.shape)*Tcorr
    10631196        fotp = FPP*Tcorr       
    … …  
    10761209        dfadfr = np.sum(fa/occ,axis=-2)        #array(2,refBlk,ntwin,nAtom) Fdata != 0 avoids /0. problem
    10771210        dfadba = np.sum(-cosp*Tcorr,axis=-2)  #array(refBlk,nAtom)
    1078         if nTwin > 1:
    1079             dfadfr = np.swapaxes(dfadfr,0,1)
    1080             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)
    1081             dfadui = np.swapaxes(np.sum(-SQfactor[nxs,:,nxs,nxs,nxs]*fa,axis=-2),0,1) #array(Ops,refBlk,nTw,nAtoms)
    1082             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)
    1083             # array(nTwin,2,refBlk,nAtom,3) & array(nTwin,2,refBlk,nAtom,6)
    1084         else:
    1085             dfadx = np.sum(twopi*Uniq[nxs,:,nxs,:,:]*np.swapaxes(fax,-2,-1)[:,:,:,:,nxs],axis=-2)
    1086             dfadui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fa,axis=-2) #array(Ops,refBlk,nTw,nAtoms)
    1087             dfadua = np.sum(-Hij[nxs,:,nxs,:,:]*np.swapaxes(fa,-2,-1)[:,:,:,:,nxs],axis=-2)
    1088             # array(2,refBlk,nAtom,3) & array(2,refBlk,nAtom,6)
     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)
    10891216        if not SGData['SGInv']:
    10901217            dfbdfr = np.sum(fb/occ,axis=-2)        #Fdata != 0 avoids /0. problem
    10911218            dfbdba = np.sum(-sinp*Tcorr,axis=-2)
    1092             if len(TwinLaw) > 1:
    1093                 dfbdfr = np.swapaxes(dfbdfr,0,1)
    1094                 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)         
    1095                 dfbdui = np.swapaxes(np.sum(-SQfactor[nxs,:,nxs,nxs,nxs]*fb,axis=-2),0,1)      #array(Ops,refBlk,nTw,nAtoms)
    1096                 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)
    1097             else:
    1098                 dfadfl = np.sum(np.sum(-FPP*Tcorr*sinp,axis=-1),axis=-1)
    1099                 dfbdfl = np.sum(np.sum(FPP*Tcorr*cosp,axis=-1),axis=-1)
    1100                 dfbdx = np.sum(twopi*Uniq[nxs,:,nxs,:,:]*np.swapaxes(fbx,-2,-1)[:,:,:,:,nxs],axis=-2)           
    1101                 dfbdui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fb,axis=-2)
    1102                 dfbdua = np.sum(-Hij[nxs,:,nxs,:,:]*np.swapaxes(fb,-2,-1)[:,:,:,:,nxs],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)
    11031223        else:
    11041224            dfbdfr = np.zeros_like(dfadfr)
    … …  
    11121232        SA = fas[0]+fas[1]
    11131233        SB = fbs[0]+fbs[1]
    1114         if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro
    1115             dFdfr[iBeg:iFin] = 2.*(fas[0,:,nxs]*dfadfr[0]+fas[1,:,nxs]*dfadfr[1])*Mdata/len(SGMT)+   \
    1116                 2.*(fbs[0,:,nxs]*dfbdfr[0]-fbs[1,:,nxs]*dfbdfr[1])*Mdata/len(SGMT)
    1117             dFdx[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadx[0]+fas[1,:,nxs,nxs]*dfadx[1])+  \
    1118                 2.*(fbs[0,:,nxs,nxs]*dfbdx[0]+fbs[1,:,nxs,nxs]*dfbdx[1])
    1119             dFdui[iBeg:iFin] = 2.*(fas[0,:,nxs]*dfadui[0]+fas[1,:,nxs]*dfadui[1])+   \
    1120                 2.*(fbs[0,:,nxs]*dfbdui[0]-fbs[1,:,nxs]*dfbdui[1])
    1121             dFdua[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadua[0]+fas[1,:,nxs,nxs]*dfadua[1])+   \
    1122                 2.*(fbs[0,:,nxs,nxs]*dfbdua[0]+fbs[1,:,nxs,nxs]*dfbdua[1])
    1123         else:
    1124             if nTwin > 1:
    1125                 dFdfr[iBeg:iFin] = np.swapaxes(np.array([2.*TwMask[:,it,nxs]*(SA[:,it,nxs]*(np.sum(dfadfr[:,:,it],axis=1))+ \
    1126                     SB[:,it,nxs]*(np.sum(dfbdfr[:,:,it],axis=1)))*Mdata/len(SGMT) for it in range(nTwin)]),0,1)
    1127                 dFdx[iBeg:iFin] = np.swapaxes(np.array([2.*TwMask[:,it,nxs,nxs]*(SA[:,it,nxs,nxs]*np.sum(dfadx[:,:,it],axis=1)+    \
    1128                     SB[:,it,nxs,nxs]*np.sum(dfbdx[:,:,it],axis=1)) for it in range(nTwin)]),0,1)
    1129                 dFdui[iBeg:iFin] = np.swapaxes(np.array([2.*TwMask[:,it,nxs]*(SA[:,it,nxs]*np.sum(dfadui[:,:,it],axis=1)+  \
    1130                     SB[:,it,nxs]*np.sum(dfbdui[:,:,it],axis=1)) for it in range(nTwin)]),0,1)
    1131                 dFdua[iBeg:iFin] = np.swapaxes(np.array([2.*TwMask[:,it,nxs,nxs]*(SA[:,it,nxs,nxs]*np.sum(dfadua[:,:,it],axis=1)+  \
    1132                     SB[:,it,nxs,nxs]*np.sum(dfbdua[:,:,it],axis=1)) for it in range(nTwin)]),0,1)
    1133                 dFdtw[iBeg:iFin] = np.sum(TwMask*fas,axis=0)**2+np.sum(TwMask*fbs,axis=0)**2
     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
    11341243               
    1135             else:   #these are good for no twin single crystals
    1136                 dFdfr[iBeg:iFin] = (2.*SA[:,nxs]*(dfadfr[0]+dfadfr[1])+2.*SB[:,nxs]*(dfbdfr[0]+dfbdfr[1]))*Mdata/len(SGMT)
    1137                 dFdx[iBeg:iFin] = 2.*SA[:,nxs,nxs]*(dfadx[0]+dfadx[1])+2.*SB[:,nxs,nxs]*(dfbdx[0]+dfbdx[1])
    1138                 dFdui[iBeg:iFin] = 2.*SA[:,nxs]*(dfadui[0]+dfadui[1])+2.*SB[:,nxs]*(dfbdui[0]+dfbdui[1])
    1139                 dFdua[iBeg:iFin] = 2.*SA[:,nxs,nxs]*(dfadua[0]+dfadua[1])+2.*SB[:,nxs,nxs]*(dfbdua[0]+dfbdua[1])
    1140                 dFdfl[iBeg:iFin] = -SA*dfadfl-SB*dfbdfl  #array(nRef,)
    1141         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)])+ \
    1142                             fbs[0,nxs]*np.array([np.sum(dfbdba.T*dBabdA,axis=0),np.sum(-dfbdba.T*parmDict[phfx+'BabA']*SQfactor*dBabdA,axis=0)])).T
     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
    11431246#        GSASIIpath.IPyBreak()
    11441247        iBeg += blkSize
    11451248    print ' %d derivative time %.4f\r'%(nRef,time.time()-time0)
    11461249        #loop over atoms - each dict entry is list of derivatives for all the reflections
    1147     if nTwin > 1:
    1148         for i in range(len(Mdata)):     #these all OK?
    1149             dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,nxs],axis=0)
    1150             dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,nxs],axis=0)
    1151             dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,nxs],axis=0)
    1152             dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,nxs],axis=0)
    1153             dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,nxs],axis=0)
    1154             dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,nxs],axis=0)
    1155             dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,nxs],axis=0)
    1156             dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,nxs],axis=0)
    1157             dFdvDict[pfx+'AU12:'+str(i)] = 2.*np.sum(dFdua.T[3][i]*TwinFr[:,nxs],axis=0)
    1158             dFdvDict[pfx+'AU13:'+str(i)] = 2.*np.sum(dFdua.T[4][i]*TwinFr[:,nxs],axis=0)
    1159             dFdvDict[pfx+'AU23:'+str(i)] = 2.*np.sum(dFdua.T[5][i]*TwinFr[:,nxs],axis=0)
    1160     else:
    1161         for i in range(len(Mdata)):
    1162             dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
    1163             dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
    1164             dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
    1165             dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
    1166             dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
    1167             dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
    1168             dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
    1169             dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
    1170             dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
    1171             dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
    1172             dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
    1173         dFdvDict[phfx+'Flack'] = 4.*dFdfl.T
     1250    for i in range(len(Mdata)):     #these all OK?
     1251        dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,nxs],axis=0)
     1252        dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,nxs],axis=0)
     1253        dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,nxs],axis=0)
     1254        dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,nxs],axis=0)
     1255        dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,nxs],axis=0)
     1256        dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,nxs],axis=0)
     1257        dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,nxs],axis=0)
     1258        dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,nxs],axis=0)
     1259        dFdvDict[pfx+'AU12:'+str(i)] = 2.*np.sum(dFdua.T[3][i]*TwinFr[:,nxs],axis=0)
     1260        dFdvDict[pfx+'AU13:'+str(i)] = 2.*np.sum(dFdua.T[4][i]*TwinFr[:,nxs],axis=0)
     1261        dFdvDict[pfx+'AU23:'+str(i)] = 2.*np.sum(dFdua.T[5][i]*TwinFr[:,nxs],axis=0)
     1262    dFdvDict[phfx+'Flack'] = 4.*dFdfl.T
    11741263    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
    11751264    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
    1176     if nTwin > 1:
    1177         for i in range(nTwin):
    1178             dFdvDict[phfx+'TwinFr:'+str(i)] = dFdtw.T[i]
     1265    for i in range(nTwin):
     1266        dFdvDict[phfx+'TwinFr:'+str(i)] = dFdtw.T[i]
    11791267    return dFdvDict
    11801268   
    … …  
    11831271    Compute super structure factors for all h,k,l,m for phase
    11841272    puts the result, F^2, in each ref[8+im] in refList
     1273    works on blocks of 100 reflections for speed
    11851274    input:
    11861275   
    11871276    :param dict refDict: where
    1188         'RefList' list where each ref = h,k,l,m,d,...
     1277        'RefList' list where each ref = h,k,l,m,it,d,...
    11891278        'FF' dict of form factors - filed in below
    11901279    :param np.array G:      reciprocal metric tensor
    … …  
    13181407
    13191408def SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
     1409    'Needs a doc string - no twins'
     1410    phfx = pfx.split(':')[0]+hfx
     1411    ast = np.sqrt(np.diag(G))
     1412    Mast = twopisq*np.multiply.outer(ast,ast)
     1413    SGInv = SGData['SGInv']
     1414    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
     1415    SGT = np.array([ops[1] for ops in SGData['SGOps']])
     1416    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
     1417    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
     1418    FFtables = calcControls['FFtables']
     1419    BLtables = calcControls['BLtables']
     1420    nRef = len(refDict['RefList'])
     1421    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
     1422    mSize = len(Mdata)  #no. atoms
     1423    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
     1424    ngl,nWaves,Fmod,Xmod,Umod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,Mast)
     1425    waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,Mast)
     1426    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
     1427    FF = np.zeros(len(Tdata))
     1428    if 'NC' in calcControls[hfx+'histType']:
     1429        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
     1430    elif 'X' in calcControls[hfx+'histType']:
     1431        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
     1432        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
     1433    Uij = np.array(G2lat.U6toUij(Uijdata)).T
     1434    bij = Mast*Uij
     1435    if not len(refDict['FF']):
     1436        if 'N' in calcControls[hfx+'histType']:
     1437            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
     1438        else:
     1439            dat = G2el.getFFvalues(FFtables,0.)       
     1440        refDict['FF']['El'] = dat.keys()
     1441        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))
     1442    dFdvDict = {}
     1443    dFdfr = np.zeros((nRef,mSize))
     1444    dFdx = np.zeros((nRef,mSize,3))
     1445    dFdui = np.zeros((nRef,mSize))
     1446    dFdua = np.zeros((nRef,mSize,6))
     1447    dFdbab = np.zeros((nRef,2))
     1448    dFdfl = np.zeros((nRef))
     1449    dFdtw = np.zeros((nRef))
     1450    dFdGf = np.zeros((nRef,mSize,FSSdata.shape[1],2))
     1451    dFdGx = np.zeros((nRef,mSize,XSSdata.shape[1],6))
     1452    dFdGz = np.zeros((nRef,mSize,5))
     1453    dFdGu = np.zeros((nRef,mSize,USSdata.shape[1],12))
     1454    Flack = 1.0
     1455    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
     1456        Flack = 1.-2.*parmDict[phfx+'Flack']
     1457    time0 = time.time()
     1458    nRef = len(refDict['RefList'])/100
     1459    for iref,refl in enumerate(refDict['RefList']):
     1460        if 'T' in calcControls[hfx+'histType']:
     1461            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12+im])
     1462        H = np.array(refl[:4])
     1463        HP = H[:3]+modQ*H[3:]            #projected hklm to hkl
     1464        SQ = 1./(2.*refl[4+im])**2             # or (sin(theta)/lambda)**2
     1465        SQfactor = 8.0*SQ*np.pi**2
     1466        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
     1467        Bab = parmDict[phfx+'BabA']*dBabdA
     1468        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
     1469        FF = refDict['FF']['FF'][iref].T[Tindx]
     1470        Uniq = np.inner(H,SSGMT)
     1471        Phi = np.inner(H,SSGT)
     1472        UniqP = np.inner(HP,SGMT)
     1473        if SGInv:   #if centro - expand HKL sets
     1474            Uniq = np.vstack((Uniq,-Uniq))
     1475            Phi = np.hstack((Phi,-Phi))
     1476            UniqP = np.vstack((UniqP,-UniqP))
     1477        phase = twopi*(np.inner(Uniq[:,:3],(dXdata+Xdata).T)+Phi[:,nxs])
     1478        sinp = np.sin(phase)
     1479        cosp = np.cos(phase)
     1480        occ = Mdata*Fdata/Uniq.shape[0]
     1481        biso = -SQfactor*Uisodata[:,nxs]
     1482        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[0],axis=1).T    #ops x atoms
     1483        HbH = -np.sum(UniqP[:,nxs,:3]*np.inner(UniqP[:,:3],bij),axis=-1)  #ops x atoms
     1484        Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in UniqP]) #atoms x 3x3
     1485        Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])                     #atoms x 6
     1486        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)     #ops x atoms
     1487        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[0]  #ops x atoms
     1488        fot = (FF+FP-Bab)*Tcorr     #ops x atoms
     1489        fotp = FPP*Tcorr            #ops x atoms
     1490        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms
     1491        dGdf,dGdx,dGdu,dGdz = G2mth.ModulationDerv(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt)
     1492        # GfpuA is 2 x ops x atoms
     1493        # derivs are: ops x atoms x waves x 2,6,12, or 5 parms as [real,imag] parts
     1494        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) # array(2,nEqv,nAtoms)
     1495        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])  #or array(2,nEqv,nAtoms)
     1496        fag = fa*GfpuA[0]-fb*GfpuA[1]
     1497        fbg = fb*GfpuA[0]+fa*GfpuA[1]
     1498       
     1499        fas = np.sum(np.sum(fag,axis=1),axis=1)     # 2 x twin
     1500        fbs = np.sum(np.sum(fbg,axis=1),axis=1)
     1501        fax = np.array([-fot*sinp,-fotp*cosp])   #positions; 2 x ops x atoms
     1502        fbx = np.array([fot*cosp,-fotp*sinp])
     1503        fax = fax*GfpuA[0]-fbx*GfpuA[1]
     1504        fbx = fbx*GfpuA[0]+fax*GfpuA[1]
     1505        #sum below is over Uniq
     1506        dfadfr = np.sum(fag/occ,axis=1)        #Fdata != 0 ever avoids /0. problem
     1507        dfbdfr = np.sum(fbg/occ,axis=1)        #Fdata != 0 avoids /0. problem
     1508        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
     1509        dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)
     1510        dfadui = np.sum(-SQfactor*fag,axis=1)
     1511        dfbdui = np.sum(-SQfactor*fbg,axis=1)
     1512        dfadx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fax,-2,-1)[:,:,:,nxs],axis=-2)  #2 x nAtom x 3xyz; sum nOps
     1513        dfbdx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fbx,-2,-1)[:,:,:,nxs],axis=-2)           
     1514        dfadua = np.sum(-Hij*np.swapaxes(fag,-2,-1)[:,:,:,nxs],axis=-2)             #2 x nAtom x 6Uij; sum nOps
     1515        dfbdua = np.sum(-Hij*np.swapaxes(fbg,-2,-1)[:,:,:,nxs],axis=-2)         #these are correct also for twins above
     1516        # array(2,nAtom,nWave,2) & array(2,nAtom,nWave,6) & array(2,nAtom,nWave,12); sum on nOps
     1517        dfadGf = np.sum(fa[:,:,:,nxs,nxs]*dGdf[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdf[1][nxs,:,:,:,:],axis=1)
     1518        dfbdGf = np.sum(fb[:,:,:,nxs,nxs]*dGdf[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdf[1][nxs,:,:,:,:],axis=1)
     1519        dfadGx = np.sum(fa[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1)
     1520        dfbdGx = np.sum(fb[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1)
     1521        dfadGz = np.sum(fa[:,:,0,nxs,nxs]*dGdz[0][nxs,:,:,:]-fb[:,:,0,nxs,nxs]*dGdz[1][nxs,:,:,:],axis=1)
     1522        dfbdGz = np.sum(fb[:,:,0,nxs,nxs]*dGdz[0][nxs,:,:,:]+fa[:,:,0,nxs,nxs]*dGdz[1][nxs,:,:,:],axis=1)
     1523        dfadGu = np.sum(fa[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1)
     1524        dfbdGu = np.sum(fb[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1)   
     1525#        GSASIIpath.IPyBreak()
     1526        if not SGData['SGInv']:   #Flack derivative
     1527            dfadfl = np.sum(-FPP*Tcorr*sinp)
     1528            dfbdfl = np.sum(FPP*Tcorr*cosp)
     1529        else:
     1530            dfadfl = 1.0
     1531            dfbdfl = 1.0
     1532        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
     1533        SA = fas[0]+fas[1]      #float = A+A' (might be array[nTwin])
     1534        SB = fbs[0]+fbs[1]      #float = B+B' (might be array[nTwin])
     1535        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro?
     1536            dFdfl[iref] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
     1537            dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)+   \
     1538                2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
     1539            dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])+  \
     1540                2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
     1541            dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])+   \
     1542                2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
     1543            dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])+   \
     1544                2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
     1545            dFdGf[iref] = 2.*(fas[0]*dfadGf[0]+fas[1]*dfadGf[1])+  \
     1546                2.*(fbs[0]*dfbdGf[0]+fbs[1]*dfbdGf[1])
     1547            dFdGx[iref] = 2.*(fas[0]*dfadGx[0]+fas[1]*dfadGx[1])+  \
     1548                2.*(fbs[0]*dfbdGx[0]-fbs[1]*dfbdGx[1])
     1549            dFdGz[iref] = 2.*(fas[0]*dfadGz[0]+fas[1]*dfadGz[1])+  \
     1550                2.*(fbs[0]*dfbdGz[0]+fbs[1]*dfbdGz[1])
     1551            dFdGu[iref] = 2.*(fas[0]*dfadGu[0]+fas[1]*dfadGu[1])+  \
     1552                2.*(fbs[0]*dfbdGu[0]+fbs[1]*dfbdGu[1])
     1553        else:                       #OK, I think
     1554            dFdfr[iref] = 2.*(SA*dfadfr[0]+SA*dfadfr[1]+SB*dfbdfr[0]+SB*dfbdfr[1])*Mdata/len(Uniq) #array(nRef,nAtom)
     1555            dFdx[iref] = 2.*(SA*dfadx[0]+SA*dfadx[1]+SB*dfbdx[0]+SB*dfbdx[1])    #array(nRef,nAtom,3)
     1556            dFdui[iref] = 2.*(SA*dfadui[0]+SA*dfadui[1]+SB*dfbdui[0]+SB*dfbdui[1])   #array(nRef,nAtom)
     1557            dFdua[iref] = 2.*(SA*dfadua[0]+SA*dfadua[1]+SB*dfbdua[0]+SB*dfbdua[1])    #array(nRef,nAtom,6)
     1558            dFdfl[iref] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
     1559                           
     1560            dFdGf[iref] = 2.*(SA*dfadGf[0]+SB*dfbdGf[1])      #array(nRef,natom,nwave,2)
     1561            dFdGx[iref] = 2.*(SA*dfadGx[0]+SB*dfbdGx[1])      #array(nRef,natom,nwave,6)
     1562            dFdGz[iref] = 2.*(SA*dfadGz[0]+SB*dfbdGz[1])      #array(nRef,natom,5)
     1563            dFdGu[iref] = 2.*(SA*dfadGu[0]+SB*dfbdGu[1])      #array(nRef,natom,nwave,12)
     1564#            GSASIIpath.IPyBreak()
     1565        dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
     1566            2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
     1567        #loop over atoms - each dict entry is list of derivatives for all the reflections
     1568        if not iref%100 :
     1569            print ' %d derivative time %.4f\r'%(iref,time.time()-time0),
     1570    for i in range(len(Mdata)):     #loop over atoms
     1571        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
     1572        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
     1573        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
     1574        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
     1575        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
     1576        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
     1577        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
     1578        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
     1579        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
     1580        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
     1581        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
     1582        for j in range(FSSdata.shape[1]):        #loop over waves Fzero & Fwid?
     1583            dFdvDict[pfx+'Fsin:'+str(i)+':'+str(j)] = dFdGf.T[0][j][i]
     1584            dFdvDict[pfx+'Fcos:'+str(i)+':'+str(j)] = dFdGf.T[1][j][i]
     1585        nx = 0
     1586        if waveTypes[i] in ['Block','ZigZag']:
     1587            nx = 1
     1588            dFdvDict[pfx+'Tmin:'+str(i)+':0'] = dFdGz.T[0][i]   #ZigZag/Block waves (if any)
     1589            dFdvDict[pfx+'Tmax:'+str(i)+':0'] = dFdGz.T[1][i]
     1590            dFdvDict[pfx+'Xmax:'+str(i)+':0'] = dFdGz.T[2][i]
     1591            dFdvDict[pfx+'Ymax:'+str(i)+':0'] = dFdGz.T[3][i]
     1592            dFdvDict[pfx+'Zmax:'+str(i)+':0'] = dFdGz.T[4][i]
     1593        for j in range(XSSdata.shape[1]-nx):       #loop over waves
     1594            dFdvDict[pfx+'Xsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[0][j][i]
     1595            dFdvDict[pfx+'Ysin:'+str(i)+':'+str(j+nx)] = dFdGx.T[1][j][i]
     1596            dFdvDict[pfx+'Zsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[2][j][i]
     1597            dFdvDict[pfx+'Xcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[3][j][i]
     1598            dFdvDict[pfx+'Ycos:'+str(i)+':'+str(j+nx)] = dFdGx.T[4][j][i]
     1599            dFdvDict[pfx+'Zcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[5][j][i]
     1600        for j in range(USSdata.shape[1]):       #loop over waves
     1601            dFdvDict[pfx+'U11sin:'+str(i)+':'+str(j)] = dFdGu.T[0][j][i]
     1602            dFdvDict[pfx+'U22sin:'+str(i)+':'+str(j)] = dFdGu.T[1][j][i]
     1603            dFdvDict[pfx+'U33sin:'+str(i)+':'+str(j)] = dFdGu.T[2][j][i]
     1604            dFdvDict[pfx+'U12sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[3][j][i]
     1605            dFdvDict[pfx+'U13sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[4][j][i]
     1606            dFdvDict[pfx+'U23sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[5][j][i]
     1607            dFdvDict[pfx+'U11cos:'+str(i)+':'+str(j)] = dFdGu.T[6][j][i]
     1608            dFdvDict[pfx+'U22cos:'+str(i)+':'+str(j)] = dFdGu.T[7][j][i]
     1609            dFdvDict[pfx+'U33cos:'+str(i)+':'+str(j)] = dFdGu.T[8][j][i]
     1610            dFdvDict[pfx+'U12cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[9][j][i]
     1611            dFdvDict[pfx+'U13cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[10][j][i]
     1612            dFdvDict[pfx+'U23cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[11][j][i]
     1613           
     1614#        GSASIIpath.IPyBreak()
     1615    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
     1616    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
     1617    return dFdvDict
     1618   
     1619def SStructureFactorDervTw(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
    13201620    'Needs a doc string'
    13211621    phfx = pfx.split(':')[0]+hfx
    … …  
    13611661        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))
    13621662    dFdvDict = {}
    1363     if nTwin > 1:
    1364         dFdfr = np.zeros((nRef,nTwin,mSize))
    1365         dFdx = np.zeros((nRef,nTwin,mSize,3))
    1366         dFdui = np.zeros((nRef,nTwin,mSize))
    1367         dFdua = np.zeros((nRef,nTwin,mSize,6))
    1368         dFdbab = np.zeros((nRef,nTwin,2))
    1369         dFdfl = np.zeros((nRef,nTwin))
    1370         dFdtw = np.zeros((nRef,nTwin))
    1371         dFdGf = np.zeros((nRef,nTwin,mSize,FSSdata.shape[1]))
    1372         dFdGx = np.zeros((nRef,nTwin,mSize,XSSdata.shape[1],3))
    1373         dFdGz = np.zeros((nRef,nTwin,mSize,5))
    1374         dFdGu = np.zeros((nRef,nTwin,mSize,USSdata.shape[1],6))
    1375     else:
    1376         dFdfr = np.zeros((nRef,mSize))
    1377         dFdx = np.zeros((nRef,mSize,3))
    1378         dFdui = np.zeros((nRef,mSize))
    1379         dFdua = np.zeros((nRef,mSize,6))
    1380         dFdbab = np.zeros((nRef,2))
    1381         dFdfl = np.zeros((nRef))
    1382         dFdtw = np.zeros((nRef))
    1383         dFdGf = np.zeros((nRef,mSize,FSSdata.shape[1],2))
    1384         dFdGx = np.zeros((nRef,mSize,XSSdata.shape[1],6))
    1385         dFdGz = np.zeros((nRef,mSize,5))
    1386         dFdGu = np.zeros((nRef,mSize,USSdata.shape[1],12))
     1663    dFdfr = np.zeros((nRef,nTwin,mSize))
     1664    dFdx = np.zeros((nRef,nTwin,mSize,3))
     1665    dFdui = np.zeros((nRef,nTwin,mSize))
     1666    dFdua = np.zeros((nRef,nTwin,mSize,6))
     1667    dFdbab = np.zeros((nRef,nTwin,2))
     1668    dFdfl = np.zeros((nRef,nTwin))
     1669    dFdtw = np.zeros((nRef,nTwin))
     1670    dFdGf = np.zeros((nRef,nTwin,mSize,FSSdata.shape[1]))
     1671    dFdGx = np.zeros((nRef,nTwin,mSize,XSSdata.shape[1],3))
     1672    dFdGz = np.zeros((nRef,nTwin,mSize,5))
     1673    dFdGu = np.zeros((nRef,nTwin,mSize,USSdata.shape[1],6))
    13871674    Flack = 1.0
    13881675    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
    … …  
    13951682        H = np.array(refl[:4])
    13961683        HP = H[:3]+modQ*H[3:]            #projected hklm to hkl
    1397 #        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(4,nTwins) or (4)
    1398 #        TwMask = np.any(H,axis=-1)
    1399 #        if TwinLaw.shape[0] > 1 and TwDict:
    1400 #            if iref in TwDict:
    1401 #                for i in TwDict[iref]:
    1402 #                    for n in range(NTL):
    1403 #                        H[i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
    1404 #            TwMask = np.any(H,axis=-1)
     1684        H = np.inner(H.T,TwinLaw)   #maybe array(4,nTwins) or (4)
     1685        TwMask = np.any(H,axis=-1)
     1686        if TwinLaw.shape[0] > 1 and TwDict:
     1687            if iref in TwDict:
     1688                for i in TwDict[iref]:
     1689                    for n in range(NTL):
     1690                        H[i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
     1691            TwMask = np.any(H,axis=-1)
    14051692        SQ = 1./(2.*refl[4+im])**2             # or (sin(theta)/lambda)**2
    14061693        SQfactor = 8.0*SQ*np.pi**2
    … …  
    14511738        dfadui = np.sum(-SQfactor*fag,axis=1)
    14521739        dfbdui = np.sum(-SQfactor*fbg,axis=1)
    1453         if nTwin > 1:
    1454             dfadx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
    1455             dfbdx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])           
    1456             dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fag,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
    1457             dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fbg,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
    1458             # array(2,nTwin,nAtom,3) & array(2,nTwin,nAtom,6) & array(2,nTwin,nAtom,12)
    1459             dfadGf = np.sum(fa[:,it,:,:,nxs,nxs]*dGdf[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdf[1][nxs,nxs,:,:,:,:],axis=1)
    1460             dfbdGf = np.sum(fb[:,it,:,:,nxs,nxs]*dGdf[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdf[1][nxs,nxs,:,:,:,:],axis=1)
    1461             dfadGx = np.sum(fa[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1)
    1462             dfbdGx = np.sum(fb[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1)
    1463             dfadGz = np.sum(fa[:,it,:,0,nxs,nxs]*dGdz[0][nxs,nxs,:,:,:]-fb[:,it,:,0,nxs,nxs]*dGdz[1][nxs,nxs,:,:,:],axis=1)
    1464             dfbdGz = np.sum(fb[:,it,:,0,nxs,nxs]*dGdz[0][nxs,nxs,:,:,:]+fa[:,it,:,0,nxs,nxs]*dGdz[1][nxs,nxs,:,:,:],axis=1)
    1465             dfadGu = np.sum(fa[:,it,:,:,nxs,nxs]*dGdu[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdu[1][nxs,nxs,:,:,:,:],axis=1)
    1466             dfbdGu = np.sum(fb[:,it,:,:,nxs,nxs]*dGdu[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdu[1][nxs,nxs,:,:,:,:],axis=1)
    1467         else:
    1468             dfadx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fax,-2,-1)[:,:,:,nxs],axis=-2)  #2 x nAtom x 3xyz; sum nOps
    1469             dfbdx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fbx,-2,-1)[:,:,:,nxs],axis=-2)           
    1470             dfadua = np.sum(-Hij*np.swapaxes(fag,-2,-1)[:,:,:,nxs],axis=-2)             #2 x nAtom x 6Uij; sum nOps
    1471             dfbdua = np.sum(-Hij*np.swapaxes(fbg,-2,-1)[:,:,:,nxs],axis=-2)         #these are correct also for twins above
    1472             # array(2,nAtom,nWave,2) & array(2,nAtom,nWave,6) & array(2,nAtom,nWave,12); sum on nOps
    1473             dfadGf = np.sum(fa[:,:,:,nxs,nxs]*dGdf[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdf[1][nxs,:,:,:,:],axis=1)
    1474             dfbdGf = np.sum(fb[:,:,:,nxs,nxs]*dGdf[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdf[1][nxs,:,:,:,:],axis=1)
    1475             dfadGx = np.sum(fa[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1)
    1476             dfbdGx = np.sum(fb[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1)
    1477             dfadGz = np.sum(fa[:,:,0,nxs,nxs]*dGdz[0][nxs,:,:,:]-fb[:,:,0,nxs,nxs]*dGdz[1][nxs,:,:,:],axis=1)
    1478             dfbdGz = np.sum(fb[:,:,0,nxs,nxs]*dGdz[0][nxs,:,:,:]+fa[:,:,0,nxs,nxs]*dGdz[1][nxs,:,:,:],axis=1)
    1479             dfadGu = np.sum(fa[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1)
    1480             dfbdGu = np.sum(fb[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1)   
     1740        dfadx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
     1741        dfbdx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])           
     1742        dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fag,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
     1743        dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fbg,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
     1744        # array(2,nTwin,nAtom,3) & array(2,nTwin,nAtom,6) & array(2,nTwin,nAtom,12)
     1745        dfadGf = np.sum(fa[:,it,:,:,nxs,nxs]*dGdf[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdf[1][nxs,nxs,:,:,:,:],axis=1)
     1746        dfbdGf = np.sum(fb[:,it,:,:,nxs,nxs]*dGdf[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdf[1][nxs,nxs,:,:,:,:],axis=1)
     1747        dfadGx = np.sum(fa[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1)
     1748        dfbdGx = np.sum(fb[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1)
     1749        dfadGz = np.sum(fa[:,it,:,0,nxs,nxs]*dGdz[0][nxs,nxs,:,:,:]-fb[:,it,:,0,nxs,nxs]*dGdz[1][nxs,nxs,:,:,:],axis=1)
     1750        dfbdGz = np.sum(fb[:,it,:,0,nxs,nxs]*dGdz[0][nxs,nxs,:,:,:]+fa[:,it,:,0,nxs,nxs]*dGdz[1][nxs,nxs,:,:,:],axis=1)
     1751        dfadGu = np.sum(fa[:,it,:,:,nxs,nxs]*dGdu[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdu[1][nxs,nxs,:,:,:,:],axis=1)
     1752        dfbdGu = np.sum(fb[:,it,:,:,nxs,nxs]*dGdu[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdu[1][nxs,nxs,:,:,:,:],axis=1)
    14811753#        GSASIIpath.IPyBreak()
    14821754        if not SGData['SGInv'] and len(TwinLaw) == 1:   #Flack derivative
    … …  
    14891761        SA = fas[0]+fas[1]      #float = A+A' (might be array[nTwin])
    14901762        SB = fbs[0]+fbs[1]      #float = B+B' (might be array[nTwin])
    1491         if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro
    1492             dFdfl[iref] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
    1493             dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)+   \
    1494                 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
    1495             dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])+  \
    1496                 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
    1497             dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])+   \
    1498                 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
    1499             dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])+   \
    1500                 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
    1501             dFdGf[iref] = 2.*(fas[0]*dfadGf[0]+fas[1]*dfadGf[1])+  \
    1502                 2.*(fbs[0]*dfbdGf[0]+fbs[1]*dfbdGf[1])
    1503             dFdGx[iref] = 2.*(fas[0]*dfadGx[0]+fas[1]*dfadGx[1])+  \
    1504                 2.*(fbs[0]*dfbdGx[0]-fbs[1]*dfbdGx[1])
    1505             dFdGz[iref] = 2.*(fas[0]*dfadGz[0]+fas[1]*dfadGz[1])+  \
    1506                 2.*(fbs[0]*dfbdGz[0]+fbs[1]*dfbdGz[1])
    1507             dFdGu[iref] = 2.*(fas[0]*dfadGu[0]+fas[1]*dfadGu[1])+  \
    1508                 2.*(fbs[0]*dfbdGu[0]+fbs[1]*dfbdGu[1])
    1509         else:
    1510             if nTwin > 1:
    1511                 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)]
    1512                 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)]
    1513                 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)]
    1514                 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)]
    1515                 dFdtw[iref] = np.sum(TwMask*fas,axis=0)**2+np.sum(TwMask*fbs,axis=0)**2
     1763        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)]
     1764        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)]
     1765        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)]
     1766        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)]
     1767        dFdtw[iref] = np.sum(TwMask*fas,axis=0)**2+np.sum(TwMask*fbs,axis=0)**2
    15161768
    1517                 dFdGf[iref] = [2.*TwMask[it]*(SA[it]*dfadGf[1]+SB[it]*dfbdGf[1]) for it in range(nTwin)]
    1518                 dFdGx[iref] = [2.*TwMask[it]*(SA[it]*dfadGx[1]+SB[it]*dfbdGx[1]) for it in range(nTwin)]
    1519                 dFdGz[iref] = [2.*TwMask[it]*(SA[it]*dfadGz[1]+SB[it]*dfbdGz[1]) for it in range(nTwin)]
    1520                 dFdGu[iref] = [2.*TwMask[it]*(SA[it]*dfadGu[1]+SB[it]*dfbdGu[1]) for it in range(nTwin)]               
    1521             else:   #these are good for no twin single crystals
    1522                 dFdfr[iref] = 2.*(SA*dfadfr[0]+SA*dfadfr[1]+SB*dfbdfr[0]+SB*dfbdfr[1])*Mdata/len(Uniq) #array(nRef,nAtom)
    1523                 dFdx[iref] = 2.*(SA*dfadx[0]+SA*dfadx[1]+SB*dfbdx[0]+SB*dfbdx[1])    #array(nRef,nAtom,3)
    1524                 dFdui[iref] = 2.*(SA*dfadui[0]+SA*dfadui[1]+SB*dfbdui[0]+SB*dfbdui[1])   #array(nRef,nAtom)
    1525                 dFdua[iref] = 2.*(SA*dfadua[0]+SA*dfadua[1]+SB*dfbdua[0]+SB*dfbdua[1])    #array(nRef,nAtom,6)
    1526                 dFdfl[iref] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
    1527                                
    1528                 dFdGf[iref] = 2.*(SA*dfadGf[0]+SB*dfbdGf[1])      #array(nRef,natom,nwave,2)
    1529                 dFdGx[iref] = 2.*(SA*dfadGx[0]+SB*dfbdGx[1])      #array(nRef,natom,nwave,6)
    1530                 dFdGz[iref] = 2.*(SA*dfadGz[0]+SB*dfbdGz[1])      #array(nRef,natom,5)
    1531                 dFdGu[iref] = 2.*(SA*dfadGu[0]+SB*dfbdGu[1])      #array(nRef,natom,nwave,12)
     1769        dFdGf[iref] = [2.*TwMask[it]*(SA[it]*dfadGf[1]+SB[it]*dfbdGf[1]) for it in range(nTwin)]
     1770        dFdGx[iref] = [2.*TwMask[it]*(SA[it]*dfadGx[1]+SB[it]*dfbdGx[1]) for it in range(nTwin)]
     1771        dFdGz[iref] = [2.*TwMask[it]*(SA[it]*dfadGz[1]+SB[it]*dfbdGz[1]) for it in range(nTwin)]
     1772        dFdGu[iref] = [2.*TwMask[it]*(SA[it]*dfadGu[1]+SB[it]*dfbdGu[1]) for it in range(nTwin)]               
    15321773#            GSASIIpath.IPyBreak()
    15331774        dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
    … …  
    29143155    if parmDict[phfx+'Scale'] < 0.:
    29153156        parmDict[phfx+'Scale'] = .001
    2916     if im:
    2917         dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
     3157    if im: # split to nontwin/twin versions
     3158        if len(TwinLaw) > 1:
     3159            dFdvDict = SStructureFactorDervTw(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
     3160        else:
     3161            dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
    29183162    else:
    29193163        if len(TwinLaw) > 1:
Note: See TracChangeset for help on using the changeset viewer.