source: trunk/GSASIIstruct.py @ 482

Last change on this file since 482 was 482, checked in by vondreele, 11 years ago

move FileDlgFixExt? to G2IO
begin to add a save selected seq results to file
move IO routines called in G2struct back there
remove wx stuff from G2pwd & G2struct
put progress bar in G2pwdGUI (not G2pwd)

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