source: trunk/GSASIIindex.py @ 303

Last change on this file since 303 was 303, checked in by vondreele, 12 years ago

Allow reading of xye (Topas style) files
Begin implementation of spherical harmonics texture
Refactor indexing; replace cell refinement from peak positions
  • Property svn:keywords set to Date Author Revision URL Id
File size: 18.7 KB
Line 
1#GSASII cell indexing program: variation on that of A. Coehlo
2#   includes cell refinement from peak positions (not zero as yet)
3########### SVN repository information ###################
4# $Date: 2011-06-20 19:37:07 +0000 (Mon, 20 Jun 2011) $
5# $Author: vondreele $
6# $Revision: 303 $
7# $URL: trunk/GSASIIindex.py $
8# $Id: GSASIIindex.py 303 2011-06-20 19:37:07Z vondreele $
9########### SVN repository information ###################
10import math
11import wx
12import time
13import numpy as np
14import numpy.linalg as nl
15import GSASIIpath
16import pypowder as pyp              #assumes path has been amended to include correctr bin directory
17import GSASIIlattice as G2lat
18import scipy.optimize as so
19
20# trig functions in degrees
21sind = lambda x: math.sin(x*math.pi/180.)
22asind = lambda x: 180.*math.asin(x)/math.pi
23tand = lambda x: math.tan(x*math.pi/180.)
24atand = lambda x: 180.*math.atan(x)/math.pi
25atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
26cosd = lambda x: math.cos(x*math.pi/180.)
27acosd = lambda x: 180.*math.acos(x)/math.pi
28rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
29#numpy versions
30npsind = lambda x: np.sin(x*np.pi/180.)
31npasind = lambda x: 180.*np.arcsin(x)/math.pi
32npcosd = lambda x: np.cos(x*math.pi/180.)
33nptand = lambda x: np.tan(x*math.pi/180.)
34npatand = lambda x: 180.*np.arctan(x)/np.pi
35npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
36   
37def scaleAbyV(A,V):
38    v = G2lat.calc_V(A)
39    scale = math.exp(math.log(v/V)/3.)**2
40    for i in range(6):
41        A[i] *= scale
42   
43def ranaxis(dmin,dmax):
44    import random as rand
45    return rand.random()*(dmax-dmin)+dmin
46   
47def ran2axis(k,N):
48    import random as rand
49    T = 1.5+0.49*k/N
50#    B = 0.99-0.49*k/N
51#    B = 0.99-0.049*k/N
52    B = 0.99-0.149*k/N
53    R = (T-B)*rand.random()+B
54    return R
55   
56#def ranNaxis(k,N):
57#    import random as rand
58#    T = 1.0+1.0*k/N
59#    B = 1.0-1.0*k/N
60#    R = (T-B)*rand.random()+B
61#    return R
62   
63def ranAbyV(Bravais,dmin,dmax,V):
64    cell = [0,0,0,0,0,0]
65    bad = True
66    while bad:
67        bad = False
68        cell = rancell(Bravais,dmin,dmax)
69        G,g = G2lat.cell2Gmat(cell)
70        A = G2lat.Gmat2A(G)
71        if G2lat.calc_rVsq(A) < 1:
72            scaleAbyV(A,V)
73            cell = G2lat.A2cell(A)
74            for i in range(3):
75                bad |= cell[i] < dmin
76    return A
77   
78def ranAbyR(Bravais,A,k,N,ranFunc):
79    R = ranFunc(k,N)
80    if Bravais in [0,1,2]:          #cubic - not used
81        A[0] = A[1] = A[2] = A[0]*R
82        A[3] = A[4] = A[5] = 0.
83    elif Bravais in [3,4]:          #hexagonal/trigonal
84        A[0] = A[1] = A[3] = A[0]*R
85        A[2] *= R
86        A[4] = A[5] = 0.       
87    elif Bravais in [5,6]:          #tetragonal
88        A[0] = A[1] = A[0]*R
89        A[2] *= R
90        A[3] = A[4] = A[5] = 0.       
91    elif Bravais in [7,8,9,10]:     #orthorhombic
92        A[0] *= R
93        A[1] *= R
94        A[2] *= R
95        A[3] = A[4] = A[5] = 0.       
96    elif Bravais in [11,12]:        #monoclinic
97        A[0] *= R
98        A[1] *= R
99        A[2] *= R
100        A[4] *= R
101        A[3] = A[5] = 0.       
102    else:                           #triclinic
103        A[0] *= R
104        A[1] *= R
105        A[2] *= R
106        A[3] *= R
107        A[4] *= R
108        A[5] *= R
109    return A
110   
111def rancell(Bravais,dmin,dmax):
112    if Bravais in [0,1,2]:          #cubic
113        a = b = c = ranaxis(dmin,dmax)
114        alp = bet = gam = 90
115    elif Bravais in [3,4]:          #hexagonal/trigonal
116        a = b = ranaxis(dmin,dmax)
117        c = ranaxis(dmin,dmax)
118        alp = bet =  90
119        gam = 120
120    elif Bravais in [5,6]:          #tetragonal
121        a = b = ranaxis(dmin,dmax)
122        c = ranaxis(dmin,dmax)
123        alp = bet = gam = 90
124    elif Bravais in [7,8,9,10]:       #orthorhombic - F,I,C,P - a<b<c convention
125        abc = [ranaxis(dmin,dmax),ranaxis(dmin,dmax),ranaxis(dmin,dmax)]
126        abc.sort()
127        a = abc[0]
128        b = abc[1]
129        c = abc[2]
130        alp = bet = gam = 90
131    elif Bravais in [11,12]:        #monoclinic - C,P - a<c convention
132        ac = [ranaxis(dmin,dmax),ranaxis(dmin,dmax)]
133        ac.sort()
134        a = ac[0]
135        b = ranaxis(dmin,dmax)
136        c = ac[1]
137        alp = gam = 90
138        bet = ranaxis(90.,130.)
139    else:                           #triclinic - a<b<c convention
140        abc = [ranaxis(dmin,dmax),ranaxis(dmin,dmax),ranaxis(dmin,dmax)]
141        abc.sort()
142        a = abc[0]
143        b = abc[1]
144        c = abc[2]
145        r = 0.5*b/c
146        alp = ranaxis(acosd(r),acosd(-r))
147        r = 0.5*a/c
148        bet = ranaxis(acosd(r),acosd(-r))
149        r = 0.5*a/b
150        gam = ranaxis(acosd(r),acosd(-r)) 
151    return [a,b,c,alp,bet,gam]
152   
153def calc_M20(peaks,HKL):
154    diff = 0
155    X20 = 0
156    for Nobs20,peak in enumerate(peaks):
157        if peak[3]:
158            Qobs = 1.0/peak[7]**2
159            Qcalc = 1.0/peak[8]**2
160            diff += abs(Qobs-Qcalc)
161        elif peak[2]:
162            X20 += 1
163        if Nobs20 == 19: 
164            d20 = peak[7]
165            break
166    else:
167        d20 = peak[7]
168        Nobs20 = len(peaks)
169    for N20,hkl in enumerate(HKL):
170        if hkl[3] < d20:
171            break               
172    eta = diff/Nobs20
173    Q20 = 1.0/d20**2
174    if diff:
175        M20 = Q20/(2.0*diff)
176    else:
177        M20 = 0
178    M20 /= (1.+X20)
179    return M20,X20
180   
181def sortM20(cells):
182    #cells is M20,X20,Bravais,a,b,c,alp,bet,gam
183    #sort highest M20 1st
184    T = []
185    for i,M in enumerate(cells):
186        T.append((M[0],i))
187    D = dict(zip(T,cells))
188    T.sort()
189    T.reverse()
190    X = []
191    for key in T:
192        X.append(D[key])
193    return X
194               
195def IndexPeaks(peaks,HKL):
196    import bisect
197    N = len(HKL)
198    if N == 0: return False
199    hklds = list(np.array(HKL).T[3])+[1000.0,0.0,]
200    hklds.sort()                                        # ascending sort - upper bound at end
201    hklmax = [0,0,0]
202    for ipk,peak in enumerate(peaks):
203        if peak[2]:
204            i = bisect.bisect_right(hklds,peak[7])          # find peak position in hkl list
205            dm = peak[7]-hklds[i-1]                         # peak to neighbor hkls in list
206            dp = hklds[i]-peak[7]
207            pos = N-i                                       # reverse the order
208            if dp > dm: pos += 1                            # closer to upper than lower
209            hkl = HKL[pos]                                 # put in hkl
210            if hkl[4] >= 0:                                 # peak already assigned - test if this one better
211                opeak = peaks[hkl[4]]
212                dold = abs(opeak[7]-hkl[3])
213                dnew = min(dm,dp)
214                if dold > dnew:                             # new better - zero out old
215                    opeak[4:7] = [0,0,0]
216                    opeak[8] = 0.
217                else:                                       # old better - do nothing
218                    continue               
219            hkl[4] = ipk
220            peak[4:7] = hkl[:3]
221            peak[8] = hkl[3]                                # fill in d-calc
222    for peak in peaks:
223        peak[3] = False
224        if peak[2]:
225            if peak[8] > 0.:
226                for j in range(3):
227                    if abs(peak[j+4]) > hklmax[j]: hklmax[j] = abs(peak[j+4])
228                peak[3] = True
229    if hklmax[0]*hklmax[1]*hklmax[2] > 0:
230        return True
231    else:
232        return False
233
234def FitHKL(ibrav,peaks,A,Pwr):
235   
236    def Values2A(ibrav,values):
237        if ibrav in [0,1,2]:
238            return [values[0],values[0],values[0],0,0,0]
239        elif ibrav in [3,4]:
240            return [values[0],values[0],values[1],values[0],0,0]
241        elif ibrav in [5,6]:
242            return [values[0],values[0],values[1],0,0,0]
243        elif ibrav in [7,8,9,10]:
244            return [values[0],values[1],values[2],0,0,0]
245        elif ibrav in [11,12]:
246            return [values[0],values[1],values[2],0,values[3],0]
247        else:
248            return values
249           
250    def A2values(ibrav,A):
251        if ibrav in [0,1,2]:
252            return [A[0],]
253        elif ibrav in [3,4,5,6]:
254            return [A[0],A[2]]
255        elif ibrav in [7,8,9,10]:
256            return [A[0],A[1],A[2]]
257        elif ibrav in [11,12]:
258            return [A[0],A[1],A[2],A[4]]
259        else:
260            return A
261   
262    def errFit(values,ibrav,d,H,Pwr):
263        A = Values2A(ibrav,values)
264        Qo = 1./d**2
265        Qc = G2lat.calc_rDsq(H,A)
266        return (Qo-Qc)*d**Pwr
267   
268    Peaks = np.array(peaks).T
269    values = A2values(ibrav,A)   
270    result = so.leastsq(errFit,values,args=(ibrav,Peaks[7],Peaks[4:7],Pwr),full_output=True)
271    A = Values2A(ibrav,result[0])
272    return True,np.sum(errFit(result[0],ibrav,Peaks[7],Peaks[4:7],Pwr)**2),A
273           
274def FitHKLZ(ibrav,peaks,Z,A):
275    return A,Z
276   
277def rotOrthoA(A):
278    return [A[1],A[2],A[0],0,0,0]
279   
280def swapMonoA(A):
281    return [A[2],A[1],A[0],0,A[4],0]
282   
283def oddPeak(indx,peaks):
284    noOdd = True
285    for peak in peaks:
286        H = peak[4:7]
287        if H[indx] % 2:
288            noOdd = False
289    return noOdd
290   
291def halfCell(ibrav,A,peaks):
292    if ibrav in [0,1,2]:
293        if oddPeak(0,peaks):
294            A[0] *= 2
295            A[1] = A[2] = A[0]
296    elif ibrav in [3,4,5,6]:
297        if oddPeak(0,peaks):
298            A[0] *= 2
299            A[1] = A[0]
300        if oddPeak(2,peaks):
301            A[2] *=2
302    else:
303        if oddPeak(0,peaks):
304            A[0] *=2
305        if oddPeak(1,peaks):
306            A[1] *=2
307        if oddPeak(2,peaks):
308            A[2] *=2
309    return A
310   
311def getDmin(peaks):
312    return peaks[-1][7]
313   
314def getDmax(peaks):
315    return peaks[0][7]
316   
317def refinePeaks(peaks,ibrav,A):
318    dmin = getDmin(peaks)
319    smin = 1.0e10
320    pwr = 3
321    maxTries = 10
322    if ibrav == 13:
323        pwr = 4
324        maxTries = 10
325    OK = False
326    tries = 0
327    HKL = G2lat.GenHBravais(dmin,ibrav,A)
328    while IndexPeaks(peaks,HKL):
329        Pwr = pwr - (tries % 2)
330        HKL = []
331        tries += 1
332        osmin = smin
333        oldA = A
334        OK,smin,A = FitHKL(ibrav,peaks,A,Pwr)
335        if min(A[:3]) <= 0:
336            A = oldA
337            OK = False
338            break
339        if OK:
340            HKL = G2lat.GenHBravais(dmin,ibrav,A)
341        if len(HKL) == 0: break                         #absurd cell obtained!
342        rat = (osmin-smin)/smin
343        if abs(rat) < 1.0e-5 or not OK: break
344        if tries > maxTries: break
345    if OK:
346        OK,smin,A = FitHKL(ibrav,peaks,A,2)
347        Peaks = np.array(peaks).T
348        H = Peaks[4:7]
349        Peaks[8] = 1./np.sqrt(G2lat.calc_rDsq(H,A))
350        peaks = Peaks.T
351       
352    M20,X20 = calc_M20(peaks,HKL)
353    return len(HKL),M20,X20,A
354       
355def findBestCell(dlg,ncMax,A,Ntries,ibrav,peaks,V1):
356# dlg & ncMax are used for wx progress bar
357# A != 0 find the best A near input A,
358# A = 0 for random cell, volume normalized to V1;
359# returns number of generated hkls, M20, X20 & A for best found
360    mHKL = [3,3,3, 5,5, 5,5, 7,7,7,7, 9,9, 10]
361    dmin = getDmin(peaks)-0.05
362    amin = 2.5
363    amax = 5.*getDmax(peaks)
364    Asave = []
365    GoOn = True
366    if A:
367        HKL = G2lat.GenHBravais(dmin,ibrav,A[:])
368        if len(HKL) > mHKL[ibrav]:
369            IndexPeaks(peaks,HKL)
370            Asave.append([calc_M20(peaks,HKL),A[:]])
371    tries = 0
372    while tries < Ntries:
373        if A:
374            Anew = ranAbyR(ibrav,A[:],tries+1,Ntries,ran2axis)
375            if ibrav in [11,12,13]:
376                Anew = ranAbyR(ibrav,A[:],tries/10+1,Ntries,ran2axis)
377        else:
378            Anew = ranAbyV(ibrav,amin,amax,V1)
379        HKL = G2lat.GenHBravais(dmin,ibrav,Anew)
380       
381        if IndexPeaks(peaks,HKL) and len(HKL) > mHKL[ibrav]:
382            Lhkl,M20,X20,Anew = refinePeaks(peaks,ibrav,Anew)
383            Asave.append([calc_M20(peaks,HKL),Anew[:]])
384            if ibrav == 9:                          #C-centered orthorhombic
385                for i in range(2):
386                    Anew = rotOrthoA(Anew[:])
387                    Lhkl,M20,X20,Anew = refinePeaks(peaks,ibrav,Anew)
388                    HKL = G2lat.GenHBravais(dmin,ibrav,Anew)
389                    IndexPeaks(peaks,HKL)
390                    Asave.append([calc_M20(peaks,HKL),Anew[:]])
391            elif ibrav == 11:                      #C-centered monoclinic
392                Anew = swapMonoA(Anew[:])
393                Lhkl,M20,X20,Anew = refinePeaks(peaks,ibrav,Anew)
394                HKL = G2lat.GenHBravais(dmin,ibrav,Anew)
395                IndexPeaks(peaks,HKL)
396                Asave.append([calc_M20(peaks,HKL),Anew[:]])
397        else:
398            break
399        Nc = len(HKL)
400        if Nc >= ncMax:
401            GoOn = False
402        elif dlg:
403            GoOn = dlg.Update(Nc)[0]
404            if not GoOn:
405                break
406        tries += 1   
407    X = sortM20(Asave)
408    if X:
409        Lhkl,M20,X20,A = refinePeaks(peaks,ibrav,X[0][1])
410        return GoOn,Lhkl,M20,X20,X[0][1]
411    else:
412        return GoOn,0,0,0,Anew
413       
414def monoCellReduce(ibrav,A):
415    a,b,c,alp,bet,gam = G2lat.A2cell(A)
416    G,g = G2lat.A2Gmat(A)
417    if ibrav in [11]:
418        u = [0,0,-1]
419        v = [1,0,2]
420        anew = math.sqrt(np.dot(np.dot(v,g),v))
421        if anew < a:
422            cang = np.dot(np.dot(u,g),v)/(anew*c)
423            beta = acosd(-abs(cang))
424            A = G2lat.cell2A([anew,b,c,90,beta,90])
425    else:
426        u = [-1,0,0]
427        v = [1,0,1]
428        cnew = math.sqrt(np.dot(np.dot(v,g),v))
429        if cnew < c:
430            cang = np.dot(np.dot(u,g),v)/(a*cnew)
431            beta = acosd(-abs(cang))
432            A = G2lat.cell2A([a,b,cnew,90,beta,90])
433    return A
434
435def DoIndexPeaks(peaks,inst,controls,bravais):
436   
437    delt = 0.005                                     #lowest d-spacing cushion - can be fixed?
438    amin = 2.5
439    amax = 5.0*getDmax(peaks)
440    dmin = getDmin(peaks)-delt
441    bravaisNames = ['Cubic-F','Cubic-I','Cubic-P','Trigonal-R','Trigonal/Hexagonal-P',
442        'Tetragonal-I','Tetragonal-P','Orthorhombic-F','Orthorhombic-I','Orthorhombic-C',
443        'Orthorhombic-P','Monoclinic-C','Monoclinic-P','Triclinic']
444    tries = ['1st','2nd','3rd','4th','5th','6th','7th','8th','9th','10th']
445    N1s = [1,1,1,   5,5,  5,5, 50,50,50,50,  50,50, 200]
446    N2s = [1,1,1,   2,2,  2,2,     2,2,2,2,   2,2,   4]
447    Nm  = [1,1,1,   1,1,  1,1,     1,1,1,1,   2,2,   4]
448    Nobs = len(peaks)
449    wave = inst[1]
450    if len(inst) > 10:
451        zero = inst[3]
452    else:
453        zero = inst[2]
454    print "%s %8.5f %6.3f" % ('wavelength, zero =',wave,zero)
455    print "%s %8.3f %8.3f" % ('lattice parameter range = ',amin,amax)
456    ifzero,maxzero,ncno = controls[:3]
457    ncMax = Nobs*ncno
458    print "%s %d %s %d %s %d" % ('change zero =',ifzero,'Nc/No max =',ncno,' Max Nc =',ncno*Nobs)
459    cells = []
460    for ibrav in range(14):
461        begin = time.time()
462        if bravais[ibrav]:
463            print 'cell search for ',bravaisNames[ibrav]
464            print '      M20  X20  Nc       a          b          c        alpha       beta      gamma     volume      V-test'
465            V1 = controls[3]
466            bestM20 = 0
467            topM20 = 0
468            cycle = 0
469            while cycle < 5:
470                dlg = wx.ProgressDialog("Generated reflections",tries[cycle]+" cell search for "+bravaisNames[ibrav],ncMax, 
471                    style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_REMAINING_TIME|wx.PD_CAN_ABORT)
472                screenSize = wx.ClientDisplayRect()
473                Size = dlg.GetSize()
474                dlg.SetPosition(wx.Point(screenSize[2]-Size[0]-305,screenSize[1]+5))
475                try:
476                    GoOn = True
477                    while GoOn:                                                 #Loop over increment of volume
478                        N2 = 0
479                        while N2 < N2s[ibrav]:                                  #Table 2 step (iii)               
480                            if ibrav > 2:
481                                if not N2:
482                                    A = []
483                                    GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,A,Nm[ibrav]*N1s[ibrav],ibrav,peaks,V1)
484                                if A:
485                                    GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,A[:],N1s[ibrav],ibrav,peaks,0)
486                            else:
487                                GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,0,Nm[ibrav]*N1s[ibrav],ibrav,peaks,V1)
488                            if Nc >= ncMax:
489                                GoOn = False
490                                break
491                            elif 3*Nc < Nobs:
492                                N2 = 10
493                                break
494                            else:
495                                if not GoOn:
496                                    break
497                                if M20 > 1.0:
498                                    bestM20 = max(bestM20,M20)
499                                    A = halfCell(ibrav,A[:],peaks)
500                                    if ibrav in [12]:
501                                        A = monoCellReduce(ibrav,A[:])
502                                    HKL = G2lat.GenHBravais(dmin,ibrav,A)
503                                    IndexPeaks(peaks,HKL)
504                                    a,b,c,alp,bet,gam = G2lat.A2cell(A)
505                                    V = G2lat.calc_V(A)
506                                    print "%10.3f %3d %3d %10.5f %10.5f %10.5f %10.3f %10.3f %10.3f %10.2f %10.2f" % (M20,X20,Nc,a,b,c,alp,bet,gam,V,V1)
507                                    if M20 >= 2.0:
508                                        cells.append([M20,X20,ibrav,a,b,c,alp,bet,gam,V,False])
509                            if not GoOn:
510                                break
511                            N2 += 1
512                        if ibrav < 11:
513                            V1 *= 1.1
514                        elif ibrav in range(11,14):
515                            V1 *= 1.05
516                        if not GoOn:
517                            if bestM20 > topM20:
518                                topM20 = bestM20
519                                if cells:
520                                    V1 = cells[0][9]
521                                else:
522                                    V1 = 25
523                                ncMax += Nobs
524                                cycle += 1
525                                print 'Restart search, new Max Nc = ',ncMax
526                            else:
527                                cycle = 10
528                finally:
529                    dlg.Destroy()
530            print '%s%s%s%s'%('finished cell search for ',bravaisNames[ibrav], \
531                ', elapsed time = ',G2lat.sec2HMS(time.time()-begin))
532           
533    if cells:
534        cells = sortM20(cells)
535        cells[0][-1] = True
536        return True,dmin,cells
537    else:
538        return False,0,0
539       
Note: See TracBrowser for help on using the repository browser.