source: trunk/GSASIIindex.py @ 102

Last change on this file since 102 was 78, checked in by vondreel, 12 years ago

split IndexPeaks?, etc. from GSASIIcomp.py & put in new GSASIIindex.py

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