Changeset 2115


Ignore:
Timestamp:
Jan 7, 2016 9:48:49 AM (6 years ago)
Author:
vondreele
Message:

remove im from SStructureFactor arguments (always = 1 & not needed)
remove Twin segments from SSructureFactor - not used
define new SStructureFactorTw for twinned super lattices (probably not working yet)
add missing Flack parm derivative for super structures

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIImath.py

    r2110 r2115  
    10091009    else:
    10101010        D = H[:,3:]*glTau[nxs,:]              #m*e*tau; refBlk x ops X ngl
     1011        HdotXD = twopi*(HdotX+D[:,nxs,:])
     1012    cosHA = np.sum(Fmod*HbH*np.cos(HdotXD)*glWt,axis=-1)       #real part; refBlk X ops x atoms; sum for G-L integration
     1013    sinHA = np.sum(Fmod*HbH*np.sin(HdotXD)*glWt,axis=-1)       #imag part; ditto
     1014    return np.array([cosHA,sinHA])             # 2 x refBlk x SGops x atoms
     1015   
     1016def ModulationTw(H,HP,nWaves,Fmod,Xmod,Umod,glTau,glWt):
     1017    '''
     1018    H: array nRefBlk x tw x ops X hklt
     1019    HP: array nRefBlk x tw x ops X hklt proj to hkl
     1020    Fmod: array 2 x atoms x waves    (sin,cos terms)
     1021    Xmod: array atoms X ngl X 3
     1022    Umod: array atoms x ngl x 3x3
     1023    glTau,glWt: arrays Gauss-Lorentzian pos & wts
     1024    '''
     1025   
     1026    if nWaves[2]:
     1027        if len(HP.shape) > 3:   #Blocks of reflections
     1028            HbH = np.exp(-np.sum(HP[:,:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
     1029        else:   #single reflections
     1030            HbH = np.exp(-np.sum(HP[:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
     1031    else:
     1032        HbH = 1.0
     1033    HdotX = np.inner(HP,Xmod)                   #refBlk x tw x ops x atoms X ngl
     1034    if len(H.shape) > 3:
     1035        D = glTau*H[:,:,:,3:,nxs]              #m*e*tau; refBlk x tw x ops X ngl
     1036        HdotXD = twopi*(HdotX+D[:,:,:,nxs,:])
     1037    else:
     1038        D = H*glTau[nxs,:]              #m*e*tau; refBlk x ops X ngl
    10111039        HdotXD = twopi*(HdotX+D[:,nxs,:])
    10121040    cosHA = np.sum(Fmod*HbH*np.cos(HdotXD)*glWt,axis=-1)       #real part; refBlk X ops x atoms; sum for G-L integration

  • trunk/GSASIIstrMath.py

    r2111 r2115  
    14251425    return dFdvDict
    14261426   
    1427 def SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
     1427def SStructureFactor(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
    14281428    '''
    1429     Compute super structure factors for all h,k,l,m for phase
     1429    Compute super structure factors for all h,k,l,m for phase - no twins
     1430    puts the result, F^2, in each ref[9] in refList
     1431    works on blocks of 32 reflections for speed
     1432    input:
     1433   
     1434    :param dict refDict: where
     1435        'RefList' list where each ref = h,k,l,m,it,d,...
     1436        'FF' dict of form factors - filed in below
     1437    :param np.array G:      reciprocal metric tensor
     1438    :param str pfx:    phase id string
     1439    :param dict SGData: space group info. dictionary output from SpcGroup
     1440    :param dict calcControls:
     1441    :param dict ParmDict:
     1442
     1443    '''
     1444    phfx = pfx.split(':')[0]+hfx
     1445    ast = np.sqrt(np.diag(G))
     1446    Mast = twopisq*np.multiply.outer(ast,ast)   
     1447    SGInv = SGData['SGInv']
     1448    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
     1449    SGT = np.array([ops[1] for ops in SGData['SGOps']])
     1450    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
     1451    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
     1452    FFtables = calcControls['FFtables']
     1453    BLtables = calcControls['BLtables']
     1454    Flack = 1.0
     1455    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
     1456        Flack = 1.-2.*parmDict[phfx+'Flack']
     1457    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
     1458    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
     1459    ngl,nWaves,Fmod,Xmod,Umod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,Mast)
     1460    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
     1461    FF = np.zeros(len(Tdata))
     1462    if 'NC' in calcControls[hfx+'histType']:
     1463        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
     1464    elif 'X' in calcControls[hfx+'histType']:
     1465        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
     1466        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
     1467    Uij = np.array(G2lat.U6toUij(Uijdata)).T
     1468    bij = Mast*Uij
     1469    blkSize = 32       #no. of reflections in a block
     1470    nRef = refDict['RefList'].shape[0]
     1471    if not len(refDict['FF']):
     1472        if 'N' in calcControls[hfx+'histType']:
     1473            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
     1474            refDict['FF']['El'] = dat.keys()
     1475            refDict['FF']['FF'] = np.ones((nRef,len(dat)))*dat.values()           
     1476        else:
     1477            dat = G2el.getFFvalues(FFtables,0.)       
     1478            refDict['FF']['El'] = dat.keys()
     1479            refDict['FF']['FF'] = np.ones((nRef,len(dat)))
     1480            for iref,ref in enumerate(refDict['RefList']):
     1481                SQ = 1./(2.*ref[5])**2
     1482                dat = G2el.getFFvalues(FFtables,SQ)
     1483                refDict['FF']['FF'][iref] *= dat.values()
     1484    time0 = time.time()
     1485#reflection processing begins here - big arrays!
     1486    iBeg = 0
     1487    while iBeg < nRef:
     1488        iFin = min(iBeg+blkSize,nRef)
     1489        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
     1490        H = refl.T[:4]                          #array(blkSize,4)
     1491        HP = H[:3]+modQ[:,nxs]*H[3:]            #projected hklm to hkl
     1492        SQ = 1./(2.*refl.T[5])**2               #array(blkSize)
     1493        SQfactor = 4.0*SQ*twopisq               #ditto prev.
     1494        Uniq = np.inner(H.T,SSGMT)
     1495        UniqP = np.inner(HP.T,SGMT)
     1496        Phi = np.inner(H.T,SSGT)
     1497        if SGInv:   #if centro - expand HKL sets
     1498            Uniq = np.hstack((Uniq,-Uniq))
     1499            Phi = np.hstack((Phi,-Phi))
     1500            UniqP = np.hstack((UniqP,-UniqP))
     1501        if 'T' in calcControls[hfx+'histType']:
     1502            if 'P' in calcControls[hfx+'histType']:
     1503                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
     1504            else:
     1505                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
     1506            FP = np.repeat(FP.T,Uniq.shape[1],axis=0)
     1507            FPP = np.repeat(FPP.T,Uniq.shape[1],axis=0)
     1508        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),Uniq.shape[1])
     1509        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
     1510        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,Uniq.shape[1],axis=0)
     1511        phase = twopi*(np.inner(Uniq[:,:,:3],(dXdata.T+Xdata.T))-Phi[:,:,nxs])
     1512        sinp = np.sin(phase)
     1513        cosp = np.cos(phase)
     1514        biso = -SQfactor*Uisodata[:,nxs]
     1515        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[1],axis=1).T
     1516        HbH = -np.sum(UniqP[:,:,nxs,:]*np.inner(UniqP[:,:,:],bij),axis=-1)  #use hklt proj to hkl
     1517        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
     1518        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[1]  #refBlk x ops x atoms
     1519        if 'T' in calcControls[hfx+'histType']:
     1520            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
     1521            fb = np.array([np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
     1522        else:
     1523            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
     1524            fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
     1525        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x refBlk x sym X atoms
     1526        fag = fa*GfpuA[0]-fb*GfpuA[1]   #real; 2 x refBlk x sym x atoms
     1527        fbg = fb*GfpuA[0]+fa*GfpuA[1]
     1528        fas = np.sum(np.sum(fag,axis=-1),axis=-1)   #2 x refBlk; sum sym & atoms
     1529        fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
     1530#        GSASIIpath.IPyBreak()
     1531        if 'P' in calcControls[hfx+'histType']:
     1532            refl.T[10] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2    #square of sums
     1533#            refl.T[10] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0)
     1534            refl.T[11] = atan2d(fbs[0],fas[0])  #ignore f' & f"
     1535        else:
     1536            refl.T[10] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2    #square of sums
     1537            refl.T[8] = np.copy(refl.T[10])               
     1538            refl.T[11] = atan2d(fbs[0],fas[0])  #ignore f' & f"
     1539        iBeg += blkSize
     1540    print 'nRef %d time %.4f\r'%(nRef,time.time()-time0)
     1541
     1542def SStructureFactorTw(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
     1543    '''
     1544    Compute super structure factors for all h,k,l,m for phase - twins only
    14301545    puts the result, F^2, in each ref[8+im] in refList
    1431     works on blocks of 100 reflections for speed
     1546    works on blocks of 32 reflections for speed
    14321547    input:
    14331548   
    … …  
    14961611        iFin = min(iBeg+blkSize,nRef)
    14971612        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
    1498         H = refl.T[:4]                          #array(blkSize,4)
    1499         HP = H[:3]+modQ[:,nxs]*H[3:]            #projected hklm to hkl
    1500 #        HP = np.squeeze(np.inner(HP.T,TwinLaw))   #maybe array(blkSize,nTwins,4) or (blkSize,4)
    1501 #        TwMask = np.any(HP,axis=-1)
    1502 #        if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]
    1503 #            for ir in range(blkSize):
    1504 #                iref = ir+iBeg
    1505 #                if iref in TwDict:
    1506 #                    for i in TwDict[iref]:
    1507 #                        for n in range(NTL):
    1508 #                            HP[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
    1509 #            TwMask = np.any(HP,axis=-1)
     1613        H = refl[:,:4]                          #array(blkSize,4)
     1614        H3 = refl[:,:3]
     1615        HP = H[:,:3]+modQ[nxs,:]*H[:,3:]        #projected hklm to hkl
     1616        HP = np.inner(HP,TwinLaw)             #array(blkSize,nTwins,4)
     1617        H3 = np.inner(H3,TwinLaw)       
     1618        TwMask = np.any(HP,axis=-1)
     1619        if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]
     1620            for ir in range(blkSize):
     1621                iref = ir+iBeg
     1622                if iref in TwDict:
     1623                    for i in TwDict[iref]:
     1624                        for n in range(NTL):
     1625                            HP[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
     1626                            H3[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
     1627            TwMask = np.any(HP,axis=-1)
    15101628        SQ = 1./(2.*refl.T[5])**2               #array(blkSize)
    15111629        SQfactor = 4.0*SQ*twopisq               #ditto prev.
    1512         Uniq = np.inner(H.T,SSGMT)
    1513         UniqP = np.inner(HP.T,SGMT)
    1514         Phi = np.inner(H.T,SSGT)
     1630        Uniq = np.inner(H,SSGMT)
     1631        Uniq3 = np.inner(H3,SGMT)
     1632        UniqP = np.inner(HP,SGMT)
     1633        Phi = np.inner(H,SSGT)
    15151634        if SGInv:   #if centro - expand HKL sets
    15161635            Uniq = np.hstack((Uniq,-Uniq))
     1636            Uniq3 = np.hstack((Uniq3,-Uniq3))
    15171637            Phi = np.hstack((Phi,-Phi))
    15181638            UniqP = np.hstack((UniqP,-UniqP))
    … …  
    15271647        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
    15281648        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,Uniq.shape[1]*len(TwinLaw),axis=0)
    1529         phase = twopi*(np.inner(Uniq[:,:,:3],(dXdata.T+Xdata.T))-Phi[:,:,nxs])
     1649        phase = twopi*(np.inner(Uniq3,(dXdata.T+Xdata.T))-Phi[:,nxs,:,nxs])
    15301650        sinp = np.sin(phase)
    15311651        cosp = np.cos(phase)
    15321652        biso = -SQfactor*Uisodata[:,nxs]
    15331653        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[1]*len(TwinLaw),axis=1).T
    1534         HbH = -np.sum(UniqP[:,:,nxs,:]*np.inner(UniqP[:,:,:],bij),axis=-1)  #use hklt proj to hkl
     1654        HbH = -np.sum(UniqP[:,:,:,nxs]*np.inner(UniqP[:,:,:],bij),axis=-1)  #use hklt proj to hkl
    15351655        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
    15361656        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[1]  #refBlk x ops x atoms
     1657#        GSASIIpath.IPyBreak()
    15371658        if 'T' in calcControls[hfx+'histType']:
    15381659            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
    … …  
    15411662            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
    15421663            fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
    1543         GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x refBlk x sym X atoms
     1664        GfpuA = G2mth.ModulationTw(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x refBlk x sym X atoms
    15441665        fag = fa*GfpuA[0]-fb*GfpuA[1]   #real; 2 x refBlk x sym x atoms
    15451666        fbg = fb*GfpuA[0]+fa*GfpuA[1]
    15461667        fas = np.sum(np.sum(fag,axis=-1),axis=-1)   #2 x refBlk; sum sym & atoms
    15471668        fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
    1548 #        GSASIIpath.IPyBreak()
    1549         if 'P' in calcControls[hfx+'histType']:
    1550             refl.T[10] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2    #square of sums
    1551 #            refl.T[10] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0)
    1552             refl.T[11] = atan2d(fbs[0],fas[0])  #ignore f' & f"
    1553         else:
    1554             if len(TwinLaw) > 1:
    1555                 refl.T[10] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2                  #FcT from primary twin element
    1556                 refl.T[8] = np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fas,axis=0)**2,axis=-1)+   \
    1557                     np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fbs,axis=0)**2,axis=-1)                 #Fc sum over twins
    1558                 refl.T[11] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f"
    1559             else:
    1560                 refl.T[10] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2    #square of sums
    1561                 refl.T[8] = np.copy(refl.T[10])               
    1562                 refl.T[11] = atan2d(fbs[0],fas[0])  #ignore f' & f"
     1669        refl.T[10] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2                  #FcT from primary twin element
     1670        refl.T[8] = np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fas,axis=0)**2,axis=-1)+   \
     1671            np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fbs,axis=0)**2,axis=-1)                 #Fc sum over twins
     1672        refl.T[11] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f"
    15631673        iBeg += blkSize
    15641674    print 'nRef %d time %.4f\r'%(nRef,time.time()-time0)
    … …  
    16981808        dfadGu = np.sum(fa[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1)
    16991809        dfbdGu = np.sum(fb[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1)   
    1700 #        GSASIIpath.IPyBreak()
    17011810        if not SGData['SGInv']:   #Flack derivative
    17021811            dfadfl = np.sum(-FPP*Tcorr*sinp)
    … …  
    17051814            dfadfl = 1.0
    17061815            dfbdfl = 1.0
     1816#        GSASIIpath.IPyBreak()
    17071817        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
    1708         SA = fas[0]+fas[1]      #float = A+A' (might be array[nTwin])
    1709         SB = fbs[0]+fbs[1]      #float = B+B' (might be array[nTwin])
     1818        SA = fas[0]+fas[1]      #float = A+A'
     1819        SB = fbs[0]+fbs[1]      #float = B+B'
    17101820        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro?
    17111821            dFdfl[iref] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
    … …  
    17881898           
    17891899#        GSASIIpath.IPyBreak()
     1900    dFdvDict[phfx+'Flack'] = 4.*dFdfl.T
    17901901    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
    17911902    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
    … …  
    31203231            time0 = time.time()
    31213232            if im:
    3122                 SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
     3233                SStructureFactor(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
    31233234            else:
    31243235                StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
    … …  
    38283939            phfx = '%d:%d:'%(Phase['pId'],hId)
    38293940            SGData = Phase['General']['SGData']
     3941            TwinLaw = calcControls[phfx+'TwinLaw']
    38303942            im = 0
    38313943            if parmDict[phfx+'Scale'] < 0.:
    … …  
    38403952            time0 = time.time()
    38413953            if im:
    3842                 SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
     3954                if len(TwinLaw) > 1:
     3955                    SStructureFactorTw(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
     3956                else:
     3957                    SStructureFactor(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
    38433958            else:
    38443959                StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
Note: See TracChangeset for help on using the changeset viewer.