source: trunk/GSASIIstruct.py @ 406

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

implement background subtraction during image integration
implement copy for instrument parms & background
continue constraint GUI development
make sure proper updates in refinement/seq refinement

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