Changeset 2184 for trunk/GSASIIpwd.py


Ignore:
Timestamp:
Mar 25, 2016 10:36:47 AM (6 years ago)
Author:
vondreele
Message:

fix problems in adding atoms, etc for 'faulted' phases in the regular phase stuff

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIpwd.py

    r2183 r2184  
    19361936def CalcStackingSADP(Layers):
    19371937   
    1938     def getXY_HK(hk,Layers):
    1939         XY_HK = {}
    1940         for layer in Layers['Layers']:
    1941             if not layer['SameAs']:     #a real layer
    1942                 XY_HK[layer['Name']] = np.inner(hk,np.array([atom[2:4] for atom in layer['Atoms']]))*2.*np.pi
    1943         return XY_HK
     1938    atTypes = Layers['AtInfo'].keys()
     1939    Adat = []
     1940    for atType in atTypes:
     1941        if atType == 'H':
     1942            blen = -.3741
     1943        else:
     1944            blen = Layers['AtInfo'][atType]['Isotopes']['Nat. Abund.']['SL'][0]
     1945        Adat.append([[Adat['fa'][i],Adat['fb'][i]] for i in range(4)]+[Adat['fc'],])
     1946    Adat = np.array(Adat).flatten()
     1947    AtomXOU = []
     1948    AtomTp = []
     1949    AtomNL = []
     1950    LayerSymm = []
     1951    LayerNum = []
     1952    layerNames = []
     1953    Natm = 1
     1954    Nuniq = 0
     1955    for layer in Layers['Layers']:
     1956        layerNames.append(layer['Name'])
     1957    for il,layer in enumerate(Layers['Layers']):
     1958        if layer['SameAs']:
     1959            LayerNum.append(layerNames.index(layer['SameAs'])+1)
     1960            continue
     1961        else:
     1962            LayerNum.append(il)
     1963            Nuniq += 1
     1964        if '-1' in layer['Symm']:
     1965            LayerSymm.append(1)
     1966        else:
     1967            LayerSymm.append(0)
     1968        AtomNL.append(Natm)
     1969        for ia,atom in enumerate(layer['Atoms']):
     1970            [name,atype,x,y,z,frac,Uiso] = atom
     1971            Natm += 1
     1972            AtomTp.append(atype)
     1973            AtomXOU.append([x,y,z,frac,Uiso*78.9568])
     1974    TransX = []
     1975    TransP = []
     1976    for Ytrans in Layers['Transitions']:
     1977        TransP.append([trans[0] for trans in Ytrans])   #get just the numbers
     1978        TransX.append([trans[1:4] for trans in Ytrans])   #get just the numbers
     1979    TransP = np.array(Trans,dtype='float')
     1980    TransX = np.array(Trans,dtype='float')
     1981    Nlayers = np.sqrt(TransP.shape[0])
     1982    Cell = Layers['Cell'][1:4]+Layers['Cell'][6:7]
    19441983       
    1945     def getDXY_HK(hk,Trans,detune):
    1946         N = Trans.shape[0]
    1947         DXY_HK = np.zeros((N,N),dtype='cfloat')
    1948         Lphi = np.zeros((N,N),dtype='cfloat')
    1949         for iY in range(N):
    1950             for iX in range(N):
    1951                 if [str(iY+1),str(iX+1)] in Layers['allowedTrans']:
    1952                     dot = 2.*np.pi*np.inner(hk,Trans[iY,iX,1:3])
    1953                     Lphi[iY,iX] = complex(np.cos(dot),np.sin(dot))
    1954                     DXY_HK[iY][iX] = detune*Trans[iY,iX,0]*Lphi[iY,iX]
    1955         return DXY_HK,Lphi
    1956                    
    1957                    
    1958    
    1959     sadpSize = 256
    1960     G,g = G2lat.cell2Gmat(Layers['Cell'][1:7])  #recip/real met. tensors
    1961     A = G2lat.Gmat2A(G)
     1984    laueId = ['-1','2/m(ab)','2/m(c)','mmm','-3','-3m','4/m','4/mmm',
     1985        '6/m','6/mmm','axial','unknown'].index(Layers['Laue'])
     1986    planeId = ['h0l','0kl','hhl','h-hl'],index(Layers['Sadp']['Plane'])+1
    19621987    lmax = float(Layers['Sadp']['Lmax'])
    1963     Smax = G2lat.calc_rDsq([0.,0.,lmax],A)
    1964     plane = Layers['Sadp']['Plane']
    1965     if plane == 'h0l':
    1966         hkmax = int(lmax*np.sqrt(A[2]/A[0]))
    1967     elif plane == '0kl':
    1968         hkmax = int(lmax*np.sqrt(A[2]/A[1]))
    1969     elif plane == 'hhl':
    1970         hkmax = int(lmax*np.sqrt(A[2]/(A[0]+A[1]+A[5])))
    1971     elif plane == 'h-hl':
    1972         hkmax = int(lmax*np.sqrt(A[2]/(A[0]+A[1]-A[5])))
    1973     dl = 2.*lmax/sadpSize
    1974     Trans = []
    1975     for Ytrans in Layers['Transitions']:
    1976         Trans.append([trans[:4] for trans in Ytrans])   #get just the numbers
    1977     Trans = np.array(Trans,dtype='float')
    1978     N = np.sqrt(Trans.shape[0])
    1979     lmax -= 0.5*dl  #shift lmax to avoid l=0
    1980     lmin = -lmax
    1981     sadpBlock = sadpSize
    1982     if Layers['Laue'] not in ['-1','2/m(c)','-3','-3m','axial']:
    1983         lmin = -0.5*dl
    1984         sadpBlock /= 2
    1985 #start
    1986     cnt = 0
    1987     Names = [layer['Name'] for layer in Layers['Layers']]
    1988     Nlayers = len(Names)
    1989     detune = 1.0-0.001
    1990     for i in range(hkmax+1):
    1991         if plane == 'h0l':
    1992             hk = np.array([i,0])
    1993         elif plane == '0kl':
    1994             hk = np.array([0,i])
    1995         elif lane == 'hhl':
    1996             hk = np.array([i,i])
    1997         else:
    1998             hk = np.array([i,-i])
    1999         print ' h = %d k = %d'%(hk[0],hk[1])
    2000         XY_HK = getXY_HK(hk,Layers)                 #good
    2001         DXY_HK,Lphi = getDXY_HK(hk,Trans,detune)    #good
    2002              
     1988    controls = [laueId,planeId,lmax,Nuniq,]
    20031989       
    2004        
    2005    
     1990    Sadp = np.array(256**2)   
     1991    Sadp = pygetsadp(controls,len(AtTypes),AtTypes,Adat,Cell,Natm,AtomTp,AtomXOU,Nlayers,TransP,TransX,stackseq,Sadp)
    20061992   
    20071993   
Note: See TracChangeset for help on using the changeset viewer.