Changeset 1919


Ignore:
Timestamp:
Jul 6, 2015 8:41:52 AM (7 years ago)
Author:
vondreele
Message:

Twin derivatives

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIstrMath.py

    r1918 r1919  
    966966        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
    967967        FF = refDict['FF']['FF'][iref].T[Tindx].T
    968         Uniq = np.inner(H,SGMT)             # array(3,nSGOp) or
     968        Uniq = np.inner(H,SGMT)             # array(nSGOp,3) or (nTwin,nSGOp,3)
    969969        Phi = np.inner(H,SGT)
    970970        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
    971971        sinp = np.sin(phase)
    972972        cosp = np.cos(phase)
    973         occ = Mdata*Fdata/len(Uniq)
     973        occ = Mdata*Fdata/len(SGT)
    974974        biso = -SQfactor*Uisodata[:,np.newaxis]
    975         Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*nTwin,axis=1).T
     975        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*nTwin,axis=1)
    976976        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
    977977        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])
    978978        Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6)))
    979         Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
    980         Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*occ
     979        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
     980        Tcorr = (np.reshape(Tiso,Tuij.shape)*Tuij).T*occ
    981981        fot = (FF+FP-Bab)*Tcorr
    982982        fotp = FPP*Tcorr       
    983         fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
     983        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) # array(2,nTwin,nEqv,nAtoms)
    984984        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
    985         fas = np.sum(np.sum(fa,axis=-1),axis=-1)      #real sum over atoms & unique hkl
     985        fas = np.sum(np.sum(fa,axis=-1),axis=-1)      #real sum over atoms & unique hkl array(2,nTwins)
    986986        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)      #imag sum over atoms & uniq hkl
    987         fax = np.array([-fot*sinp,-fotp*cosp])   #positions
     987        fax = np.array([-fot*sinp,-fotp*cosp])   #positions array(2,ntwi,nEqv,nAtoms)
    988988        fbx = np.array([fot*cosp,-fotp*sinp])
    989989        #sum below is over Uniq
    990         dfadfr = np.sum(fa/occ,axis=-2)        #Fdata != 0 avoids /0. problem
     990        dfadfr = np.sum(fa/occ,axis=-2)        #array(2,ntwin,nAtom) Fdata != 0 avoids /0. problem
    991991        dfadba = np.sum(-cosp*Tcorr[:,np.newaxis],axis=1)
    992992        dfadui = np.sum(-SQfactor*fa,axis=-2)
    993         if len(TwinLaw) > 1:
     993        if nTwin > 1:
    994994            dfadx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fax,-2,-1)[:,it,:,:,np.newaxis],axis=-2) for it in range(nTwin)])
    995995            dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fa,-2,-1)[:,it,:,:,np.newaxis],axis=-2) for it in range(nTwin)])
     996            # array(nTwin,2,nAtom,3) & array(nTwin,2,nAtom,6)
    996997        else:
    997998            dfadx = np.sum(twopi*Uniq*np.swapaxes(fax,-2,-1)[:,:,:,np.newaxis],axis=-2)
    998999            dfadua = np.sum(-Hij*np.swapaxes(fa,-2,-1)[:,:,:,np.newaxis],axis=-2)
     1000            # array(2,nAtom,3) & array(2,nAtom,6)
    9991001        if not SGData['SGInv']:
    10001002            dfbdfr = np.sum(fb/occ,axis=-2)        #Fdata != 0 avoids /0. problem
     
    10311033            SB = fbs[0]+fas[1]
    10321034            if nTwin > 1:
    1033                 dFdfr[iref] = TwMask[:,np.newaxis]*[SA[it]*(dfadfr[0][it]+dfbdfr[1][it])*Mdata/len(Uniq[it])+ \
    1034                     SB[it]*(dfbdfr[0][it]+dfadfr[1][it])*Mdata/len(Uniq[it]) for it in range(nTwin)]
    1035                 dFdx[iref] = TwMask[:,np.newaxis,np.newaxis]*[SA[it]*(dfadx[it][0]+dfbdx[it][1])+SB[it]*(dfbdx[it][0]+dfadx[it][1]) for it in range(nTwin)]
    1036                 dFdui[iref] = TwMask[:,np.newaxis]*[SA[it]*(dfadui[0][it]+dfbdui[1][it])+SB[it]*(dfbdui[0][it]+dfadui[1][it]) for it in range(nTwin)]
    1037                 dFdua[iref] = TwMask[:,np.newaxis,np.newaxis]*[SA[it]*(dfadua[it][0]+dfbdua[it][1])+SB[it]*(dfbdua[it][0]+dfadua[it][1]) for it in range(nTwin)]
    1038                 dFdtw[iref] = TwMask*np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2
     1035                dFdfr[iref] = [2.*TwMask[it]*SA[it]*(dfadfr[0][it]+dfbdfr[1][it])*Mdata/len(Uniq[it])+ \
     1036                    2.*TwMask[it]*SB[it]*(dfbdfr[0][it]+dfadfr[1][it])*Mdata/len(Uniq[it]) for it in range(nTwin)]
     1037                dFdx[iref] = [2.*TwMask[it]*SA[it]*(dfadx[it][0]+dfbdx[it][1])+2.*TwMask[it]*SB[it]*(dfbdx[it][0]+dfadx[it][1]) for it in range(nTwin)]
     1038                dFdui[iref] = [2.*TwMask[it]*SA[it]*(dfadui[0][it]+dfbdui[1][it])+2.*TwMask[it]*SB[it]*(dfbdui[0][it]+dfadui[1][it]) for it in range(nTwin)]
     1039                dFdua[iref] = [2.*TwMask[it]*SA[it]*(dfadua[it][0]+dfbdua[it][1])+2.*TwMask[it]*SB[it]*(dfbdua[it][0]+dfadua[it][1]) for it in range(nTwin)]
     1040                dFdtw[iref] = 2.*TwMask*np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2
    10391041            else:   #these are good for no twin single crystals
    10401042                dFdfr[iref] = 2.*SA*(dfadfr[0]+dfbdfr[1])*Mdata/len(Uniq)+ \
    1041                     2.*SB*(dfbdfr[0]+dfadfr[1])*Mdata/len(Uniq)
    1042                 dFdx[iref] = 2.*SA*(dfadx[0]+dfbdx[1])+2.*SB*(dfbdx[0]+dfadx[1])
    1043                 dFdui[iref] = 2.*SA*(dfadui[0]+dfbdui[1])+2.*SB*(dfbdui[0]+dfadui[1])
    1044                 dFdua[iref] = 2.*SA*(dfadua[0]+dfbdua[1])+2.*SB*(dfbdua[0]+dfadua[1])
    1045                 dFdfl[iref] = -SA*dfadfl-SB*dfbdfl
     1043                    2.*SB*(dfbdfr[0]+dfadfr[1])*Mdata/len(Uniq) #array(nRef,nAtom)
     1044                dFdx[iref] = 2.*SA*(dfadx[0]+dfbdx[1])+2.*SB*(dfbdx[0]+dfadx[1])    #array(nRef,nAtom,3)
     1045                dFdui[iref] = 2.*SA*(dfadui[0]+dfbdui[1])+2.*SB*(dfbdui[0]+dfadui[1])   #array(nRef,nAtom)
     1046                dFdua[iref] = 2.*SA*(dfadua[0]+dfbdua[1])+2.*SB*(dfbdua[0]+dfadua[1])    #array(nRef,nAtom,6)
     1047                dFdfl[iref] = -SA*dfadfl-SB*dfbdfl  #array(nRef,)
    10461048                dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
    10471049                    2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
    1048            
    10491050        #loop over atoms - each dict entry is list of derivatives for all the reflections
    10501051    if nTwin > 1:
    10511052        for i in range(len(Mdata)):     #these all OK?
    10521053            dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,np.newaxis],axis=0)
    1053             dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,np.newaxis],axis=0)/2.
    1054             dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,np.newaxis],axis=0)/2.
    1055             dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,np.newaxis],axis=0)/2.
    1056             dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,np.newaxis],axis=0)/2.
    1057             dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,np.newaxis],axis=0)/2.
    1058             dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,np.newaxis],axis=0)/2.
    1059             dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,np.newaxis],axis=0)/2.
    1060             dFdvDict[pfx+'AU12:'+str(i)] = np.sum(0.5*dFdua.T[3][i]*TwinFr[:,np.newaxis],axis=0)/2.
    1061             dFdvDict[pfx+'AU13:'+str(i)] = np.sum(0.5*dFdua.T[4][i]*TwinFr[:,np.newaxis],axis=0)/2.
    1062             dFdvDict[pfx+'AU23:'+str(i)] = np.sum(0.5*dFdua.T[5][i]*TwinFr[:,np.newaxis],axis=0)/2.
     1054            dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,np.newaxis],axis=0)
     1055            dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,np.newaxis],axis=0)
     1056            dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,np.newaxis],axis=0)
     1057            dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,np.newaxis],axis=0)
     1058            dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,np.newaxis],axis=0)
     1059            dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,np.newaxis],axis=0)
     1060            dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,np.newaxis],axis=0)
     1061            dFdvDict[pfx+'AU12:'+str(i)] = np.sum(0.5*dFdua.T[3][i]*TwinFr[:,np.newaxis],axis=0)
     1062            dFdvDict[pfx+'AU13:'+str(i)] = np.sum(0.5*dFdua.T[4][i]*TwinFr[:,np.newaxis],axis=0)
     1063            dFdvDict[pfx+'AU23:'+str(i)] = np.sum(0.5*dFdua.T[5][i]*TwinFr[:,np.newaxis],axis=0)
    10631064    else:
    10641065        for i in range(len(Mdata)):
Note: See TracChangeset for help on using the changeset viewer.