Changeset 2038


Ignore:
Timestamp:
Oct 30, 2015 4:02:02 PM (7 years ago)
Author:
vondreele
Message:

add betaij2Uij to G2lattice
revisions to SS structure factor calcs. Use hklt proj to hkl is several places
fxn works for thiourea derivs OK except X,Y,Zcos modulations; no Uijsin/cos derivatives yet
adj scaling of 4D charge flip maps
convert betaij vals from Jana2K files to Uij
start on SS read phase from cif
added a hklt F import (might vanish)

Location:
trunk
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIindex.py

    r2025 r2038  
    813813            if not GoOn:
    814814                break
    815         tries += 1   
     815        tries += 1
    816816    X = sortM20(Asave)
    817817    if X:

  • trunk/GSASIIlattice.py

    r2022 r2038  
    5050SQ2 = np.sqrt(2.)
    5151RSQPI = 1./np.sqrt(np.pi)
     52R2pisq = 1./(2.*np.pi**2)
    5253
    5354def sec2HMS(sec):
    … …  
    285286    return U6
    286287
     288def betaij2Uij(betaij,G):
     289    """
     290    Convert beta-ij to Uij tensors
     291   
     292    :param beta-ij - numpy array [beta-ij]
     293    :param G: reciprocal metric tensor
     294    :returns: Uij: numpy array [Uij]
     295    """
     296    ast = np.sqrt(np.diag(G))   #a*, b*, c*
     297    Mast = np.multiply.outer(ast,ast)   
     298    return R2pisq*UijtoU6(U6toUij(betaij)/Mast)
     299   
    287300def Uij2betaij(Uij,G):
    288301    """

  • trunk/GSASIImath.py

    r2025 r2038  
    927927################################################################################
    928928   
    929 def Modulation(waveTypes,H,FSSdata,XSSdata,USSdata,Mast):
     929def Modulation(waveTypes,H,HP,FSSdata,XSSdata,USSdata,Mast):
    930930    '''
    931931    H: array nRefBlk x ops X hklt
     932    HP: array nRefBlk x ops X hklt proj to hkl
    932933    FSSdata: array 2 x atoms x waves    (sin,cos terms)
    933934    XSSdata: array 2x3 x atoms X waves (sin,cos terms)
    … …  
    970971    XmodB = Bx[:,nx:,:,nxs]*np.cos(twopi*tauX[nxs,:,nxs,:]) #ditto
    971972    Xmod = np.sum(XmodA+XmodB+XmodZ,axis=1)                #atoms X 3 X 32; sum waves
    972     Tau = np.reshape(np.tile(glTau,Xmod.shape[0]),(Xmod.shape[0],-1))   #atoms x 32
    973     XmodT = np.swapaxes(np.concatenate((Xmod,Tau[:,nxs,:]),axis=1),1,2)      # atoms x 32 x 4 (x,y,z,t)
     973    Xmod = np.swapaxes(Xmod,1,2)                            #agree with J2K ParSup
    974974    if Au.shape[1]:
    975975        tauU = np.arange(1.,Au.shape[1]+1)[:,nxs]*glTau     #Uwaves x 32
    … …  
    977977        UmodB = Bu[:,:,:,:,nxs]*np.cos(twopi*tauU[nxs,:,nxs,nxs,:]) #ditto
    978978        Umod = np.swapaxes(np.sum(UmodA+UmodB,axis=1),1,3)      #atoms x 3x3 x 32; sum waves
    979         HbH = np.exp(-np.sum(H[:,:,nxs,nxs,:3]*np.inner(H[:,:,:3],Umod),axis=-1)) # refBlk x ops x atoms x 32 add Overhauser corr.?
     979        HbH = np.exp(-np.sum(HP[:,:,nxs,nxs]*np.inner(HP[:,:],Umod),axis=-1)) # refBlk x ops x atoms x 32 add Overhauser corr.?
    980980    else:
    981981        HbH = 1.0
    982     HdotX = np.inner(H,XmodT)     #refBlk X ops x atoms X 32
     982    D = H[:,:,3:]*glTau[nxs,nxs,:]              #m*e*tau; refBlk x  ops X 32
     983    HdotX = np.inner(HP,Xmod)+D[:,:,nxs,:]     #refBlk x ops x atoms X 32
    983984    cosHA = np.sum(Fmod*HbH*np.cos(twopi*HdotX)*glWt,axis=-1)       #real part; refBlk X ops x atoms; sum for G-L integration
    984985    sinHA = np.sum(Fmod*HbH*np.sin(twopi*HdotX)*glWt,axis=-1)       #imag part; ditto
    985986#    GSASIIpath.IPyBreak()
    986 #    if np.any(cosHA<0.) or np.any(cosHA>1.) or np.any(sinHA<0.) or np.any(sinHA>1.):
    987 #        GSASIIpath.IPyBreak()
    988987    return np.array([cosHA,sinHA])             # 2 x refBlk x SGops x atoms
    989988   
    990 def ModulationDerv(waveTypes,H,Hij,FSSdata,XSSdata,USSdata,Mast):
    991     '''
    992     H: array ops X hklt
     989def ModulationDerv(waveTypes,H,HP,Hij,FSSdata,XSSdata,USSdata,Mast):
     990    '''
     991    H: array ops X hklt proj to hkl
    993992    FSSdata: array 2 x atoms x waves    (sin,cos terms)
    994993    XSSdata: array 2x3 x atoms X waves (sin,cos terms)
    … …  
    999998    nxs = np.newaxis
    1000999    numeric = True   
    1001     cosHA,sinHA = Modulation(waveTypes,np.array([H,]),FSSdata,XSSdata,USSdata,Mast)
     1000    cosHA,sinHA = Modulation(waveTypes,np.array([H,]),np.array([HP,]),FSSdata,XSSdata,USSdata,Mast)
    10021001    Mf = [H.shape[0],]+list(FSSdata.T.shape)    #ops x atoms x waves x 2 (sin+cos frac mods)
    10031002    dGdMfC = np.zeros(Mf)
    … …  
    10371036    XmodB = Bx[:,nx:,:,nxs]*CtauX #ditto
    10381037    Xmod = np.sum(XmodA+XmodB+XmodZ,axis=1)                #atoms X pos X 32; sum waves
    1039     Tau = np.reshape(np.tile(glTau,Xmod.shape[0]),(Xmod.shape[0],-1))   #atoms x 32
    1040     XmodT = np.swapaxes(np.concatenate((Xmod,Tau[:,nxs,:]),axis=1),1,2)      # atoms x 32 x 4 (x,y,z,t)
    1041     D = H[:,3][:,nxs]*glTau[nxs,:]              #m*tau; ops X 32
    1042 #    HdotX = np.inner(H[:,:3],np.swapaxes(Xmod,1,2))+D[:,nxs,:]     #ops x atoms X 32
    1043     HdotX = np.inner(H,XmodT)     #refBlk X ops x atoms X 32
     1038    Xmod = np.swapaxes(Xmod,1,2)
     1039    D = H[:,3][:,nxs]*glTau[nxs,:]              #m*e*tau; ops X 32
     1040    HdotX = np.inner(HP,Xmod)+D[:,nxs,:]     #refBlk X ops x atoms X 32
    10441041    if Af.shape[1]:
    10451042        tauF = np.arange(1.,Af.shape[1]+1-nf)[:,nxs]*glTau  #Fwaves x 32
    … …  
    10581055        UmodB = Bu[:,:,:,:,nxs]*CtauU #ditto
    10591056        Umod = np.swapaxes(np.sum(UmodA+UmodB,axis=1),1,3)      #atoms x 3x3 x 32; sum waves
    1060         HbH = np.exp(-np.sum(H[:,nxs,nxs,:3]*np.inner(H[:,:3],Umod),axis=-1)) # ops x atoms x 32 add Overhauser corr.?
    1061        
     1057        HbH = np.exp(-np.sum(HP[:,nxs,nxs]*np.inner(HP[:],Umod),axis=-1)) # ops x atoms x 32 add Overhauser corr.?
     1058        #derivatives??
    10621059    else:
    10631060        HbH = np.ones_like(HdotX)
    1064     HdotXA = H[:,nxs,nxs,nxs,:3]*np.swapaxes(XmodA,-1,-2)[nxs,:,:,:,:]+D[:,nxs,nxs,:,nxs]  #ops x atoms x waves x 32 x xyz
    1065     HdotXB = H[:,nxs,nxs,nxs,:3]*np.swapaxes(XmodB,-1,-2)[nxs,:,:,:,:]+D[:,nxs,nxs,:,nxs]
    1066     dHdXA = H[:,nxs,nxs,nxs,:3]*np.swapaxes(StauX,-1,-2)[nxs,:,:,:,:]    #ops x atoms x waves x 32 x xyz
    1067     dHdXB = H[:,nxs,nxs,nxs,:3]*np.swapaxes(CtauX,-1,-2)[nxs,:,:,:,:]              #ditto
     1061    HdotXA = HP[:,nxs,nxs,nxs,:]*np.swapaxes(XmodA,-1,-2)[nxs,:,:,:,:]+D[:,nxs,nxs,:,nxs]  #ops x atoms x waves x 32 x xyz
     1062    HdotXB = HP[:,nxs,nxs,nxs,:]*np.swapaxes(XmodB,-1,-2)[nxs,:,:,:,:]+D[:,nxs,nxs,:,nxs]
     1063    dHdXA = HP[:,nxs,nxs,nxs,:]*np.swapaxes(StauX,-1,-2)[nxs,:,:,:,:]    #ops x atoms x waves x 32 x xyz
     1064    dHdXB = HP[:,nxs,nxs,nxs,:]*np.swapaxes(CtauX,-1,-2)[nxs,:,:,:,:]              #ditto
    10681065    dGdMxCa = np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(-twopi*dHdXA*np.sin(twopi*HdotXA))*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
    10691066    dGdMxCb = np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(-twopi*dHdXB*np.sin(twopi*HdotXB))*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
    … …  
    22522249    return mapData
    22532250   
    2254 def Fourier4DMap(data,reflDict):
    2255     '''default doc string
    2256    
    2257     :param type name: description
    2258    
    2259     :returns: type name: description
    2260    
    2261     '''
    2262     generalData = data['General']
    2263     map4DData = generalData['4DmapData']
    2264     mapData = generalData['Map']
    2265     dmin = mapData['Resolution']
    2266     SGData = generalData['SGData']
    2267     SSGData = generalData['SSGData']
    2268     SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
    2269     SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
    2270     cell = generalData['Cell'][1:8]       
    2271     A = G2lat.cell2A(cell[:6])
    2272     maxM = generalData['SuperVec'][2]
    2273     Hmax = G2lat.getHKLmax(dmin,SGData,A)+[maxM,]
    2274     adjHKLmax(SGData,Hmax)
    2275     Hmax = np.asarray(Hmax,dtype='i')+1
    2276     Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
    2277     time0 = time.time()
    2278     for iref,ref in enumerate(reflDict['RefList']):
    2279         if ref[5] > dmin:
    2280             Fosq,Fcsq,ph = ref[9:12]
    2281             Fosq = np.where(Fosq>0.,Fosq,0.)    #can't use Fo^2 < 0
    2282             Uniq = np.inner(ref[:4],SSGMT)
    2283             Phi = np.inner(ref[:4],SSGT)
    2284             for i,hkl in enumerate(Uniq):        #uses uniq
    2285                 hkl = np.asarray(hkl,dtype='i')
    2286                 dp = 360.*Phi[i]                #and phi
    2287                 a = cosd(ph+dp)
    2288                 b = sind(ph+dp)
    2289                 phasep = complex(a,b)
    2290                 phasem = complex(a,-b)
    2291                 if 'Fobs' in mapData['MapType']:
    2292                     F = np.sqrt(Fosq)
    2293                     h,k,l,m = hkl+Hmax
    2294                     Fhkl[h,k,l,m] = F*phasep
    2295                     h,k,l,m = -hkl+Hmax
    2296                     Fhkl[h,k,l,m] = F*phasem
    2297                 elif 'Fcalc' in mapData['MapType']:
    2298                     F = np.sqrt(Fcsq)
    2299                     h,k,l,m = hkl+Hmax
    2300                     Fhkl[h,k,l,m] = F*phasep
    2301                     h,k,l,m = -hkl+Hmax
    2302                     Fhkl[h,k,l,m] = F*phasem                   
    2303                 elif 'delt-F' in mapData['MapType']:
    2304                     dF = np.sqrt(Fosq)-np.sqrt(Fcsq)
    2305                     h,k,l,m = hkl+Hmax
    2306                     Fhkl[h,k,l,m] = dF*phasep
    2307                     h,k,l,m = -hkl+Hmax
    2308                     Fhkl[h,k,l,m] = dF*phasem
    2309     SSrho = fft.fftn(fft.fftshift(Fhkl))/cell[6]          #4D map
    2310     rho = fft.fftn(fft.fftshift(Fhkl[:,:,:,maxM+1]))/cell[6]    #3D map
    2311     map4DData['rho'] = np.real(SSrho)
    2312     map4DData['rhoMax'] = max(np.max(map4DData['rho']),-np.min(map4DData['rho']))
    2313     map4DData['minmax'] = [np.max(map4DData['rho']),np.min(map4DData['rho'])]
    2314     map4DData['Type'] = reflDict['Type']
    2315     mapData['Type'] = reflDict['Type']
    2316     mapData['rho'] = np.real(rho)
    2317     mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
    2318     mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
    2319     print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
    2320 
    23212251def FourierMap(data,reflDict):
    23222252    '''default doc string
    … …  
    23882318    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
    23892319   
     2320def Fourier4DMap(data,reflDict):
     2321    '''default doc string
     2322   
     2323    :param type name: description
     2324   
     2325    :returns: type name: description
     2326   
     2327    '''
     2328    generalData = data['General']
     2329    map4DData = generalData['4DmapData']
     2330    mapData = generalData['Map']
     2331    dmin = mapData['Resolution']
     2332    SGData = generalData['SGData']
     2333    SSGData = generalData['SSGData']
     2334    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
     2335    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
     2336    cell = generalData['Cell'][1:8]       
     2337    A = G2lat.cell2A(cell[:6])
     2338    maxM = generalData['SuperVec'][2]
     2339    Hmax = G2lat.getHKLmax(dmin,SGData,A)+[maxM,]
     2340    adjHKLmax(SGData,Hmax)
     2341    Hmax = np.asarray(Hmax,dtype='i')+1
     2342    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
     2343    time0 = time.time()
     2344    for iref,ref in enumerate(reflDict['RefList']):
     2345        if ref[5] > dmin:
     2346            Fosq,Fcsq,ph = ref[9:12]
     2347            Fosq = np.where(Fosq>0.,Fosq,0.)    #can't use Fo^2 < 0
     2348            Uniq = np.inner(ref[:4],SSGMT)
     2349            Phi = np.inner(ref[:4],SSGT)
     2350            for i,hkl in enumerate(Uniq):        #uses uniq
     2351                hkl = np.asarray(hkl,dtype='i')
     2352                dp = 360.*Phi[i]                #and phi
     2353                a = cosd(ph+dp)
     2354                b = sind(ph+dp)
     2355                phasep = complex(a,b)
     2356                phasem = complex(a,-b)
     2357                if 'Fobs' in mapData['MapType']:
     2358                    F = np.sqrt(Fosq)
     2359                    h,k,l,m = hkl+Hmax
     2360                    Fhkl[h,k,l,m] = F*phasep
     2361                    h,k,l,m = -hkl+Hmax
     2362                    Fhkl[h,k,l,m] = F*phasem
     2363                elif 'Fcalc' in mapData['MapType']:
     2364                    F = np.sqrt(Fcsq)
     2365                    h,k,l,m = hkl+Hmax
     2366                    Fhkl[h,k,l,m] = F*phasep
     2367                    h,k,l,m = -hkl+Hmax
     2368                    Fhkl[h,k,l,m] = F*phasem                   
     2369                elif 'delt-F' in mapData['MapType']:
     2370                    dF = np.sqrt(Fosq)-np.sqrt(Fcsq)
     2371                    h,k,l,m = hkl+Hmax
     2372                    Fhkl[h,k,l,m] = dF*phasep
     2373                    h,k,l,m = -hkl+Hmax
     2374                    Fhkl[h,k,l,m] = dF*phasem
     2375    SSrho = fft.fftn(fft.fftshift(Fhkl))/cell[6]          #4D map
     2376    rho = fft.fftn(fft.fftshift(Fhkl[:,:,:,maxM+1]))/cell[6]    #3D map
     2377    map4DData['rho'] = np.real(SSrho)
     2378    map4DData['rhoMax'] = max(np.max(map4DData['rho']),-np.min(map4DData['rho']))
     2379    map4DData['minmax'] = [np.max(map4DData['rho']),np.min(map4DData['rho'])]
     2380    map4DData['Type'] = reflDict['Type']
     2381    mapData['Type'] = reflDict['Type']
     2382    mapData['rho'] = np.real(rho)
     2383    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
     2384    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
     2385    print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
     2386
    23902387# map printing for testing purposes
    23912388def printRho(SGLaue,rho,rhoMax):                         
    … …  
    25742571    np.seterr(**old)
    25752572    print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
    2576     CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))
     2573    CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))/10.  #? to get on same scale as e-map
    25772574    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
    25782575    roll = findOffset(SGData,A,CEhkl)               #CEhkl needs to be just the observed set, not the full set!
    … …  
    27362733    np.seterr(**old)
    27372734    print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
    2738     CErho = np.real(fft.fftn(fft.fftshift(CEhkl[:,:,:,maxM+1])))
    2739     SSrho = np.real(fft.fftn(fft.fftshift(CEhkl)))
     2735    CErho = np.real(fft.fftn(fft.fftshift(CEhkl[:,:,:,maxM+1])))/10.    #? to get on same scale as e-map
     2736    SSrho = np.real(fft.fftn(fft.fftshift(CEhkl)))/10.                  #? ditto
    27402737    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
    27412738    roll = findSSOffset(SGData,SSGData,A,CEhkl)               #CEhkl needs to be just the observed set, not the full set!

  • trunk/GSASIIstrMath.py

    r2024 r2038  
    974974    Mast = twopisq*np.multiply.outer(ast,ast)   
    975975    SGInv = SGData['SGInv']
     976    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
     977    SGT = np.array([ops[1] for ops in SGData['SGOps']])
    976978    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
    977979    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
    … …  
    10041006        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
    10051007        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
    1006     Uij = np.array(G2lat.U6toUij(Uijdata))
    1007     bij = Mast*Uij.T
     1008    Uij = np.array(G2lat.U6toUij(Uijdata)).T
     1009    bij = Mast*Uij
    10081010    blkSize = 100       #no. of reflections in a block
    10091011    nRef = refDict['RefList'].shape[0]
    … …  
    10281030        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
    10291031        H = refl.T[:4]                          #array(blkSize,4)
     1032        HP = H[:3]+modQ[:,nxs]*H[3:]            #projected hklm to hkl
    10301033#        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(blkSize,nTwins,4) or (blkSize,4)
    10311034#        TwMask = np.any(H,axis=-1)
    … …  
    10411044        SQfactor = 4.0*SQ*twopisq               #ditto prev.
    10421045        Uniq = np.inner(H.T,SSGMT)
     1046        UniqP = np.inner(HP.T,SGMT)
    10431047        Phi = np.inner(H.T,SSGT)
     1048        PhiP = np.inner(HP.T,SGT)
    10441049        if SGInv:   #if centro - expand HKL sets
    10451050            Uniq = np.hstack((Uniq,-Uniq))
    10461051            Phi = np.hstack((Phi,-Phi))
     1052            UniqP = np.hstack((UniqP,-UniqP))
     1053            PhiP = np.hstack((PhiP,-PhiP))
    10471054        if 'T' in calcControls[hfx+'histType']:
    10481055            if 'P' in calcControls[hfx+'histType']:
    … …  
    10551062        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
    10561063        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,Uniq.shape[1]*len(TwinLaw),axis=0)
    1057         phase = twopi*(np.inner(Uniq[:,:,:3],(dXdata.T+Xdata.T))+Phi[:,:,nxs]) #+TauT.T[nxs,:,:]*eps[nxs,:,nxs]
     1064        phase = twopi*(np.inner(Uniq[:,:,:3],(dXdata.T+Xdata.T))+Phi[:,:,nxs])
    10581065        sinp = np.sin(phase)
    10591066        cosp = np.cos(phase)
    10601067        biso = -SQfactor*Uisodata[:,nxs]
    10611068        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[1]*len(TwinLaw),axis=1).T
    1062         HbH = -np.sum(Uniq[:,:,nxs,:3]*np.inner(Uniq[:,:,:3],bij),axis=-1)
     1069        HbH = -np.sum(UniqP[:,:,nxs,:]*np.inner(UniqP[:,:,:],bij),axis=-1)  #use hklt proj to hkl
    10631070        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
    1064         Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[1]
     1071        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[1]  #refBlk x ops x atoms
    10651072        if 'T' in calcControls[hfx+'histType']:
    10661073            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
    … …  
    10691076            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
    10701077            fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
    1071         GfpuA = G2mth.Modulation(waveTypes,Uniq,FSSdata,XSSdata,USSdata,Mast) #2 x refBlk x sym X atoms
     1078        GfpuA = G2mth.Modulation(waveTypes,Uniq,UniqP,FSSdata,XSSdata,USSdata,Mast) #2 x refBlk x sym X atoms
    10721079        fag = fa*GfpuA[0]-fb*GfpuA[1]   #real; 2 x refBlk x sym x atoms
    10731080        fbg = fb*GfpuA[0]+fa*GfpuA[1]
    10741081        fas = np.sum(np.sum(fag,axis=-1),axis=-1)   #2 x refBlk; sum sym & atoms
    10751082        fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
    1076    #     GSASIIpath.IPyBreak()
     1083#        GSASIIpath.IPyBreak()
    10771084        if 'P' in calcControls[hfx+'histType']:
    10781085            refl.T[10] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0)
    … …  
    10981105    Mast = twopisq*np.multiply.outer(ast,ast)
    10991106    SGInv = SGData['SGInv']
     1107    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
     1108    SGT = np.array([ops[1] for ops in SGData['SGOps']])
    11001109    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
    11011110    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
    … …  
    11181127    if SGInv:
    11191128        TauT = np.hstack((TauT,-TauT))
     1129    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
    11201130    FF = np.zeros(len(Tdata))
    11211131    if 'NC' in calcControls[hfx+'histType']:
    … …  
    11241134        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
    11251135        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
    1126     Uij = np.array(G2lat.U6toUij(Uijdata))
    1127     bij = Mast*Uij.T
     1136    Uij = np.array(G2lat.U6toUij(Uijdata)).T
     1137    bij = Mast*Uij
    11281138    if not len(refDict['FF']):
    11291139        if 'N' in calcControls[hfx+'histType']:
    … …  
    11681178            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12+im])
    11691179        H = np.array(refl[:4])
     1180        HP = H[:3]+modQ*H[3:]            #projected hklm to hkl
    11701181#        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(4,nTwins) or (4)
    11711182#        TwMask = np.any(H,axis=-1)
    … …  
    11841195        Uniq = np.inner(H,SSGMT)
    11851196        Phi = np.inner(H,SSGT)
     1197        UniqP = np.inner(HP,SGMT)
     1198        PhiP = np.inner(HP,SGT)
    11861199        if SGInv:   #if centro - expand HKL sets
    11871200            Uniq = np.vstack((Uniq,-Uniq))
    11881201            Phi = np.hstack((Phi,-Phi))
     1202            UniqP = np.vstack((UniqP,-UniqP))
     1203            PhiP = np.hstack((PhiP,-PhiP))
     1204#        GSASIIpath.IPyBreak()                 
    11891205        phase = twopi*(np.inner(Uniq[:,:3],(dXdata+Xdata).T)+Phi[:,nxs])
    11901206        sinp = np.sin(phase)
    … …  
    11931209        biso = -SQfactor*Uisodata[:,nxs]
    11941210        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[0]*len(TwinLaw),axis=1).T    #ops x atoms
    1195         HbH = -np.sum(Uniq[:,nxs,:3]*np.inner(Uniq[:,:3],bij),axis=-1)  #ops x atoms
    1196         Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in Uniq]) #atoms x 3x3
     1211        HbH = -np.sum(UniqP[:,nxs,:3]*np.inner(UniqP[:,:3],bij),axis=-1)  #ops x atoms
     1212        Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in UniqP]) #atoms x 3x3
    11971213        Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6)))
    11981214        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)     #ops x atoms
    … …  
    12001216        fot = (FF+FP-Bab)*Tcorr     #ops x atoms
    12011217        fotp = FPP*Tcorr            #ops x atoms
    1202         GfpuA,dGdf,dGdx,dGdu = G2mth.ModulationDerv(waveTypes,Uniq,Hij,FSSdata,XSSdata,USSdata,Mast)
     1218        GfpuA,dGdf,dGdx,dGdu = G2mth.ModulationDerv(waveTypes,Uniq,UniqP,Hij,FSSdata,XSSdata,USSdata,Mast)
    12031219        # derivs are: ops x atoms x waves x 1,3,or 6 parms as [real,imag] parts
    12041220        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) # array(2,nTwin,nEqv,nAtoms)
    … …  
    12201236        dfadui = np.sum(-SQfactor*fag,axis=1)
    12211237        dfbdui = np.sum(-SQfactor*fbg,axis=1)
    1222 #        GSASIIpath.IPyBreak()                 
    12231238        if nTwin > 1:
    12241239            dfadx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
    … …  
    12291244            dfadGx = np.sum(fa[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1)
    12301245            dfbdGx = np.sum(fb[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1)
     1246            dfadGx = np.sum(fa[:,it,:,:,nxs,nxs]*dGdu[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdu[1][nxs,nxs,:,:,:,:],axis=1)
     1247            dfbdGx = np.sum(fb[:,it,:,:,nxs,nxs]*dGdu[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdu[1][nxs,nxs,:,:,:,:],axis=1)
    12311248            # array (2,nAtom,wave,3)
    12321249        else:
    … …  
    12381255            dfadGx = np.sum(fa[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1)
    12391256            dfbdGx = np.sum(fb[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1)
     1257            dfadGu = np.sum(fa[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1)
     1258            dfbdGu = np.sum(fb[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1)   
    12401259            # array (2,nAtom,wave,3)
    12411260        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for al2O3!   
    … …  
    12571276            dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])+   \
    12581277                2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
     1278            dFdGx[iref] = 2.*(fas[0]*dfadGx[0]+fas[1]*dfadGx[1])+  \
     1279                2.*(fbs[0]*dfbdGx[0]+fbs[1]*dfbdGx[1])
     1280            dFdGu[iref] = 2.*(fas[0]*dfadGu[0]+fas[1]*dfadGu[1])+  \
     1281                2.*(fbs[0]*dfbdGu[0]+fbs[1]*dfbdGu[1])
    12591282        else:
    12601283            SA = fas[0]-fbs[1]
    … …  
    12661289                dFdua[iref] = [2.*TwMask[it]*(fas[0][it]*dfadua[it][0]+fas[1][it]*dfadua[it][1]+fbs[0][it]*dfbdua[it][0]+fbs[1][it]*dfbdua[it][1]) for it in range(nTwin)]
    12671290                dFdGx[iref] = [2.*TwMask[it]*(fas[0][it]*dfadGx[0]+fas[1][it]*dfadGx[1]+fbs[0][it]*dfbdGx[0]+fbs[1][it]*dfbdGx[1]) for it in range(nTwin)]
     1291                dFdGu[iref] = [2.*TwMask[it]*(fas[0][it]*dfadGu[0]+fas[1][it]*dfadGu[1]+fbs[0][it]*dfbdGu[0]+fbs[1][it]*dfbdGu[1]) for it in range(nTwin)]
    12681292                dFdtw[iref] = np.sum(TwMask*fas,axis=0)**2+np.sum(TwMask*fbs,axis=0)**2
    12691293               
    … …  
    12751299                dFdfl[iref] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
    12761300                dFdGx[iref] = 2.*(fas[0]*dfadGx[0]+fas[1]*dfadGx[1]+fbs[0]*dfbdGx[0]+fbs[1]*dfbdGx[1])      #array(nRef,natom,nwave,6)
     1301                dFdGu[iref] = 2.*(fas[0]*dfadGu[0]+fas[1]*dfadGu[1]+fbs[0]*dfbdGu[0]+fbs[1]*dfbdGu[1])      #array(nRef,natom,nwave,6)
    12771302        dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
    12781303            2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
    … …  
    12991324            dFdvDict[pfx+'Ycos:'+str(i)+':'+str(j)] = dFdGx.T[4][j][i]
    13001325            dFdvDict[pfx+'Zcos:'+str(i)+':'+str(j)] = dFdGx.T[5][j][i]
     1326        for j in range(USSdata.shape[1]):
     1327            dFdvDict[pfx+'U11sin:'+str(i)+':'+str(j)] = dFdGu.T[0][j][i]
     1328            dFdvDict[pfx+'U22sin:'+str(i)+':'+str(j)] = dFdGu.T[1][j][i]
     1329            dFdvDict[pfx+'U33sin:'+str(i)+':'+str(j)] = dFdGu.T[2][j][i]
     1330            dFdvDict[pfx+'U12sin:'+str(i)+':'+str(j)] = dFdGu.T[3][j][i]
     1331            dFdvDict[pfx+'U13sin:'+str(i)+':'+str(j)] = dFdGu.T[4][j][i]
     1332            dFdvDict[pfx+'U23sin:'+str(i)+':'+str(j)] = dFdGu.T[5][j][i]
     1333            dFdvDict[pfx+'U11cos:'+str(i)+':'+str(j)] = dFdGu.T[6][j][i]
     1334            dFdvDict[pfx+'U22cos:'+str(i)+':'+str(j)] = dFdGu.T[7][j][i]
     1335            dFdvDict[pfx+'U33cos:'+str(i)+':'+str(j)] = dFdGu.T[8][j][i]
     1336            dFdvDict[pfx+'U12cos:'+str(i)+':'+str(j)] = dFdGu.T[9][j][i]
     1337            dFdvDict[pfx+'U13cos:'+str(i)+':'+str(j)] = dFdGu.T[10][j][i]
     1338            dFdvDict[pfx+'U23cos:'+str(i)+':'+str(j)] = dFdGu.T[11][j][i]
     1339           
    13011340    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
    13021341    dFdvDict[phfx+'BabU'] = dFdbab.T[1]

  • trunk/imports/G2phase.py

    r1762 r2038  
    2929import GSASIIpath
    3030GSASIIpath.SetVersionNumber("$Revision$")
     31R2pisq = 1./(2.*np.pi**2)
    3132
    3233class PDB_ReaderClass(G2IO.ImportPhase):
    … …  
    402403                    float(cell[3]),float(cell[4]),float(cell[5])]
    403404                Volume = G2lat.calc_V(G2lat.cell2A(cell))
     405                G,g = G2lat.cell2Gmat(cell)
     406                ast = np.sqrt(np.diag(G))
     407                Mast = np.multiply.outer(ast,ast)   
     408               
    404409            elif 'spgroup' in S:
    405410                if 'X' in S:
    … …  
    479484                Uiso = 0.
    480485                Uij = [float(S2[:9]),float(S2[9:18]),float(S2[18:27]),
    481                     float(S2[27:36]),float(S2[36:45]),float(S2[45:54])]
     486                    float(S2[27:36]),float(S2[36:45]),float(S2[45:54])] #these things are betaij! need to convert to Uij
     487                Uij = R2pisq*G2lat.UijtoU6(G2lat.U6toUij(Uij)/Mast)
    482488            for i in range(S1N[0]):
    483489                if not i:
    … …  
    516522                Spos[i] = [vals,False]
    517523            for i,it in enumerate(Sadp):
     524                #these are betaij modulations! need to convert to Uij modulations
    518525                vals = [float(it[:9]),float(it[9:18]),float(it[18:27]),float(it[27:36]),float(it[36:45]),float(it[45:54]),
    519526                    float(it[54:63]),float(it[63:72]),float(it[72:81]),float(it[81:90]),float(it[90:99]),float(it[99:108])]               
     527                vals[:6] = R2pisq*G2lat.UijtoU6(G2lat.U6toUij(vals[:6])/Mast)    #convert sin bij to Uij
     528                vals[6:] = R2pisq*G2lat.UijtoU6(G2lat.U6toUij(vals[6:])/Mast)    #convert cos bij to Uij
    520529                Sadp[i] = [vals,False]
    521530            Atom = [Name,aType,'',XYZ[0],XYZ[1],XYZ[2],1.0,SytSym,Mult,IA,Uiso]

  • trunk/imports/G2phase_CIF.py

    r1985 r2038  
    196196                                '_atom_site_aniso_u_13' : 15,
    197197                                '_atom_site_aniso_u_23' : 16, }
     198                G2SSAtomDict = {'_atom_site_occ_Fourier_wave_vector_seq_id' : 1,
     199                                '_atom_site_occ_Fourier_param_cos' : 2,
     200                                '_atom_site_occ_Fourier_param_sin' : 3,
     201                                '_atom_site_displace_Fourier_axis' : 2,
     202                                '_atom_site_displace_Fourier_wave_vector_seq_id' : 3,
     203                                '_atom_site_displace_Fourier_param_cos' : 5,
     204                                '_atom_site_displace_Fourier_param_sin' : 4,
     205                                '_atom_site_U_Fourier_tens_elem' : 7,
     206                                '_atom_site_U_Fourier_wave_vector_seq_id' : 8,
     207                                '_atom_site_U_Fourier_param_cos' : 9,
     208                                '_atom_site_U_Fourier_param_sin' : 10,  }
     209
     210                                   
     211
    198212                ranIdlookup = {}
    199213                for aitem in atomloop:

  • trunk/imports/G2sfact.py

    r1924 r2038  
    7676            self.RefDict['Type'] = 'SXC'
    7777            self.RefDict['Super'] = 0
     78            self.UpdateParameters(Type='SXC',Wave=None) # histogram type
     79            return True
     80        except Exception as detail:
     81            self.errors += '\n  '+str(detail)
     82            print '\n\n'+self.formatName+' read error: '+str(detail) # for testing
     83            import traceback
     84            traceback.print_exc(file=sys.stdout)
     85            return False
     86
     87class HKLMF_ReaderClass(G2IO.ImportStructFactor):
     88    'Routines to import F, reflections from a REMOS HKLMF file'
     89    def __init__(self):
     90        super(self.__class__,self).__init__( # fancy way to self-reference
     91            extensionlist=('.fo','.FO'),
     92            strictExtension=False,
     93            formatName = 'HKLM F',
     94            longFormatName = 'REMOS [hklm, Fo] Structure factor text file'
     95            )
     96
     97    def ContentsValidator(self, filepointer):
     98        'Make sure file contains the expected columns on numbers'
     99        return ColumnValidator(self, filepointer)
     100
     101    def Reader(self,filename,filepointer, ParentFrame=None, **unused):
     102        'Read the file'
     103        try:
     104            for line,S in enumerate(filepointer):
     105                self.errors = '  Error reading line '+str(line+1)
     106                if S[0] == '#': continue       #ignore comments, if any
     107                h,k,l,m,Fo= S.split()
     108                h,k,l,m = [int(h),int(k),int(l),int(m)]
     109                if h == 999 or not any([h,k,l]):
     110                    break
     111                Fo = float(Fo)
     112                sigFo2 = Fo
     113                if Fo < 1.0:
     114                    sigFo2 = 1.0
     115               # h,k,l,m,tw,dsp,Fo2,sig,Fc2,Fot2,Fct2,phase,...
     116                self.RefDict['RefList'].append([h,k,l,m,1,0,Fo**2,sigFo2,0,Fo**2,0,0,0])
     117            self.errors = 'Error after reading reflections (unexpected!)'
     118            self.RefDict['RefList'] = np.array(self.RefDict['RefList'])
     119            self.RefDict['Type'] = 'SXC'
     120            self.RefDict['Super'] = 1
    78121            self.UpdateParameters(Type='SXC',Wave=None) # histogram type
    79122            return True
Note: See TracChangeset for help on using the changeset viewer.