source: trunk/GSASIIstruct.py @ 484

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

change authorship
split GSASIIelem into a GUI & non-GUI parts
replace MsgDialogs? in GSASIIElem.py with simple prints
cleanup & refactor distance/angle/torsion calcs.

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