source: trunk/GSASIIstruct.py @ 630

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

major mods for HKLF data
remove some dead code
mark more code as dead (#)
implement cif data style as 'val(esd)' for f & f_squared
continue implementation of HKLF data in refinement
HKLF now OK in Fourier & charge flip calcs.

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