Ignore:
Timestamp:
Oct 7, 2015 1:06:38 PM (7 years ago)
Author:
vondreele
Message:

Structure factor derivs for SS continued...

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIstrMath.py

    r1988 r1990  
    12571257        Uniq = np.inner(H,SSGMT)
    12581258        Phi = np.inner(H,SSGT)
     1259        if SGInv:   #if centro - expand HKL sets
     1260            Uniq = np.vstack((Uniq,-Uniq))
     1261            Phi = np.hstack((Phi,-Phi))
    12591262        phase = twopi*(np.inner(Uniq[:,:3],(dXdata+Xdata).T)+Phi[:,nxs])
    12601263        sinp = np.sin(phase)
    12611264        cosp = np.cos(phase)
    1262         occ = Mdata*Fdata/len(Uniq)
     1265        occ = Mdata*Fdata/Uniq.shape[0]
    12631266        biso = -SQfactor*Uisodata[:,nxs]
    12641267        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[0]*len(TwinLaw),axis=1).T    #ops x atoms
     
    12681271        Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6)))
    12691272        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)     #ops x atoms
    1270         Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[1]  #ops x atoms
     1273        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[0]  #ops x atoms
    12711274        fot = (FF+FP-Bab)*occ*Tcorr     #ops x atoms
    12721275        fotp = FPP*occ*Tcorr            #ops x atoms
     
    12861289        #sum below is over Uniq
    12871290        dfadfr = np.sum(fag/occ,axis=1)        #Fdata != 0 ever avoids /0. problem
     1291        dfbdfr = np.sum(fbg/occ,axis=1)        #Fdata != 0 avoids /0. problem
    12881292        dfadba = np.sum(-cosp*(occ*Tcorr)[:,nxs],axis=1)
     1293        dfbdba = np.sum(-sinp*(occ*Tcorr)[:,nxs],axis=1)
    12891294        dfadui = np.sum(-SQfactor*fag,axis=1)
     1295        dfbdui = np.sum(-SQfactor*fbg,axis=1)
    12901296        if nTwin > 1:
    12911297            dfadx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fax,-2,-1)[:,it,:,:,np.newaxis],axis=-2) for it in range(nTwin)])
     1298            dfbdx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fbx,-2,-1)[:,it,:,:,np.newaxis],axis=-2) for it in range(nTwin)])           
    12921299            dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fag,-2,-1)[:,it,:,:,np.newaxis],axis=-2) for it in range(nTwin)])
     1300            dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fbg,-2,-1)[:,it,:,:,np.newaxis],axis=-2) for it in range(nTwin)])
    12931301            # array(nTwin,2,nAtom,3) & array(nTwin,2,nAtom,6)
    12941302        else:
    12951303            dfadx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fax,-2,-1)[:,:,:,np.newaxis],axis=-2)
     1304            dfbdx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fbx,-2,-1)[:,:,:,np.newaxis],axis=-2)           
    12961305            dfadua = np.sum(-Hij*np.swapaxes(fag,-2,-1)[:,:,:,np.newaxis],axis=-2)
     1306            dfbdua = np.sum(-Hij*np.swapaxes(fbg,-2,-1)[:,:,:,np.newaxis],axis=-2)
    12971307            # array(2,nAtom,3) & array(2,nAtom,6)
    12981308        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for al2O3!   
    1299         if not SGData['SGInv']:
    1300             dfbdfr = np.sum(fbg/occ,axis=1)        #Fdata != 0 avoids /0. problem
    1301             dfbdba = np.sum(-sinp*Tcorr[:,np.newaxis],axis=1)
    1302             dfbdui = np.sum(-SQfactor*fb,axis=1)
    1303             if len(TwinLaw) > 1:
    1304                 dfbdx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fbx,-2,-1)[:,it,:,:,np.newaxis],axis=2) for it in range(nTwin)])           
    1305                 dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fbg,-2,-1)[:,it,:,:,np.newaxis],axis=2) for it in range(nTwin)])
    1306             else:
    1307                 dfadfl = np.sum(-FPP*Tcorr*sinp)
    1308                 dfbdfl = np.sum(FPP*Tcorr*cosp)
    1309                 dfbdx = np.sum(twopi*Uniq*np.swapaxes(fbx,-2,-1)[:,:,:,np.newaxis],axis=2)           
    1310                 dfbdua = np.sum(-Hij*np.swapaxes(fbg,-2,-1)[:,:,:,np.newaxis],axis=2)
     1309        if not SGData['SGInv'] and len(TwinLaw) == 1:   #Flack derivative
     1310            dfadfl = np.sum(-FPP*Tcorr*sinp)
     1311            dfbdfl = np.sum(FPP*Tcorr*cosp)
    13111312        else:
    1312             dfbdfr = np.zeros_like(dfadfr)
    1313             dfbdx = np.zeros_like(dfadx)
    1314             dfbdui = np.zeros_like(dfadui)
    1315             dfbdua = np.zeros_like(dfadua)
    1316             dfbdba = np.zeros_like(dfadba)
    1317             dfadfl = 0.0
    1318             dfbdfl = 0.0
     1313            dfadfl = 1.0
     1314            dfbdfl = 1.0
    13191315        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
    13201316        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro
Note: See TracChangeset for help on using the changeset viewer.