Changeset 550


Ignore:
Timestamp:
Apr 18, 2012 2:51:22 PM (10 years ago)
Author:
vondreele
Message:

charge flipping

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r537 r550  
    1010import os
    1111import os.path as ospath
     12import random as rn
    1213import numpy as np
    1314import numpy.linalg as nl
     
    1819import copy
    1920import GSASIIpath
     21import GSASIIElem as G2el
    2022import GSASIIlattice as G2lat
    2123import GSASIIspc as G2spc
     
    535537    for ref in reflData:
    536538        if ref[4] >= dmin:
     539            Fosq,Fcsq,ph = ref[8:11]
    537540            for i,hkl in enumerate(ref[11]):
    538541                hkl = np.asarray(hkl,dtype='i')
    539                 Fosq,Fcsq,ph = ref[8:11]
    540542                dp = 360.*ref[12][i]
    541543                a = cosd(ph+dp)
     
    576578    print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
    577579    mapData['rho'] = np.real(rho)
     580    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
     581    return mapData
     582
     583def ChargeFlip(data,reflData,pgbar):
     584    print 'charge flip'
     585    import scipy.fftpack as fft
     586#    import numpy.fft as fft
     587    generalData = data['General']
     588    mapData = generalData['Map']
     589    flipData = generalData['Flip']
     590    normElem = flipData['Norm element'].upper()
     591    FFtable = {}
     592    FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
     593    for ff in FFs:
     594        if ff['Symbol'] == normElem:
     595            FFtable.update(ff)
     596    dmin = flipData['Resolution']
     597    SGData = generalData['SGData']
     598    cell = generalData['Cell'][1:8]       
     599    A = G2lat.cell2A(cell[:6])
     600    Vol = cell[6]
     601    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
     602    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')
     603    time0 = time.time()
     604    for ref in reflData:
     605        dsp = ref[4]
     606        if dsp >= dmin:
     607            SQ = 0.25/dsp**2
     608            ff = 0.1*Vol*G2el.ScatFac(FFtable,SQ)[0]
     609            E = np.sqrt(ref[8])/ff
     610            ph = ref[10]
     611            if SGData['SGInv']:
     612                ph = rn.randint(0,1)*180.
     613            else:
     614                ph = rn.uniform(0.,360.)
     615            for i,hkl in enumerate(ref[11]):
     616                hkl = np.asarray(hkl,dtype='i')
     617                dp = 360.*ref[12][i]
     618                a = cosd(ph+dp)
     619                b = sind(ph+dp)
     620                phasep = complex(a,b)
     621                phasem = complex(a,-b)
     622                h,k,l = hkl+Hmax
     623                Ehkl[h,k,l] = E*phasep
     624                h,k,l = -hkl+Hmax       #Friedel pair
     625                Ehkl[h,k,l] = E*phasem
     626    CEhkl = copy.copy(Ehkl)
     627    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
     628    dmask = ma.getmask(MEhkl)
     629    sumE2 = np.sum(ma.array(np.absolute(CEhkl),mask=dmask)**2)
     630    Ncyc = 0
     631    while True:       
     632        CFrho = np.real(fft.fftn(fft.fftshift(CEhkl)))
     633        CFsig = np.std(CFrho)
     634        CErho = np.where(CFrho >= flipData['k-factor']*CFsig,CFrho,-CFrho)
     635        CFhkl = fft.ifftshift(fft.ifftn(CErho))
     636        CEhkl = np.absolute(Ehkl)*CFhkl/np.absolute(CFhkl)
     637        Ncyc += 1
     638        DEhkl = np.absolute(Ehkl)-np.absolute(CFhkl)
     639        sumDE2 = np.sum(ma.array(DEhkl,mask=dmask)**2)
     640        Rcf = min(100.,np.sqrt(sumDE2/sumE2)*100.)
     641        print '%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc)
     642        if Rcf < 5.:
     643            break
     644        GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
     645        if not GoOn:
     646            break
     647    print 'Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
     648    mapData['rho'] = np.real(fft.fftn(fft.fftshift(CEhkl)))
    578649    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
    579650    return mapData
Note: See TracChangeset for help on using the changeset viewer.