source: trunk/GSASIIstruct.py @ 399

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

Implement sequential refinement
remove print "load" & "save" for each item in Tree
revise application of azimuth offset - azimuths are now all "true" with correction

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