Changeset 2495


Ignore:
Timestamp:
Oct 18, 2016 3:58:13 PM (5 years ago)
Author:
vondreele
Message:

some cleanup of Data tab
fix (!) crash after Refine for Costraints & Rigid Bodies - DELETED windows problem
remove Babinet & Flack from magnetic structures
remove magnetic from StructureFactorDeriv2
make new StructureFactorDerivMag? & work in that

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIconstrGUI.py

    r2491 r2495  
    11001100            G2frame.dataFrame.ConstraintEdit.Enable(G2gd.wxID_EQUIVALANCEATOMS,True)
    11011101#            G2frame.dataFrame.ConstraintEdit.Enable(G2gd.wxID_ADDRIDING,True)
    1102             UpdateConstraintPanel(PhaseConstr,'Phase')
     1102            if not 'DELETED' in str(PhaseConstr):
     1103                UpdateConstraintPanel(PhaseConstr,'Phase')
    11031104        elif text == 'Global constraints':
    11041105            G2frame.Page = [page,'glb']
     
    14771478        AtInfo = data['Vector']['AtInfo']
    14781479        refChoice = {}
     1480        if 'DELETED' in str(Status):
     1481            return
    14791482        SetStatusLine(' You may use e.g. "c60" or "s60" for a vector entry')
    14801483        def rbNameSizer(rbId,rbData):
  • trunk/GSASIIddataGUI.py

    r2493 r2495  
    419419        Axis.Bind(wx.EVT_TEXT_ENTER,OnAxis)
    420420        Axis.Bind(wx.EVT_KILL_FOCUS,OnAxis)
    421         uniSizer.Add(Axis,0,WACV)
     421        uniSizer.Add(Axis,0,WACV|wx.LEFT,5)
    422422        return uniSizer
    423423       
     
    432432            Indx[sizeRef.GetId()] = [G2frame.hist,id]
    433433            sizeRef.Bind(wx.EVT_CHECKBOX, OnRef)
    434             dataSizer.Add(sizeRef,0,WACV)
     434            dataSizer.Add(sizeRef,0,WACV|wx.LEFT,5)
    435435            sizeVal = wx.TextCtrl(DData,wx.ID_ANY,fmt%(val),style=wx.TE_PROCESS_ENTER)
    436436            Indx[sizeVal.GetId()] = [G2frame.hist,id]
     
    524524        poRef.SetValue(POData[2])
    525525        poRef.Bind(wx.EVT_CHECKBOX,OnPORef)
    526         poSizer.Add(poRef,0,WACV)
     526        poSizer.Add(poRef,0,WACV|wx.LEFT,5)
    527527        poVal = wx.TextCtrl(DData,wx.ID_ANY,
    528528            '%.3f'%(POData[1]),style=wx.TE_PROCESS_ENTER)
     
    595595        shPenalty.Add(wx.StaticText(DData,wx.ID_ANY,' Negative MRD penalty list: '),0,WACV)
    596596        shPenalty.Add(wx.ComboBox(DData,value=POData[6][0],choices=POData[6],
    597             style=wx.CB_DROPDOWN),0,WACV)
     597            style=wx.CB_DROPDOWN),0,WACV|wx.LEFT,5)
    598598        hklList = wx.Button(DData,label='Select penalty hkls')
    599599        hklList.Bind(wx.EVT_BUTTON,OnHKLList)
     
    668668            Obj.SetValue("%.2f"%(UseList[G2frame.hist]['Extinction'][0]))
    669669
    670         extSizer = wx.BoxSizer(wx.VERTICAL)
    671670        if Type == 'HKLF':
     671            extSizer = wx.BoxSizer(wx.VERTICAL)
    672672            typeSizer = wx.BoxSizer(wx.HORIZONTAL)           
    673673            typeSizer.Add(wx.StaticText(DData,wx.ID_ANY,' Extinction type: '),0,WACV)
     
    727727                extSizer.Add(val2Sizer,0,WACV)
    728728        else:   #PWDR
     729            extSizer = wx.BoxSizer(wx.HORIZONTAL)
    729730            extRef = wx.CheckBox(DData,wx.ID_ANY,label=' Extinction: ')
    730731            extRef.SetValue(UseList[G2frame.hist]['Extinction'][1])
     
    876877        addtwin = wx.CheckBox(DData,wx.ID_ANY,label=' Add Twin Law')
    877878        addtwin.Bind(wx.EVT_CHECKBOX, OnAddTwin)
    878         topsizer.Add(addtwin,0,WACV)
     879        topsizer.Add(addtwin,0,WACV|wx.LEFT,5)
    879880        twinsizer.Add(topsizer)
    880881        Indx = {}
     
    890891                    TwVal = Twin[1][0]
    891892                if 'bool' not in str(type(Twin[0])):
    892                     matSizer.Add(wx.StaticText(DData,-1,' Twin Law: '),0,WACV)
     893                    matSizer.Add(wx.StaticText(DData,-1,' Twin Law: '),0,WACV|wx.LEFT,5)
    893894                    for im,Mat in enumerate(twinMat):
    894895                        mat = wx.TextCtrl(DData,wx.ID_ANY,'%3d %3d %3d'%(Mat[0],Mat[1],Mat[2]),
     
    902903                        matSizer.Add(mat,0,WACV|wx.LEFT,5)
    903904                else:
    904                     matSizer.Add(wx.StaticText(DData,-1,' Nonmerohedral twin component %d: '%(it)),0,WACV)
     905                    matSizer.Add(wx.StaticText(DData,-1,' Nonmerohedral twin component %d: '%(it)),0,WACV|wx.LEFT,5)
    905906                    if not SGData['SGInv']:
    906907                        twinv = wx.CheckBox(DData,wx.ID_ANY,label=' Use enantiomorph?')
     
    10871088            bottomSizer.Add(poSizer,0,WACV|wx.TOP|wx.BOTTOM,5)
    10881089            bottomSizer.Add(ExtSizer('PWDR'),0,WACV|wx.TOP|wx.BOTTOM,5)
    1089             bottomSizer.Add(BabSizer(),0,WACV|wx.BOTTOM,5)
     1090            if generalData['Type'] != 'magnetic':
     1091                bottomSizer.Add(BabSizer(),0,WACV|wx.BOTTOM,5)
    10901092        elif G2frame.hist[:4] == 'HKLF':
    10911093#patch
     
    11091111            useList.append(name)
    11101112    mainSizer = wx.BoxSizer(wx.VERTICAL)
    1111     mainSizer.Add(wx.StaticText(DData,wx.ID_ANY,' Histogram data for '+PhaseName+':'),0,WACV)
     1113    mainSizer.Add(wx.StaticText(DData,wx.ID_ANY,' Histogram data for '+PhaseName+':'),0,WACV|wx.LEFT,5)
    11121114    if G2frame.hist != '':
    11131115        topSizer = wx.FlexGridSizer(1,2,5,5)
  • trunk/GSASIIstrIO.py

    r2483 r2495  
    12931293            if Print:
    12941294                print >>pFile,'\n Phase name: ',General['Name']
    1295                 print >>pFile,135*'-'
     1295                print >>pFile,135*'='
    12961296                PrintFFtable(FFtable)
    12971297                PrintBLtable(BLtable)
     
    13371337            if Print:
    13381338                print >>pFile,'\n Phase name: ',General['Name']
    1339                 print >>pFile,135*'-'
     1339                print >>pFile,135*'='
    13401340                print >>pFile,''
    13411341                if len(SSGtext):    #if superstructure
     
    18841884    for phase in Phases:
    18851885        print >>pFile,' Result for phase: ',phase
     1886        print >>pFile,135*'='
    18861887        Phase = Phases[phase]
    18871888        General = Phase['General']
     
    22742275                            if hapData[item][5][i]:
    22752276                                hapVary.append(pfx+item+sfx)
    2276                 for bab in ['BabA','BabU']:
    2277                     hapDict[pfx+bab] = hapData['Babinet'][bab][0]
    2278                     if hapData['Babinet'][bab][1]:
    2279                         hapVary.append(pfx+bab)
     2277                if Phases[phase]['General']['Type'] != 'magnetic':
     2278                    for bab in ['BabA','BabU']:
     2279                        hapDict[pfx+bab] = hapData['Babinet'][bab][0]
     2280                        if hapData['Babinet'][bab][1]:
     2281                            hapVary.append(pfx+bab)
    22802282                               
    22812283                if Print:
    22822284                    print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
    2283                     print >>pFile,135*'-'
     2285                    print >>pFile,135*'='
    22842286                    print >>pFile,' Phase fraction  : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
    22852287                    print >>pFile,' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
     
    22942296                    PrintMuStrain(hapData['Mustrain'],SGData)
    22952297                    PrintHStrain(hapData['HStrain'],SGData)
    2296                     if hapData['Babinet']['BabA'][0]:
    2297                         PrintBabinet(hapData['Babinet'])                       
     2298                    if Phases[phase]['General']['Type'] != 'magnetic':
     2299                        if hapData['Babinet']['BabA'][0]:
     2300                            PrintBabinet(hapData['Babinet'])                       
    22982301                if resetRefList:
    22992302                    refList = []
     
    24172420                if Print:
    24182421                    print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
    2419                     print >>pFile,135*'-'
     2422                    print >>pFile,135*'='
    24202423                    print >>pFile,' Scale factor     : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
    24212424                    if extType != 'None':
     
    26802683                    if pfx+name in sigDict:
    26812684                        SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
    2682                 for name in ['BabA','BabU']:
    2683                     hapData['Babinet'][name][0] = parmDict[pfx+name]
    2684                     if pfx+name in sigDict:
    2685                         BabSig[pfx+name] = sigDict[pfx+name]               
     2685                if Phases[phase]['General']['Type'] != 'magnetic':
     2686                    for name in ['BabA','BabU']:
     2687                        hapData['Babinet'][name][0] = parmDict[pfx+name]
     2688                        if pfx+name in sigDict:
     2689                            BabSig[pfx+name] = sigDict[pfx+name]               
    26862690               
    26872691            elif 'HKLF' in histogram:
     
    27302734                    continue
    27312735                print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
    2732                 print >>pFile,130*'-'
     2736                print >>pFile,135*'='
    27332737                hapData = HistoPhase[histogram]
    27342738                hId = Histogram['hId']
     
    27562760                    PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
    27572761                    PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
    2758                     if len(BabSig):
    2759                         PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
     2762                    if Phases[phase]['General']['Type'] != 'magnetic':
     2763                        if len(BabSig):
     2764                            PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
    27602765                   
    27612766                elif 'HKLF' in histogram:
     
    29902995            if Print:
    29912996                print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
    2992                 print >>pFile,135*'-'
     2997                print >>pFile,135*'='
    29932998                Units = {'C':' deg','T':' msec'}
    29942999                units = Units[controlDict[pfx+'histType'][2]]
     
    32103215
    32113216            print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
    3212             print >>pFile,135*'-'
     3217            print >>pFile,135*'='
    32133218            print >>pFile,' PWDR histogram weight factor = '+'%.3f'%(Histogram['wtFactor'])
    32143219            print >>pFile,' Final refinement wR = %.2f%% on %d observations in this histogram'%(Histogram['Residuals']['wR'],Histogram['Residuals']['Nobs'])
  • trunk/GSASIIstrMath.py

    r2493 r2495  
    753753            FP = np.repeat(FP.T,len(SGT)*len(TwinLaw),axis=0)
    754754            FPP = np.repeat(FPP.T,len(SGT)*len(TwinLaw),axis=0)
    755         Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT)*len(TwinLaw))
    756755        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
    757756        Uniq = np.inner(H,SGMT)
     
    787786            fbms = np.sum(np.sum(fbm,axis=-1),axis=-1)                          #ditto
    788787#            GSASIIpath.IPyBreak()
    789         if 'T' in calcControls[hfx+'histType']: #fa,fb are 2 X blkSize X nTwin X nOps x nAtoms
    790             fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
    791             fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr])
     788            refl.T[9] = np.sum(fams**2,axis=0)+np.sum(fbms**2,axis=0)
     789            refl.T[7] = np.copy(refl.T[9])               
     790            refl.T[10] = 0.0    #atan2d(fbs[0],fas[0]) - what is phase for mag refl?
    792791        else:
    793             fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
    794             fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
    795         fas = np.sum(np.sum(fa,axis=-1),axis=-1)  #real 2 x blkSize x nTwin; sum over atoms & uniq hkl
    796         fbs = np.sum(np.sum(fb,axis=-1),axis=-1)  #imag
    797         if SGData['SGInv']: #centrosymmetric; B=0
    798             fbs[0] *= 0.
    799             fas[1] *= 0.
    800         if 'P' in calcControls[hfx+'histType']:     #PXC, PNC & PNT: F^2 = A[0]^2 + A[1]^2 + B[0]^2 + B[1]^2
    801             refl.T[9] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0) #add fam**2 & fbm**2 here   
    802             refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
    803             if 'N' in calcControls[hfx+'histType'] and parmDict[pfx+'isMag']:
    804                 refl.T[9] = np.sum(fams**2,axis=0)+np.sum(fbms**2,axis=0)
    805         else:                                       #HKLF: F^2 = (A[0]+A[1])^2 + (B[0]+B[1])^2
    806             if len(TwinLaw) > 1:
    807                 refl.T[9] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2   #FcT from primary twin element
    808                 refl.T[7] = np.sum(TwinFr*TwMask*np.sum(fas,axis=0)**2,axis=-1)+   \
    809                     np.sum(TwinFr*TwMask*np.sum(fbs,axis=0)**2,axis=-1)                        #Fc sum over twins
    810                 refl.T[10] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f" & use primary twin
    811             else:   # checked correct!!
    812                 refl.T[9] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2
    813                 refl.T[7] = np.copy(refl.T[9])               
     792            Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT)*len(TwinLaw))
     793            if 'T' in calcControls[hfx+'histType']: #fa,fb are 2 X blkSize X nTwin X nOps x nAtoms
     794                fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
     795                fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr])
     796            else:
     797                fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
     798                fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
     799            fas = np.sum(np.sum(fa,axis=-1),axis=-1)  #real 2 x blkSize x nTwin; sum over atoms & uniq hkl
     800            fbs = np.sum(np.sum(fb,axis=-1),axis=-1)  #imag
     801            if SGData['SGInv']: #centrosymmetric; B=0
     802                fbs[0] *= 0.
     803                fas[1] *= 0.
     804            if 'P' in calcControls[hfx+'histType']:     #PXC, PNC & PNT: F^2 = A[0]^2 + A[1]^2 + B[0]^2 + B[1]^2
     805                refl.T[9] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0) #add fam**2 & fbm**2 here   
    814806                refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
    815807                if 'N' in calcControls[hfx+'histType'] and parmDict[pfx+'isMag']:
    816808                    refl.T[9] = np.sum(fams**2,axis=0)+np.sum(fbms**2,axis=0)
     809            else:                                       #HKLF: F^2 = (A[0]+A[1])^2 + (B[0]+B[1])^2
     810                if len(TwinLaw) > 1:
     811                    refl.T[9] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2   #FcT from primary twin element
     812                    refl.T[7] = np.sum(TwinFr*TwMask*np.sum(fas,axis=0)**2,axis=-1)+   \
     813                        np.sum(TwinFr*TwMask*np.sum(fbs,axis=0)**2,axis=-1)                        #Fc sum over twins
     814                    refl.T[10] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f" & use primary twin
     815                else:   # checked correct!!
     816                    refl.T[9] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2
     817                    refl.T[7] = np.copy(refl.T[9])               
     818                    refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
    817819#        GSASIIpath.IPyBreak()
    818820#                refl.T[10] = atan2d(np.sum(fbs,axis=0),np.sum(fas,axis=0)) #include f' & f"
     
    10581060        GetAtomFXU(pfx,calcControls,parmDict)
    10591061    mSize = len(Mdata)
    1060     if parmDict[pfx+'isMag']:
    1061         Mag = np.sqrt(np.sum(Gdata**2,axis=0))      #magnitude of moments for uniq atoms
    1062         Gdata = np.where(Mag>0.,Gdata/Mag,0.)       #normalze mag. moments
    1063         Gdata = np.inner(Bmat,Gdata.T)              #convert to crystal space
    1064         Gdata = np.inner(Gdata.T,SGMT).T            #apply sym. ops.
    1065         if SGData['SGInv']:
    1066             Gdata = np.hstack((Gdata,-Gdata))       #inversion if any
    1067         Gdata = np.hstack([Gdata for icen in range(Ncen)])        #dup over cell centering
    1068         Gdata = SGData['MagMom'][nxs,:,nxs]*Gdata   #flip vectors according to spin flip
    1069         Gdata = np.inner(Amat,Gdata.T)              #convert back to cart. space MXYZ, Natoms, NOps*Inv*Ncen
    1070         Gdata = np.swapaxes(Gdata,1,2)              # put Natoms last
    1071 #        GSASIIpath.IPyBreak()
    1072         Mag = np.tile(Mag[:,nxs],len(SGMT)*Ncen).T  #Mag same length as Gdata
    1073         if SGData['SGInv']:
    1074             Mag = np.repeat(Mag,2,axis=0)
    1075         dFdMx = np.zeros((nRef,mSize,3))
    10761062    FF = np.zeros(len(Tdata))
    10771063    if 'NC' in calcControls[hfx+'histType']:
     
    10851071    dFdfr = np.zeros((nRef,mSize))
    10861072    dFdx = np.zeros((nRef,mSize,3))
    1087     dFdmx = np.zeros((nRef,mSize,3))
    10881073    dFdui = np.zeros((nRef,mSize))
    10891074    dFdua = np.zeros((nRef,mSize,6))
     
    11341119        else:
    11351120            fotp = FPP*Tcorr     
    1136         if 'N' in calcControls[hfx+'histType'] and parmDict[pfx+'isMag']:
    1137             MF = refDict['FF']['MF'][iBeg:iFin].T[Tindx].T   #Nref,Natm
    1138             TMcorr = 0.539*(np.reshape(Tiso,Tuij.shape)*Tuij)[:,0,:]*Mdata*MF/(2*Nops*Ncen)     #Nref,Natm                                  #Nref,Natm
    1139             if SGData['SGInv']:
    1140                 mphase = np.hstack((phase,-phase))
    1141                 TMcorr = TMcorr/2.
    1142             else:
    1143                 mphase = phase
    1144             mphase = np.array([mphase+twopi*np.inner(cen,H)[:,nxs,nxs] for cen in SGData['SGCen']])
    1145             mphase = np.concatenate(mphase,axis=1)              #Nref,Nop,Natm
    1146             sinm = np.sin(mphase)                               #ditto - match magstrfc.for
    1147             cosm = np.cos(mphase)                               #ditto
    1148             HM = np.inner(Bmat.T,H)                             #put into cartesian space
    1149             HM = HM/np.sqrt(np.sum(HM**2,axis=0))               #Gdata = MAGS & HM = UVEC in magstrfc.for both OK
    1150             eDotK = np.sum(HM[:,:,nxs,nxs]*Gdata[:,nxs,:,:],axis=0)
    1151             Q = HM[:,:,nxs,nxs]*eDotK[nxs,:,:,:]-Gdata[:,nxs,:,:] #xyz,Nref,Nop,Natm = BPM in magstrfc.for OK
    1152             fam = Q*TMcorr[nxs,:,nxs,:]*cosm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #ditto
    1153             fbm = Q*TMcorr[nxs,:,nxs,:]*sinm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #ditto
    1154             fams = np.sum(np.sum(fam,axis=-1),axis=-1)                          #xyz,Nref
    1155             fbms = np.sum(np.sum(fbm,axis=-1),axis=-1)                          #ditto
    1156             Hij = np.hstack((Hij,Hij))
    11571121#            GSASIIpath.IPyBreak()
    11581122        if 'T' in calcControls[hfx+'histType']:
     
    11671131        fbx = np.array([fot*cosp,-fotp*sinp])
    11681132        #sum below is over Uniq
    1169         if 'N' in calcControls[hfx+'histType'] and parmDict[pfx+'isMag']:
    1170             dfadfr = np.sum(fams/occ,axis=-2)        #array(refBlk,nAtom) Fdata != 0 avoids /0. problem
    1171             dfadba = np.sum(-cosm*TMcorr[:,nxs,:],axis=-2)    #array(refBlk,nAtom)
    1172             dfadx = np.sum(twopi*Uniq[:,nxs,:,:]*np.swapaxes(fax[0],-2,-1)[:,:,:,nxs],axis=-2)
    1173             dfadui = np.sum(-SQfactor[:,nxs,nxs]*fams,axis=-2) #array(Ops,refBlk,nAtoms)
    1174             dfadua = np.sum(-Hij[nxs,:,:,nxs,:]*fams[:,:,:,:,nxs],axis=-3)
    1175         else:
    1176             dfadfr = np.sum(fa/occ,axis=-2)        #array(2,refBlk,nAtom) Fdata != 0 avoids /0. problem
    1177             dfadba = np.sum(-cosp*Tcorr,axis=-2)  #array(refBlk,nAtom)
    1178             dfadx = np.sum(twopi*Uniq[nxs,:,nxs,:,:]*np.swapaxes(fax,-2,-1)[:,:,:,:,nxs],axis=-2)
    1179             dfadui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fa,axis=-2) #array(Ops,refBlk,nAtoms)
    1180             dfadua = np.sum(-Hij[nxs,:,nxs,:,:]*np.swapaxes(fa,-2,-1)[:,:,:,:,nxs],axis=-2)
     1133        dfadfr = np.sum(fa/occ,axis=-2)        #array(2,refBlk,nAtom) Fdata != 0 avoids /0. problem
     1134        dfadba = np.sum(-cosp*Tcorr,axis=-2)  #array(refBlk,nAtom)
     1135        dfadx = np.sum(twopi*Uniq[nxs,:,nxs,:,:]*np.swapaxes(fax,-2,-1)[:,:,:,:,nxs],axis=-2)
     1136        dfadui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fa,axis=-2) #array(Ops,refBlk,nAtoms)
     1137        dfadua = np.sum(-Hij[nxs,:,nxs,:,:]*np.swapaxes(fa,-2,-1)[:,:,:,:,nxs],axis=-2)
    11811138        # array(2,refBlk,nAtom,3) & array(2,refBlk,nAtom,6)
    11821139        if not SGData['SGInv']:
    1183             if 'N' in calcControls[hfx+'histType'] and parmDict[pfx+'isMag']:
    1184                 dfbdfr = np.sum(fbm/occ,axis=-2)        #array(refBlk,nAtom) Fdata != 0 avoids /0. problem
    1185                 dfbdba = np.sum(-sinm*TMcorr[:,nxs,:],axis=-2)    #array(refBlk,nAtom)
    1186                 dfbdx = np.sum(twopi*Uniq[:,nxs,:,:]*np.swapaxes(fbx[0],-2,-1)[:,:,:,nxs],axis=-2)
    1187                 dfbdui = np.sum(-SQfactor[:,nxs,nxs]*fbm,axis=-2) #array(Ops,refBlk,nAtoms)
    1188                 dfbdua = np.sum(-Hij[nxs,:,:,nxs,:]*fbm[:,:,:,:,nxs],axis=-3)
    1189             else:
    1190                 dfbdfr = np.sum(fb/occ,axis=-2)        #Fdata != 0 avoids /0. problem
    1191                 dfbdba = np.sum(-sinp*Tcorr,axis=-2)
    1192                 dfadfl = np.sum(np.sum(-fotp*sinp,axis=-1),axis=-1)
    1193                 dfbdfl = np.sum(np.sum(fotp*cosp,axis=-1),axis=-1)
    1194                 dfbdx = np.sum(twopi*Uniq[nxs,:,nxs,:,:]*np.swapaxes(fbx,-2,-1)[:,:,:,:,nxs],axis=-2)           
    1195                 dfbdui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fb,axis=-2)
    1196                 dfbdua = np.sum(-Hij[nxs,:,nxs,:,:]*np.swapaxes(fb,-2,-1)[:,:,:,:,nxs],axis=-2)
     1140            dfbdfr = np.sum(fb/occ,axis=-2)        #Fdata != 0 avoids /0. problem
     1141            dfbdba = np.sum(-sinp*Tcorr,axis=-2)
     1142            dfadfl = np.sum(np.sum(-fotp*sinp,axis=-1),axis=-1)
     1143            dfbdfl = np.sum(np.sum(fotp*cosp,axis=-1),axis=-1)
     1144            dfbdx = np.sum(twopi*Uniq[nxs,:,nxs,:,:]*np.swapaxes(fbx,-2,-1)[:,:,:,:,nxs],axis=-2)           
     1145            dfbdui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fb,axis=-2)
     1146            dfbdua = np.sum(-Hij[nxs,:,nxs,:,:]*np.swapaxes(fb,-2,-1)[:,:,:,:,nxs],axis=-2)
    11971147        else:
    11981148            dfbdfr = np.zeros_like(dfadfr)
     
    12071157        SB = fbs[0]+fbs[1]
    12081158        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro
    1209             if 'N' in calcControls[hfx+'histType'] and parmDict[pfx+'isMag']:
    1210                 dFdfr[iBeg:iFin] = 2.*(fams[:,nxs]*dfadfr+fbms[:,nxs]*dfbdfr)*Mdata/len(SGMT)
    1211                 dFdx[iBeg:iFin] = 2.*(fams[:,nxs,nxs]*dfadx+fbms[:,nxs,nxs]*dfbdx)
    1212                 dFdui[iBeg:iFin] = 2.*(fams[:,nxs]*dfadui+fbms[:,nxs]*dfbdui)
    1213                 dFdua[iBeg:iFin] = 2.*(fams[:,nxs,nxs]*dfadua+fbms[:,nxs,nxs]*dfbdua)
    1214             else:
    1215                 dFdfr[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs]*dfadfr+fbs[:,:,nxs]*dfbdfr,axis=0)*Mdata/len(SGMT)
    1216                 dFdx[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs,nxs]*dfadx+fbs[:,:,nxs,nxs]*dfbdx,axis=0)
    1217                 dFdui[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs]*dfadui+fbs[:,:,nxs]*dfbdui,axis=0)
    1218                 dFdua[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs,nxs]*dfadua+fbs[:,:,nxs,nxs]*dfbdua,axis=0)
     1159            dFdfr[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs]*dfadfr+fbs[:,:,nxs]*dfbdfr,axis=0)*Mdata/len(SGMT)
     1160            dFdx[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs,nxs]*dfadx+fbs[:,:,nxs,nxs]*dfbdx,axis=0)
     1161            dFdui[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs]*dfadui+fbs[:,:,nxs]*dfbdui,axis=0)
     1162            dFdua[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs,nxs]*dfadua+fbs[:,:,nxs,nxs]*dfbdua,axis=0)
    12191163        else:
    12201164            dFdfr[iBeg:iFin] = (2.*SA[:,nxs]*(dfadfr[0]+dfadfr[1])+2.*SB[:,nxs]*(dfbdfr[0]+dfbdfr[1]))*Mdata/len(SGMT)
     
    12441188    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
    12451189    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
     1190    return dFdvDict
     1191   
     1192def StructureFactorDervMag(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
     1193    '''Compute structure factor derivatives on blocks of reflections - for powders/nontwins only
     1194    faster than StructureFactorDerv - correct for powders/nontwins!!
     1195    input:
     1196   
     1197    :param dict refDict: where
     1198        'RefList' list where each ref = h,k,l,it,d,...
     1199        'FF' dict of form factors - filled in below
     1200    :param np.array G:      reciprocal metric tensor
     1201    :param str hfx:    histogram id string
     1202    :param str pfx:    phase id string
     1203    :param dict SGData: space group info. dictionary output from SpcGroup
     1204    :param dict calcControls:
     1205    :param dict parmDict:
     1206   
     1207    :returns: dict dFdvDict: dictionary of derivatives
     1208    '''
     1209    phfx = pfx.split(':')[0]+hfx
     1210    ast = np.sqrt(np.diag(G))
     1211    Mast = twopisq*np.multiply.outer(ast,ast)
     1212    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
     1213    SGT = np.array([ops[1] for ops in SGData['SGOps']])
     1214    Ncen = len(SGData['SGCen'])
     1215    Nops = len(SGMT)
     1216    MFtables = calcControls['MFtables']
     1217    Amat,Bmat = G2lat.Gmat2AB(G)
     1218    nRef = len(refDict['RefList'])
     1219    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
     1220        GetAtomFXU(pfx,calcControls,parmDict)
     1221    mSize = len(Mdata)
     1222    Mag = np.sqrt(np.sum(Gdata**2,axis=0))      #magnitude of moments for uniq atoms
     1223    Gdata = np.where(Mag>0.,Gdata/Mag,0.)       #normalze mag. moments
     1224    Gdata = np.inner(Bmat,Gdata.T)              #convert to crystal space
     1225    Gdata = np.inner(Gdata.T,SGMT).T            #apply sym. ops.
     1226    if SGData['SGInv']:
     1227        Gdata = np.hstack((Gdata,-Gdata))       #inversion if any
     1228    Gdata = np.hstack([Gdata for icen in range(Ncen)])        #dup over cell centering
     1229    Gdata = SGData['MagMom'][nxs,:,nxs]*Gdata   #flip vectors according to spin flip
     1230    Gdata = np.inner(Amat,Gdata.T)              #convert back to cart. space MXYZ, Natoms, NOps*Inv*Ncen
     1231    Gdata = np.swapaxes(Gdata,1,2)              # put Natoms last
     1232#    GSASIIpath.IPyBreak()
     1233    Mag = np.tile(Mag[:,nxs],len(SGMT)*Ncen).T  #Mag same length as Gdata
     1234    if SGData['SGInv']:
     1235        Mag = np.repeat(Mag,2,axis=0)
     1236    dFdMx = np.zeros((nRef,mSize,3))
     1237    Uij = np.array(G2lat.U6toUij(Uijdata))
     1238    bij = Mast*Uij.T
     1239    dFdvDict = {}
     1240    dFdfr = np.zeros((nRef,mSize))
     1241    dFdx = np.zeros((nRef,mSize,3))
     1242    dFdmx = np.zeros((nRef,mSize,3))
     1243    dFdui = np.zeros((nRef,mSize))
     1244    dFdua = np.zeros((nRef,mSize,6))
     1245    time0 = time.time()
     1246    nref = len(refDict['RefList'])/100   
     1247#reflection processing begins here - big arrays!
     1248    iBeg = 0
     1249    blkSize = 32       #no. of reflections in a block - optimized for speed
     1250    while iBeg < nRef:
     1251        iFin = min(iBeg+blkSize,nRef)
     1252        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
     1253        H = refl.T[:3].T
     1254        SQ = 1./(2.*refl.T[4])**2             # or (sin(theta)/lambda)**2
     1255        SQfactor = 8.0*SQ*np.pi**2
     1256#        GSASIIpath.IPyBreak()
     1257        Uniq = np.inner(H,SGMT)             # array(nSGOp,3)
     1258        Phi = np.inner(H,SGT)
     1259        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
     1260        occ = Mdata*Fdata/len(SGT)
     1261        biso = -SQfactor*Uisodata[:,nxs]
     1262        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT),axis=1).T
     1263        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
     1264        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
     1265        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])
     1266        Hij = np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(-1,len(SGT),6))
     1267        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
     1268        MF = refDict['FF']['MF'][iBeg:iFin].T[Tindx].T   #Nref,Natm
     1269        TMcorr = 0.539*(np.reshape(Tiso,Tuij.shape)*Tuij)[:,0,:]*Mdata*MF/(2*Nops*Ncen)     #Nref,Natm                                  #Nref,Natm
     1270        if SGData['SGInv']:
     1271            mphase = np.hstack((phase,-phase))
     1272            TMcorr = TMcorr/2.
     1273            Uniq = np.hstack((Uniq,-Uniq))      #Nref,Nops,hkl
     1274            Hij = np.hstack((Hij,Hij))
     1275        else:
     1276            mphase = phase
     1277        Hij = np.concatenate(np.array([Hij for cen in SGData['SGCen']]),axis=1)
     1278        Uniq = np.hstack([Uniq for cen in SGData['SGCen']])
     1279        mphase = np.array([mphase+twopi*np.inner(cen,H)[:,nxs,nxs] for cen in SGData['SGCen']])
     1280        mphase = np.concatenate(mphase,axis=1)              #Nref,Nop,Natm
     1281        sinm = np.sin(mphase)                               #ditto - match magstrfc.for
     1282        cosm = np.cos(mphase)                               #ditto
     1283        HM = np.inner(Bmat.T,H)                             #put into cartesian space
     1284        HM = HM/np.sqrt(np.sum(HM**2,axis=0))               #Gdata = MAGS & HM = UVEC in magstrfc.for both OK
     1285        eDotK = np.sum(HM[:,:,nxs,nxs]*Gdata[:,nxs,:,:],axis=0)
     1286        Q = HM[:,:,nxs,nxs]*eDotK[nxs,:,:,:]-Gdata[:,nxs,:,:] #xyz,Nref,Nop,Natm = BPM in magstrfc.for OK
     1287        fam = Q*TMcorr[nxs,:,nxs,:]*cosm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #ditto
     1288        fbm = Q*TMcorr[nxs,:,nxs,:]*sinm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #ditto
     1289        fams = np.sum(np.sum(fam,axis=-1),axis=-1)                          #xyz,Nref
     1290        fbms = np.sum(np.sum(fbm,axis=-1),axis=-1)                          #ditto
     1291        famx = Q*TMcorr[nxs,:,nxs,:]*Mag[nxs,nxs,:,:]*sinm[nxs,:,:,:]   #Mxyz,Nref,Nops,Natom
     1292        fbmx = Q*TMcorr[nxs,:,nxs,:]*Mag[nxs,nxs,:,:]*cosm[nxs,:,:,:]
     1293        #sum below is over Uniq
     1294        dfadfr = np.sum(fam/occ,axis=-2)        #array(mxyz,refBlk,nAtom) Fdata != 0 avoids /0. problem
     1295#        GSASIIpath.IPyBreak()
     1296#        dfadx = np.sum(twopi*Uniq.T[:,:,:,nxs]*famx,axis=2)
     1297#        dfadmx = twopi*Uniq.T[:,:,:,nxs]*famx) #?
     1298        dfadui = np.sum(-SQfactor[:,nxs,nxs]*fam,axis=2) #array(Ops,refBlk,nAtoms)
     1299        dfadua = np.sum(-Hij[nxs,:,:,nxs,:]*fam[:,:,:,:,nxs],axis=2)
     1300        # array(2,refBlk,nAtom,3) & array(2,refBlk,nAtom,6)
     1301        if not SGData['SGInv']:
     1302            dfbdfr = np.sum(fbm/occ,axis=-2)        #array(mxyz,refBlk,nAtom) Fdata != 0 avoids /0. problem
     1303#            dfbdx = np.sum(twopi*Uniq.T[:,:,:,nxs]*fbmx,axis=2)
     1304#            dfbdmx = twopi*Uniq[:,nxs,:,:]*np.transpose(fbmx,(1,2,3,0))
     1305            dfbdui = np.sum(-SQfactor[:,nxs,nxs]*fbm,axis=2) #array(Ops,refBlk,nAtoms)
     1306            dfbdua = np.sum(-Hij[nxs,:,:,nxs,:]*fbm[:,:,:,:,nxs],axis=2)
     1307        else:
     1308            dfbdfr = np.zeros_like(dfadfr)
     1309#            dfbdx = np.zeros_like(dfadx)
     1310#            dfbdmx = np.zeros_like(dfadmx)
     1311            dfbdui = np.zeros_like(dfadui)
     1312            dfbdua = np.zeros_like(dfadua)
     1313        dFdfr[iBeg:iFin] = np.sum(2.*(fams[:,:,nxs]*dfadfr+fbms[:,:,nxs]*dfbdfr)*Mdata/(2*Nops*Ncen),axis=0)
     1314#        dFdx[iBeg:iFin] = 2.*(fams[:,:,nxs]*dfadx+fbms[:,:,nxs]*dfbdx)
     1315#        dFdMx[iBeg:iFin] = np.sum(2.*(fams[:,nxs,nxs]*dfadmx+fbms[:,nxs,nxs]*dfbdmx),axis=0)
     1316        dFdui[iBeg:iFin] = np.sum(2.*(fams[:,:,nxs]*dfadui+fbms[:,:,nxs]*dfbdui),axis=0)
     1317        dFdua[iBeg:iFin] = np.sum(2.*(fams[:,:,nxs,nxs]*dfadua+fbms[:,:,nxs,nxs]*dfbdua),axis=0)
     1318#        GSASIIpath.IPyBreak()
     1319        iBeg += blkSize
     1320    print ' %d derivative time %.4f\r'%(nRef,time.time()-time0)
     1321        #loop over atoms - each dict entry is list of derivatives for all the reflections
     1322    for i in range(len(Mdata)):
     1323        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
     1324        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
     1325        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
     1326        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
     1327        dFdvDict[pfx+'dAMx:'+str(i)] = dFdMx.T[0][i]
     1328        dFdvDict[pfx+'dAMy:'+str(i)] = dFdMx.T[1][i]
     1329        dFdvDict[pfx+'dAMz:'+str(i)] = dFdMx.T[2][i]
     1330        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
     1331        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
     1332        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
     1333        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
     1334        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
     1335        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
     1336        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
    12461337    return dFdvDict
    12471338   
     
    35663657                dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
    35673658            else:
    3568                 dFdvDict = StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
     3659                if Phase['General']['Type'] == 'magnetic':
     3660                    dFdvDict = StructureFactorDervMag(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
     3661                else:
     3662                    dFdvDict = StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
    35693663#            print 'sf-derv time %.3fs'%(time.time()-time0)
    35703664            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
     
    38263920            dFdvDict = StructureFactorDervTw2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
    38273921        else:   #correct!!
    3828             dFdvDict = StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
     3922            if Phase['General']['Type'] == 'magnetic':  #is this going to work for single crystal mag data?
     3923                dFdvDict = StructureFactorDervMag(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
     3924            else:
     3925                dFdvDict = StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
    38293926    ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
    38303927    dMdvh = np.zeros((len(varylist),len(refDict['RefList'])))
Note: See TracChangeset for help on using the changeset viewer.