source: trunk/GSASIIstruct.py @ 342

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

major modifications to get 1st "working" version of Refine
Add GSASIImapvars.py
fix cell refinement in indexing window
single cycle for peak refinement
  • Property svn:keywords set to Date Author Revision URL Id
File size: 54.9 KB
Line 
1#GSASIIstructure - structure computation routines
2########### SVN repository information ###################
3# $Date: 2011-08-05 19:35:43 +0000 (Fri, 05 Aug 2011) $
4# $Author: vondreele $
5# $Revision: 342 $
6# $URL: trunk/GSASIIstruct.py $
7# $Id: GSASIIstruct.py 342 2011-08-05 19:35:43Z vondreele $
8########### SVN repository information ###################
9import sys
10import os
11import os.path as ospath
12import numpy as np
13import cPickle
14import time
15import math
16import wx
17import GSASIIpath
18import GSASIIElem as G2el
19import GSASIIlattice as G2lat
20import GSASIIspc as G2spc
21import GSASIIpwd as G2pwd
22import GSASIImapvars as G2mv
23import scipy.optimize as so
24
25sind = lambda x: np.sin(x*np.pi/180.)
26cosd = lambda x: np.cos(x*np.pi/180.)
27tand = lambda x: np.tan(x*np.pi/180.)
28asind = lambda x: 180.*np.arcsin(x)/np.pi
29
30
31def ShowBanner():
32    print 80*'*'
33    print '   General Structure Analysis System-II Crystal Structure Refinement'
34    print '     by Robert B. Von Dreele, Argonne National Laboratory(C), 2010'
35    print ' This product includes software developed by the UChicago Argonne, LLC,' 
36    print '            as Operator of Argonne National Laboratory.'
37    print 80*'*','\n'
38
39def GetControls(GPXfile):
40    ''' Returns dictionary of control items found in GSASII gpx file
41    input:
42        GPXfile = .gpx full file name
43    return:
44        Controls = dictionary of control items
45    '''
46    file = open(GPXfile,'rb')
47    while True:
48        try:
49            data = cPickle.load(file)
50        except EOFError:
51            break
52        datum = data[0]
53        if datum[0] == 'Controls':
54            Controls = datum[1]
55    file.close()
56    return Controls
57   
58def ShowControls(Controls):
59    print ' Controls:'
60    if Controls['bandWidth']:
61        print ' Band approximated least squares matrix refinement, width: ',Controls['bandWidth']
62    else:
63        print ' Full matrix least squares refinement'
64    print ' Maximum number of refinement cycles: ',Controls['Ncycles']
65    print ' Minimum sum(shift/esd)^2 for convergence: ','%.2f'%(Controls['minSumShftESD'])
66    print ' The Marquardt damping factor: ','%.2f'%(Controls['Marquardt'])
67    print ' Maximum allowed atom shift: ','%.2f'%(Controls['maxShift']),'A'
68    if Controls['restraintWeight']:
69        print ' The weights of the restraint histograms will be modified to normalize their contribution','\n'
70    else:
71        print ' The restraint histogram weights are fixed','\n'
72   
73def GetPhaseNames(GPXfile):
74    ''' Returns a list of phase names found under 'Phases' in GSASII gpx file
75    input:
76        GPXfile = gpx full file name
77    return:
78        PhaseNames = list of phase names
79    '''
80    file = open(GPXfile,'rb')
81    PhaseNames = []
82    while True:
83        try:
84            data = cPickle.load(file)
85        except EOFError:
86            break
87        datum = data[0]
88        if 'Phases' == datum[0]:
89            for datus in data[1:]:
90                PhaseNames.append(datus[0])
91    file.close()
92    return PhaseNames
93
94def GetAllPhaseData(GPXfile,PhaseName):
95    ''' Returns the entire dictionary for PhaseName from GSASII gpx file
96    input:
97        GPXfile = gpx full file name
98        PhaseName = phase name
99    return:
100        phase dictionary
101    '''       
102    file = open(GPXfile,'rb')
103    General = {}
104    Atoms = []
105    while True:
106        try:
107            data = cPickle.load(file)
108        except EOFError:
109            break
110        datum = data[0]
111        if 'Phases' == datum[0]:
112            for datus in data[1:]:
113                if datus[0] == PhaseName:
114                    break
115    file.close()
116    return datus[1]
117   
118def GetHistogramNames(GPXfile):
119    ''' Returns a list of histogram names found in GSASII gpx file
120    input:
121        GPXfile = .gpx full file name
122    return:
123        HistogramNames = list of histogram names (types = PWDR & HKLF)
124    '''
125    file = open(GPXfile,'rb')
126    HistogramNames = []
127    while True:
128        try:
129            data = cPickle.load(file)
130        except EOFError:
131            break
132        datum = data[0]
133        if datum[0][:4] in ['PWDR','HKLF']:
134            HistogramNames.append(datum[0])
135    file.close()
136    return HistogramNames
137   
138def GetUsedHistogramsAndPhases(GPXfile):
139    ''' Returns all histograms that are found in any phase
140    and any phase that uses a histogram
141    input:
142        GPXfile = .gpx full file name
143    return:
144        Histograms = dictionary of histograms as {name:data,...}
145        Phases = dictionary of phases that use histograms
146    '''
147    phaseNames = GetPhaseNames(GPXfile)
148    phaseData = {}
149    for name in phaseNames: 
150        phaseData[name] =  GetAllPhaseData(GPXfile,name)
151    Histograms = {}
152    Phases = {}
153    pId = 0
154    hId = 0
155    for phase in phaseData:
156        Phase = phaseData[phase]
157        if Phase['Histograms']:
158            if phase not in Phases:
159                Phase['pId'] = pId
160                pId += 1
161                Phases[phase] = Phase
162            for hist in Phase['Histograms']:
163                if hist not in Histograms:
164                    if 'PWDR' in hist[:4]: 
165                        Histograms[hist] = GetPWDRdata(GPXfile,hist)
166                    elif 'HKLF' in hist[:4]:
167                        Histograms[hist] = GetHKLFdata(GPXfile,hist)
168                    #future restraint, etc. histograms here           
169                    Histograms[hist]['hId'] = hId
170                    hId += 1
171    return Histograms,Phases
172   
173def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases):
174    ''' Updates gpxfile from all histograms that are found in any phase
175    and any phase that used a histogram
176    input:
177        GPXfile = .gpx full file name
178        Histograms = dictionary of histograms as {name:data,...}
179        Phases = dictionary of phases that use histograms
180    '''
181                       
182    def GPXBackup(GPXfile):
183        import distutils.file_util as dfu
184        GPXpath,GPXname = ospath.split(GPXfile)
185        Name = ospath.splitext(GPXname)[0]
186        files = os.listdir(GPXpath)
187        last = 0
188        for name in files:
189            name = name.split('.')
190            if len(name) == 3 and name[0] == Name and 'bak' in name[1]:
191                last = max(last,int(name[1].strip('bak'))+1)
192        GPXback = ospath.join(GPXpath,ospath.splitext(GPXname)[0]+'.bak'+str(last)+'.gpx')
193        dfu.copy_file(GPXfile,GPXback)
194        return GPXback
195       
196    GPXback = GPXBackup(GPXfile)
197    print '\n',130*'-'
198    print 'Read from file:',GPXback
199    print 'Save to file  :',GPXfile
200    infile = open(GPXback,'rb')
201    outfile = open(GPXfile,'wb')
202    while True:
203        try:
204            data = cPickle.load(infile)
205        except EOFError:
206            break
207        datum = data[0]
208        print 'read: ',datum[0]
209        if datum[0] == 'Phases':
210            for iphase in range(len(data)):
211                if data[iphase][0] in Phases:
212                    phaseName = data[iphase][0]
213                    data[iphase][1] = Phases[phaseName]
214        try:
215            histogram = Histograms[datum[0]]
216            print 'found ',datum[0]
217            data[0][1][1] = histogram['Data']
218            for datus in data[1:]:
219                print '    read: ',datus[0]
220                if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
221                    datus[1] = histogram[datus[0]]
222        except KeyError:
223            pass
224                               
225        cPickle.dump(data,outfile,1)
226    infile.close()
227    outfile.close()
228    print 'refinement save successful'
229                   
230def GetPWDRdata(GPXfile,PWDRname):
231    ''' Returns powder data from GSASII gpx file
232    input:
233        GPXfile = .gpx full file name
234        PWDRname = powder histogram name as obtained from GetHistogramNames
235    return:
236        PWDRdata = powder data dictionary with:
237            Data - powder data arrays, Limits, Instrument Parameters, Sample Parameters
238       
239    '''
240    file = open(GPXfile,'rb')
241    PWDRdata = {}
242    while True:
243        try:
244            data = cPickle.load(file)
245        except EOFError:
246            break
247        datum = data[0]
248        if datum[0] == PWDRname:
249            PWDRdata['Data'] = datum[1][1]          #powder data arrays
250            PWDRdata[data[2][0]] = data[2][1]       #Limits
251            PWDRdata[data[3][0]] = data[3][1]       #Background
252            PWDRdata[data[4][0]] = data[4][1]       #Instrument parameters
253            PWDRdata[data[5][0]] = data[5][1]       #Sample parameters
254            try:
255                PWDRdata[data[9][0]] = data[9][1]       #Reflection lists might be missing
256            except IndexError:
257                PWDRdata['Reflection lists'] = {}
258    file.close()
259    return PWDRdata
260   
261def GetHKLFdata(GPXfile,HKLFname):
262    ''' Returns single crystal data from GSASII gpx file
263    input:
264        GPXfile = .gpx full file name
265        HKLFname = single crystal histogram name as obtained from GetHistogramNames
266    return:
267        HKLFdata = single crystal data list of reflections: for each reflection:
268            HKLF = [np.array([h,k,l]),FoSq,sigFoSq,FcSq,Fcp,Fcpp,phase]
269    '''
270    file = open(GPXfile,'rb')
271    HKLFdata = []
272    while True:
273        try:
274            data = cPickle.load(file)
275        except EOFError:
276            break
277        datum = data[0]
278        if datum[0] == HKLFname:
279            HKLFdata = datum[1:][0]
280    file.close()
281    return HKLFdata
282   
283def GetFFtable(General):
284    ''' returns a dictionary of form factor data for atom types found in General
285    input:
286        General = dictionary of phase info.; includes AtomTypes
287    return:
288        FFtable = dictionary of form factor data; key is atom type
289    '''
290    atomTypes = General['AtomTypes']
291    FFtable = {}
292    for El in atomTypes:
293        FFs = G2el.GetFormFactorCoeff(El.split('+')[0].split('-')[0])
294        for item in FFs:
295            if item['Symbol'] == El:
296                FFtable[El] = item
297    return FFtable
298   
299def GetPawleyConstr(SGLaue,PawleyRef,pawleyVary):
300    if SGLaue in ['-1','2/m','mmm']:
301        return                      #no Pawley symmetry required constraints
302    for varyI in pawleyVary:
303        refI = int(varyI.split(':')[-1])
304        ih,ik,il = PawleyRef[refI][:3]
305        for varyJ in pawleyVary[:refI]:
306            refJ = int(varyJ.split(':')[-1])
307            jh,jk,jl = PawleyRef[refJ][:3]
308            if SGLaue in ['4/m','4/mmm']:
309                isum = ih**2+ik**2
310                jsum = jh**2+jk**2
311                if abs(il) == abs(jl) and isum == jsum:
312                    G2mv.StoreEquivalence(varyJ,(varyI,))
313            elif SGLaue in ['3R','3mR']:
314                isum = ih**2+ik**2+il**2
315                jsum = jh**2+jk**2*jl**2
316                isum2 = ih*ik+ih*il+ik*il
317                jsum2 = jh*jk+jh*jl+jk*jl
318                if isum == jsum and isum2 == jsum2:
319                    G2mv.StoreEquivalence(varyJ,(varyI,))
320            elif SGLaue in ['3','3m1','31m','6/m','6/mmm']:
321                isum = ih**2+ik**2+ih*ik
322                jsum = jh**2+jk**2+jh*jk
323                if abs(il) == abs(jl) and isum == jsum:
324                    G2mv.StoreEquivalence(varyJ,(varyI,))
325            elif SGLaue in ['m3','m3m']:
326                isum = ih**2+ik**2+il**2
327                jsum = jh**2+jk**2+jl**2
328                if isum == jsum:
329                    G2mv.StoreEquivalence(varyJ,(varyI,))
330       
331def GetPhaseData(PhaseData):
332   
333    def cellVary(pfx,SGData): 
334        if SGData['SGLaue'] in ['-1',]:
335            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
336        elif SGData['SGLaue'] in ['2/m',]:
337            if SGData['SGUniq'] == 'a':
338                return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3']
339            elif SGData['SGUniq'] == 'b':
340                return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A4']
341            else:
342                return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A5']
343        elif SGData['SGLaue'] in ['mmm',]:
344            return [pfx+'A0',pfx+'A1',pfx+'A2']
345        elif SGData['SGLaue'] in ['4/m','4/mmm']:
346            return [pfx+'A0',pfx+'A2']
347        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
348            return [pfx+'A0',pfx+'A2']
349        elif SGData['SGLaue'] in ['3R', '3mR']:
350            return [pfx+'A0',pfx+'A3']                       
351        elif SGData['SGLaue'] in ['m3m','m3']:
352            return [pfx+'A0']
353       
354       
355    print ' Phases:'
356    phaseVary = []
357    phaseDict = {}
358    phaseConstr = {}
359    pawleyLookup = {}
360    for name in PhaseData:
361        General = PhaseData[name]['General']
362        pId = PhaseData[name]['pId']
363        pfx = str(pId)+'::'
364        Atoms = PhaseData[name]['Atoms']
365        try:
366            PawleyRef = PhaseData[name]['Pawley ref']
367        except KeyError:
368            PawleyRef = []
369        print '\n Phase name: ',General['Name']
370        SGData = General['SGData']
371        SGtext = G2spc.SGPrint(SGData)
372        for line in SGtext: print line
373        cell = General['Cell']
374        A = G2lat.cell2A(cell[1:7])
375        phaseDict.update({pfx+'A0':A[0],pfx+'A1':A[1],pfx+'A2':A[2],pfx+'A3':A[3],pfx+'A4':A[4],pfx+'A5':A[5]})
376        if cell[0]:
377            phaseVary = cellVary(pfx,SGData)
378        print '\n Unit cell: a =','%.5f'%(cell[1]),' b =','%.5f'%(cell[2]),' c =','%.5f'%(cell[3]), \
379            ' alpha =','%.3f'%(cell[4]),' beta =','%.3f'%(cell[5]),' gamma =', \
380            '%.3f'%(cell[6]),' volume =','%.3f'%(cell[7]),' Refine?',cell[0]
381        if Atoms:
382            print '\n Atoms:'
383            line = '   name    type  refine?   x         y         z    '+ \
384                '  frac site sym  mult I/A   Uiso     U11     U22     U33     U12     U13     U23'
385            if General['Type'] == 'magnetic':
386                line += '   Mx     My     Mz'
387            elif General['Type'] == 'macromolecular':
388                line = ' res no  residue  chain '+line
389            print line
390            if General['Type'] == 'nuclear':
391                print 135*'-'
392                for at in Atoms:
393                    line = '%7s'%(at[0])+'%7s'%(at[1])+'%7s'%(at[2])+'%10.5f'%(at[3])+'%10.5f'%(at[4])+ \
394                        '%10.5f'%(at[5])+'%8.3f'%(at[6])+'%7s'%(at[7])+'%5d'%(at[8])+'%5s'%(at[9])
395                    if at[9] == 'I':
396                        line += '%8.4f'%(at[10])+48*' '
397                    else:
398                        line += 8*' '
399                        for i in range(6):
400                            line += '%8.4f'%(at[11+i])
401                    print line
402                    if 'X' in at[2]:
403                        xId,xCoef = G2spc.GetCSxinel(at[7])
404                    if 'U' in at[2]:
405                        uId,uCoef = G2spc.GetCSuinel(at[7])[:2]
406#        elif General['Type'] == 'magnetic':
407#        elif General['Type'] == 'macromolecular':
408#       PWDR: harmonics texture, unit cell, atom parameters
409        elif PawleyRef:
410            pawleyVary = []
411            for i,refl in enumerate(PawleyRef):
412                phaseDict[pfx+'PWLref:'+str(i)] = refl[6]
413                pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
414                if refl[5]:
415                    pawleyVary.append(pfx+'PWLref:'+str(i))
416            GetPawleyConstr(SGData['SGLaue'],PawleyRef,pawleyVary)      #does G2mv.StoreEquivalence
417            phaseVary += pawleyVary   
418               
419    return phaseVary,phaseDict,pawleyLookup
420   
421def SetPhaseData(parmDict,sigDict,Phases):
422   
423    def cellFill(pfx,SGData,parmDict,sigDict): 
424        if SGData['SGLaue'] in ['-1',]:
425            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
426                parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
427            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
428                sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
429        elif SGData['SGLaue'] in ['2/m',]:
430            if SGData['SGUniq'] == 'a':
431                A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
432                    parmDict[pfx+'A3'],0,0]
433                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
434                    sigDict[pfx+'A3'],0,0]
435            elif SGData['SGUniq'] == 'b':
436                A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
437                    0,parmDict[pfx+'A4'],0]
438                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
439                    0,sigDict[pfx+'A4'],0]
440            else:
441                A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
442                    0,0,parmDict[pfx+'A5']]
443                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
444                    0,0,sigDict[pfx+'A5']]
445        elif SGData['SGLaue'] in ['mmm',]:
446            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
447            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
448        elif SGData['SGLaue'] in ['4/m','4/mmm']:
449            A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
450            sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
451        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
452            A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
453                parmDict[pfx+'A0'],0,0]
454            sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
455        elif SGData['SGLaue'] in ['3R', '3mR']:
456            A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
457                parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
458            sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
459        elif SGData['SGLaue'] in ['m3m','m3']:
460            A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
461            sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
462        return A,sigA
463           
464    print '\n Phases:'
465    for phase in Phases:
466        print ' Result for phase: ',phase
467        Phase = Phases[phase]
468        General = Phase['General']
469        SGData = General['SGData']
470        cell = General['Cell']
471        pId = Phase['pId']
472        pfx = str(pId)+'::'
473        if cell[0]:
474            A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
475            print ' Reciprocal metric tensor: '
476            ptfmt = "%15.9f"
477            names = ['A11','A22','A33','A12','A13','A23']
478            namstr = '  names :'
479            ptstr =  '  values:'
480            sigstr = '  esds  :'
481            for name,a,siga in zip(names,A,sigA):
482                namstr += '%15s'%(name)
483                ptstr += ptfmt%(a)
484                if siga:
485                    sigstr += ptfmt%(siga)
486                else:
487                    sigstr += 15*' '
488            print namstr
489            print ptstr
490            print sigstr
491            cell[1:7] = G2lat.A2cell(A)
492            cell[7] = G2lat.calc_V(A)
493            print ' New unit cell: a =%.5f'%(cell[1]),' b =%.5f'%(cell[2]), \
494                ' c =%.5f'%(cell[3]),' alpha =%.3f'%(cell[4]),' beta =%.3f'%(cell[5]), \
495                ' gamma =%.3f'%(cell[6]),' volume =%.3f'%(cell[7])
496        if 'Pawley' in Phase['General']['Type']:
497            pawleyRef = Phase['Pawley ref']
498            for i,refl in enumerate(pawleyRef):
499                key = pfx+'PWLref:'+str(i)
500                refl[6] = parmDict[key]
501                if key in sigDict:
502                    refl[7] = sigDict[key]
503                else:
504                    refl[7] = 0
505
506def GetHistogramPhaseData(Phases,Histograms):
507   
508    def PrintSize(hapData):
509        line = '\n Size model    : '+hapData[0]
510        if hapData[0] in ['isotropic','uniaxial']:
511            line += ' equatorial:'+'%12.2f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
512            if hapData[0] == 'uniaxial':
513                line += ' axial:'+'%12.2f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
514            print line
515        else:
516            print line
517            Snames = ['S11','S22','S33','S12','S13','S23']
518            ptlbls = ' names :'
519            ptstr =  ' values:'
520            varstr = ' refine:'
521            for i,name in enumerate(Snames):
522                ptlbls += '%12s' % (name)
523                ptstr += '%12.6f' % (hapData[4][i])
524                varstr += '%12s' % (str(hapData[5][i]))
525            print ptlbls
526            print ptstr
527            print varstr
528       
529    def PrintMuStrain(hapData,SGData):
530        line = '\n Mustrain model: '+hapData[0]
531        if hapData[0] in ['isotropic','uniaxial']:
532            line += ' equatorial:'+'%12.4f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
533            if hapData[0] == 'uniaxial':
534                line += ' axial:'+'%12.4f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
535            print line
536        else:
537            print line
538            Snames = G2spc.MustrainNames(SGData)
539            ptlbls = ' names :'
540            ptstr =  ' values:'
541            varstr = ' refine:'
542            for i,name in enumerate(Snames):
543                ptlbls += '%12s' % (name)
544                ptstr += '%12.6f' % (hapData[4][i])
545                varstr += '%12s' % (str(hapData[5][i]))
546            print ptlbls
547            print ptstr
548            print varstr
549       
550   
551    hapDict = {}
552    hapVary = []
553    controlDict = {}
554    poType = {}
555    poAxes = {}
556    spAxes = {}
557    spType = {}
558   
559    for phase in Phases:
560        HistoPhase = Phases[phase]['Histograms']
561        SGData = Phases[phase]['General']['SGData']
562        cell = Phases[phase]['General']['Cell'][1:7]
563        A = G2lat.cell2A(cell)
564        pId = Phases[phase]['pId']
565        for histogram in HistoPhase:
566            print '\n Phase: ',phase,' in histogram: ',histogram
567            hapData = HistoPhase[histogram]
568            Histogram = Histograms[histogram]
569            hId = Histogram['hId']
570            limits = Histogram['Limits'][1]
571            inst = Histogram['Instrument Parameters']
572            inst = dict(zip(inst[3],inst[1]))
573            Zero = inst['Zero']
574            if 'C' in inst['Type']:
575                try:
576                    wave = inst['Lam']
577                except KeyError:
578                    wave = inst['Lam1']
579                dmin = wave/(2.0*sind(limits[1]/2.0))
580            pfx = str(pId)+':'+str(hId)+':'
581            for item in ['Scale','Extinction']:
582                hapDict[pfx+item] = hapData[item][0]
583                if hapData[item][1]:
584                    hapVary.append(pfx+item)
585            controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
586            if hapData['Pref.Ori.'][0] == 'MD':
587                hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
588                controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
589                if hapData['Pref.Ori.'][2]:
590                    hapVary.append(pfx+'MD')
591            else:                           #'SH' spherical harmonics
592                for item in hapData['Pref.Ori.'][5]:
593                    hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
594                    if hapData['Pref.Ori.'][2]:
595                        hapVary.append(pfx+item)                   
596            for item in ['Mustrain','Size']:
597                controlDict[pfx+item+'Type'] = hapData[item][0]
598                if hapData[item][0] in ['isotropic','uniaxial']:
599                    hapDict[pfx+item+':0'] = hapData[item][1][0]
600                    if hapData[item][2][0]:
601                        hapVary.append(pfx+item+':0')
602                    if hapData[item][0] == 'uniaxial':
603                        controlDict[pfx+item+'Axis'] = hapData[item][3]
604                        hapDict[pfx+item+':1'] = hapData[item][1][1]
605                        if hapData[item][2][1]:
606                            hapVary.append(pfx+item+':1')
607                else:       #generalized for mustrain or ellipsoidal for size
608                    if item == 'Mustrain':
609                        names = G2spc.MustrainNames(SGData)
610                        pwrs = []
611                        for name in names:
612                            h,k,l = name[1:]
613                            pwrs.append([int(h),int(k),int(l)])
614                        controlDict[pfx+'MuPwrs'] = pwrs
615                    for i in range(len(hapData[item][4])):
616                        sfx = ':'+str(i)
617                        hapDict[pfx+item+sfx] = hapData[item][4][i]
618                        if hapData[item][5][i]:
619                            hapVary.append(pfx+item+sfx)
620                           
621            print '\n Phase fraction  : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
622            print ' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
623            if hapData['Pref.Ori.'][0] == 'MD':
624                Ax = hapData['Pref.Ori.'][3]
625                print ' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
626                    ' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])               
627            PrintSize(hapData['Size'])
628            PrintMuStrain(hapData['Mustrain'],SGData)
629            HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
630            refList = []
631            for h,k,l,d in HKLd:
632                ext,mul,Uniq = G2spc.GenHKLf([h,k,l],SGData)
633                if ext:
634                    continue
635                if 'C' in inst['Type']:
636                    pos = 2.0*asind(wave/(2.0*d))
637                    if limits[0] < pos < limits[1]:
638                        refList.append([h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,Uniq])
639                else:
640                    raise ValueError 
641            Histogram['Reflection Lists'][phase] = refList
642    return hapVary,hapDict,controlDict
643   
644def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms):
645   
646    def PrintSizeAndSig(hapData,sizeSig):
647        line = '\n Size model: '+hapData[0]
648        if hapData[0] in ['isotropic','uniaxial']:
649            line += ' equatorial:%12.2f'%(hapData[1][0])
650            if sizeSig[0][0]:
651                line += ', sig: %8.2f'%(sizeSig[0][0])
652            if hapData[0] == 'uniaxial':
653                line += ' axial:%12.2f'%(hapData[1][1])
654                if sizeSig[0][1]:
655                    line += ', sig: %8.2f'%(sizeSig[0][1])
656            print line
657        else:
658            print line
659            Snames = ['S11','S22','S33','S12','S13','S23']
660            ptlbls = ' name  :'
661            ptstr =  ' value :'
662            sigstr = ' sig   :'
663            for i,name in enumerate(Snames):
664                ptlbls += '%12s' % (name)
665                ptstr += '%12.6f' % (hapData[4][i])
666                if sizeSig[1][i]:
667                    sigstr += '%12.6f' % (sizeSig[1][i])
668                else:
669                    sigstr += 12*' '
670            print ptlbls
671            print ptstr
672            print sigstr
673       
674    def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
675        line = '\n Mustrain model: '+hapData[0]
676        if hapData[0] in ['isotropic','uniaxial']:
677            line += ' equatorial:%12.4f'%(hapData[1][0])
678            if mustrainSig[0][0]:
679                line += ',sig: %12.4f'%(mustrainSig[0][0])
680            if hapData[0] == 'uniaxial':
681                line += ' axial:%12.4f'%(hapData[1][1])
682                if mustrainSig[0][1]:
683                     line += ',sig: %12.4f'%(mustrainSig[0][1])
684            print line
685        else:
686            print line
687            Snames = G2spc.MustrainNames(SGData)
688            ptlbls = ' name  :'
689            ptstr =  ' value :'
690            sigstr = ' sig   :'
691            for i,name in enumerate(Snames):
692                ptlbls += '%12s' % (name)
693                ptstr += '%12.6f' % (hapData[4][i])
694                sigstr += '%12.6f' % (mustrainSig[1][i])
695            print ptlbls
696            print ptstr
697            print sigstr
698       
699    for phase in Phases:
700        HistoPhase = Phases[phase]['Histograms']
701        SGData = Phases[phase]['General']['SGData']
702        pId = Phases[phase]['pId']
703        for histogram in HistoPhase:
704            print '\n Phase: ',phase,' in histogram: ',histogram
705            hapData = HistoPhase[histogram]
706            Histogram = Histograms[histogram]
707            hId = Histogram['hId']
708            pfx = str(pId)+':'+str(hId)+':'
709           
710# Add this stuff for Rietveld refinement!!!!           
711#            for item in ['Scale','Extinction']:
712#                hapDict[pfx+item] = hapData[item][0]
713#                if hapData[item][1]:
714#                    hapVary.append(pfx+item)
715#            controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
716#            if hapData['Pref.Ori.'][0] == 'MD':
717#                hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
718#                controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
719#                if hapData['Pref.Ori.'][2]:
720#                    hapVary.append(pfx+'MD')
721#            else:                           #'SH' spherical harmonics
722#                for item in hapData['Pref.Ori.'][5]:
723#                    hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
724#                    if hapData['Pref.Ori.'][2]:
725#                        hapVary.append(pfx+item)
726#            print '\n Phase fraction  : %10.4f'%(hapData['Scale'][0])
727#            print ' Extinction coeff: %10.4f'%(hapData['Extinction'][0])
728#            if hapData['Pref.Ori.'][0] == 'MD':
729#                Ax = hapData['Pref.Ori.'][3]
730#                print ' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
731#                    ' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
732               
733            SizeMuStrSig = {'Mustrain':[[0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
734                'Size':[[0,0],[0 for i in range(len(hapData['Size'][4]))]]}                 
735            for item in ['Mustrain','Size']:
736                if hapData[item][0] in ['isotropic','uniaxial']:
737                    hapData[item][1][0] = parmDict[pfx+item+':0']
738                    if hapData[item][2][0]: 
739                        SizeMuStrSig[item][0][0] = sigDict[pfx+item+':0']
740                    if hapData[item][0] == 'uniaxial':
741                        hapData[item][1][1] = parmDict[pfx+item+':1']
742                        if hapData[item][2][1]:
743                            SizeMuStrSig[item][0][1] = sigDict[pfx+item+':1']
744                else:       #generalized for mustrain or ellipsoidal for size
745                    for i in range(len(hapData[item][4])):
746                        sfx = ':'+str(i)
747                        hapData[item][4][i] = parmDict[pfx+item+sfx]
748                        if hapData[item][5][i]:
749                            SizeMuStrSig[item][1][i] = sigDict[pfx+item+sfx]
750                           
751            PrintSizeAndSig(hapData['Size'],SizeMuStrSig['Size'])
752            PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig['Mustrain'],SGData)
753   
754def GetHistogramData(Histograms):
755   
756    def GetBackgroundParms(hId,Background):
757        bakType,bakFlag = Background[:2]
758        backVals = Background[3:]
759        backNames = [':'+str(hId)+':Back:'+str(i) for i in range(len(backVals))]
760        if bakFlag:                                 #returns backNames as varyList = backNames
761            return bakType,dict(zip(backNames,backVals)),backNames
762        else:                                       #no background varied; varyList = []
763            return bakType,dict(zip(backNames,backVals)),[]
764       
765    def GetInstParms(hId,Inst):
766        insVals,insFlags,insNames = Inst[1:4]
767        dataType = insVals[0]
768        instDict = {}
769        insVary = []
770        pfx = ':'+str(hId)+':'
771        for i,flag in enumerate(insFlags):
772            insName = pfx+insNames[i]
773            instDict[insName] = insVals[i]
774            if flag:
775                insVary.append(insName)
776        instDict[pfx+'X'] = max(instDict[pfx+'X'],0.1)
777        instDict[pfx+'Y'] = max(instDict[pfx+'Y'],0.1)
778        instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.001)
779        return dataType,instDict,insVary
780       
781    def GetSampleParms(hId,Sample):
782        sampVary = []
783        pfx = ':'+str(hId)+':'       
784        sampDict = {pfx+'Gonio. radius':Sample['Gonio. radius']}
785        Type = Sample['Type']
786        if 'Bragg' in Type:             #Bragg-Brentano
787            for item in ['Scale','Shift','Transparency']:       #surface roughness?, diffuse scattering?
788                sampDict[pfx+item] = Sample[item][0]
789                if Sample[item][1]:
790                    sampVary.append(pfx+item)
791        elif 'Debye' in Type:        #Debye-Scherrer
792            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
793                sampDict[pfx+item] = Sample[item][0]
794                if Sample[item][1]:
795                    sampVary.append(pfx+item)
796        return Type,sampDict,sampVary
797       
798    def PrintBackground(Background):
799        print '\n Background function: ',Background[0],' Refine?',bool(Background[1])
800        line = ' Coefficients: '
801        for back in Background[3:]:
802            line += '%10.3f'%(back)
803        print line
804       
805    def PrintInstParms(Inst):
806        print '\n Instrument Parameters:'
807        ptlbls = ' name  :'
808        ptstr =  ' value :'
809        varstr = ' refine:'
810        instNames = Inst[3][1:]
811        for i,name in enumerate(instNames):
812            ptlbls += '%12s' % (name)
813            ptstr += '%12.6f' % (Inst[1][i+1])
814            if name in ['Lam1','Lam2','Azimuth']:
815                varstr += 12*' '
816            else:
817                varstr += '%12s' % (str(bool(Inst[2][i+1])))
818        print ptlbls
819        print ptstr
820        print varstr
821       
822    def PrintSampleParms(Sample):
823        print '\n Sample Parameters:'
824        ptlbls = ' name  :'
825        ptstr =  ' value :'
826        varstr = ' refine:'
827        if 'Bragg' in Sample['Type']:
828            for item in ['Scale','Shift','Transparency']:
829                ptlbls += '%14s'%(item)
830                ptstr += '%14.4f'%(Sample[item][0])
831                varstr += '%14s'%(str(bool(Sample[item][1])))
832           
833        elif 'Debye' in Type:        #Debye-Scherrer
834            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
835                ptlbls += '%14s'%(item)
836                ptstr += '%14.4f'%(Sample[item][0])
837                varstr += '%14s'%(str(bool(Sample[item][1])))
838
839        print ptlbls
840        print ptstr
841        print varstr
842       
843
844    histDict = {}
845    histVary = []
846    controlDict = {}
847    for histogram in Histograms:
848        Histogram = Histograms[histogram]
849        hId = Histogram['hId']
850        pfx = ':'+str(hId)+':'
851        controlDict[pfx+'Limits'] = Histogram['Limits'][1]
852       
853        Background = Histogram['Background'][0]
854        Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
855        controlDict[pfx+'bakType'] = Type
856        histDict.update(bakDict)
857        histVary += bakVary
858       
859        Inst = Histogram['Instrument Parameters']
860        Type,instDict,insVary = GetInstParms(hId,Inst)
861        controlDict[pfx+'histType'] = Type
862        histDict.update(instDict)
863        histVary += insVary
864       
865        Sample = Histogram['Sample Parameters']
866        Type,sampDict,sampVary = GetSampleParms(hId,Sample)
867        controlDict[pfx+'instType'] = Type
868        histDict.update(sampDict)
869        histVary += sampVary
870
871        print '\n Histogram: ',histogram,' histogram Id: ',hId
872        print 135*'-'
873        Units = {'C':' deg','T':' msec'}
874        units = Units[controlDict[pfx+'histType'][2]]
875        Limits = controlDict[pfx+'Limits']
876        print ' Instrument type: ',Sample['Type']
877        print ' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)     
878        PrintSampleParms(Sample)
879        PrintInstParms(Inst)
880        PrintBackground(Background)
881       
882    return histVary,histDict,controlDict
883   
884def SetHistogramData(parmDict,sigDict,Histograms):
885   
886    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
887        backVals = Background[3:]
888        backSig = [0 for i in range(len(Background[3:]))]
889        for i in range(len(backVals)):
890            backVals[i] = parmDict[pfx+'Back:'+str(i)]
891            if Background[1]:
892                backSig[i] = sigDict[pfx+'Back:'+str(i)]
893        return backSig
894       
895    def SetInstParms(pfx,Inst,parmDict,sigDict):
896        insVals,insFlags,insNames = Inst[1:4]
897        instSig = [0 for i in range(len(insVals))]
898        for i,flag in enumerate(insFlags):
899            insName = pfx+insNames[i]
900            insVals[i] = parmDict[insName]
901            if flag:
902                instSig[i] = sigDict[insName]
903        return instSig
904       
905    def SetSampleParms(pfx,Sample,parmDict,sigDict):
906        if 'Bragg' in Sample['Type']:             #Bragg-Brentano
907            sampSig = [0 for i in range(3)]
908            for i,item in enumerate(['Scale','Shift','Transparency']):       #surface roughness?, diffuse scattering?
909                Sample[item][0] = parmDict[pfx+item]
910                if Sample[item][1]:
911                    sampSig[i] = sigDict[pfx+item]
912        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
913            sampSig = [0 for i in range(4)]
914            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
915                Sample[item][0] = parmDict[pfx+item]
916                if Sample[item][1]:
917                    sampSig[i] = sigDict[pfx+item]
918        return sampSig
919       
920    def PrintBackgroundSig(Background,backSig):
921        print '\n Background function: ',Background[0]
922        valstr = ' value : '
923        sigstr = ' sig   : '
924        for i,back in enumerate(Background[3:]):
925            valstr += '%10.4f'%(back)
926            if Background[1]:
927                sigstr += '%10.4f'%(backSig[i])
928            else:
929                sigstr += 10*' '
930        print valstr
931        print sigstr
932       
933    def PrintInstParmsSig(Inst,instSig):
934        print '\n Instrument Parameters:'
935        ptlbls = ' names :'
936        ptstr =  ' value :'
937        sigstr = ' sig   :'
938        instNames = Inst[3][1:]
939        for i,name in enumerate(instNames):
940            ptlbls += '%12s' % (name)
941            ptstr += '%12.6f' % (Inst[1][i+1])
942            if instSig[i+1]:
943                sigstr += '%12.6f' % (instSig[i+1])
944            else:
945                sigstr += 12*' '
946        print ptlbls
947        print ptstr
948        print sigstr
949       
950    def PrintSampleParmsSig(Sample,sampleSig):
951        print '\n Sample Parameters:'
952        ptlbls = ' names :'
953        ptstr =  ' values:'
954        sigstr = ' sig   :'
955        if 'Bragg' in Sample['Type']:
956            for i,item in enumerate(['Scale','Shift','Transparency']):
957                ptlbls += '%14s'%(item)
958                ptstr += '%14.4f'%(Sample[item][0])
959                sigstr += '%14.4f'%(sampleSig[i])
960           
961        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
962            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
963                ptlbls += '%14s'%(item)
964                ptstr += '%14.4f'%(Sample[item][0])
965                sigstr += '%14.4f'%(sampleSig[i])
966
967        print ptlbls
968        print ptstr
969        print sigstr
970       
971    for histogram in Histograms:
972        if 'PWDR' in histogram:
973            Histogram = Histograms[histogram]
974            hId = Histogram['hId']
975            pfx = ':'+str(hId)+':'
976           
977            Background = Histogram['Background'][0]
978            backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
979           
980            Inst = Histogram['Instrument Parameters']
981            instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
982       
983            Sample = Histogram['Sample Parameters']
984            sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
985
986            print '\n Histogram: ',histogram,' histogram Id: ',hId
987            print 135*'-'
988            print ' Instrument type: ',Sample['Type']
989            PrintSampleParmsSig(Sample,sampSig)
990            PrintInstParmsSig(Inst,instSig)
991            PrintBackgroundSig(Background,backSig)
992
993
994def SetupSFcalc(General,Atoms):
995    ''' setup data for use in StructureFactor; mostly rearranging arrays
996    input:
997        General = dictionary of general phase info.
998        Atoms = list of atom parameters
999    returns:
1000        G = reciprocal metric tensor
1001        StrData = list of arrays; one entry per atom:
1002            T = atom types
1003            R = refinement flags, e.g. 'FXU'
1004            F = site fractions
1005            X = atom coordinates as numpy array
1006            M = atom multiplicities
1007            IA = isotropic/anisotropic thermal flags
1008            Uiso = isotropic thermal parameters if IA = 'I'; else = 0
1009            Uij = numpy array of 6 anisotropic thermal parameters if IA='A'; else = zeros
1010    '''           
1011    SGData = General['SGData']
1012    Cell = General['Cell']
1013    G,g = G2lat.cell2Gmat(Cell[1:7])        #skip refine & volume; get recip & real metric tensors
1014    cx,ct,cs,cia = General['AtomPtrs']
1015    X = [];F = [];T = [];IA = [];Uiso = [];Uij = [];R = [];M = []
1016    for atom in Atoms:
1017        T.append(atom[ct])
1018        R.append(atom[ct+1])
1019        F.append(atom[cx+3])
1020        X.append(np.array(atom[cx:cx+3]))
1021        M.append(atom[cia-1])
1022        IA.append(atom[cia])
1023        Uiso.append(atom[cia+1])
1024        Uij.append(np.array(atom[cia+2:cia+8]))
1025    return G,[T,R,np.array(F),np.array(X),np.array(M),IA,np.array(Uiso),np.array(Uij)]
1026           
1027def StructureFactor(H,G,SGData,StrData,FFtable):
1028    ''' Compute structure factor for a single h,k,l
1029    H = np.array(h,k,l)
1030    G = reciprocal metric tensor
1031    SGData = space group info. dictionary output from SpcGroup
1032    StrData = [
1033        [atom types],
1034        [refinement flags],
1035        [site fractions],
1036        np.array(coordinates),
1037        [multiplicities],
1038        [I/A flag],
1039        [Uiso],
1040        np.array(Uij)]
1041    FFtable = dictionary of form factor coeff. for atom types used in StrData
1042    '''
1043    twopi = 2.0*math.pi
1044    twopisq = 2.0*math.pi**2
1045    SQ = G2lat.calc_rDsq2(H,G)/4.0          # SQ = (sin(theta)/lambda)**2
1046    SQfactor = 8.0*SQ*math.pi**2
1047    print 'SQ',SQfactor
1048    FF = {}
1049    for type in FFtable.keys():
1050        FF[type] = G2el.ScatFac(FFtable[type],SQ)
1051    print 'FF',FF
1052    iabsnt,mulp,Uniq,Phs = G2spc.GenHKL(H,SGData,False)
1053    fa = [0,0,0]        #real
1054    fb = [0,0,0]        #imaginary
1055    if not iabsnt:
1056        phase = twopi*np.inner(Uniq,StrData[3])       #2pi[hx+ky+lz] for each atom in each equiv. hkl
1057        sinp = np.sin(phase)
1058        cosp = np.cos(phase)
1059        occ = StrData[2]*StrData[4]/mulp
1060        Tiso = np.multiply(StrData[6],SQfactor)
1061        Tuij = np.multiply(-SQfactor,np.inner(H,np.inner(G2spc.Uij2U(StrData[7]),H)))
1062        print 'sinp',sinp
1063        print 'cosp',cosp
1064        print 'occ',occ
1065        print 'Tiso',Tiso
1066        print 'Tuij',Tuij
1067    else:
1068        print 'Absent'
1069       
1070def Dict2Values(parmdict, varylist):
1071    '''Use before call to leastsq to setup list of values for the parameters
1072    in parmdict, as selected by key in varylist'''
1073    return [parmdict[key] for key in varylist] 
1074   
1075def Values2Dict(parmdict, varylist, values):
1076    ''' Use after call to leastsq to update the parameter dictionary with
1077    values corresponding to keys in varylist'''
1078    parmdict.update(zip(varylist,values))
1079   
1080def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
1081   
1082    def GetSampleGam(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse):
1083        if calcControls[phfx+'SizeType'] == 'isotropic':
1084            gam = 1.8*wave/(np.pi*parmDict[phfx+'Size:0']*cosd(refl[5]/2.))
1085        elif calcControls[phfx+'SizeType'] == 'uniaxial':
1086            H = np.array(refl[:3])
1087            P = np.array(calcControls[phfx+'SizeAxis'])
1088            cosP,sinP = G2lat.CosSinAngle(H,V,P)
1089            gam = (1.8*wave/np.pi)*parmDict[phfx+'Size:0']*parmDict[phfx+'Size:1']
1090            gam /= np.sqrt((cosP*parmDict[phfx+'Size:1'])**2+(sinP*parmDict[phfx+'Size:0'])**2)*cosd(refl[5]/2.)
1091        else:           #ellipsoidal crystallites
1092            H = np.array(refl[:3])
1093            gam += 1.8*wave/(np.pi*cosd(refl[5])*np.inner(H,np.inner(sizeEllipse,H)))           
1094        if calcControls[phfx+'MustrainType'] == 'isotropic':
1095            gam += 0.018*parmDict[phfx+'Mustrain:0']*tand(refl[5]/2.)/np.pi
1096        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1097            H = np.array(refl[:3])
1098            P = np.array(calcControls[phfx+'MustrainAxis'])
1099            cosP,sinP = G2lat.CosSinAngle(H,V,P)
1100            Si = parmDict[phfx+'Mustrain:0']
1101            Sa = parmDict[phfx+'Mustrain:1']
1102            gam += 0.018*Si*Sa*tand(refl[5]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
1103        else:       #generalized - P.W. Stephens model
1104            pwrs = calcControls[phfx+'MuPwrs']
1105            sum = 0
1106            for i,pwr in enumerate(pwrs):
1107                sum += parmDict[phfx+'Mustrain:'+str(i)]*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
1108            gam += 0.018*refl[4]**2*tand(refl[5]/2.)*sum
1109           
1110        return gam 
1111       
1112       
1113    hId = Histogram['hId']
1114    hfx = ':%d:'%(hId)
1115    bakType = calcControls[hfx+'bakType']
1116    yb = G2pwd.getBackground(hfx,parmDict,bakType,x)
1117    yc = np.zeros_like(yb)
1118       
1119    if 'C' in calcControls[hfx+'histType']:   
1120        U = parmDict[hfx+'U']
1121        V = parmDict[hfx+'V']
1122        W = parmDict[hfx+'W']
1123        X = parmDict[hfx+'X']
1124        Y = parmDict[hfx+'Y']
1125        shl = max(parmDict[hfx+'SH/L'],0.001)
1126        Zero = parmDict[hfx+'Zero']
1127        pola = parmDict[hfx+'Polariz.']
1128        Ka2 = False
1129        if hfx+'Lam1' in parmDict.keys():
1130            wave = parmDict[hfx+'Lam1']
1131            Ka2 = True
1132            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1133            kRatio = parmDict[hfx+'I(L2)/I(L1)']
1134        else:
1135            wave = parmDict[hfx+'Lam']
1136    else:
1137        print 'TOF Undefined at present'
1138        raise ValueError
1139    for phase in Histogram['Reflection Lists']:
1140        refList = Histogram['Reflection Lists'][phase]
1141        Phase = Phases[phase]
1142        pId = Phase['pId']
1143        pfx = '%d::'%(pId)
1144        phfx = '%d:%d:'%(pId,hId)
1145        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1146        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1147        sizeEllipse = []
1148        if calcControls[phfx+'SizeType'] == 'ellipsoidal':
1149            sizeEllipse = G2lat.U6toUij([parmDIct[phfx+'Size:%d'%(i)] for i in range(6)])
1150        for refl in refList:
1151            if 'C' in calcControls[hfx+'histType']:
1152                h,k,l,m = refl[:4]
1153                dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G)
1154                d = np.sqrt(dsq)
1155                pos = 2.0*asind(wave/(2.0*d))+Zero
1156                const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1157                if 'Bragg' in calcControls[hfx+'instType']:
1158                    pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
1159                        1.e-7*parmDict[hfx+'Transparency']*sind(pos))            #trans(=1/mueff) in Angstroms
1160                else:               #Debye-Scherrer - simple but maybe not right
1161                    pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
1162                refl[5] = pos
1163                tanPos = tand(pos/2.0)
1164                refl[6] = U*tanPos**2+V*tanPos+W        #save peak sigma
1165                refl[7] = X/cosd(pos/2.0)+Y*tanPos+GetSampleGam(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse) #save peak gamma
1166                if 'Pawley' in Phase['General']['Type']:
1167                    try:
1168                        refl[8] = parmDict[hfx+'Scale']*m*parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])]
1169                    except KeyError:
1170                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
1171                        continue
1172                else:
1173                    raise ValueError       #wants strctrfacr calc here
1174                Wd,fmin,fmax = G2pwd.getWidths(pos,refl[6],refl[7],shl)
1175                if pos > 90:
1176                    fmin,fmax = [fmax,fmin]         
1177                iBeg = np.searchsorted(x,pos-fmin)
1178                iFin = np.searchsorted(x,pos+fmax)
1179                if not iBeg+iFin:       #peak below low limit - skip peak
1180                    continue
1181                elif not iBeg-iFin:     #peak above high limit - done
1182                    return yc,yb
1183                yc[iBeg:iFin] += G2pwd.getFCJVoigt(pos,refl[8],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
1184                if Ka2:
1185                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
1186                    Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
1187                    if pos > 90:
1188                        fmin,fmax = [fmax,fmin]         
1189                    iBeg = np.searchsorted(x,pos2-fmin)
1190                    iFin = np.searchsorted(x,pos2+fmax)
1191                    yc[iBeg:iFin] += G2pwd.getFCJVoigt(pos2,refl[8]*kRatio,refl[6],refl[7],shl,x[iBeg:iFin])        #and here
1192            else:
1193                raise ValueError
1194    return yc,yb   
1195           
1196
1197   
1198def Refine(GPXfile):
1199   
1200    def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):       
1201        parmdict.update(zip(varylist,values))
1202        G2mv.Dict2Map(parmDict)
1203        Histograms,Phases = HistoPhases
1204        M = np.empty(0)
1205        sumwYo = 0
1206        Nobs = 0
1207        for histogram in Histograms:
1208            if 'PWDR' in histogram[:4]:
1209                Histogram = Histograms[histogram]
1210                hId = Histogram['hId']
1211                hfx = ':%d:'%(hId)
1212                Limits = calcControls[hfx+'Limits']
1213                x,y,w,yc,yb,yd = Histogram['Data']
1214                yc *= 0.0                           #zero full calcd profiles
1215                yb *= 0.0
1216                yd *= 0.0
1217                xB = np.searchsorted(x,Limits[0])
1218                xF = np.searchsorted(x,Limits[1])
1219                Nobs += xF-xB
1220                Histogram['sumwYo'] = np.sum(w[xB:xF]*y[xB:xF]**2)
1221                sumwYo += Histogram['sumwYo']
1222                yc[xB:xF],yb[xB:xF] = getPowderProfile(parmdict,x[xB:xF],
1223                    varylist,Histogram,Phases,calcControls,pawleyLookup)
1224                yc[xB:xF] *= parmDict[hfx+'Scale']
1225                yc[xB:xF] += yb[xB:xF]
1226                yd[xB:xF] = y[xB:xF]-yc[xB:xF]
1227                Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF]))
1228                M = np.concatenate((M,np.sqrt(w[xB:xF])*(yd[xB:xF])))
1229#                print 'sum M^2,y,yb,yc',np.sum(M**2),np.sum(y[xB:xF]),np.sum(yb[xB:xF]),np.sum(yc[xB:xF])
1230        Histograms['sumwYo'] = sumwYo
1231        Histograms['Nobs'] = Nobs
1232        Rwp = min(100.,np.sqrt(np.sum(M**2)/sumwYo)*100.)
1233        if dlg:
1234            GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Powder profile Rwp =',Rwp,'%'))[0]
1235            if not GoOn:
1236                return -M           #abort!!
1237        return M
1238   
1239    ShowBanner()
1240    varyList = []
1241    parmDict = {}
1242    calcControls = {}   
1243#    Controls = GetControls(GPXfile)
1244#    ShowControls(Controls)           
1245    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
1246    if not Phases:
1247        print ' *** ERROR - you have no phases to refine! ***'
1248        print ' *** Refine aborted ***'
1249        raise Exception
1250    if not Histograms:
1251        print ' *** ERROR - you have no data to refine with! ***'
1252        print ' *** Refine aborted ***'
1253        raise Exception
1254    phaseVary,phaseDict,pawleyLookup = GetPhaseData(Phases)
1255    hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms)
1256    calcControls.update(controlDict)
1257    histVary,histDict,controlDict = GetHistogramData(Histograms)
1258    calcControls.update(controlDict)
1259    varyList = phaseVary+histVary+hapVary
1260    parmDict.update(phaseDict)
1261    parmDict.update(hapDict)
1262    parmDict.update(histDict)
1263    constrDict,constrFlag,fixedList = G2mv.InputParse([])        #constraints go here?
1264    groups,parmlist = G2mv.GroupConstraints(constrDict)
1265    G2mv.GenerateConstraints(groups,parmlist,constrDict,constrFlag,fixedList)
1266    G2mv.Map2Dict(parmDict,varyList)
1267
1268    while True:
1269        begin = time.time()
1270        values =  np.array(Dict2Values(parmDict, varyList))
1271        dlg = wx.ProgressDialog('Residual','Powder profile Rwp =',101.0, 
1272            style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_REMAINING_TIME|wx.PD_CAN_ABORT)
1273        screenSize = wx.ClientDisplayRect()
1274        Size = dlg.GetSize()
1275        dlg.SetPosition(wx.Point(screenSize[2]-Size[0]-305,screenSize[1]+5))
1276        try:
1277            result = so.leastsq(errRefine,values,full_output=True,  #factor=1.,epsfcn=0.00001,ftol=0.0001,
1278                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
1279            print len(result)
1280        finally:
1281            dlg.Destroy()
1282        runtime = time.time()-begin   
1283        chisq = np.sum(result[2]['fvec']**2)
1284        ncyc = int(result[2]['nfev']/len(varyList))
1285        Values2Dict(parmDict, varyList, result[0])
1286        G2mv.Dict2Map(parmDict)
1287        Rwp = np.sqrt(chisq/Histograms['sumwYo'])*100.      #to %
1288        GOF = chisq/(Histograms['Nobs']-len(varyList))
1289        print '\n Refinement results:'
1290        print 135*'-'
1291        print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',Histograms['Nobs'],' Number of parameters: ',len(varyList)
1292        print 'Refinement time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
1293        print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
1294        try:
1295            sig = np.sqrt(np.diag(result[1])*GOF)
1296            if np.any(np.isnan(sig)):
1297                print '*** Least squares aborted - some invalid esds possible ***'
1298            table = dict(zip(varyList,zip(values,result[0],(result[0]-values)/sig)))
1299#            for item in table: print item,table[item]               #useful debug - are things shifting?
1300            break                   #refinement succeeded - finish up!
1301        except ValueError:          #result[1] is None on singular matrix
1302            print '**** Refinement failed - singular matrix ****'
1303            Ipvt = result[2]['ipvt']
1304            for i,ipvt in enumerate(Ipvt):
1305                if not np.sum(result[2]['fjac'],axis=1)[i]:
1306                    print 'Removing parameter: ',varyList[ipvt-1]
1307                    del(varyList[ipvt-1])
1308                    break
1309
1310    sigDict = dict(zip(varyList,sig))
1311    SetPhaseData(parmDict,sigDict,Phases)
1312    SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms)
1313    SetHistogramData(parmDict,sigDict,Histograms)
1314    SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases)
1315
1316def main():
1317    arg = sys.argv
1318    if len(arg) > 1:
1319        GPXfile = arg[1]
1320        if not ospath.exists(GPXfile):
1321            print 'ERROR - ',GPXfile," doesn't exist!"
1322            exit()
1323        GPXpath = ospath.dirname(arg[1])
1324        Refine(GPXfile)
1325    else:
1326        print 'ERROR - missing filename'
1327        exit()
1328    print "Done"
1329         
1330if __name__ == '__main__':
1331    main()
Note: See TracBrowser for help on using the repository browser.