Changeset 550
 Timestamp:
 Apr 18, 2012 2:51:22 PM (11 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIImath.py
r537 r550 10 10 import os 11 11 import os.path as ospath 12 import random as rn 12 13 import numpy as np 13 14 import numpy.linalg as nl … … 18 19 import copy 19 20 import GSASIIpath 21 import GSASIIElem as G2el 20 22 import GSASIIlattice as G2lat 21 23 import GSASIIspc as G2spc … … 535 537 for ref in reflData: 536 538 if ref[4] >= dmin: 539 Fosq,Fcsq,ph = ref[8:11] 537 540 for i,hkl in enumerate(ref[11]): 538 541 hkl = np.asarray(hkl,dtype='i') 539 Fosq,Fcsq,ph = ref[8:11]540 542 dp = 360.*ref[12][i] 541 543 a = cosd(ph+dp) … … 576 578 print 'Fourier map time: %.4f'%(time.time()time0),'no. elements: %d'%(Fhkl.size) 577 579 mapData['rho'] = np.real(rho) 580 mapData['rhoMax'] = max(np.max(mapData['rho']),np.min(mapData['rho'])) 581 return mapData 582 583 def 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['kfactor']*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))) 578 649 mapData['rhoMax'] = max(np.max(mapData['rho']),np.min(mapData['rho'])) 579 650 return mapData
Note: See TracChangeset
for help on using the changeset viewer.