Changeset 1992


Ignore:
Timestamp:
Oct 7, 2015 3:35:11 PM (6 years ago)
Author:
vondreele
Message:

fix SS position derivative error

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r1990 r1992  
    10621062    dGdBx = np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(twopi*dHdXB*np.cos(twopi*HdotXB)-np.sin(twopi*HdotXB))*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
    10631063#    GSASIIpath.IPyBreak()                 
    1064     return np.array([cosHA,sinHA]),dGdAf,dGdBf,dGdAx,dGdBx,dGdAu,dGdBu      #ops X atoms
     1064    return np.array([cosHA,sinHA]),np.concatenate((dGdAf,dGdBf),-1), \
     1065        np.concatenate((dGdAx,dGdBx),-1),np.concatenate((dGdAu,dGdBu),-1)      #ops X atoms
    10651066   
    10661067def posFourier(tau,psin,pcos,smul):
  • trunk/GSASIIstrMath.py

    r1990 r1992  
    10861086    print 'nRef %d time %.4f\r'%(nRef,time.time()-time0)
    10871087
    1088 #def SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
    1089 #    '''
    1090 #    Compute super structure factors for all h,k,l,m for phase
    1091 #    puts the result, F^2, in each ref[8+im] in refList
    1092 #    input:
    1093 #   
    1094 #    :param dict refDict: where
    1095 #        'RefList' list where each ref = h,k,l,m,d,...
    1096 #        'FF' dict of form factors - filed in below
    1097 #    :param np.array G:      reciprocal metric tensor
    1098 #    :param str pfx:    phase id string
    1099 #    :param dict SGData: space group info. dictionary output from SpcGroup
    1100 #    :param dict calcControls:
    1101 #    :param dict ParmDict:
    1102 #
    1103 #    '''       
    1104 #    phfx = pfx.split(':')[0]+hfx
    1105 #    ast = np.sqrt(np.diag(G))
    1106 #    Mast = twopisq*np.multiply.outer(ast,ast)
    1107 #   
    1108 #    SGInv = SGData['SGInv']
    1109 #    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
    1110 #    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
    1111 #    FFtables = calcControls['FFtables']
    1112 #    BLtables = calcControls['BLtables']
    1113 #    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
    1114 #    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
    1115 #    FF = np.zeros(len(Tdata))
    1116 #    if 'NC' in calcControls[hfx+'histType']:
    1117 #        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
    1118 #    else:
    1119 #        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
    1120 #        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
    1121 #    Uij = np.array(G2lat.U6toUij(Uijdata))
    1122 #    bij = Mast*Uij.T
    1123 #    if not len(refDict['FF']):
    1124 #        if 'N' in calcControls[hfx+'histType']:
    1125 #            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
    1126 #        else:
    1127 #            dat = G2el.getFFvalues(FFtables,0.)       
    1128 #        refDict['FF']['El'] = dat.keys()
    1129 #        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))
    1130 #    time0 = time.time()
    1131 #    nref = len(refDict['RefList'])/100   
    1132 #    for iref,refl in enumerate(refDict['RefList']):
    1133 #        if 'NT' in calcControls[hfx+'histType']:
    1134 #            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl[14+im])
    1135 #        fbs = np.array([0,0])
    1136 #        H = refl[:4]
    1137 #        SQ = 1./(2.*refl[4+im])**2
    1138 #        SQfactor = 4.0*SQ*twopisq
    1139 #        Bab = parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor)
    1140 #        if not np.any(refDict['FF']['FF'][iref]):                #no form factors - 1st time thru StructureFactor
    1141 #            if 'N' in calcControls[hfx+'histType']:
    1142 #                dat = G2el.getBLvalues(BLtables)
    1143 #                refDict['FF']['FF'][iref] = dat.values()
    1144 #            else:       #'X'
    1145 #                dat = G2el.getFFvalues(FFtables,SQ)
    1146 #                refDict['FF']['FF'][iref] = dat.values()
    1147 #        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
    1148 #        FF = refDict['FF']['FF'][iref][Tindx]
    1149 #        SSUniq = np.inner(H,SSGMT)
    1150 #        SSPhi = np.inner(H,SSGT)
    1151 #        if SGInv:   #if centro - expand HKL sets
    1152 #            SSUniq = np.vstack((SSUniq,-SSUniq))
    1153 #            SSPhi = np.hstack((SSPhi,-SSPhi))
    1154 #        phase = twopi*(np.inner(SSUniq[:,:3],(dXdata.T+Xdata.T))+SSPhi[:,np.newaxis])
    1155 #        sinp = np.sin(phase)
    1156 #        cosp = np.cos(phase)
    1157 #        biso = -SQfactor*Uisodata
    1158 #        Tiso = np.where(biso<1.,np.exp(biso),1.0)
    1159 #        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in SSUniq[:,:3]])
    1160 #        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
    1161 #        Tcorr = Tiso*Tuij*Mdata*Fdata/len(SSUniq)
    1162 #        fa = np.array([(FF+FP-Bab)*cosp*Tcorr,-FPP*sinp*Tcorr])     #2 x sym x atoms
    1163 #        fb = np.array([FPP*cosp*Tcorr,(FF+FP-Bab)*sinp*Tcorr])      #swapped around - better?
    1164 #        GfpuA = G2mth.Modulation(waveTypes,SSUniq,FSSdata,XSSdata,USSdata,Mast) #2 x sym X atoms
    1165 #        GSASIIpath.IPyBreak()
    1166 #        fa *= GfpuA
    1167 #        fb *= GfpuA       
    1168 #        fas = np.sum(np.sum(fa,axis=1),axis=1)
    1169 #        fbs = np.sum(np.sum(fb,axis=1),axis=1)
    1170 #        fasq = fas**2
    1171 #        fbsq = fbs**2        #imaginary
    1172 #        refl[9+im] = np.sum(fasq)+np.sum(fbsq)
    1173 #        refl[7+im] = np.sum(fasq)+np.sum(fbsq)
    1174 #        refl[10+im] = atan2d(fbs[0],fas[0])
    1175 #        if not iref%nref: print 'ref no. %d time %.4f\r'%(iref,time.time()-time0),
    1176 #    print '\nref no. %d time %.4f\r'%(iref,time.time()-time0)
    1177 #
    11781088def SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
    11791089    'Needs a doc string'
     
    11981108    nRef = len(refDict['RefList'])
    11991109    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
    1200     mSize = len(Mdata)
     1110    mSize = len(Mdata)  #no. atoms
    12011111    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
    12021112    FF = np.zeros(len(Tdata))
     
    12241134        dFdfl = np.zeros((nRef,nTwin))
    12251135        dFdtw = np.zeros((nRef,nTwin))
     1136        dFdGf = np.zeros((nRef,nTwin,mSize,FSSdata.shape[1],2))
     1137        dFdGx = np.zeros((nRef,nTwin,mSize,XSSdata.shape[1],6))
     1138        dFdGu = np.zeros((nRef,nTwin,mSize,USSdata.shape[1],12))
    12261139    else:
    12271140        dFdfr = np.zeros((nRef,mSize))
     
    12321145        dFdfl = np.zeros((nRef))
    12331146        dFdtw = np.zeros((nRef))
     1147        dFdGf = np.zeros((nRef,mSize,FSSdata.shape[1],2))
     1148        dFdGx = np.zeros((nRef,mSize,XSSdata.shape[1],6))
     1149        dFdGu = np.zeros((nRef,mSize,USSdata.shape[1],12))
    12341150    Flack = 1.0
    12351151    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
     
    12671183        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[0]*len(TwinLaw),axis=1).T    #ops x atoms
    12681184        HbH = -np.sum(Uniq[:,nxs,:3]*np.inner(Uniq[:,:3],bij),axis=-1)  #ops x atoms
    1269 #        GSASIIpath.IPyBreak()
    12701185        Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in Uniq]) #atoms x 3x3
    12711186        Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6)))
    12721187        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)     #ops x atoms
    12731188        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[0]  #ops x atoms
    1274         fot = (FF+FP-Bab)*occ*Tcorr     #ops x atoms
    1275         fotp = FPP*occ*Tcorr            #ops x atoms
    1276         GfpuA,dGdAf,dGdBf,dGdAx,dGdBx,dGdAu,dGdBu =     \
    1277             G2mth.ModulationDerv(waveTypes,Uniq,FSSdata,XSSdata,USSdata,Mast)
     1189        fot = (FF+FP-Bab)*Tcorr     #ops x atoms
     1190        fotp = FPP*Tcorr            #ops x atoms
     1191        GfpuA,dGdf,dGdx,dGdu = G2mth.ModulationDerv(waveTypes,Uniq,FSSdata,XSSdata,USSdata,Mast)
    12781192        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) # array(2,nTwin,nEqv,nAtoms)
    12791193        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
     
    12901204        dfadfr = np.sum(fag/occ,axis=1)        #Fdata != 0 ever avoids /0. problem
    12911205        dfbdfr = np.sum(fbg/occ,axis=1)        #Fdata != 0 avoids /0. problem
    1292         dfadba = np.sum(-cosp*(occ*Tcorr)[:,nxs],axis=1)
    1293         dfbdba = np.sum(-sinp*(occ*Tcorr)[:,nxs],axis=1)
     1206        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
     1207        dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)
    12941208        dfadui = np.sum(-SQfactor*fag,axis=1)
    12951209        dfbdui = np.sum(-SQfactor*fbg,axis=1)
     
    13071221            # array(2,nAtom,3) & array(2,nAtom,6)
    13081222        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for al2O3!   
     1223#        GSASIIpath.IPyBreak()
    13091224        if not SGData['SGInv'] and len(TwinLaw) == 1:   #Flack derivative
    13101225            dfadfl = np.sum(-FPP*Tcorr*sinp)
Note: See TracChangeset for help on using the changeset viewer.