Changeset 1622 for trunk/GSASIImath.py


Ignore:
Timestamp:
Jan 5, 2015 1:02:33 PM (7 years ago)
Author:
vondreele
Message:

add 4D charge flipping - works to give 3D structure; only enabled for 4D data
fix problems with hkl limits in Fourier calcs. - see adjHKLmax
fix reflection table/2D & 3D display issues
disable MC/SA for 4D data

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r1604 r1622  
    14521452    '''
    14531453    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
    1454         Hmax[0] = ((Hmax[0]+3)/6)*6
    1455         Hmax[1] = ((Hmax[1]+3)/6)*6
    1456         Hmax[2] = ((Hmax[2]+1)/4)*4
     1454        Hmax[0] = int(math.ceil(Hmax[0]/6.))*6
     1455        Hmax[1] = int(math.ceil(Hmax[1]/6.))*6
     1456        Hmax[2] = int(math.ceil(Hmax[2]/4.))*4
    14571457    else:
    1458         Hmax[0] = ((Hmax[0]+2)/4)*4
    1459         Hmax[1] = ((Hmax[1]+2)/4)*4
    1460         Hmax[2] = ((Hmax[2]+2)/4)*4
     1458        Hmax[0] = int(math.ceil(Hmax[0]/4.))*4
     1459        Hmax[1] = int(math.ceil(Hmax[1]/4.))*4
     1460        Hmax[2] = int(math.ceil(Hmax[2]/4.))*4
    14611461
    14621462def OmitMap(data,reflDict,pgbar=None):
     
    15541554    time0 = time.time()
    15551555    for iref,ref in enumerate(reflDict['RefList']):
    1556         if ref[5] >= dmin:
     1556        if ref[5] > dmin:
    15571557            Fosq,Fcsq,ph = ref[9:12]
    15581558            Uniq = np.inner(ref[:4],SSGMT)
     
    16081608    time0 = time.time()
    16091609    for iref,ref in enumerate(reflDict['RefList']):
    1610         if ref[4] >= dmin:
     1610        if ref[4] > dmin:
    16111611            Fosq,Fcsq,ph = ref[8:11]
    16121612            Uniq = np.inner(ref[:3],SGMT)
     
    17851785    for iref,ref in enumerate(reflDict['RefList']):
    17861786        dsp = ref[4+im]
    1787         if im and ref[3]:   #skip super lattice reflections
     1787        if im and ref[3]:   #skip super lattice reflections - result is 3D projection
    17881788            continue
    1789         if dsp >= dmin:
     1789        if dsp > dmin:
    17901790            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
    17911791            if FFtable:
     
    18481848    return mapData
    18491849   
    1850 def SearchMap(data):
     1850def SSChargeFlip(data,reflDict,pgbar):
     1851    '''default doc string
     1852   
     1853    :param type name: description
     1854   
     1855    :returns: type name: description
     1856   
     1857    '''
     1858    generalData = data['General']
     1859    mapData = generalData['Map']
     1860    map4DData = {}
     1861    flipData = generalData['Flip']
     1862    FFtable = {}
     1863    if 'None' not in flipData['Norm element']:
     1864        normElem = flipData['Norm element'].upper()
     1865        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
     1866        for ff in FFs:
     1867            if ff['Symbol'] == normElem:
     1868                FFtable.update(ff)
     1869    dmin = flipData['Resolution']
     1870    SGData = generalData['SGData']
     1871    SSGData = generalData['SSGData']
     1872    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
     1873    SGT = np.array([ops[1] for ops in SGData['SGOps']])
     1874    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
     1875    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
     1876    cell = generalData['Cell'][1:8]       
     1877    A = G2lat.cell2A(cell[:6])
     1878    Vol = cell[6]
     1879    maxM = generalData['SuperVec'][2]
     1880    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A)+[maxM,],dtype='i')+1
     1881    adjHKLmax(SGData,Hmax)
     1882    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
     1883    time0 = time.time()
     1884    for iref,ref in enumerate(reflDict['RefList']):
     1885        dsp = ref[5]
     1886        if dsp > dmin:
     1887            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
     1888            if FFtable:
     1889                SQ = 0.25/dsp**2
     1890                ff *= G2el.ScatFac(FFtable,SQ)[0]
     1891            if ref[9] > 0.:         #use only +ve Fobs**2
     1892                E = np.sqrt(ref[9])/ff
     1893            else:
     1894                E = 0.
     1895            ph = ref[11]
     1896            ph = rn.uniform(0.,360.)
     1897            Uniq = np.inner(ref[:4],SSGMT)
     1898            Phi = np.inner(ref[:4],SSGT)
     1899            for i,hklm in enumerate(Uniq):        #uses uniq
     1900                hklm = np.asarray(hklm,dtype='i')
     1901                dp = 360.*Phi[i]                #and phi
     1902                a = cosd(ph+dp)
     1903                b = sind(ph+dp)
     1904                phasep = complex(a,b)
     1905                phasem = complex(a,-b)
     1906                h,k,l,m = hklm+Hmax
     1907                Ehkl[h,k,l,m] = E*phasep
     1908                h,k,l,m = -hklm+Hmax       #Friedel pair refl.
     1909                Ehkl[h,k,l,m] = E*phasem
     1910#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
     1911    CEhkl = copy.copy(Ehkl)
     1912    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
     1913    Emask = ma.getmask(MEhkl)
     1914    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
     1915    Ncyc = 0
     1916    old = np.seterr(all='raise')
     1917    while True:       
     1918        CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
     1919        CEsig = np.std(CErho)
     1920        CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
     1921        CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho)      #solves U atom problem!
     1922        CFhkl = fft.ifftshift(fft.ifftn(CFrho))
     1923        CFhkl = np.where(CFhkl,CFhkl,1.0)           #avoid divide by zero
     1924        phase = CFhkl/np.absolute(CFhkl)
     1925        CEhkl = np.absolute(Ehkl)*phase
     1926        Ncyc += 1
     1927        sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
     1928        DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
     1929        Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
     1930        if Rcf < 5.:
     1931            break
     1932        GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
     1933        if not GoOn or Ncyc > 10000:
     1934            break
     1935    np.seterr(**old)
     1936    print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
     1937    CErho = np.real(fft.fftn(fft.fftshift(CEhkl[:,:,:,maxM+1])))
     1938    SSrho = np.real(fft.fftn(fft.fftshift(CEhkl)))
     1939    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
     1940    roll = findOffset(SGData,A,CEhkl[:,:,:,maxM+1])               #CEhkl needs to be just the observed set, not the full set!
     1941
     1942    mapData['Rcf'] = Rcf
     1943    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
     1944    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
     1945    mapData['rollMap'] = [0,0,0]
     1946
     1947    map4DData['Rcf'] = Rcf
     1948    map4DData['rho'] = np.real(np.roll(np.roll(np.roll(SSrho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2))
     1949    map4DData['rhoMax'] = max(np.max(map4DData['rho']),-np.min(map4DData['rho']))
     1950    return mapData,map4DData
     1951   
     1952def SearchMap(generalData,drawingData):
    18511953    '''Does a search of a density map for peaks meeting the criterion of peak
    18521954    height is greater than mapData['cutOff']/100 of mapData['rhoMax'] where
     
    19172019        return Vec,Hess
    19182020       
    1919     generalData = data['General']
    19202021    phaseName = generalData['Name']
    19212022    SGData = generalData['SGData']
    19222023    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
    1923     drawingData = data['Drawing']
    19242024    peaks = []
    19252025    mags = []
Note: See TracChangeset for help on using the changeset viewer.