source: trunk/GSASIIstruct.py @ 407

Last change on this file since 407 was 407, checked in by vondreele, 10 years ago

constraints round 2

  • Property svn:keywords set to Date Author Revision URL Id
File size: 114.0 KB
Line 
1#GSASIIstructure - structure computation routines
2########### SVN repository information ###################
3# $Date: 2011-11-05 22:30:24 +0000 (Sat, 05 Nov 2011) $
4# $Author: vondreele $
5# $Revision: 407 $
6# $URL: trunk/GSASIIstruct.py $
7# $Id: GSASIIstruct.py 407 2011-11-05 22:30:24Z vondreele $
8########### SVN repository information ###################
9import sys
10import os
11import os.path as ospath
12import numpy as np
13import numpy.linalg as nl
14import cPickle
15import time
16import math
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
29atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
30
31
32def ShowBanner():
33    print 80*'*'
34    print '   General Structure Analysis System-II Crystal Structure Refinement'
35    print '     by Robert B. Von Dreele, Argonne National Laboratory(C), 2010'
36    print ' This product includes software developed by the UChicago Argonne, LLC,' 
37    print '            as Operator of Argonne National Laboratory.'
38    print 80*'*','\n'
39
40def GetControls(GPXfile):
41    ''' Returns dictionary of control items found in GSASII gpx file
42    input:
43        GPXfile = .gpx full file name
44    return:
45        Controls = dictionary of control items
46    '''
47    Controls = {'deriv type':'analytical','min dM/M':0.0001,'shift factor':1.}
48    file = open(GPXfile,'rb')
49    while True:
50        try:
51            data = cPickle.load(file)
52        except EOFError:
53            break
54        datum = data[0]
55        if datum[0] == 'Controls':
56            Controls.update(datum[1])
57    file.close()
58    return Controls
59   
60def ShowControls(Controls):
61    print ' Least squares controls:'
62    print ' Derivative type: ',Controls['deriv type']
63    print ' Minimum delta-M/M for convergence: ','%.2g'%(Controls['min dM/M'])
64    print ' Initial shift factor: ','%.3f'%(Controls['shift factor'])
65   
66def GetConstraints(GPXfile):
67    constrList = []
68    file = open(GPXfile,'rb')
69    while True:
70        try:
71            data = cPickle.load(file)
72        except EOFError:
73            break
74        datum = data[0]
75        if datum[0] == 'Constraints':
76            constDict = datum[1]
77            for item in constDict:
78                constrList += constDict[item]
79    file.close()
80    return constrList
81   
82def GetPhaseNames(GPXfile):
83    ''' Returns a list of phase names found under 'Phases' in GSASII gpx file
84    input:
85        GPXfile = gpx full file name
86    return:
87        PhaseNames = list of phase names
88    '''
89    file = open(GPXfile,'rb')
90    PhaseNames = []
91    while True:
92        try:
93            data = cPickle.load(file)
94        except EOFError:
95            break
96        datum = data[0]
97        if 'Phases' == datum[0]:
98            for datus in data[1:]:
99                PhaseNames.append(datus[0])
100    file.close()
101    return PhaseNames
102
103def GetAllPhaseData(GPXfile,PhaseName):
104    ''' Returns the entire dictionary for PhaseName from GSASII gpx file
105    input:
106        GPXfile = gpx full file name
107        PhaseName = phase name
108    return:
109        phase dictionary
110    '''       
111    file = open(GPXfile,'rb')
112    General = {}
113    Atoms = []
114    while True:
115        try:
116            data = cPickle.load(file)
117        except EOFError:
118            break
119        datum = data[0]
120        if 'Phases' == datum[0]:
121            for datus in data[1:]:
122                if datus[0] == PhaseName:
123                    break
124    file.close()
125    return datus[1]
126   
127def GetHistograms(GPXfile,hNames):
128    """ Returns a dictionary of histograms found in GSASII gpx file
129    input:
130        GPXfile = .gpx full file name
131        hNames = list of histogram names
132    return:
133        Histograms = dictionary of histograms (types = PWDR & HKLF)
134    """
135    file = open(GPXfile,'rb')
136    Histograms = {}
137    while True:
138        try:
139            data = cPickle.load(file)
140        except EOFError:
141            break
142        datum = data[0]
143        hist = datum[0]
144        if hist in hNames:
145            if 'PWDR' in hist[:4]:
146                PWDRdata = {}
147                PWDRdata['Data'] = datum[1][1]          #powder data arrays
148                PWDRdata[data[2][0]] = data[2][1]       #Limits
149                PWDRdata[data[3][0]] = data[3][1]       #Background
150                PWDRdata[data[4][0]] = data[4][1]       #Instrument parameters
151                PWDRdata[data[5][0]] = data[5][1]       #Sample parameters
152                try:
153                    PWDRdata[data[9][0]] = data[9][1]       #Reflection lists might be missing
154                except IndexError:
155                    PWDRdata['Reflection lists'] = {}
156   
157                Histograms[hist] = PWDRdata
158            elif 'HKLF' in hist[:4]:
159                HKLFdata = []
160                datum = data[0]
161                HKLFdata = datum[1:][0]
162                Histograms[hist] = HKLFdata           
163    file.close()
164    return Histograms
165   
166def GetHistogramNames(GPXfile,hType):
167    """ Returns a list of histogram names found in GSASII gpx file
168    input:
169        GPXfile = .gpx full file name
170        hType = list ['PWDR','HKLF']
171    return:
172        HistogramNames = list of histogram names (types = PWDR & HKLF)
173    """
174    file = open(GPXfile,'rb')
175    HistogramNames = []
176    while True:
177        try:
178            data = cPickle.load(file)
179        except EOFError:
180            break
181        datum = data[0]
182        if datum[0][:4] in hType:
183            HistogramNames.append(datum[0])
184    file.close()
185    return HistogramNames
186   
187def GetUsedHistogramsAndPhases(GPXfile):
188    ''' Returns all histograms that are found in any phase
189    and any phase that uses a histogram
190    input:
191        GPXfile = .gpx full file name
192    return:
193        Histograms = dictionary of histograms as {name:data,...}
194        Phases = dictionary of phases that use histograms
195    '''
196    phaseNames = GetPhaseNames(GPXfile)
197    histoList = GetHistogramNames(GPXfile,['PWDR','HKLF'])
198    allHistograms = GetHistograms(GPXfile,histoList)
199    phaseData = {}
200    for name in phaseNames: 
201        phaseData[name] =  GetAllPhaseData(GPXfile,name)
202    Histograms = {}
203    Phases = {}
204    for phase in phaseData:
205        Phase = phaseData[phase]
206        if Phase['Histograms']:
207            if phase not in Phases:
208                pId = phaseNames.index(phase)
209                Phase['pId'] = pId
210                Phases[phase] = Phase
211            for hist in Phase['Histograms']:
212                if hist not in Histograms:
213                    Histograms[hist] = allHistograms[hist]
214                    #future restraint, etc. histograms here           
215                    hId = histoList.index(hist)
216                    Histograms[hist]['hId'] = hId
217    return Histograms,Phases
218   
219def GPXBackup(GPXfile,makeBack=True):
220    import distutils.file_util as dfu
221    GPXpath,GPXname = ospath.split(GPXfile)
222    Name = ospath.splitext(GPXname)[0]
223    files = os.listdir(GPXpath)
224    last = 0
225    for name in files:
226        name = name.split('.')
227        if len(name) == 3 and name[0] == Name and 'bak' in name[1]:
228            if makeBack:
229                last = max(last,int(name[1].strip('bak'))+1)
230            else:
231                last = max(last,int(name[1].strip('bak')))
232    GPXback = ospath.join(GPXpath,ospath.splitext(GPXname)[0]+'.bak'+str(last)+'.gpx')
233    dfu.copy_file(GPXfile,GPXback)
234    return GPXback
235       
236def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,CovData,makeBack=True):
237    ''' Updates gpxfile from all histograms that are found in any phase
238    and any phase that used a histogram
239    input:
240        GPXfile = .gpx full file name
241        Histograms = dictionary of histograms as {name:data,...}
242        Phases = dictionary of phases that use histograms
243        CovData = dictionary of refined variables, varyList, & covariance matrix
244        makeBack = True if new backup of .gpx file is to be made; else use the last one made
245    '''
246                       
247    GPXback = GPXBackup(GPXfile,makeBack)
248    print '\n',135*'-'
249    print 'Read from file:',GPXback
250    print 'Save to file  :',GPXfile
251    infile = open(GPXback,'rb')
252    outfile = open(GPXfile,'wb')
253    while True:
254        try:
255            data = cPickle.load(infile)
256        except EOFError:
257            break
258        datum = data[0]
259#        print 'read: ',datum[0]
260        if datum[0] == 'Phases':
261            for iphase in range(len(data)):
262                if data[iphase][0] in Phases:
263                    phaseName = data[iphase][0]
264                    data[iphase][1].update(Phases[phaseName])
265        elif datum[0] == 'Covariance':
266            data[0][1] = CovData
267        try:
268            histogram = Histograms[datum[0]]
269#            print 'found ',datum[0]
270            data[0][1][1] = histogram['Data']
271            for datus in data[1:]:
272#                print '    read: ',datus[0]
273                if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
274                    datus[1] = histogram[datus[0]]
275        except KeyError:
276            pass
277                               
278        cPickle.dump(data,outfile,1)
279    infile.close()
280    outfile.close()
281    print 'refinement save successful'
282   
283def SetSeqResult(GPXfile,Histograms,SeqResult):
284    GPXback = GPXBackup(GPXfile)
285    print '\n',135*'-'
286    print 'Read from file:',GPXback
287    print 'Save to file  :',GPXfile
288    infile = open(GPXback,'rb')
289    outfile = open(GPXfile,'wb')
290    while True:
291        try:
292            data = cPickle.load(infile)
293        except EOFError:
294            break
295        datum = data[0]
296        if datum[0] == 'Sequental results':
297            data[0][1] = SeqResult
298        try:
299            histogram = Histograms[datum[0]]
300            data[0][1][1] = histogram['Data']
301            for datus in data[1:]:
302                if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
303                    datus[1] = histogram[datus[0]]
304        except KeyError:
305            pass
306                               
307        cPickle.dump(data,outfile,1)
308    infile.close()
309    outfile.close()
310    print 'refinement save successful'
311   
312                   
313def GetPWDRdata(GPXfile,PWDRname):
314    ''' Returns powder data from GSASII gpx file
315    input:
316        GPXfile = .gpx full file name
317        PWDRname = powder histogram name as obtained from GetHistogramNames
318    return:
319        PWDRdata = powder data dictionary with:
320            Data - powder data arrays, Limits, Instrument Parameters, Sample Parameters
321       
322    '''
323    file = open(GPXfile,'rb')
324    PWDRdata = {}
325    while True:
326        try:
327            data = cPickle.load(file)
328        except EOFError:
329            break
330        datum = data[0]
331        if datum[0] == PWDRname:
332            PWDRdata['Data'] = datum[1][1]          #powder data arrays
333            PWDRdata[data[2][0]] = data[2][1]       #Limits
334            PWDRdata[data[3][0]] = data[3][1]       #Background
335            PWDRdata[data[4][0]] = data[4][1]       #Instrument parameters
336            PWDRdata[data[5][0]] = data[5][1]       #Sample parameters
337            try:
338                PWDRdata[data[9][0]] = data[9][1]       #Reflection lists might be missing
339            except IndexError:
340                PWDRdata['Reflection lists'] = {}
341    file.close()
342    return PWDRdata
343   
344def GetHKLFdata(GPXfile,HKLFname):
345    ''' Returns single crystal data from GSASII gpx file
346    input:
347        GPXfile = .gpx full file name
348        HKLFname = single crystal histogram name as obtained from GetHistogramNames
349    return:
350        HKLFdata = single crystal data list of reflections: for each reflection:
351            HKLF = [np.array([h,k,l]),FoSq,sigFoSq,FcSq,Fcp,Fcpp,phase]
352    '''
353    file = open(GPXfile,'rb')
354    HKLFdata = []
355    while True:
356        try:
357            data = cPickle.load(file)
358        except EOFError:
359            break
360        datum = data[0]
361        if datum[0] == HKLFname:
362            HKLFdata = datum[1:][0]
363    file.close()
364    return HKLFdata
365   
366def GetFFtable(General):
367    ''' returns a dictionary of form factor data for atom types found in General
368    input:
369        General = dictionary of phase info.; includes AtomTypes
370    return:
371        FFtable = dictionary of form factor data; key is atom type
372    '''
373    atomTypes = General['AtomTypes']
374    FFtable = {}
375    for El in atomTypes:
376        FFs = G2el.GetFormFactorCoeff(El.split('+')[0].split('-')[0])
377        for item in FFs:
378            if item['Symbol'] == El.upper():
379                FFtable[El] = item
380    return FFtable
381   
382def GetPawleyConstr(SGLaue,PawleyRef,pawleyVary):
383    if SGLaue in ['-1','2/m','mmm']:
384        return                      #no Pawley symmetry required constraints
385    for i,varyI in enumerate(pawleyVary):
386        refI = int(varyI.split(':')[-1])
387        ih,ik,il = PawleyRef[refI][:3]
388        for varyJ in pawleyVary[0:i]:
389            refJ = int(varyJ.split(':')[-1])
390            jh,jk,jl = PawleyRef[refJ][:3]
391            if SGLaue in ['4/m','4/mmm']:
392                isum = ih**2+ik**2
393                jsum = jh**2+jk**2
394                if abs(il) == abs(jl) and isum == jsum:
395                    G2mv.StoreEquivalence(varyJ,(varyI,))
396            elif SGLaue in ['3R','3mR']:
397                isum = ih**2+ik**2+il**2
398                jsum = jh**2+jk**2*jl**2
399                isum2 = ih*ik+ih*il+ik*il
400                jsum2 = jh*jk+jh*jl+jk*jl
401                if isum == jsum and isum2 == jsum2:
402                    G2mv.StoreEquivalence(varyJ,(varyI,))
403            elif SGLaue in ['3','3m1','31m','6/m','6/mmm']:
404                isum = ih**2+ik**2+ih*ik
405                jsum = jh**2+jk**2+jh*jk
406                if abs(il) == abs(jl) and isum == jsum:
407                    G2mv.StoreEquivalence(varyJ,(varyI,))
408            elif SGLaue in ['m3','m3m']:
409                isum = ih**2+ik**2+il**2
410                jsum = jh**2+jk**2+jl**2
411                if isum == jsum:
412                    G2mv.StoreEquivalence(varyJ,(varyI,))
413                   
414def cellVary(pfx,SGData): 
415    if SGData['SGLaue'] in ['-1',]:
416        return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
417    elif SGData['SGLaue'] in ['2/m',]:
418        if SGData['SGUniq'] == 'a':
419            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3']
420        elif SGData['SGUniq'] == 'b':
421            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A4']
422        else:
423            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A5']
424    elif SGData['SGLaue'] in ['mmm',]:
425        return [pfx+'A0',pfx+'A1',pfx+'A2']
426    elif SGData['SGLaue'] in ['4/m','4/mmm']:
427        return [pfx+'A0',pfx+'A2']
428    elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
429        return [pfx+'A0',pfx+'A2']
430    elif SGData['SGLaue'] in ['3R', '3mR']:
431        return [pfx+'A0',pfx+'A3']                       
432    elif SGData['SGLaue'] in ['m3m','m3']:
433        return [pfx+'A0',]
434       
435           
436def GetPhaseData(PhaseData,Print=True):
437           
438    def PrintFFtable(FFtable):
439        print '\n Scattering factors:'
440        print '   Symbol     fa                                      fb                                      fc'
441        print 99*'-'
442        for Ename in FFtable:
443            ffdata = FFtable[Ename]
444            fa = ffdata['fa']
445            fb = ffdata['fb']
446            print ' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f' %  \
447                (Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],ffdata['fc'])
448               
449    def PrintAtoms(General,Atoms):
450        print '\n Atoms:'
451        line = '   name    type  refine?   x         y         z    '+ \
452            '  frac site sym  mult I/A   Uiso     U11     U22     U33     U12     U13     U23'
453        if General['Type'] == 'magnetic':
454            line += '   Mx     My     Mz'
455        elif General['Type'] == 'macromolecular':
456            line = ' res no  residue  chain '+line
457        print line
458        if General['Type'] == 'nuclear':
459            print 135*'-'
460            for i,at in enumerate(Atoms):
461                line = '%7s'%(at[0])+'%7s'%(at[1])+'%7s'%(at[2])+'%10.5f'%(at[3])+'%10.5f'%(at[4])+ \
462                    '%10.5f'%(at[5])+'%8.3f'%(at[6])+'%7s'%(at[7])+'%5d'%(at[8])+'%5s'%(at[9])
463                if at[9] == 'I':
464                    line += '%8.4f'%(at[10])+48*' '
465                else:
466                    line += 8*' '
467                    for j in range(6):
468                        line += '%8.4f'%(at[11+j])
469                print line
470       
471    def PrintTexture(textureData):
472        print '\n Spherical harmonics texture: Order:' + \
473            str(textureData['Order'])+' Refine? '+str(textureData['SH Coeff'][0])
474        names = ['omega','chi','phi']
475        line = '\n'
476        for name in names:
477            line += ' SH '+name+':'+'%12.4f'%(textureData['Sample '+name][1])+' Refine? '+str(textureData['Sample '+name][0])
478        print line
479        print '\n Texture coefficients:'
480        ptlbls = ' names :'
481        ptstr =  ' values:'
482        SHcoeff = textureData['SH Coeff'][1]
483        for item in SHcoeff:
484            ptlbls += '%12s'%(item)
485            ptstr += '%12.4f'%(SHcoeff[item]) 
486        print ptlbls
487        print ptstr   
488       
489    if Print: print ' Phases:'
490    phaseVary = []
491    phaseDict = {}
492    phaseConstr = {}
493    pawleyLookup = {}
494    FFtables = {}
495    Natoms = {}
496    AtMults = {}
497    AtIA = {}
498    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
499    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
500    for name in PhaseData:
501        General = PhaseData[name]['General']
502        pId = PhaseData[name]['pId']
503        pfx = str(pId)+'::'
504        FFtable = GetFFtable(General)
505        FFtables.update(FFtable)
506        Atoms = PhaseData[name]['Atoms']
507        try:
508            PawleyRef = PhaseData[name]['Pawley ref']
509        except KeyError:
510            PawleyRef = []
511        if Print: print '\n Phase name: ',General['Name']
512        SGData = General['SGData']
513        SGtext = G2spc.SGPrint(SGData)
514        if Print: 
515            for line in SGtext: print line
516        cell = General['Cell']
517        A = G2lat.cell2A(cell[1:7])
518        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]})
519        if cell[0]:
520            phaseVary += cellVary(pfx,SGData)
521        if Print: print '\n Unit cell: a =','%.5f'%(cell[1]),' b =','%.5f'%(cell[2]),' c =','%.5f'%(cell[3]), \
522            ' alpha =','%.3f'%(cell[4]),' beta =','%.3f'%(cell[5]),' gamma =', \
523            '%.3f'%(cell[6]),' volume =','%.3f'%(cell[7]),' Refine?',cell[0]
524        if Print and FFtable:
525            PrintFFtable(FFtable)
526        Natoms[pfx] = 0
527        if Atoms:
528            if General['Type'] == 'nuclear':
529                Natoms[pfx] = len(Atoms)
530                for i,at in enumerate(Atoms):
531                    phaseDict.update({pfx+'Atype:'+str(i):at[1],pfx+'Afrac:'+str(i):at[6],pfx+'Amul:'+str(i):at[8],
532                        pfx+'Ax:'+str(i):at[3],pfx+'Ay:'+str(i):at[4],pfx+'Az:'+str(i):at[5],
533                        pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx+'dAz:'+str(i):0.,         #refined shifts for x,y,z
534                        pfx+'AI/A:'+str(i):at[9],})
535                    if at[9] == 'I':
536                        phaseDict[pfx+'AUiso:'+str(i)] = at[10]
537                    else:
538                        phaseDict.update({pfx+'AU11:'+str(i):at[11],pfx+'AU22:'+str(i):at[12],pfx+'AU33:'+str(i):at[13],
539                            pfx+'AU12:'+str(i):at[14],pfx+'AU13:'+str(i):at[15],pfx+'AU23:'+str(i):at[16]})
540                    if 'F' in at[2]:
541                        phaseVary.append(pfx+'Afrac:'+str(i))
542                    if 'X' in at[2]:
543                        xId,xCoef = G2spc.GetCSxinel(at[7])
544                        delnames = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)]
545                        for j in range(3):
546                            if xId[j] > 0:                               
547                                phaseVary.append(delnames[j])
548                                for k in range(j):
549                                    if xId[j] == xId[k]:
550                                        G2mv.StoreEquivalence(delnames[k],((delnames[j],xCoef[j]),)) 
551                    if 'U' in at[2]:
552                        if at[9] == 'I':
553                            phaseVary.append(pfx+'AUiso:'+str(i))
554                        else:
555                            uId,uCoef = G2spc.GetCSuinel(at[7])[:2]
556                            names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i),
557                                pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)]
558                            for j in range(6):
559                                if uId[j] > 0:                               
560                                    phaseVary.append(names[j])
561                                    for k in range(j):
562                                        if uId[j] == uId[k]:
563                                            G2mv.StoreEquivalence(names[k],((names[j],uCoef[j]),))
564#            elif General['Type'] == 'magnetic':
565#            elif General['Type'] == 'macromolecular':
566
567            if 'SH Texture' in General:
568                textureData = General['SH Texture']
569                phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']]
570                phaseDict[pfx+'SHorder'] = textureData['Order']
571                for name in ['omega','chi','phi']:
572                    phaseDict[pfx+'SH '+name] = textureData['Sample '+name][1]
573                    if textureData['Sample '+name][0]:
574                        phaseVary.append(pfx+'SH '+name)
575                for name in textureData['SH Coeff'][1]:
576                    phaseDict[pfx+name] = textureData['SH Coeff'][1][name]
577                    if textureData['SH Coeff'][0]:
578                        phaseVary.append(pfx+name)
579               
580            if Print:
581                PrintAtoms(General,Atoms)
582                if 'SH Texture' in General:
583                    PrintTexture(textureData)
584                   
585        elif PawleyRef:
586            pawleyVary = []
587            for i,refl in enumerate(PawleyRef):
588                phaseDict[pfx+'PWLref:'+str(i)] = refl[6]
589                pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
590                if refl[5]:
591                    pawleyVary.append(pfx+'PWLref:'+str(i))
592            GetPawleyConstr(SGData['SGLaue'],PawleyRef,pawleyVary)      #does G2mv.StoreEquivalence
593            phaseVary += pawleyVary
594               
595    return Natoms,phaseVary,phaseDict,pawleyLookup,FFtables
596   
597def getVCov(varyNames,varyList,covMatrix):
598    vcov = np.zeros((len(varyNames),len(varyNames)))
599    for i1,name1 in enumerate(varyNames):
600        for i2,name2 in enumerate(varyNames):
601            try:
602                vcov[i1][i2] = covMatrix[varyList.index(name1)][varyList.index(name2)]
603            except ValueError:
604                vcov[i1][i2] = 0.0
605    return vcov
606   
607def getCellEsd(pfx,SGData,A,covData):
608    dpr = 180./np.pi
609    rVsq = G2lat.calc_rVsq(A)
610    G,g = G2lat.A2Gmat(A)       #get recip. & real metric tensors
611    cell = np.array(G2lat.Gmat2cell(g))   #real cell
612    cellst = np.array(G2lat.Gmat2cell(G)) #recip. cell
613    scos = cosd(cellst[3:6])
614    ssin = sind(cellst[3:6])
615    scot = scos/ssin
616    rcos = cosd(cell[3:6])
617    rsin = sind(cell[3:6])
618    rcot = rcos/rsin
619    RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
620    varyList = covData['varyList']
621    covMatrix = covData['covMatrix']
622    vcov = getVCov(RMnames,varyList,covMatrix)
623    Ax = np.array(A)
624    Ax[3:] /= 2.
625    drVdA = np.array([Ax[1]*Ax[2]-Ax[5]**2,Ax[0]*Ax[2]-Ax[4]**2,Ax[0]*Ax[1]-Ax[3]**2,
626        Ax[4]*Ax[5]-Ax[2]*Ax[3],Ax[3]*Ax[5]-Ax[1]*Ax[4],Ax[3]*Ax[4]-Ax[0]*Ax[5]])
627    srcvlsq = np.inner(drVdA,np.inner(vcov,drVdA.T))
628    Vol = 1/np.sqrt(rVsq)
629    sigVol = Vol**3*np.sqrt(srcvlsq)/2.
630    R123 = Ax[0]*Ax[1]*Ax[2]
631    dsasdg = np.zeros((3,6))
632    dadg = np.zeros((6,6))
633    for i0 in range(3):         #0  1   2
634        i1 = (i0+1)%3           #1  2   0
635        i2 = (i1+1)%3           #2  0   1
636        i3 = 5-i2               #3  5   4
637        i4 = 5-i1               #4  3   5
638        i5 = 5-i0               #5  4   3
639        dsasdg[i0][i1] = 0.5*scot[i0]*scos[i0]/Ax[i1]
640        dsasdg[i0][i2] = 0.5*scot[i0]*scos[i0]/Ax[i2]
641        dsasdg[i0][i5] = -scot[i0]/np.sqrt(Ax[i1]*Ax[i2])
642        denmsq = Ax[i0]*(R123-Ax[i1]*Ax[i4]**2-Ax[i2]*Ax[i3]**2+(Ax[i4]*Ax[i3])**2)
643        denom = np.sqrt(denmsq)
644        dadg[i5][i0] = -Ax[i5]/denom-rcos[i0]/denmsq*(R123-0.5*Ax[i1]*Ax[i4]**2-0.5*Ax[i2]*Ax[i3]**2)
645        dadg[i5][i1] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i2]-Ax[i0]*Ax[i4]**2)
646        dadg[i5][i2] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i1]-Ax[i0]*Ax[i3]**2)
647        dadg[i5][i3] = Ax[i4]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i2]*Ax[i3]-Ax[i3]*Ax[i4]**2)
648        dadg[i5][i4] = Ax[i3]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i1]*Ax[i4]-Ax[i3]**2*Ax[i4])
649        dadg[i5][i5] = -Ax[i0]/denom
650    for i0 in range(3):
651        i1 = (i0+1)%3
652        i2 = (i1+1)%3
653        i3 = 5-i2
654        for ij in range(6):
655            dadg[i0][ij] = cell[i0]*(rcot[i2]*dadg[i3][ij]/rsin[i2]-dsasdg[i1][ij]/ssin[i1])
656            if ij == i0:
657                dadg[i0][ij] = dadg[i0][ij]-0.5*cell[i0]/Ax[i0]
658            dadg[i3][ij] = -dadg[i3][ij]*rsin[2-i0]*dpr
659    sigMat = np.inner(dadg,np.inner(vcov,dadg.T))
660    var = np.diag(sigMat)
661    CS = np.where(var>0.,np.sqrt(var),0.)
662    cellSig = [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol]  #exchange sig(alp) & sig(gam) to get in right order
663    return cellSig           
664   
665def SetPhaseData(parmDict,sigDict,Phases,covData):
666   
667    def cellFill(pfx,SGData,parmDict,sigDict): 
668        if SGData['SGLaue'] in ['-1',]:
669            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
670                parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
671            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
672                sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
673        elif SGData['SGLaue'] in ['2/m',]:
674            if SGData['SGUniq'] == 'a':
675                A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
676                    parmDict[pfx+'A3'],0,0]
677                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
678                    sigDict[pfx+'A3'],0,0]
679            elif SGData['SGUniq'] == 'b':
680                A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
681                    0,parmDict[pfx+'A4'],0]
682                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
683                    0,sigDict[pfx+'A4'],0]
684            else:
685                A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
686                    0,0,parmDict[pfx+'A5']]
687                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
688                    0,0,sigDict[pfx+'A5']]
689        elif SGData['SGLaue'] in ['mmm',]:
690            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
691            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
692        elif SGData['SGLaue'] in ['4/m','4/mmm']:
693            A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
694            sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
695        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
696            A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
697                parmDict[pfx+'A0'],0,0]
698            sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
699        elif SGData['SGLaue'] in ['3R', '3mR']:
700            A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
701                parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
702            sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
703        elif SGData['SGLaue'] in ['m3m','m3']:
704            A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
705            sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
706        return A,sigA
707       
708    def PrintAtomsAndSig(General,Atoms,atomsSig):
709        print '\n Atoms:'
710        line = '   name      x         y         z      frac   Uiso     U11     U22     U33     U12     U13     U23'
711        if General['Type'] == 'magnetic':
712            line += '   Mx     My     Mz'
713        elif General['Type'] == 'macromolecular':
714            line = ' res no  residue  chain '+line
715        print line
716        if General['Type'] == 'nuclear':
717            print 135*'-'
718            fmt = {0:'%7s',1:'%7s',3:'%10.5f',4:'%10.5f',5:'%10.5f',6:'%8.3f',10:'%8.5f',
719                11:'%8.5f',12:'%8.5f',13:'%8.5f',14:'%8.5f',15:'%8.5f',16:'%8.5f'}
720            noFXsig = {3:[10*' ','%10s'],4:[10*' ','%10s'],5:[10*' ','%10s'],6:[8*' ','%8s']}
721            for i,at in enumerate(Atoms):
722                name = fmt[0]%(at[0])+fmt[1]%(at[1])+':'
723                valstr = ' values:'
724                sigstr = ' sig   :'
725                for ind in [3,4,5,6]:
726                    sigind = str(i)+':'+str(ind)
727                    valstr += fmt[ind]%(at[ind])                   
728                    if sigind in atomsSig:
729                        sigstr += fmt[ind]%(atomsSig[sigind])
730                    else:
731                        sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
732                if at[9] == 'I':
733                    valstr += fmt[10]%(at[10])
734                    if str(i)+':10' in atomsSig:
735                        sigstr += fmt[10]%(atomsSig[str(i)+':10'])
736                    else:
737                        sigstr += 8*' '
738                else:
739                    valstr += 8*' '
740                    sigstr += 8*' '
741                    for ind in [11,12,13,14,15,16]:
742                        sigind = str(i)+':'+str(ind)
743                        valstr += fmt[ind]%(at[ind])
744                        if sigind in atomsSig:                       
745                            sigstr += fmt[ind]%(atomsSig[sigind])
746                        else:
747                            sigstr += 8*' '
748                print name
749                print valstr
750                print sigstr
751               
752    def PrintSHtextureAndSig(textureData,SHtextureSig):
753        print '\n Spherical harmonics texture: Order:' + str(textureData['Order'])
754        names = ['omega','chi','phi']
755        namstr = '  names :'
756        ptstr =  '  values:'
757        sigstr = '  esds  :'
758        for name in names:
759            namstr += '%12s'%(name)
760            ptstr += '%12.3f'%(textureData['Sample '+name][1])
761            if 'Sample '+name in SHtextureSig:
762                sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
763            else:
764                sigstr += 12*' '
765        print namstr
766        print ptstr
767        print sigstr
768        print '\n Texture coefficients:'
769        namstr = '  names :'
770        ptstr =  '  values:'
771        sigstr = '  esds  :'
772        SHcoeff = textureData['SH Coeff'][1]
773        for name in SHcoeff:
774            namstr += '%12s'%(name)
775            ptstr += '%12.3f'%(SHcoeff[name])
776            if name in SHtextureSig:
777                sigstr += '%12.3f'%(SHtextureSig[name])
778            else:
779                sigstr += 12*' '
780        print namstr
781        print ptstr
782        print sigstr
783       
784           
785    print '\n Phases:'
786    for phase in Phases:
787        print ' Result for phase: ',phase
788        Phase = Phases[phase]
789        General = Phase['General']
790        SGData = General['SGData']
791        Atoms = Phase['Atoms']
792        cell = General['Cell']
793        textureData = General['SH Texture']
794        pId = Phase['pId']
795        pfx = str(pId)+'::'
796        if cell[0]:
797            A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
798            cellSig = getCellEsd(pfx,SGData,A,covData)  #includes sigVol
799            print ' Reciprocal metric tensor: '
800            ptfmt = "%15.9f"
801            names = ['A11','A22','A33','A12','A13','A23']
802            namstr = '  names :'
803            ptstr =  '  values:'
804            sigstr = '  esds  :'
805            for name,a,siga in zip(names,A,sigA):
806                namstr += '%15s'%(name)
807                ptstr += ptfmt%(a)
808                if siga:
809                    sigstr += ptfmt%(siga)
810                else:
811                    sigstr += 15*' '
812            print namstr
813            print ptstr
814            print sigstr
815            cell[1:7] = G2lat.A2cell(A)
816            cell[7] = G2lat.calc_V(A)
817            print ' New unit cell:'
818            ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
819            names = ['a','b','c','alpha','beta','gamma','Volume']
820            namstr = '  names :'
821            ptstr =  '  values:'
822            sigstr = '  esds  :'
823            for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
824                namstr += '%12s'%(name)
825                ptstr += fmt%(a)
826                if siga:
827                    sigstr += fmt%(siga)
828                else:
829                    sigstr += 12*' '
830            print namstr
831            print ptstr
832            print sigstr
833           
834        if 'Pawley' in Phase['General']['Type']:
835            pawleyRef = Phase['Pawley ref']
836            for i,refl in enumerate(pawleyRef):
837                key = pfx+'PWLref:'+str(i)
838                refl[6] = abs(parmDict[key])        #suppress negative Fsq
839                if key in sigDict:
840                    refl[7] = sigDict[key]
841                else:
842                    refl[7] = 0
843        else:
844            atomsSig = {}
845            if General['Type'] == 'nuclear':
846                for i,at in enumerate(Atoms):
847                    names = {3:pfx+'Ax:'+str(i),4:pfx+'Ay:'+str(i),5:pfx+'Az:'+str(i),6:pfx+'Afrac:'+str(i),
848                        10:pfx+'AUiso:'+str(i),11:pfx+'AU11:'+str(i),12:pfx+'AU22:'+str(i),13:pfx+'AU33:'+str(i),
849                        14:pfx+'AU12:'+str(i),15:pfx+'AU13:'+str(i),16:pfx+'AU23:'+str(i)}
850                    for ind in [3,4,5,6]:
851                        at[ind] = parmDict[names[ind]]
852                        if ind in [3,4,5]:
853                            name = names[ind].replace('A','dA')
854                        else:
855                            name = names[ind]
856                        if name in sigDict:
857                            atomsSig[str(i)+':'+str(ind)] = sigDict[name]
858                    if at[9] == 'I':
859                        at[10] = parmDict[names[10]]
860                        if names[10] in sigDict:
861                            atomsSig[str(i)+':10'] = sigDict[names[10]]
862                    else:
863                        for ind in [11,12,13,14,15,16]:
864                            at[ind] = parmDict[names[ind]]
865                            if names[ind] in sigDict:
866                                atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
867            PrintAtomsAndSig(General,Atoms,atomsSig)
868           
869        if textureData['Order']:
870            SHtextureSig = {}
871            for name in ['omega','chi','phi']:
872                aname = pfx+'SH '+name
873                textureData['Sample '+name][1] = parmDict[aname]
874                if aname in sigDict:
875                    SHtextureSig['Sample '+name] = sigDict[aname]
876            for name in textureData['SH Coeff'][1]:
877                aname = pfx+name
878                textureData['SH Coeff'][1][name] = parmDict[aname]
879                if aname in sigDict:
880                    SHtextureSig[name] = sigDict[aname]
881            PrintSHtextureAndSig(textureData,SHtextureSig)
882
883def GetHistogramPhaseData(Phases,Histograms,Print=True):
884   
885    def PrintSize(hapData):
886        if hapData[0] in ['isotropic','uniaxial']:
887            line = '\n Size model    : %9s'%(hapData[0])
888            line += ' equatorial:'+'%12.3f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
889            if hapData[0] == 'uniaxial':
890                line += ' axial:'+'%12.3f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
891            print line
892        else:
893            print '\n Size model    : %s'%(hapData[0])
894            Snames = ['S11','S22','S33','S12','S13','S23']
895            ptlbls = ' names :'
896            ptstr =  ' values:'
897            varstr = ' refine:'
898            for i,name in enumerate(Snames):
899                ptlbls += '%12s' % (name)
900                ptstr += '%12.6f' % (hapData[4][i])
901                varstr += '%12s' % (str(hapData[5][i]))
902            print ptlbls
903            print ptstr
904            print varstr
905       
906    def PrintMuStrain(hapData,SGData):
907        if hapData[0] in ['isotropic','uniaxial']:
908            line = '\n Mustrain model: %9s'%(hapData[0])
909            line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
910            if hapData[0] == 'uniaxial':
911                line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
912            print line
913        else:
914            print '\n Mustrain model: %s'%(hapData[0])
915            Snames = G2spc.MustrainNames(SGData)
916            ptlbls = ' names :'
917            ptstr =  ' values:'
918            varstr = ' refine:'
919            for i,name in enumerate(Snames):
920                ptlbls += '%12s' % (name)
921                ptstr += '%12.6f' % (hapData[4][i])
922                varstr += '%12s' % (str(hapData[5][i]))
923            print ptlbls
924            print ptstr
925            print varstr
926
927    def PrintHStrain(hapData,SGData):
928        print '\n Hydrostatic/elastic strain: '
929        Hsnames = G2spc.HStrainNames(SGData)
930        ptlbls = ' names :'
931        ptstr =  ' values:'
932        varstr = ' refine:'
933        for i,name in enumerate(Hsnames):
934            ptlbls += '%12s' % (name)
935            ptstr += '%12.6f' % (hapData[0][i])
936            varstr += '%12s' % (str(hapData[1][i]))
937        print ptlbls
938        print ptstr
939        print varstr
940
941    def PrintSHPO(hapData):
942        print '\n Spherical harmonics preferred orientation: Order:' + \
943            str(hapData[4])+' Refine? '+str(hapData[2])
944        ptlbls = ' names :'
945        ptstr =  ' values:'
946        for item in hapData[5]:
947            ptlbls += '%12s'%(item)
948            ptstr += '%12.3f'%(hapData[5][item]) 
949        print ptlbls
950        print ptstr
951   
952    hapDict = {}
953    hapVary = []
954    controlDict = {}
955    poType = {}
956    poAxes = {}
957    spAxes = {}
958    spType = {}
959   
960    for phase in Phases:
961        HistoPhase = Phases[phase]['Histograms']
962        SGData = Phases[phase]['General']['SGData']
963        cell = Phases[phase]['General']['Cell'][1:7]
964        A = G2lat.cell2A(cell)
965        pId = Phases[phase]['pId']
966        for histogram in HistoPhase:
967            try:
968                Histogram = Histograms[histogram]
969            except KeyError:                       
970                #skip if histogram not included e.g. in a sequential refinement
971                continue
972            hapData = HistoPhase[histogram]
973            hId = Histogram['hId']
974            limits = Histogram['Limits'][1]
975            inst = Histogram['Instrument Parameters']
976            inst = dict(zip(inst[3],inst[1]))
977            Zero = inst['Zero']
978            if 'C' in inst['Type']:
979                try:
980                    wave = inst['Lam']
981                except KeyError:
982                    wave = inst['Lam1']
983                dmin = wave/(2.0*sind(limits[1]/2.0))
984            pfx = str(pId)+':'+str(hId)+':'
985            for item in ['Scale','Extinction']:
986                hapDict[pfx+item] = hapData[item][0]
987                if hapData[item][1]:
988                    hapVary.append(pfx+item)
989            names = G2spc.HStrainNames(SGData)
990            for i,name in enumerate(names):
991                hapDict[pfx+name] = hapData['HStrain'][0][i]
992                if hapData['HStrain'][1][i]:
993                    hapVary.append(pfx+name)
994            controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
995            if hapData['Pref.Ori.'][0] == 'MD':
996                hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
997                controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
998                if hapData['Pref.Ori.'][2]:
999                    hapVary.append(pfx+'MD')
1000            else:                           #'SH' spherical harmonics
1001                controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
1002                controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
1003                for item in hapData['Pref.Ori.'][5]:
1004                    hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
1005                    if hapData['Pref.Ori.'][2]:
1006                        hapVary.append(pfx+item)
1007            for item in ['Mustrain','Size']:
1008                controlDict[pfx+item+'Type'] = hapData[item][0]
1009                if hapData[item][0] in ['isotropic','uniaxial']:
1010                    hapDict[pfx+item+':0'] = hapData[item][1][0]
1011                    if hapData[item][2][0]:
1012                        hapVary.append(pfx+item+':0')
1013                    if hapData[item][0] == 'uniaxial':
1014                        controlDict[pfx+item+'Axis'] = hapData[item][3]
1015                        hapDict[pfx+item+':1'] = hapData[item][1][1]
1016                        if hapData[item][2][1]:
1017                            hapVary.append(pfx+item+':1')
1018                else:       #generalized for mustrain or ellipsoidal for size
1019                    if item == 'Mustrain':
1020                        names = G2spc.MustrainNames(SGData)
1021                        pwrs = []
1022                        for name in names:
1023                            h,k,l = name[1:]
1024                            pwrs.append([int(h),int(k),int(l)])
1025                        controlDict[pfx+'MuPwrs'] = pwrs
1026                    for i in range(len(hapData[item][4])):
1027                        sfx = ':'+str(i)
1028                        hapDict[pfx+item+sfx] = hapData[item][4][i]
1029                        if hapData[item][5][i]:
1030                            hapVary.append(pfx+item+sfx)
1031                           
1032            if Print: 
1033                print '\n Phase: ',phase,' in histogram: ',histogram
1034                print 135*'-'
1035                print ' Phase fraction  : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1036                print ' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
1037                if hapData['Pref.Ori.'][0] == 'MD':
1038                    Ax = hapData['Pref.Ori.'][3]
1039                    print ' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
1040                        ' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
1041                else: #'SH' for spherical harmonics
1042                    PrintSHPO(hapData['Pref.Ori.'])
1043                PrintSize(hapData['Size'])
1044                PrintMuStrain(hapData['Mustrain'],SGData)
1045                PrintHStrain(hapData['HStrain'],SGData)
1046            HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
1047            refList = []
1048            for h,k,l,d in HKLd:
1049                ext,mul,Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
1050                if ext:
1051                    continue
1052                if 'C' in inst['Type']:
1053                    pos = 2.0*asind(wave/(2.0*d))
1054                    if limits[0] < pos < limits[1]:
1055                        refList.append([h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,0.0,Uniq,phi,0.0])
1056                else:
1057                    raise ValueError 
1058            Histogram['Reflection Lists'][phase] = refList
1059    return hapVary,hapDict,controlDict
1060   
1061def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,Print=True):
1062   
1063    def PrintSizeAndSig(hapData,sizeSig):
1064        line = '\n Size model:     %9s'%(hapData[0])
1065        if hapData[0] in ['isotropic','uniaxial']:
1066            line += ' equatorial:%12.3f'%(hapData[1][0])
1067            if sizeSig[0][0]:
1068                line += ', sig: %8.3f'%(sizeSig[0][0])
1069            if hapData[0] == 'uniaxial':
1070                line += ' axial:%12.3f'%(hapData[1][1])
1071                if sizeSig[0][1]:
1072                    line += ', sig: %8.3f'%(sizeSig[0][1])
1073            print line
1074        else:
1075            print line
1076            Snames = ['S11','S22','S33','S12','S13','S23']
1077            ptlbls = ' name  :'
1078            ptstr =  ' value :'
1079            sigstr = ' sig   :'
1080            for i,name in enumerate(Snames):
1081                ptlbls += '%12s' % (name)
1082                ptstr += '%12.6f' % (hapData[4][i])
1083                if sizeSig[1][i]:
1084                    sigstr += '%12.6f' % (sizeSig[1][i])
1085                else:
1086                    sigstr += 12*' '
1087            print ptlbls
1088            print ptstr
1089            print sigstr
1090       
1091    def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
1092        line = '\n Mustrain model: %9s'%(hapData[0])
1093        if hapData[0] in ['isotropic','uniaxial']:
1094            line += ' equatorial:%12.1f'%(hapData[1][0])
1095            if mustrainSig[0][0]:
1096                line += ', sig: %8.1f'%(mustrainSig[0][0])
1097            if hapData[0] == 'uniaxial':
1098                line += ' axial:%12.1f'%(hapData[1][1])
1099                if mustrainSig[0][1]:
1100                     line += ', sig: %8.1f'%(mustrainSig[0][1])
1101            print line
1102        else:
1103            print line
1104            Snames = G2spc.MustrainNames(SGData)
1105            ptlbls = ' name  :'
1106            ptstr =  ' value :'
1107            sigstr = ' sig   :'
1108            for i,name in enumerate(Snames):
1109                ptlbls += '%12s' % (name)
1110                ptstr += '%12.6f' % (hapData[4][i])
1111                if mustrainSig[1][i]:
1112                    sigstr += '%12.6f' % (mustrainSig[1][i])
1113                else:
1114                    sigstr += 12*' '
1115            print ptlbls
1116            print ptstr
1117            print sigstr
1118           
1119    def PrintHStrainAndSig(hapData,strainSig,SGData):
1120        print '\n Hydrostatic/elastic strain: '
1121        Hsnames = G2spc.HStrainNames(SGData)
1122        ptlbls = ' name  :'
1123        ptstr =  ' value :'
1124        sigstr = ' sig   :'
1125        for i,name in enumerate(Hsnames):
1126            ptlbls += '%12s' % (name)
1127            ptstr += '%12.6g' % (hapData[0][i])
1128            if name in strainSig:
1129                sigstr += '%12.6g' % (strainSig[name])
1130            else:
1131                sigstr += 12*' '
1132        print ptlbls
1133        print ptstr
1134        print sigstr
1135       
1136    def PrintSHPOAndSig(hapData,POsig):
1137        print '\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
1138        ptlbls = ' names :'
1139        ptstr =  ' values:'
1140        sigstr = ' sig   :'
1141        for item in hapData[5]:
1142            ptlbls += '%12s'%(item)
1143            ptstr += '%12.3f'%(hapData[5][item])
1144            if item in POsig:
1145                sigstr += '%12.3f'%(POsig[item])
1146            else:
1147                sigstr += 12*' ' 
1148        print ptlbls
1149        print ptstr
1150        print sigstr
1151   
1152    for phase in Phases:
1153        HistoPhase = Phases[phase]['Histograms']
1154        SGData = Phases[phase]['General']['SGData']
1155        pId = Phases[phase]['pId']
1156        for histogram in HistoPhase:
1157            try:
1158                Histogram = Histograms[histogram]
1159            except KeyError:                       
1160                #skip if histogram not included e.g. in a sequential refinement
1161                continue
1162            print '\n Phase: ',phase,' in histogram: ',histogram
1163            print 130*'-'
1164            hapData = HistoPhase[histogram]
1165            hId = Histogram['hId']
1166            pfx = str(pId)+':'+str(hId)+':'
1167            print ' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections'   \
1168                %(Histogram[pfx+'Rf'],Histogram[pfx+'Rf^2'],Histogram[pfx+'Nref'])
1169           
1170            PhFrExtPOSig = {}
1171            for item in ['Scale','Extinction']:
1172                hapData[item][0] = parmDict[pfx+item]
1173                if pfx+item in sigDict:
1174                    PhFrExtPOSig[item] = sigDict[pfx+item]
1175            if hapData['Pref.Ori.'][0] == 'MD':
1176                hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
1177                if pfx+'MD' in sigDict:
1178                    PhFrExtPOSig['MD'] = sigDict[pfx+'MD']
1179            else:                           #'SH' spherical harmonics
1180                for item in hapData['Pref.Ori.'][5]:
1181                    hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
1182                    if pfx+item in sigDict:
1183                        PhFrExtPOSig[item] = sigDict[pfx+item]
1184            if Print:
1185                if 'Scale' in PhFrExtPOSig:
1186                    print ' Phase fraction  : %10.4f, sig %10.4f'%(hapData['Scale'][0],PhFrExtPOSig['Scale'])
1187                if 'Extinction' in PhFrExtPOSig:
1188                    print ' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig['Extinction'])
1189                if hapData['Pref.Ori.'][0] == 'MD':
1190                    if 'MD' in PhFrExtPOSig:
1191                        print ' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig['MD'])
1192                else:
1193                    PrintSHPOAndSig(hapData['Pref.Ori.'],PhFrExtPOSig)
1194            SizeMuStrSig = {'Mustrain':[[0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
1195                'Size':[[0,0],[0 for i in range(len(hapData['Size'][4]))]],
1196                'HStrain':{}}                 
1197            for item in ['Mustrain','Size']:
1198                if hapData[item][0] in ['isotropic','uniaxial']:                   
1199                    hapData[item][1][0] = parmDict[pfx+item+':0']
1200                    if item == 'Size':
1201                        hapData[item][1][0] = min(10.,max(0.01,hapData[item][1][0]))
1202                    if pfx+item+':0' in sigDict: 
1203                        SizeMuStrSig[item][0][0] = sigDict[pfx+item+':0']
1204                    if hapData[item][0] == 'uniaxial':
1205                        hapData[item][1][1] = parmDict[pfx+item+':1']
1206                        if item == 'Size':
1207                            hapData[item][1][1] = min(10.,max(0.01,hapData[item][1][1]))                       
1208                        if pfx+item+':1' in sigDict:
1209                            SizeMuStrSig[item][0][1] = sigDict[pfx+item+':1']
1210                else:       #generalized for mustrain or ellipsoidal for size
1211                    for i in range(len(hapData[item][4])):
1212                        sfx = ':'+str(i)
1213                        hapData[item][4][i] = parmDict[pfx+item+sfx]
1214                        if pfx+item+sfx in sigDict:
1215                            SizeMuStrSig[item][1][i] = sigDict[pfx+item+sfx]
1216            names = G2spc.HStrainNames(SGData)
1217            for i,name in enumerate(names):
1218                hapData['HStrain'][0][i] = parmDict[pfx+name]
1219                if pfx+name in sigDict:
1220                    SizeMuStrSig['HStrain'][name] = sigDict[pfx+name]
1221            if Print:
1222                PrintSizeAndSig(hapData['Size'],SizeMuStrSig['Size'])
1223                PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig['Mustrain'],SGData)
1224                PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig['HStrain'],SGData)
1225   
1226def GetHistogramData(Histograms,Print=True):
1227   
1228    def GetBackgroundParms(hId,Background):
1229        bakType,bakFlag = Background[:2]
1230        backVals = Background[3:]
1231        backNames = [':'+str(hId)+':Back:'+str(i) for i in range(len(backVals))]
1232        if bakFlag:                                 #returns backNames as varyList = backNames
1233            return bakType,dict(zip(backNames,backVals)),backNames
1234        else:                                       #no background varied; varyList = []
1235            return bakType,dict(zip(backNames,backVals)),[]
1236       
1237    def GetInstParms(hId,Inst):
1238        insVals,insFlags,insNames = Inst[1:4]
1239        dataType = insVals[0]
1240        instDict = {}
1241        insVary = []
1242        pfx = ':'+str(hId)+':'
1243        for i,flag in enumerate(insFlags):
1244            insName = pfx+insNames[i]
1245            instDict[insName] = insVals[i]
1246            if flag:
1247                insVary.append(insName)
1248        instDict[pfx+'X'] = max(instDict[pfx+'X'],0.01)
1249        instDict[pfx+'Y'] = max(instDict[pfx+'Y'],0.01)
1250        instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
1251        return dataType,instDict,insVary
1252       
1253    def GetSampleParms(hId,Sample):
1254        sampVary = []
1255        hfx = ':'+str(hId)+':'       
1256        sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
1257            hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
1258        Type = Sample['Type']
1259        if 'Bragg' in Type:             #Bragg-Brentano
1260            for item in ['Scale','Shift','Transparency']:       #surface roughness?, diffuse scattering?
1261                sampDict[hfx+item] = Sample[item][0]
1262                if Sample[item][1]:
1263                    sampVary.append(hfx+item)
1264        elif 'Debye' in Type:        #Debye-Scherrer
1265            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
1266                sampDict[hfx+item] = Sample[item][0]
1267                if Sample[item][1]:
1268                    sampVary.append(hfx+item)
1269        return Type,sampDict,sampVary
1270       
1271    def PrintBackground(Background):
1272        print '\n Background function: ',Background[0],' Refine?',bool(Background[1])
1273        line = ' Coefficients: '
1274        for i,back in enumerate(Background[3:]):
1275            line += '%10.3f'%(back)
1276            if i and not i%10:
1277                line += '\n'+15*' '
1278        print line
1279       
1280    def PrintInstParms(Inst):
1281        print '\n Instrument Parameters:'
1282        ptlbls = ' name  :'
1283        ptstr =  ' value :'
1284        varstr = ' refine:'
1285        instNames = Inst[3][1:]
1286        for i,name in enumerate(instNames):
1287            ptlbls += '%12s' % (name)
1288            ptstr += '%12.6f' % (Inst[1][i+1])
1289            if name in ['Lam1','Lam2','Azimuth']:
1290                varstr += 12*' '
1291            else:
1292                varstr += '%12s' % (str(bool(Inst[2][i+1])))
1293        print ptlbls
1294        print ptstr
1295        print varstr
1296       
1297    def PrintSampleParms(Sample):
1298        print '\n Sample Parameters:'
1299        print ' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \
1300            (Sample['Omega'],Sample['Chi'],Sample['Phi'])
1301        ptlbls = ' name  :'
1302        ptstr =  ' value :'
1303        varstr = ' refine:'
1304        if 'Bragg' in Sample['Type']:
1305            for item in ['Scale','Shift','Transparency']:
1306                ptlbls += '%14s'%(item)
1307                ptstr += '%14.4f'%(Sample[item][0])
1308                varstr += '%14s'%(str(bool(Sample[item][1])))
1309           
1310        elif 'Debye' in Type:        #Debye-Scherrer
1311            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
1312                ptlbls += '%14s'%(item)
1313                ptstr += '%14.4f'%(Sample[item][0])
1314                varstr += '%14s'%(str(bool(Sample[item][1])))
1315
1316        print ptlbls
1317        print ptstr
1318        print varstr
1319       
1320
1321    histDict = {}
1322    histVary = []
1323    controlDict = {}
1324    for histogram in Histograms:
1325        Histogram = Histograms[histogram]
1326        hId = Histogram['hId']
1327        pfx = ':'+str(hId)+':'
1328        controlDict[pfx+'Limits'] = Histogram['Limits'][1]
1329       
1330        Background = Histogram['Background'][0]
1331        Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
1332        controlDict[pfx+'bakType'] = Type
1333        histDict.update(bakDict)
1334        histVary += bakVary
1335       
1336        Inst = Histogram['Instrument Parameters']
1337        Type,instDict,insVary = GetInstParms(hId,Inst)
1338        controlDict[pfx+'histType'] = Type
1339        if pfx+'Lam1' in instDict:
1340            controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
1341        else:
1342            controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']           
1343        histDict.update(instDict)
1344        histVary += insVary
1345       
1346        Sample = Histogram['Sample Parameters']
1347        Type,sampDict,sampVary = GetSampleParms(hId,Sample)
1348        controlDict[pfx+'instType'] = Type
1349        histDict.update(sampDict)
1350        histVary += sampVary
1351
1352        if Print: 
1353            print '\n Histogram: ',histogram,' histogram Id: ',hId
1354            print 135*'-'
1355            Units = {'C':' deg','T':' msec'}
1356            units = Units[controlDict[pfx+'histType'][2]]
1357            Limits = controlDict[pfx+'Limits']
1358            print ' Instrument type: ',Sample['Type']
1359            print ' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)     
1360            PrintSampleParms(Sample)
1361            PrintInstParms(Inst)
1362            PrintBackground(Background)
1363       
1364    return histVary,histDict,controlDict
1365   
1366def SetHistogramData(parmDict,sigDict,Histograms,Print=True):
1367   
1368    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
1369        lenBack = len(Background[3:])
1370        backSig = [0 for i in range(lenBack)]
1371        for i in range(lenBack):
1372            Background[3+i] = parmDict[pfx+'Back:'+str(i)]
1373            if pfx+'Back:'+str(i) in sigDict:
1374                backSig[i] = sigDict[pfx+'Back:'+str(i)]
1375        return backSig
1376       
1377    def SetInstParms(pfx,Inst,parmDict,sigDict):
1378        insVals,insFlags,insNames = Inst[1:4]
1379        instSig = [0 for i in range(len(insVals))]
1380        for i,flag in enumerate(insFlags):
1381            insName = pfx+insNames[i]
1382            insVals[i] = parmDict[insName]
1383            if insName in sigDict:
1384                instSig[i] = sigDict[insName]
1385        return instSig
1386       
1387    def SetSampleParms(pfx,Sample,parmDict,sigDict):
1388        if 'Bragg' in Sample['Type']:             #Bragg-Brentano
1389            sampSig = [0 for i in range(3)]
1390            for i,item in enumerate(['Scale','Shift','Transparency']):       #surface roughness?, diffuse scattering?
1391                Sample[item][0] = parmDict[pfx+item]
1392                if pfx+item in sigDict:
1393                    sampSig[i] = sigDict[pfx+item]
1394        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
1395            sampSig = [0 for i in range(4)]
1396            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
1397                Sample[item][0] = parmDict[pfx+item]
1398                if pfx+item in sigDict:
1399                    sampSig[i] = sigDict[pfx+item]
1400        return sampSig
1401       
1402    def PrintBackgroundSig(Background,backSig):
1403        print '\n Background function: ',Background[0]
1404        valstr = ' value : '
1405        sigstr = ' sig   : '
1406        for i,back in enumerate(Background[3:]):
1407            valstr += '%10.4f'%(back)
1408            if Background[1]:
1409                sigstr += '%10.4f'%(backSig[i])
1410            else:
1411                sigstr += 10*' '
1412        print valstr
1413        print sigstr
1414       
1415    def PrintInstParmsSig(Inst,instSig):
1416        print '\n Instrument Parameters:'
1417        ptlbls = ' names :'
1418        ptstr =  ' value :'
1419        sigstr = ' sig   :'
1420        instNames = Inst[3][1:]
1421        for i,name in enumerate(instNames):
1422            ptlbls += '%12s' % (name)
1423            ptstr += '%12.6f' % (Inst[1][i+1])
1424            if instSig[i+1]:
1425                sigstr += '%12.6f' % (instSig[i+1])
1426            else:
1427                sigstr += 12*' '
1428        print ptlbls
1429        print ptstr
1430        print sigstr
1431       
1432    def PrintSampleParmsSig(Sample,sampleSig):
1433        print '\n Sample Parameters:'
1434        ptlbls = ' names :'
1435        ptstr =  ' values:'
1436        sigstr = ' sig   :'
1437        if 'Bragg' in Sample['Type']:
1438            for i,item in enumerate(['Scale','Shift','Transparency']):
1439                ptlbls += '%14s'%(item)
1440                ptstr += '%14.4f'%(Sample[item][0])
1441                if sampleSig[i]:
1442                    sigstr += '%14.4f'%(sampleSig[i])
1443                else:
1444                    sigstr += 14*' '
1445           
1446        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
1447            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
1448                ptlbls += '%14s'%(item)
1449                ptstr += '%14.4f'%(Sample[item][0])
1450                if sampleSig[i]:
1451                    sigstr += '%14.4f'%(sampleSig[i])
1452                else:
1453                    sigstr += 14*' '
1454
1455        print ptlbls
1456        print ptstr
1457        print sigstr
1458       
1459    for histogram in Histograms:
1460        if 'PWDR' in histogram:
1461            Histogram = Histograms[histogram]
1462            hId = Histogram['hId']
1463            pfx = ':'+str(hId)+':'
1464            Background = Histogram['Background'][0]
1465            backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
1466           
1467            Inst = Histogram['Instrument Parameters']
1468            instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
1469       
1470            Sample = Histogram['Sample Parameters']
1471            sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
1472
1473            print '\n Histogram: ',histogram,' histogram Id: ',hId
1474            print 135*'-'
1475            print ' Final refinement wRp = %.2f%% on %d observations in this histogram'%(Histogram['wRp'],Histogram['Nobs'])
1476            if Print:
1477                print ' Instrument type: ',Sample['Type']
1478                PrintSampleParmsSig(Sample,sampSig)
1479                PrintInstParmsSig(Inst,instSig)
1480                PrintBackgroundSig(Background,backSig)
1481
1482def GetAtomFXU(pfx,FFtables,calcControls,parmDict):
1483    Natoms = calcControls['Natoms'][pfx]
1484    Tdata = Natoms*[' ',]
1485    Mdata = np.zeros(Natoms)
1486    IAdata = Natoms*[' ',]
1487    Fdata = np.zeros(Natoms)
1488    FFdata = []
1489    Xdata = np.zeros((3,Natoms))
1490    dXdata = np.zeros((3,Natoms))
1491    Uisodata = np.zeros(Natoms)
1492    Uijdata = np.zeros((6,Natoms))
1493    keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata,
1494        'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2],
1495        'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata,
1496        'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2],
1497        'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5]}
1498    for iatm in range(Natoms):
1499        for key in keys:
1500            parm = pfx+key+str(iatm)
1501            if parm in parmDict:
1502                keys[key][iatm] = parmDict[parm]
1503        FFdata.append(FFtables[Tdata[iatm]])
1504    return FFdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata
1505   
1506def StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict):
1507    ''' Compute structure factors for all h,k,l for phase
1508    input:
1509        refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv]
1510        G:      reciprocal metric tensor
1511        pfx:    phase id string
1512        SGData: space group info. dictionary output from SpcGroup
1513        calcControls:
1514        ParmDict:
1515    puts result F^2 in each ref[8] in refList
1516    '''       
1517    twopi = 2.0*np.pi
1518    twopisq = 2.0*np.pi**2
1519    ast = np.sqrt(np.diag(G))
1520    Mast = twopisq*np.multiply.outer(ast,ast)
1521    FFtables = calcControls['FFtables']
1522    FFdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,FFtables,calcControls,parmDict)
1523    FP = np.array([El[hfx+'FP'] for El in FFdata])
1524    FPP = np.array([El[hfx+'FPP'] for El in FFdata])
1525    maxPos = len(SGData['SGOps'])
1526    Uij = np.array(G2lat.U6toUij(Uijdata))
1527    bij = Mast*Uij.T
1528    for refl in refList:
1529        fbs = np.array([0,0])
1530        H = refl[:3]
1531        SQ = 1./(2.*refl[4])**2
1532        FF = np.array([G2el.ScatFac(El,SQ)[0] for El in FFdata])
1533        SQfactor = 4.0*SQ*twopisq
1534        Uniq = refl[11]
1535        phi = refl[12]
1536        phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+phi[:,np.newaxis])
1537        sinp = np.sin(phase)
1538        cosp = np.cos(phase)
1539        occ = Mdata*Fdata/len(Uniq)
1540        biso = -SQfactor*Uisodata
1541        Tiso = np.where(biso<1.,np.exp(biso),1.0)
1542        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
1543        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1544        Tcorr = Tiso*Tuij
1545        fa = np.array([(FF+FP)*occ*cosp*Tcorr,-FPP*occ*sinp*Tcorr])
1546        fas = np.sum(np.sum(fa,axis=1),axis=1)        #real
1547        if not SGData['SGInv']:
1548            fb = np.array([(FF+FP)*occ*sinp*Tcorr,FPP*occ*cosp*Tcorr])
1549            fbs = np.sum(np.sum(fb,axis=1),axis=1)
1550        fasq = fas**2
1551        fbsq = fbs**2        #imaginary
1552        refl[9] = np.sum(fasq)+np.sum(fbsq)
1553        refl[10] = atan2d(fbs[0],fas[0])
1554    return refList
1555   
1556def StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict):
1557    twopi = 2.0*np.pi
1558    twopisq = 2.0*np.pi**2
1559    ast = np.sqrt(np.diag(G))
1560    Mast = twopisq*np.multiply.outer(ast,ast)
1561    FFtables = calcControls['FFtables']
1562    FFdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,FFtables,calcControls,parmDict)
1563    FP = np.array([El[hfx+'FP'] for El in FFdata])
1564    FPP = np.array([El[hfx+'FPP'] for El in FFdata])
1565    maxPos = len(SGData['SGOps'])       
1566    Uij = np.array(G2lat.U6toUij(Uijdata))
1567    bij = Mast*Uij.T
1568    dFdvDict = {}
1569    dFdfr = np.zeros((len(refList),len(Mdata)))
1570    dFdx = np.zeros((len(refList),len(Mdata),3))
1571    dFdui = np.zeros((len(refList),len(Mdata)))
1572    dFdua = np.zeros((len(refList),len(Mdata),6))
1573    for iref,refl in enumerate(refList):
1574        H = np.array(refl[:3])
1575        SQ = 1./(2.*refl[4])**2          # or (sin(theta)/lambda)**2
1576        FF = np.array([G2el.ScatFac(El,SQ)[0] for El in FFdata])
1577        SQfactor = 8.0*SQ*np.pi**2
1578        Uniq = refl[11]
1579        phi = refl[12]
1580        phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+phi[np.newaxis,:])
1581        sinp = np.sin(phase)
1582        cosp = np.cos(phase)
1583        occ = Mdata*Fdata/len(Uniq)
1584        biso = -SQfactor*Uisodata
1585        Tiso = np.where(biso<1.,np.exp(biso),1.0)
1586#        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
1587        HbH = -np.inner(H,np.inner(bij,H))
1588        Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
1589        Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
1590        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1591        Tcorr = Tiso*Tuij
1592        fot = (FF+FP)*occ*Tcorr
1593        fotp = FPP*occ*Tcorr
1594        fa = np.array([fot[:,np.newaxis]*cosp,fotp[:,np.newaxis]*cosp])       #non positions
1595        fb = np.array([fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])
1596       
1597        fas = np.sum(np.sum(fa,axis=1),axis=1)
1598        fbs = np.sum(np.sum(fb,axis=1),axis=1)
1599        fax = np.array([-fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])   #positions
1600        fbx = np.array([fot[:,np.newaxis]*cosp,-fot[:,np.newaxis]*cosp])
1601        #sum below is over Uniq
1602        dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)
1603        dfadx = np.sum(twopi*Uniq*fax[:,:,:,np.newaxis],axis=2)
1604        dfadui = np.sum(-SQfactor*fa,axis=2)
1605        dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
1606        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for     
1607        dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)
1608        dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])
1609        dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])
1610        dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])
1611        if not SGData['SGInv']:
1612            dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2)
1613            dfbdx = np.sum(twopi*Uniq*fbx[:,:,:,np.newaxis],axis=2)         
1614            dfbdui = np.sum(-SQfactor*fb,axis=2)
1615            dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
1616            dFdfr[iref] += 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
1617            dFdx[iref] += 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
1618            dFdui[iref] += 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
1619            dFdua[iref] += 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
1620        #loop over atoms - each dict entry is list of derivatives for all the reflections
1621    for i in range(len(Mdata)):     
1622        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
1623        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
1624        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
1625        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
1626        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
1627        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
1628        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
1629        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
1630        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
1631        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
1632        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
1633    return dFdvDict
1634       
1635def Dict2Values(parmdict, varylist):
1636    '''Use before call to leastsq to setup list of values for the parameters
1637    in parmdict, as selected by key in varylist'''
1638    return [parmdict[key] for key in varylist] 
1639   
1640def Values2Dict(parmdict, varylist, values):
1641    ''' Use after call to leastsq to update the parameter dictionary with
1642    values corresponding to keys in varylist'''
1643    parmdict.update(zip(varylist,values))
1644   
1645def ApplyXYZshifts(parmDict,varyList):
1646    ''' takes atom x,y,z shift and applies it to corresponding atom x,y,z value
1647        input:
1648            parmDict - parameter dictionary
1649            varyList - list of variables
1650        returns:
1651            newAtomDict - dictionary of new atomic coordinate names & values;
1652                key is parameter shift name
1653    '''
1654    newAtomDict = {}
1655    for item in parmDict:
1656        if 'dA' in item:
1657            parm = ''.join(item.split('d'))
1658            parmDict[parm] += parmDict[item]
1659            newAtomDict[item] = [parm,parmDict[parm]]
1660    return newAtomDict
1661   
1662def SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict):
1663    IFCoup = 'Bragg' in calcControls[hfx+'instType']
1664    odfCor = 1.0
1665    H = refl[:3]
1666    cell = G2lat.Gmat2cell(g)
1667    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
1668    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
1669    phi,beta = G2lat.CrsAng(H,cell,SGData)
1670    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
1671    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
1672    for item in SHnames:
1673        L,M,N = eval(item.strip('C'))
1674        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1675        Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
1676        Lnorm = G2lat.Lnorm(L)
1677        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
1678    return odfCor
1679   
1680def SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict):
1681    FORPI = 12.5663706143592
1682    IFCoup = 'Bragg' in calcControls[hfx+'instType']
1683    odfCor = 1.0
1684    dFdODF = {}
1685    dFdSA = [0,0,0]
1686    H = refl[:3]
1687    cell = G2lat.Gmat2cell(g)
1688    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
1689    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
1690    phi,beta = G2lat.CrsAng(H,cell,SGData)
1691    psi,gam,dPSdA,dGMdA = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup)
1692    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
1693    for item in SHnames:
1694        L,M,N = eval(item.strip('C'))
1695        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1696        Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
1697        Lnorm = G2lat.Lnorm(L)
1698        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
1699        dFdODF[pfx+item] = Lnorm*Kcl*Ksl
1700        for i in range(3):
1701            dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i])
1702    return odfCor,dFdODF,dFdSA
1703   
1704def SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict):
1705    odfCor = 1.0
1706    H = refl[:3]
1707    cell = G2lat.Gmat2cell(g)
1708    Sangl = [0.,0.,0.]
1709    if 'Bragg' in calcControls[hfx+'instType']:
1710        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
1711        IFCoup = True
1712    else:
1713        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
1714        IFCoup = False
1715    phi,beta = G2lat.CrsAng(H,cell,SGData)
1716    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
1717    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
1718    for item in SHnames:
1719        L,N = eval(item.strip('C'))
1720        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) 
1721        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
1722    return odfCor
1723   
1724def SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict):
1725    FORPI = 12.5663706143592
1726    odfCor = 1.0
1727    dFdODF = {}
1728    H = refl[:3]
1729    cell = G2lat.Gmat2cell(g)
1730    Sangl = [0.,0.,0.]
1731    if 'Bragg' in calcControls[hfx+'instType']:
1732        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
1733        IFCoup = True
1734    else:
1735        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
1736        IFCoup = False
1737    phi,beta = G2lat.CrsAng(H,cell,SGData)
1738    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
1739    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
1740    for item in SHnames:
1741        L,N = eval(item.strip('C'))
1742        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) 
1743        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
1744        dFdODF[phfx+item] = Kcsl*Lnorm
1745    return odfCor,dFdODF
1746   
1747def GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
1748    if calcControls[phfx+'poType'] == 'MD':
1749        MD = parmDict[phfx+'MD']
1750        MDAxis = calcControls[phfx+'MDAxis']
1751        sumMD = 0
1752        for H in refl[11]:           
1753            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1754            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
1755            sumMD += A**3
1756        POcorr = sumMD/len(refl[11])
1757    else:   #spherical harmonics
1758        POcorr = SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict)
1759    return POcorr
1760   
1761def GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
1762    POderv = {}
1763    if calcControls[phfx+'poType'] == 'MD':
1764        MD = parmDict[phfx+'MD']
1765        MDAxis = calcControls[phfx+'MDAxis']
1766        sumMD = 0
1767        sumdMD = 0
1768        for H in refl[11]:           
1769            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1770            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
1771            sumMD += A**3
1772            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
1773        POcorr = sumMD/len(refl[11])
1774        POderv[phfx+'MD'] = sumdMD/len(refl[11])
1775    else:   #spherical harmonics
1776        POcorr,POderv = SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict)
1777    return POcorr,POderv
1778   
1779def GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1780    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3]               #scale*multiplicity
1781    Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
1782    Icorr *= GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
1783    if pfx+'SHorder' in parmDict:
1784        Icorr *= SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict)
1785    refl[13] = Icorr       
1786   
1787def GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1788    dIdsh = 1./parmDict[hfx+'Scale']
1789    dIdsp = 1./parmDict[phfx+'Scale']
1790    pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
1791    dIdPola /= pola
1792    POcorr,dIdPO = GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
1793    for iPO in dIdPO:
1794        dIdPO[iPO] /= POcorr
1795    dFdODF = {}
1796    dFdSA = [0,0,0]
1797    if pfx+'SHorder' in parmDict:
1798        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict)
1799        for iSH in dFdODF:
1800            dFdODF[iSH] /= odfCor
1801        for i in range(3):
1802            dFdSA[i] /= odfCor
1803    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA
1804       
1805def GetSampleGam(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse):
1806    costh = cosd(refl[5]/2.)
1807    #crystallite size
1808    if calcControls[phfx+'SizeType'] == 'isotropic':
1809        gam = 1.8*wave/(np.pi*parmDict[phfx+'Size:0']*costh)
1810    elif calcControls[phfx+'SizeType'] == 'uniaxial':
1811        H = np.array(refl[:3])
1812        P = np.array(calcControls[phfx+'SizeAxis'])
1813        cosP,sinP = G2lat.CosSinAngle(H,P,G)
1814        gam = (1.8*wave/np.pi)/(parmDict[phfx+'Size:0']*parmDict[phfx+'Size:1']*costh)
1815        gam *= np.sqrt((cosP*parmDict[phfx+'Size:1'])**2+(sinP*parmDict[phfx+'Size:0'])**2)
1816    else:           #ellipsoidal crystallites - wrong not make sense
1817        H = np.array(refl[:3])
1818        gam += 1.8*wave/(np.pi*costh*np.inner(H,np.inner(sizeEllipse,H)))
1819    #microstrain               
1820    if calcControls[phfx+'MustrainType'] == 'isotropic':
1821        gam += 0.018*parmDict[phfx+'Mustrain:0']*tand(refl[5]/2.)/np.pi
1822    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1823        H = np.array(refl[:3])
1824        P = np.array(calcControls[phfx+'MustrainAxis'])
1825        cosP,sinP = G2lat.CosSinAngle(H,P,G)
1826        Si = parmDict[phfx+'Mustrain:0']
1827        Sa = parmDict[phfx+'Mustrain:1']
1828        gam += 0.018*Si*Sa*tand(refl[5]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
1829    else:       #generalized - P.W. Stephens model
1830        pwrs = calcControls[phfx+'MuPwrs']
1831        sum = 0
1832        for i,pwr in enumerate(pwrs):
1833            sum += parmDict[phfx+'Mustrain:'+str(i)]*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
1834        gam += 0.018*refl[4]**2*tand(refl[5]/2.)*sum           
1835    return gam
1836       
1837def GetSampleGamDerv(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse):
1838    gamDict = {}
1839    costh = cosd(refl[5]/2.)
1840    tanth = tand(refl[5]/2.)
1841    #crystallite size derivatives
1842    if calcControls[phfx+'SizeType'] == 'isotropic':
1843        gamDict[phfx+'Size:0'] = -1.80*wave/(np.pi*costh)
1844    elif calcControls[phfx+'SizeType'] == 'uniaxial':
1845        H = np.array(refl[:3])
1846        P = np.array(calcControls[phfx+'SizeAxis'])
1847        cosP,sinP = G2lat.CosSinAngle(H,P,G)
1848        Si = parmDict[phfx+'Size:0']
1849        Sa = parmDict[phfx+'Size:1']
1850        gami = (1.80*wave/np.pi)/(Si*Sa)
1851        sqtrm = np.sqrt((cosP*Sa)**2+(sinP*Si)**2)
1852        gam = gami*sqtrm/costh           
1853        gamDict[phfx+'Size:0'] = gami*Si*sinP**2/(sqtrm*costh)-gam/Si
1854        gamDict[phfx+'Size:1'] = gami*Sa*cosP**2/(sqtrm*costh)-gam/Sa         
1855    else:           #ellipsoidal crystallites - do numerically? - not right not make sense
1856        H = np.array(refl[:3])
1857        gam = 1.8*wave/(np.pi*costh*np.inner(H,np.inner(sizeEllipse,H)))
1858    #microstrain derivatives               
1859    if calcControls[phfx+'MustrainType'] == 'isotropic':
1860        gamDict[phfx+'Mustrain:0'] =  0.018*tanth/np.pi           
1861    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1862        H = np.array(refl[:3])
1863        P = np.array(calcControls[phfx+'MustrainAxis'])
1864        cosP,sinP = G2lat.CosSinAngle(H,P,G)
1865        Si = parmDict[phfx+'Mustrain:0']
1866        Sa = parmDict[phfx+'Mustrain:1']
1867        gami = 0.018*Si*Sa*tanth/np.pi
1868        sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1869        gam = gami/sqtrm
1870        gamDict[phfx+'Mustrain:0'] = gam/Si-gami*Si*cosP**2/sqtrm**3
1871        gamDict[phfx+'Mustrain:1'] = gam/Sa-gami*Sa*sinP**2/sqtrm**3
1872    else:       #generalized - P.W. Stephens model
1873        pwrs = calcControls[phfx+'MuPwrs']
1874        const = 0.018*refl[4]**2*tanth
1875        for i,pwr in enumerate(pwrs):
1876            gamDict[phfx+'Mustrain:'+str(i)] = const*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
1877    return gamDict
1878       
1879def GetReflPos(refl,wave,G,hfx,calcControls,parmDict):
1880    h,k,l = refl[:3]
1881    dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G)
1882    d = np.sqrt(dsq)
1883    refl[4] = d
1884    pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
1885    const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1886    if 'Bragg' in calcControls[hfx+'instType']:
1887        pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
1888            parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
1889    else:               #Debye-Scherrer - simple but maybe not right
1890        pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
1891    return pos
1892
1893def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict):
1894    dpr = 180./np.pi
1895    h,k,l = refl[:3]
1896    dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
1897    dst = np.sqrt(dstsq)
1898    pos = refl[5]
1899    const = dpr/np.sqrt(1.0-wave*dst/4.0)
1900    dpdw = const*dst
1901    dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])
1902    dpdA *= const*wave/(2.0*dst)
1903    dpdZ = 1.0
1904    const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1905    if 'Bragg' in calcControls[hfx+'instType']:
1906        dpdSh = -4.*const*cosd(pos/2.0)
1907        dpdTr = -const*sind(pos)*100.0
1908        return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.
1909    else:               #Debye-Scherrer - simple but maybe not right
1910        dpdXd = -const*cosd(pos)
1911        dpdYd = -const*sind(pos)
1912        return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd
1913           
1914def GetHStrainShift(refl,SGData,phfx,parmDict):
1915    laue = SGData['SGLaue']
1916    uniq = SGData['SGUniq']
1917    h,k,l = refl[:3]
1918    if laue in ['m3','m3m']:
1919        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
1920            refl[4]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
1921    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1922        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
1923    elif laue in ['3R','3mR']:
1924        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
1925    elif laue in ['4/m','4/mmm']:
1926        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
1927    elif laue in ['mmm']:
1928        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1929    elif laue in ['2/m']:
1930        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1931        if uniq == 'a':
1932            Dij += parmDict[phfx+'D23']*k*l
1933        elif uniq == 'b':
1934            Dij += parmDict[phfx+'D13']*h*l
1935        elif uniq == 'c':
1936            Dij += parmDict[phfx+'D12']*h*k
1937    else:
1938        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
1939            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
1940    return Dij*refl[4]**2*tand(refl[5]/2.0)
1941           
1942def GetHStrainShiftDerv(refl,SGData,phfx):
1943    laue = SGData['SGLaue']
1944    uniq = SGData['SGUniq']
1945    h,k,l = refl[:3]
1946    if laue in ['m3','m3m']:
1947        dDijDict = {phfx+'D11':h**2+k**2+l**2,
1948            phfx+'eA':((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
1949    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1950        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
1951    elif laue in ['3R','3mR']:
1952        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
1953    elif laue in ['4/m','4/mmm']:
1954        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
1955    elif laue in ['mmm']:
1956        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1957    elif laue in ['2/m']:
1958        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1959        if uniq == 'a':
1960            dDijDict[phfx+'D23'] = k*l
1961        elif uniq == 'b':
1962            dDijDict[phfx+'D13'] = h*l
1963        elif uniq == 'c':
1964            dDijDict[phfx+'D12'] = h*k
1965            names.append()
1966    else:
1967        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
1968            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
1969    for item in dDijDict:
1970        dDijDict[item] *= refl[4]**2*tand(refl[5]/2.0)
1971    return dDijDict
1972   
1973def GetFprime(controlDict,Histograms):
1974    FFtables = controlDict['FFtables']
1975    if not FFtables:
1976        return
1977    for histogram in Histograms:
1978        if 'PWDR' in histogram[:4]:
1979            Histogram = Histograms[histogram]
1980            hId = Histogram['hId']
1981            hfx = ':%d:'%(hId)
1982            keV = controlDict[hfx+'keV']
1983            for El in FFtables:
1984                Orbs = G2el.GetXsectionCoeff(El.split('+')[0].split('-')[0])
1985                FP,FPP,Mu = G2el.FPcalc(Orbs, keV)
1986                FFtables[El][hfx+'FP'] = FP
1987                FFtables[El][hfx+'FPP'] = FPP               
1988           
1989def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
1990   
1991    def GetReflSIgGam(refl,wave,G,hfx,phfx,calcControls,parmDict,sizeEllipse):
1992        U = parmDict[hfx+'U']
1993        V = parmDict[hfx+'V']
1994        W = parmDict[hfx+'W']
1995        X = parmDict[hfx+'X']
1996        Y = parmDict[hfx+'Y']
1997        tanPos = tand(refl[5]/2.0)
1998        sig = U*tanPos**2+V*tanPos+W        #save peak sigma
1999        sig = max(0.001,sig)
2000        gam = X/cosd(refl[5]/2.0)+Y*tanPos+GetSampleGam(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse) #save peak gamma
2001        gam = max(0.001,gam)
2002        return sig,gam
2003               
2004    hId = Histogram['hId']
2005    hfx = ':%d:'%(hId)
2006    bakType = calcControls[hfx+'bakType']
2007    yb = G2pwd.getBackground(hfx,parmDict,bakType,x)
2008    yc = np.zeros_like(yb)
2009       
2010    if 'C' in calcControls[hfx+'histType']:   
2011        shl = max(parmDict[hfx+'SH/L'],0.002)
2012        Ka2 = False
2013        if hfx+'Lam1' in parmDict.keys():
2014            wave = parmDict[hfx+'Lam1']
2015            Ka2 = True
2016            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2017            kRatio = parmDict[hfx+'I(L2)/I(L1)']
2018        else:
2019            wave = parmDict[hfx+'Lam']
2020    else:
2021        print 'TOF Undefined at present'
2022        raise ValueError
2023    for phase in Histogram['Reflection Lists']:
2024        refList = Histogram['Reflection Lists'][phase]
2025        Phase = Phases[phase]
2026        pId = Phase['pId']
2027        pfx = '%d::'%(pId)
2028        phfx = '%d:%d:'%(pId,hId)
2029        hfx = ':%d:'%(hId)
2030        SGData = Phase['General']['SGData']
2031        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2032        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2033        Vst = np.sqrt(nl.det(G))    #V*
2034        if 'Pawley' not in Phase['General']['Type']:
2035            refList = StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict)
2036        sizeEllipse = []
2037        if calcControls[phfx+'SizeType'] == 'ellipsoidal':
2038            sizeEllipse = G2lat.U6toUij([parmDIct[phfx+'Size:%d'%(i)] for i in range(6)])
2039        for refl in refList:
2040            if 'C' in calcControls[hfx+'histType']:
2041                h,k,l = refl[:3]
2042                refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict)         #corrected reflection position
2043                Lorenz = 1./(2.*sind(refl[5]/2.)**2*cosd(refl[5]/2.))           #Lorentz correction
2044                refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict)               #apply hydrostatic strain shift
2045                refl[6:8] = GetReflSIgGam(refl,wave,G,hfx,phfx,calcControls,parmDict,sizeEllipse)    #peak sig & gam
2046                GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)    #puts corrections in refl[13]
2047                refl[13] *= Vst*Lorenz
2048                if 'Pawley' in Phase['General']['Type']:
2049                    try:
2050                        refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])])
2051                    except KeyError:
2052#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
2053                        continue
2054                Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
2055                iBeg = np.searchsorted(x,refl[5]-fmin)
2056                iFin = np.searchsorted(x,refl[5]+fmax)
2057                if not iBeg+iFin:       #peak below low limit - skip peak
2058                    continue
2059                elif not iBeg-iFin:     #peak above high limit - done
2060                    break
2061                yc[iBeg:iFin] += refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
2062                if Ka2:
2063                    pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
2064                    Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
2065                    iBeg = np.searchsorted(x,pos2-fmin)
2066                    iFin = np.searchsorted(x,pos2+fmax)
2067                    if not iBeg+iFin:       #peak below low limit - skip peak
2068                        continue
2069                    elif not iBeg-iFin:     #peak above high limit - done
2070                        return yc,yb
2071                    yc[iBeg:iFin] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin])        #and here
2072            elif 'T' in calcControls[hfx+'histType']:
2073                print 'TOF Undefined at present'
2074                raise Exception    #no TOF yet
2075    return yc,yb
2076   
2077def GetFobsSq(Histograms,Phases,parmDict,calcControls):
2078    for histogram in Histograms:
2079        if 'PWDR' in histogram[:4]:
2080            Histogram = Histograms[histogram]
2081            hId = Histogram['hId']
2082            hfx = ':%d:'%(hId)
2083            Limits = calcControls[hfx+'Limits']
2084            shl = max(parmDict[hfx+'SH/L'],0.002)
2085            Ka2 = False
2086            kRatio = 0.0
2087            if hfx+'Lam1' in parmDict.keys():
2088                Ka2 = True
2089                lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2090                kRatio = parmDict[hfx+'I(L2)/I(L1)']
2091            x,y,w,yc,yb,yd = Histogram['Data']
2092            ymb = np.array(y-yb)
2093            ycmb = np.array(yc-yb)
2094            ratio = np.where(ycmb!=0.,ymb/ycmb,0.0)           
2095            xB = np.searchsorted(x,Limits[0])
2096            xF = np.searchsorted(x,Limits[1])
2097            refLists = Histogram['Reflection Lists']
2098            for phase in refLists:
2099                Phase = Phases[phase]
2100                pId = Phase['pId']
2101                phfx = '%d:%d:'%(pId,hId)
2102                refList = refLists[phase]
2103                sumFo = 0.0
2104                sumdF = 0.0
2105                sumFosq = 0.0
2106                sumdFsq = 0.0
2107                for refl in refList:
2108                    if 'C' in calcControls[hfx+'histType']:
2109                        yp = np.zeros_like(yb)
2110                        Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
2111                        iBeg = np.searchsorted(x[xB:xF],refl[5]-fmin)
2112                        iFin = np.searchsorted(x[xB:xF],refl[5]+fmax)
2113                        iFin2 = iFin
2114                        yp[iBeg:iFin] = refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here                           
2115                        if Ka2:
2116                            pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
2117                            Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
2118                            iBeg2 = np.searchsorted(x,pos2-fmin)
2119                            iFin2 = np.searchsorted(x,pos2+fmax)
2120                            yp[iBeg2:iFin2] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])        #and here
2121                        refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[13]*(1.+kRatio)),0.0))
2122                    elif 'T' in calcControls[hfx+'histType']:
2123                        print 'TOF Undefined at present'
2124                        raise Exception    #no TOF yet
2125                    Fo = np.sqrt(np.abs(refl[8]))
2126                    Fc = np.sqrt(np.abs(refl[9]))
2127                    sumFo += Fo
2128                    sumFosq += refl[8]**2
2129                    sumdF += np.abs(Fo-Fc)
2130                    sumdFsq += (refl[8]-refl[9])**2
2131                Histogram[phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
2132                Histogram[phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
2133                Histogram[phfx+'Nref'] = len(refList)
2134               
2135def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
2136   
2137    def cellVaryDerv(pfx,SGData,dpdA): 
2138        if SGData['SGLaue'] in ['-1',]:
2139            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
2140                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
2141        elif SGData['SGLaue'] in ['2/m',]:
2142            if SGData['SGUniq'] == 'a':
2143                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
2144            elif SGData['SGUniq'] == 'b':
2145                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
2146            else:
2147                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
2148        elif SGData['SGLaue'] in ['mmm',]:
2149            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
2150        elif SGData['SGLaue'] in ['4/m','4/mmm']:
2151            return [[pfx+'A0',dpdA[0]+dpdA[1]],[pfx+'A2',dpdA[2]]]
2152        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
2153            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[3]],[pfx+'A2',dpdA[2]]]
2154        elif SGData['SGLaue'] in ['3R', '3mR']:
2155            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
2156        elif SGData['SGLaue'] in ['m3m','m3']:
2157            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]]]
2158   
2159    lenX = len(x)               
2160    hId = Histogram['hId']
2161    hfx = ':%d:'%(hId)
2162    bakType = calcControls[hfx+'bakType']
2163    dMdv = np.zeros(shape=(len(varylist),len(x)))
2164    if hfx+'Back:0' in varylist:
2165        dMdb = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x)
2166        bBpos =varylist.index(hfx+'Back:0')
2167        dMdv[bBpos:bBpos+len(dMdb)] = dMdb
2168       
2169    if 'C' in calcControls[hfx+'histType']:   
2170        dx = x[1]-x[0]
2171        shl = max(parmDict[hfx+'SH/L'],0.002)
2172        Ka2 = False
2173        if hfx+'Lam1' in parmDict.keys():
2174            wave = parmDict[hfx+'Lam1']
2175            Ka2 = True
2176            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2177            kRatio = parmDict[hfx+'I(L2)/I(L1)']
2178        else:
2179            wave = parmDict[hfx+'Lam']
2180    else:
2181        print 'TOF Undefined at present'
2182        raise ValueError
2183    for phase in Histogram['Reflection Lists']:
2184        refList = Histogram['Reflection Lists'][phase]
2185        Phase = Phases[phase]
2186        SGData = Phase['General']['SGData']
2187        pId = Phase['pId']
2188        pfx = '%d::'%(pId)
2189        phfx = '%d:%d:'%(pId,hId)
2190        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2191        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2192        if 'Pawley' not in Phase['General']['Type']:
2193            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)
2194        sizeEllipse = []
2195        if calcControls[phfx+'SizeType'] == 'ellipsoidal':
2196            sizeEllipse = G2lat.U6toUij([parmDIct[phfx+'Size:%d'%(i)] for i in range(6)])
2197        for iref,refl in enumerate(refList):
2198            if 'C' in calcControls[hfx+'histType']:        #CW powder
2199                h,k,l = refl[:3]
2200                dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA = GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2201                if 'Pawley' in Phase['General']['Type']:
2202                    try:
2203                        refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])])
2204                    except KeyError:
2205#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
2206                        continue
2207                Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
2208                iBeg = np.searchsorted(x,refl[5]-fmin)
2209                iFin = np.searchsorted(x,refl[5]+fmax)
2210                if not iBeg+iFin:       #peak below low limit - skip peak
2211                    continue
2212                elif not iBeg-iFin:     #peak above high limit - done
2213                    break
2214                pos = refl[5]
2215                tanth = tand(pos/2.0)
2216                costh = cosd(pos/2.0)
2217                dMdpk = np.zeros(shape=(6,len(x)))
2218                dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])
2219                for i in range(1,5):
2220                    dMdpk[i][iBeg:iFin] += 100.*dx*refl[13]*refl[9]*dMdipk[i]
2221                dMdpk[0][iBeg:iFin] += 100.*dx*refl[13]*refl[9]*dMdipk[0]
2222                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
2223                if Ka2:
2224                    pos2 = refl[5]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
2225                    kdelt = int((pos2-refl[5])/dx)               
2226                    iBeg2 = min(lenX,iBeg+kdelt)
2227                    iFin2 = min(lenX,iFin+kdelt)
2228                    if iBeg2-iFin2:
2229                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])
2230                        for i in range(1,5):
2231                            dMdpk[i][iBeg2:iFin2] += 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[i]
2232                        dMdpk[0][iBeg2:iFin2] += 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[0]
2233                        dMdpk[5][iBeg2:iFin2] += 100.*dx*refl[13]*dMdipk2[0]
2234                        dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*refl[9]}
2235                if 'Pawley' in Phase['General']['Type']:
2236                    try:
2237                        idx = varylist.index(pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)]))
2238                        dMdv[idx] = dervDict['int']/refl[9]
2239                    except ValueError:
2240                        pass
2241                dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
2242                names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
2243                    hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
2244                    hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
2245                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
2246                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
2247                    hfx+'DisplaceY':[dpdY,'pos'],}
2248                for name in names:
2249                    if name in varylist:
2250                        item = names[name]
2251                        dMdv[varylist.index(name)] += item[0]*dervDict[item[1]]
2252                for iPO in dIdPO:
2253                    if iPO in varylist:
2254                        dMdv[varylist.index(iPO)] += dIdPO[iPO]*dervDict['int']
2255                for i,name in enumerate(['omega','chi','phi']):
2256                    aname = pfx+'SH '+name
2257                    if aname in varylist:
2258                        dMdv[varylist.index(aname)] += dFdSA[i]*dervDict['int']
2259                for iSH in dFdODF:
2260                    if iSH in varylist:
2261                        dMdv[varylist.index(iSH)] += dFdODF[iSH]*dervDict['int']
2262                cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
2263                for name,dpdA in cellDervNames:
2264                    if name in varylist:
2265                        dMdv[varylist.index(name)] += dpdA*dervDict['pos']
2266                dDijDict = GetHStrainShiftDerv(refl,SGData,phfx)
2267                for name in dDijDict:
2268                    if name in varylist:
2269                        dMdv[varylist.index(name)] += dDijDict[name]*dervDict['pos']
2270                gamDict = GetSampleGamDerv(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse)
2271                for name in gamDict:
2272                    if name in varylist:
2273                        dMdv[varylist.index(name)] += gamDict[name]*dervDict['gam']
2274                                               
2275            elif 'T' in calcControls[hfx+'histType']:
2276                print 'TOF Undefined at present'
2277                raise Exception    #no TOF yet
2278#do atom derivatives -  for F,X & U so far             
2279            corr = dervDict['int']/refl[9]
2280            for name in varylist:
2281                try:
2282                    aname = name.split(pfx)[1][:2]
2283                    if aname in ['Af','dA','AU']:
2284                         dMdv[varylist.index(name)] += dFdvDict[name][iref]*corr
2285                except IndexError:
2286                    pass
2287    return dMdv
2288
2289def dervRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
2290    parmdict.update(zip(varylist,values))
2291    G2mv.Dict2Map(parmdict)
2292    Histograms,Phases = HistoPhases
2293    dMdv = np.empty(0)
2294    for histogram in Histograms:
2295        if 'PWDR' in histogram[:4]:
2296            Histogram = Histograms[histogram]
2297            hId = Histogram['hId']
2298            hfx = ':%d:'%(hId)
2299            Limits = calcControls[hfx+'Limits']
2300            x,y,w,yc,yb,yd = Histogram['Data']
2301            xB = np.searchsorted(x,Limits[0])
2302            xF = np.searchsorted(x,Limits[1])
2303            dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF],
2304                varylist,Histogram,Phases,calcControls,pawleyLookup)
2305            if len(dMdv):
2306                dMdv = np.concatenate((dMdv.T,dMdvh.T)).T
2307            else:
2308                dMdv = dMdvh
2309    return dMdv
2310
2311def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):       
2312    parmdict.update(zip(varylist,values))
2313    Values2Dict(parmdict, varylist, values)
2314    G2mv.Dict2Map(parmdict)
2315    Histograms,Phases = HistoPhases
2316    M = np.empty(0)
2317    sumwYo = 0
2318    Nobs = 0
2319    for histogram in Histograms:
2320        if 'PWDR' in histogram[:4]:
2321            Histogram = Histograms[histogram]
2322            hId = Histogram['hId']
2323            hfx = ':%d:'%(hId)
2324            Limits = calcControls[hfx+'Limits']
2325            x,y,w,yc,yb,yd = Histogram['Data']
2326            yc *= 0.0                           #zero full calcd profiles
2327            yb *= 0.0
2328            yd *= 0.0
2329            xB = np.searchsorted(x,Limits[0])
2330            xF = np.searchsorted(x,Limits[1])
2331            Histogram['Nobs'] = xF-xB
2332            Nobs += Histogram['Nobs']
2333            Histogram['sumwYo'] = np.sum(w[xB:xF]*y[xB:xF]**2)
2334            sumwYo += Histogram['sumwYo']
2335            yc[xB:xF],yb[xB:xF] = getPowderProfile(parmdict,x[xB:xF],
2336                varylist,Histogram,Phases,calcControls,pawleyLookup)
2337            yc[xB:xF] += yb[xB:xF]
2338            yd[xB:xF] = yc[xB:xF]-y[xB:xF]          #yc-yo then all dydv have no '-' needed
2339            Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF]))
2340            wdy = np.sqrt(w[xB:xF])*(yd[xB:xF])
2341            Histogram['wRp'] = min(100.,np.sqrt(np.sum(wdy**2)/Histogram['sumwYo'])*100.)
2342            M = np.concatenate((M,wdy))
2343    Histograms['sumwYo'] = sumwYo
2344    Histograms['Nobs'] = Nobs
2345    Rwp = min(100.,np.sqrt(np.sum(M**2)/sumwYo)*100.)
2346    if dlg:
2347        GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Powder profile wRp =',Rwp,'%'))[0]
2348        if not GoOn:
2349            parmDict['saved values'] = values
2350            raise Exception         #Abort!!
2351    return M
2352   
2353                   
2354def Refine(GPXfile,dlg):
2355    import cPickle
2356    import pytexture as ptx
2357    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
2358   
2359    ShowBanner()
2360    varyList = []
2361    parmDict = {}
2362    calcControls = {}
2363    G2mv.InitVars()   
2364    Controls = GetControls(GPXfile)
2365    ShowControls(Controls)           
2366    constrList = GetConstraints(GPXfile)
2367    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
2368    if not Phases:
2369        print ' *** ERROR - you have no histograms to refine! ***'
2370        print ' *** Refine aborted ***'
2371        raise Exception
2372    if not Histograms:
2373        print ' *** ERROR - you have no data to refine with! ***'
2374        print ' *** Refine aborted ***'
2375        raise Exception
2376    Natoms,phaseVary,phaseDict,pawleyLookup,FFtables = GetPhaseData(Phases)
2377    calcControls['Natoms'] = Natoms
2378    calcControls['FFtables'] = FFtables
2379    hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms)
2380    calcControls.update(controlDict)
2381    histVary,histDict,controlDict = GetHistogramData(Histograms)
2382    calcControls.update(controlDict)
2383    varyList = phaseVary+hapVary+histVary
2384    parmDict.update(phaseDict)
2385    parmDict.update(hapDict)
2386    parmDict.update(histDict)
2387    GetFprime(calcControls,Histograms)
2388    for item in constrList: print item
2389    constrDict,constrFlag,fixedList = G2mv.InputParse([])        #constraints go here?
2390    groups,parmlist = G2mv.GroupConstraints(constrDict)
2391    G2mv.GenerateConstraints(groups,parmlist,constrDict,constrFlag,fixedList)
2392    G2mv.Map2Dict(parmDict,varyList)
2393
2394    while True:
2395        begin = time.time()
2396        values =  np.array(Dict2Values(parmDict, varyList))
2397        Ftol = Controls['min dM/M']
2398        Factor = Controls['shift factor']
2399        if Controls['deriv type'] == 'analytic':
2400            result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
2401                ftol=Ftol,col_deriv=True,factor=Factor,
2402                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
2403            ncyc = int(result[2]['nfev']/2)               
2404        else:           #'numeric'
2405            result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
2406                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
2407            ncyc = int(result[2]['nfev']/len(varyList))
2408#        table = dict(zip(varyList,zip(values,result[0],(result[0]-values))))
2409#        for item in table: print item,table[item]               #useful debug - are things shifting?
2410        runtime = time.time()-begin
2411        chisq = np.sum(result[2]['fvec']**2)
2412        Values2Dict(parmDict, varyList, result[0])
2413        G2mv.Dict2Map(parmDict)
2414        newAtomDict = ApplyXYZshifts(parmDict,varyList)
2415       
2416        Rwp = np.sqrt(chisq/Histograms['sumwYo'])*100.      #to %
2417        GOF = chisq/(Histograms['Nobs']-len(varyList))
2418        print '\n Refinement results:'
2419        print 135*'-'
2420        print ' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histograms['Nobs'],' Number of parameters: ',len(varyList)
2421        print ' Refinement time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
2422        print ' wRp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
2423        try:
2424            covMatrix = result[1]*GOF
2425            sig = np.sqrt(np.diag(covMatrix))
2426            if np.any(np.isnan(sig)):
2427                print '*** Least squares aborted - some invalid esds possible ***'
2428#            table = dict(zip(varyList,zip(values,result[0],(result[0]-values)/sig)))
2429#            for item in table: print item,table[item]               #useful debug - are things shifting?
2430            break                   #refinement succeeded - finish up!
2431        except TypeError:          #result[1] is None on singular matrix
2432            print '**** Refinement failed - singular matrix ****'
2433            Ipvt = result[2]['ipvt']
2434            for i,ipvt in enumerate(Ipvt):
2435                if not np.sum(result[2]['fjac'],axis=1)[i]:
2436                    print 'Removing parameter: ',varyList[ipvt-1]
2437                    del(varyList[ipvt-1])
2438                    break
2439
2440#    print 'dependentParmList: ',G2mv.dependentParmList
2441#    print 'arrayList: ',G2mv.arrayList
2442#    print 'invarrayList: ',G2mv.invarrayList
2443#    print 'indParmList: ',G2mv.indParmList
2444#    print 'fixedDict: ',G2mv.fixedDict
2445    GetFobsSq(Histograms,Phases,parmDict,calcControls)
2446    sigDict = dict(zip(varyList,sig))
2447    covData = {'variables':result[0],'varyList':varyList,'sig':sig,
2448        'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict}
2449    SetPhaseData(parmDict,sigDict,Phases,covData)
2450    SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms)
2451    SetHistogramData(parmDict,sigDict,Histograms)
2452    SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,covData)
2453#for testing purposes!!!
2454#    file = open('structTestdata.dat','wb')
2455#    cPickle.dump(parmDict,file,1)
2456#    cPickle.dump(varyList,file,1)
2457#    for histogram in Histograms:
2458#        if 'PWDR' in histogram[:4]:
2459#            Histogram = Histograms[histogram]
2460#    cPickle.dump(Histogram,file,1)
2461#    cPickle.dump(Phases,file,1)
2462#    cPickle.dump(calcControls,file,1)
2463#    cPickle.dump(pawleyLookup,file,1)
2464#    file.close()
2465
2466def SeqRefine(GPXfile,dlg):
2467    import cPickle
2468    import pytexture as ptx
2469    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
2470   
2471    ShowBanner()
2472    print ' Sequential Refinement'
2473    G2mv.InitVars()   
2474    Controls = GetControls(GPXfile)
2475    ShowControls(Controls)           
2476    constrList = GetConstraints(GPXfile)
2477    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
2478    if not Phases:
2479        print ' *** ERROR - you have no histograms to refine! ***'
2480        print ' *** Refine aborted ***'
2481        raise Exception
2482    if not Histograms:
2483        print ' *** ERROR - you have no data to refine with! ***'
2484        print ' *** Refine aborted ***'
2485        raise Exception
2486    Natoms,phaseVary,phaseDict,pawleyLookup,FFtables = GetPhaseData(Phases,False)
2487    if 'Seq Data' in Controls:
2488        histNames = Controls['Seq Data']
2489    else:
2490        histNames = GetHistogramNames(GPXfile,['PWDR',])
2491    if 'Reverse Seq' in Controls:
2492        if Controls['Reverse Seq']:
2493            histNames.reverse()
2494    SeqResult = {'histNames':histNames}
2495    makeBack = True
2496    for ihst,histogram in enumerate(histNames):
2497        ifPrint = False
2498        if dlg:
2499            dlg.SetTitle('Residual for histogram '+str(ihst))
2500        calcControls = {}
2501        calcControls['Natoms'] = Natoms
2502        calcControls['FFtables'] = FFtables
2503        varyList = []
2504        parmDict = {}
2505        Histo = {histogram:Histograms[histogram],}
2506        hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histo,False)
2507        calcControls.update(controlDict)
2508        histVary,histDict,controlDict = GetHistogramData(Histo,False)
2509        calcControls.update(controlDict)
2510        varyList = phaseVary+hapVary+histVary
2511        if not ihst:
2512            saveVaryList = varyList[:]
2513            for i,item in enumerate(saveVaryList):
2514                items = item.split(':')
2515                if items[1]:
2516                    items[1] = ''
2517                item = ':'.join(items)
2518                saveVaryList[i] = item
2519            SeqResult['varyList'] = saveVaryList
2520        else:
2521            newVaryList = varyList[:]
2522            for i,item in enumerate(newVaryList):
2523                items = item.split(':')
2524                if items[1]:
2525                    items[1] = ''
2526                item = ':'.join(items)
2527                newVaryList[i] = item
2528            if newVaryList != SeqResult['varyList']:
2529                print newVaryList
2530                print SeqResult['varyList']
2531                print '**** ERROR - variable list for this histogram does not match previous'
2532                raise Exception
2533        parmDict.update(phaseDict)
2534        parmDict.update(hapDict)
2535        parmDict.update(histDict)
2536        GetFprime(calcControls,Histo)
2537        constrDict,constrFlag,fixedList = G2mv.InputParse([])        #constraints go here?
2538        groups,parmlist = G2mv.GroupConstraints(constrDict)
2539        G2mv.GenerateConstraints(groups,parmlist,constrDict,constrFlag,fixedList)
2540        G2mv.Map2Dict(parmDict,varyList)
2541   
2542        while True:
2543            begin = time.time()
2544            values =  np.array(Dict2Values(parmDict, varyList))
2545            Ftol = Controls['min dM/M']
2546            Factor = Controls['shift factor']
2547            if Controls['deriv type'] == 'analytic':
2548                result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
2549                    ftol=Ftol,col_deriv=True,factor=Factor,
2550                    args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
2551                ncyc = int(result[2]['nfev']/2)               
2552            else:           #'numeric'
2553                result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
2554                    args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
2555                ncyc = int(result[2]['nfev']/len(varyList))
2556            runtime = time.time()-begin
2557            chisq = np.sum(result[2]['fvec']**2)
2558            Values2Dict(parmDict, varyList, result[0])
2559            G2mv.Dict2Map(parmDict)
2560            newAtomDict = ApplyXYZshifts(parmDict,varyList)
2561           
2562            Rwp = np.sqrt(chisq/Histo['sumwYo'])*100.      #to %
2563            GOF = chisq/(Histo['Nobs']-len(varyList))
2564            print '\n Refinement results for histogram: v'+histogram
2565            print 135*'-'
2566            print ' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histo['Nobs'],' Number of parameters: ',len(varyList)
2567            print ' Refinement time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
2568            print ' wRp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
2569            try:
2570                covMatrix = result[1]*GOF
2571                sig = np.sqrt(np.diag(covMatrix))
2572                if np.any(np.isnan(sig)):
2573                    print '*** Least squares aborted - some invalid esds possible ***'
2574                    ifPrint = True
2575                break                   #refinement succeeded - finish up!
2576            except TypeError:          #result[1] is None on singular matrix
2577                print '**** Refinement failed - singular matrix ****'
2578                Ipvt = result[2]['ipvt']
2579                for i,ipvt in enumerate(Ipvt):
2580                    if not np.sum(result[2]['fjac'],axis=1)[i]:
2581                        print 'Removing parameter: ',varyList[ipvt-1]
2582                        del(varyList[ipvt-1])
2583                        break
2584   
2585        GetFobsSq(Histo,Phases,parmDict,calcControls)
2586        sigDict = dict(zip(varyList,sig))
2587        covData = {'variables':result[0],'varyList':varyList,'sig':sig,
2588            'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict}
2589        SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,ifPrint)
2590        SetHistogramData(parmDict,sigDict,Histo,ifPrint)
2591        SeqResult[histogram] = covData
2592        SetUsedHistogramsAndPhases(GPXfile,Histo,Phases,covData,makeBack)
2593        makeBack = False
2594    SetSeqResult(GPXfile,Histograms,SeqResult)
2595
2596
2597def main():
2598    arg = sys.argv
2599    if len(arg) > 1:
2600        GPXfile = arg[1]
2601        if not ospath.exists(GPXfile):
2602            print 'ERROR - ',GPXfile," doesn't exist!"
2603            exit()
2604        GPXpath = ospath.dirname(arg[1])
2605        Refine(GPXfile,None)
2606    else:
2607        print 'ERROR - missing filename'
2608        exit()
2609    print "Done"
2610         
2611if __name__ == '__main__':
2612    main()
Note: See TracBrowser for help on using the repository browser.