source: trunk/GSASIIstruct.py @ 408

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

add copy to limits menu
more constraint GUI stuff

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