Changeset 1887


Ignore:
Timestamp:
Jun 12, 2015 3:13:05 PM (7 years ago)
Author:
vondreele
Message:

default new twin to inversion twin - usual case
fix twin results output
twin functions now work - derivatives work but incorrect.

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIddataGUI.py

    r1885 r1887  
    780780       
    781781        def OnAddTwin(event):
    782             twinMat = np.array([[1,0,0],[0,1,0],[0,0,1]])
    783             twinVal = [1.0,False]
     782            twinMat = np.array([[-1,0,0],[0,-1,0],[0,0,-1]])    #inversion by default
     783            twinVal = [0.0,False]
    784784            UseList[G2frame.hist]['Twins'].append([twinMat,twinVal])
    785785            addtwin.SetValue(False)
    … …  
    802802            try:
    803803                val = float(Obj.GetValue())
    804                 if 0. > val > 1.:\
     804                if 0. > val > 1.:
    805805                    raise ValueError
    806806            except ValueError:
    807807                val = UseList[G2frame.hist]['Twins'][it][1][0]
    808808            UseList[G2frame.hist]['Twins'][it][1][0] = val
    809             Obj.SetValue('%.3f'%(val))
     809            sumTw = 0.
     810            for it,twin in enumerate(UseList[G2frame.hist]['Twins']):
     811                if it:
     812                    sumTw += twin[1][0]
     813            UseList[G2frame.hist]['Twins'][0][1][0] = 1.-sumTw
     814            wx.CallLater(100,RepaintHistogramInfo)
    810815           
    811816        def OnTwinRef(event):
    … …  
    818823            it = Indx[Obj.GetId()]
    819824            del UseList[G2frame.hist]['Twins'][it]
     825            sumTw = 0.
     826            for it,twin in enumerate(UseList[G2frame.hist]['Twins']):
     827                if it:
     828                    sumTw += twin[1][0]
     829            UseList[G2frame.hist]['Twins'][0][1][0] = 1.-sumTw
    820830            wx.CallLater(100,RepaintHistogramInfo)           
    821831           

  • trunk/GSASIIstrIO.py

    r1886 r1887  
    24262426        print >>pFile,sigstr
    24272427       
    2428     def PrintTwinsAndSig(pfx,hapData,TwinSig):
     2428    def PrintTwinsAndSig(pfx,twinData,TwinSig):
    24292429        print >>pFile,'\n Twin Law fractions : '
    24302430        ptlbls = ' names :'
    24312431        ptstr =  ' values:'
    24322432        sigstr = ' sig   :'
    2433         for item in hapData:
    2434             ptlbls += '%12s'%(item)
    2435             ptstr += '%12.3f'%(hapData[item][0])
    2436             if pfx+item in TwinSig:
    2437                 sigstr += '%12.3f'%(TwinSig[pfx+item])
     2433        for it,item in enumerate(twinData):
     2434            ptlbls += '     twin #%d'%(it)
     2435            ptstr += '%12.3f'%(item[1][0])
     2436            if pfx+'TwinFr;'+str(it) in TwinSig:
     2437                sigstr += '%12.3f'%(TwinSig[pfx+'TwinFr;'+str(it)])
    24382438            else:
    24392439                sigstr += 12*' '
    … …  
    25342534                    if pfx+item in sigDict:
    25352535                        BabSig[pfx+item] = sigDict[pfx+item]
    2536                 item = 'TwinFr'
    2537                 it = 0
    2538                 sumTwFr = 0.
    2539                 while True:
    2540                     try:
    2541                         hapData['TwinFr'][1][0] = parmDict[pfx+'TwinFr;'+str(it)]
    2542                         if pfx+'TwinFr;'+str(it) in sigDict:
    2543                             TwinFrSig[pfx+'TwinFr;'+str(it)] = sigDict[pfx+'TwinFr;'+str(it)]
    2544                         if it:
    2545                             sumTwFr += hapData['TwinFr'][1][0]
    2546                         it += 1
    2547                     except KeyError:
    2548                         break
    2549 #                hapData['TwinFr'][1][0]
     2536                if 'Twins' in hapData:
     2537                    it = 0
     2538                    sumTwFr = 0.
     2539                    while True:
     2540                        try:
     2541                            hapData['Twins'][it][1][0] = parmDict[pfx+'TwinFr;'+str(it)]
     2542                            if pfx+'TwinFr;'+str(it) in sigDict:
     2543                                TwinFrSig[pfx+'TwinFr;'+str(it)] = sigDict[pfx+'TwinFr;'+str(it)]
     2544                            if it:
     2545                                sumTwFr += hapData['Twins'][it][1][0]
     2546                            it += 1
     2547                        except KeyError:
     2548                            break
     2549                    hapData['Twins'][0][1][0] = 1.-sumTwFr
    25502550
    25512551    if Print:
    … …  
    26092609                    if pfx+'Flack' in ScalExtSig:
    26102610                        print >>pFile,' Flack parameter : %10.3f, sig %10.3f'%(hapData['Flack'][0],ScalExtSig[pfx+'Flack'])
    2611                     if pfx+'TwinFr;1' in TwinFrSig:
    2612                         PrintTwinFrAndSig(pfx,hapData['TwinFr'],TwinFrSig)
     2611                    if len(TwinFrSig):
     2612                        PrintTwinsAndSig(pfx,hapData['Twins'],TwinFrSig)
    26132613
    26142614################################################################################

  • trunk/GSASIIstrMath.py

    r1886 r1887  
    918918        TwinFr = np.array([parmDict[phfx+'TwinFr;'+str(i)] for i in range(len(TwinLaw))])
    919919        if len(TwinLaw) > 1:
    920             TwinFr[0] = 1.-np.sum(TwinFr[1:])       
     920            TwinFr[0] = 1.-np.sum(TwinFr[1:])
     921    nTwin = len(TwinLaw)       
    921922    nRef = len(refDict['RefList'])
    922923    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
    … …  
    931932    bij = Mast*Uij.T
    932933    dFdvDict = {}
    933     dFdfr = np.zeros((nRef,mSize))
    934     dFdx = np.zeros((nRef,mSize,3))
    935     dFdui = np.zeros((nRef,mSize))
    936     dFdua = np.zeros((nRef,mSize,6))
    937     dFdbab = np.zeros((nRef,2))
    938     dFdfl = np.zeros(nRef)
     934    dFdfr = np.squeeze(np.zeros((nRef,nTwin,mSize)))
     935    dFdx = np.squeeze(np.zeros((nRef,nTwin,mSize,3)))
     936    dFdui = np.squeeze(np.zeros((nRef,nTwin,mSize)))
     937    dFdua = np.squeeze(np.zeros((nRef,nTwin,mSize,6)))
     938    dFdbab = np.squeeze(np.zeros((nRef,nTwin,2)))
     939    dFdfl = np.squeeze(np.zeros((nRef,nTwin)))
     940    dFdtw = np.zeros((nRef,nTwin))
    939941    Flack = 1.0
    940942    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType']:
    … …  
    958960        occ = Mdata*Fdata/len(Uniq)
    959961        biso = -SQfactor*Uisodata[:,np.newaxis]
    960         Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*len(TwinLaw),axis=1).T
     962        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*nTwin,axis=1).T
    961963        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
    962964        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])
    963         Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(len(TwinLaw),-1,6)))
     965        Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6)))
    964966        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
    965967        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/len(SGMT)
    966968        fot = (FF+FP-Bab)*occ*Tcorr
    967969        fotp = FPP*occ*Tcorr       
    968         if 'T' in calcControls[hfx+'histType']:
    969             fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
    970             fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
    971         else:
    972             fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
    973             fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
     970        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
     971        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
    974972        fas = np.sum(np.sum(fa,axis=-1),axis=-1)      #real sum over atoms & unique hkl
    975973        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)      #imag sum over atoms & uniq hkl
    976974        fax = np.array([-fot*sinp,-fotp*cosp])   #positions
    977975        fbx = np.array([fot*cosp,-fotp*sinp])
    978         #sum below is over Uniq - twin effects??
     976        #sum below is over Uniq
    979977        dfadfr = np.sum(fa/occ,axis=-2)        #Fdata != 0 ever avoids /0. problem
    980         dfadx = np.sum(twopi*Uniq*np.swapaxes(fax,-2,-1)[:,:,:,np.newaxis],axis=-2)
     978        dfadba = np.sum(-cosp*(occ*Tcorr)[:,np.newaxis],axis=1)
    981979        dfadui = np.sum(-SQfactor*fa,axis=-2)
    982         dfadua = np.sum(-Hij*np.swapaxes(fa,-2,-1)[:,:,:,np.newaxis],axis=-2)
    983         dfadba = np.sum(-cosp*(occ*Tcorr)[:,np.newaxis],axis=1)
     980        if len(TwinLaw) > 1:
     981            dfadx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fax[it],-2,-1)[:,:,:,np.newaxis],axis=-2) for it in range(nTwin)])
     982            dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fa[it],-2,-1)[:,:,:,np.newaxis],axis=-2) for it in range(nTwin)])
     983        else:
     984            dfadx = np.sum(twopi*Uniq*np.swapaxes(fax,-2,-1)[:,:,:,np.newaxis],axis=-2)
     985            dfadua = np.sum(-Hij*np.swapaxes(fa,-2,-1)[:,:,:,np.newaxis],axis=-2)
    984986        if not SGData['SGInv']:
    985987            dfbdfr = np.sum(fb/occ,axis=-2)        #Fdata != 0 ever avoids /0. problem
    986             dfbdx = np.sum(twopi*Uniq*np.swapaxes(fbx,-2,-1)[:,:,:,np.newaxis],axis=2)           
    987             dfbdui = np.sum(-SQfactor*fb,axis=-2)
    988             dfbdua = np.sum(-Hij*np.swapaxes(fb,-2,-1)[:,:,:,np.newaxis],axis=2)
    989988            dfbdba = np.sum(-sinp*(occ*Tcorr)[:,np.newaxis],axis=1)
    990989            dfadfl = np.sum(-fotp[:,np.newaxis]*sinp)
    991990            dfbdfl = np.sum(fotp[:,np.newaxis]*cosp)
     991            dfbdui = np.sum(-SQfactor*fb,axis=-2)
     992            if len(TwinLaw) > 1:
     993                dfbdx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fbx[it],-2,-1)[:,:,:,np.newaxis],axis=2) for it in range(nTwin)])           
     994                dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fb[it],-2,-1)[:,:,:,np.newaxis],axis=2) for it in range(nTwin)])
     995            else:
     996                dfbdx = np.sum(twopi*Uniq*np.swapaxes(fbx,-2,-1)[:,:,:,np.newaxis],axis=2)           
     997                dfbdua = np.sum(-Hij*np.swapaxes(fb,-2,-1)[:,:,:,np.newaxis],axis=2)
    992998        else:
    993999            dfbdfr = np.zeros_like(dfadfr)
    … …  
    10111017            SA = fas[0]-fbs[1]
    10121018            SB = fbs[0]+fas[1]
    1013             dFdfr[iref] = 2.*SA*(dfadfr[0]+dfbdfr[1])*Mdata/len(Uniq)+ \
    1014                 2.*SB*(dfbdfr[0]+dfadfr[1])*Mdata/len(Uniq)
    1015             dFdx[iref] = 2.*SA*(dfadx[0]+dfbdx[1])+2.*SB*(dfbdx[0]+dfadx[1])
    1016             dFdui[iref] = 2.*SA*(dfadui[0]+dfbdui[1])+2.*SB*(dfbdui[0]+dfadui[1])
    1017             dFdua[iref] = 2.*SA*(dfadua[0]+dfbdua[1])+2.*SB*(dfbdua[0]+dfadua[1])
    1018             dFdfl[iref] = -4.*SA*dfadfl-4.*SB*dfbdfl
     1019            if nTwin > 1:
     1020                dFdfr[iref] = [2.*SA[it]*(dfadfr[it][0]+dfbdfr[it][1])*Mdata/len(Uniq[it])+ \
     1021                    2.*SB[it]*(dfbdfr[it][0]+dfadfr[it][1])*Mdata/len(Uniq[it]) for it in range(nTwin)]
     1022                dFdx[iref] = [2.*SA[it]*(dfadx[it][0]+dfbdx[it][1])+2.*SB[it]*(dfbdx[it][0]+dfadx[it][1]) for it in range(nTwin)]
     1023                dFdui[iref] = [2.*SA[it]*(dfadui[it][0]+dfbdui[it][1])+2.*SB[it]*(dfbdui[it][0]+dfadui[it][1]) for it in range(nTwin)]
     1024                dFdua[iref] = [2.*SA[it]*(dfadua[it][0]+dfbdua[it][1])+2.*SB[it]*(dfbdua[it][0]+dfadua[it][1]) for it in range(nTwin)]
     1025                dFdfl[iref] = -SA*dfadfl-SB*dfbdfl
     1026                dFdtw[iref] = 2.*SA+2.*SB
     1027            else:
     1028                dFdfr[iref] = 2.*SA*(dfadfr[0]+dfbdfr[1])*Mdata/len(Uniq)+ \
     1029                    2.*SB*(dfbdfr[0]+dfadfr[1])*Mdata/len(Uniq)
     1030                dFdx[iref] = 2.*SA*(dfadx[0]+dfbdx[1])+2.*SB*(dfbdx[0]+dfadx[1])
     1031                dFdui[iref] = 2.*SA*(dfadui[0]+dfbdui[1])+2.*SB*(dfbdui[0]+dfadui[1])
     1032                dFdua[iref] = 2.*SA*(dfadua[0]+dfbdua[1])+2.*SB*(dfbdua[0]+dfadua[1])
     1033                dFdfl[iref] = -SA*dfadfl-SB*dfbdfl
    10191034        dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
    10201035            2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
    10211036           
    10221037        #loop over atoms - each dict entry is list of derivatives for all the reflections
    1023     for i in range(len(Mdata)):
    1024         dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
    1025         dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
    1026         dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
    1027         dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
    1028         dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
    1029         dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
    1030         dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
    1031         dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
    1032         dFdvDict[pfx+'AU12:'+str(i)] = 0.5*dFdua.T[3][i]
    1033         dFdvDict[pfx+'AU13:'+str(i)] = 0.5*dFdua.T[4][i]
    1034         dFdvDict[pfx+'AU23:'+str(i)] = 0.5*dFdua.T[5][i]
     1038    if nTwin > 1:
     1039        for i in range(len(Mdata)):
     1040            dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,np.newaxis],axis=0)
     1041            dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,np.newaxis],axis=0)
     1042            dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,np.newaxis],axis=0)
     1043            dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,np.newaxis],axis=0)
     1044            dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,np.newaxis],axis=0)
     1045            dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,np.newaxis],axis=0)
     1046            dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,np.newaxis],axis=0)
     1047            dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,np.newaxis],axis=0)
     1048            dFdvDict[pfx+'AU12:'+str(i)] = np.sum(0.5*dFdua.T[3][i]*TwinFr[:,np.newaxis],axis=0)
     1049            dFdvDict[pfx+'AU13:'+str(i)] = np.sum(0.5*dFdua.T[4][i]*TwinFr[:,np.newaxis],axis=0)
     1050            dFdvDict[pfx+'AU23:'+str(i)] = np.sum(0.5*dFdua.T[5][i]*TwinFr[:,np.newaxis],axis=0)
     1051        dFdvDict[phfx+'Flack'] = np.sum(dFdfl.T*TwinFr[:,np.newaxis],axis=0)
     1052    else:
     1053        for i in range(len(Mdata)):
     1054            dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
     1055            dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
     1056            dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
     1057            dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
     1058            dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
     1059            dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
     1060            dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
     1061            dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
     1062            dFdvDict[pfx+'AU12:'+str(i)] = 0.5*dFdua.T[3][i]
     1063            dFdvDict[pfx+'AU13:'+str(i)] = 0.5*dFdua.T[4][i]
     1064            dFdvDict[pfx+'AU23:'+str(i)] = 0.5*dFdua.T[5][i]
     1065        dFdvDict[phfx+'Flack'] = dFdfl.T
    10351066    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
    10361067    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
    1037     dFdvDict[phfx+'Flack'] = dFdfl.T
     1068    if nTwin > 1:
     1069        for i in range(nTwin-1):    #skip the base twin element
     1070            dFdvDict[phfx+'TwinFr;'+str(i+1)] = dFdtw.T[i+1]
    10381071    return dFdvDict
    10391072   
Note: See TracChangeset for help on using the changeset viewer.