source: trunk/GSASIIstruct.py @ 656

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

fix Debye background model

  • Property svn:keywords set to Date Author Revision URL Id
File size: 163.1 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASIIstructure - structure computation routines
3########### SVN repository information ###################
4# $Date: 2012-06-26 14:55:52 +0000 (Tue, 26 Jun 2012) $
5# $Author: vondreele $
6# $Revision: 656 $
7# $URL: trunk/GSASIIstruct.py $
8# $Id: GSASIIstruct.py 656 2012-06-26 14:55:52Z 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,Print=False)
130    hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms,Print=False)
131    histVary,histDict,controlDict = GetHistogramData(Histograms,Print=False)
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            elif abs(dspI-dspJ) < 1.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                        names = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)]
687                        equivs = [[],[],[]]
688                        for j in range(3):
689                            if xId[j] > 0:                               
690                                phaseVary.append(names[j])
691                                equivs[xId[j]-1].append([names[j],xCoef[j]])
692                        for equiv in equivs:
693                            if len(equiv) > 1:
694                                name = equiv[0][0]
695                                for eqv in equiv[1:]:
696                                    G2mv.StoreEquivalence(name,(eqv,))
697                    if 'U' in at[2]:
698                        if at[9] == 'I':
699                            phaseVary.append(pfx+'AUiso:'+str(i))
700                        else:
701                            uId,uCoef = G2spc.GetCSuinel(at[7])[:2]
702                            names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i),
703                                pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)]
704                            equivs = [[],[],[],[],[],[]]
705                            for j in range(6):
706                                if uId[j] > 0:                               
707                                    phaseVary.append(names[j])
708                                    equivs[uId[j]-1].append([names[j],uCoef[j]])
709                            for equiv in equivs:
710                                if len(equiv) > 1:
711                                    name = equiv[0][0]
712                                    for eqv in equiv[1:]:
713                                        G2mv.StoreEquivalence(name,(eqv,))
714#            elif General['Type'] == 'magnetic':
715#            elif General['Type'] == 'macromolecular':
716
717            textureData = General['SH Texture']
718            if textureData['Order']:
719                phaseDict[pfx+'SHorder'] = textureData['Order']
720                phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']]
721                for name in ['omega','chi','phi']:
722                    phaseDict[pfx+'SH '+name] = textureData['Sample '+name][1]
723                    if textureData['Sample '+name][0]:
724                        phaseVary.append(pfx+'SH '+name)
725                for name in textureData['SH Coeff'][1]:
726                    phaseDict[pfx+name] = textureData['SH Coeff'][1][name]
727                    if textureData['SH Coeff'][0]:
728                        phaseVary.append(pfx+name)
729               
730            if Print:
731                print '\n Phase name: ',General['Name']
732                print 135*'-'
733                PrintFFtable(FFtable)
734                PrintBLtable(BLtable)
735                print ''
736                for line in SGtext: print line
737                PrintAtoms(General,Atoms)
738                print '\n Unit cell: a =','%.5f'%(cell[1]),' b =','%.5f'%(cell[2]),' c =','%.5f'%(cell[3]), \
739                    ' alpha =','%.3f'%(cell[4]),' beta =','%.3f'%(cell[5]),' gamma =', \
740                    '%.3f'%(cell[6]),' volume =','%.3f'%(cell[7]),' Refine?',cell[0]
741                PrintTexture(textureData)
742                   
743        elif PawleyRef:
744            pawleyVary = []
745            for i,refl in enumerate(PawleyRef):
746                phaseDict[pfx+'PWLref:'+str(i)] = refl[6]
747                pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
748                if refl[5]:
749                    pawleyVary.append(pfx+'PWLref:'+str(i))
750            GetPawleyConstr(SGData['SGLaue'],PawleyRef,pawleyVary)      #does G2mv.StoreEquivalence
751            phaseVary += pawleyVary
752               
753    return Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables
754   
755def cellFill(pfx,SGData,parmDict,sigDict): 
756    if SGData['SGLaue'] in ['-1',]:
757        A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
758            parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
759        sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
760            sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
761    elif SGData['SGLaue'] in ['2/m',]:
762        if SGData['SGUniq'] == 'a':
763            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
764                parmDict[pfx+'A3'],0,0]
765            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
766                sigDict[pfx+'A3'],0,0]
767        elif SGData['SGUniq'] == 'b':
768            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
769                0,parmDict[pfx+'A4'],0]
770            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
771                0,sigDict[pfx+'A4'],0]
772        else:
773            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
774                0,0,parmDict[pfx+'A5']]
775            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
776                0,0,sigDict[pfx+'A5']]
777    elif SGData['SGLaue'] in ['mmm',]:
778        A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
779        sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
780    elif SGData['SGLaue'] in ['4/m','4/mmm']:
781        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
782        sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
783    elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
784        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
785            parmDict[pfx+'A0'],0,0]
786        sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
787    elif SGData['SGLaue'] in ['3R', '3mR']:
788        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
789            parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
790        sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
791    elif SGData['SGLaue'] in ['m3m','m3']:
792        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
793        sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
794    return A,sigA
795       
796def getCellEsd(pfx,SGData,A,covData):
797    dpr = 180./np.pi
798    rVsq = G2lat.calc_rVsq(A)
799    G,g = G2lat.A2Gmat(A)       #get recip. & real metric tensors
800    cell = np.array(G2lat.Gmat2cell(g))   #real cell
801    cellst = np.array(G2lat.Gmat2cell(G)) #recip. cell
802    scos = cosd(cellst[3:6])
803    ssin = sind(cellst[3:6])
804    scot = scos/ssin
805    rcos = cosd(cell[3:6])
806    rsin = sind(cell[3:6])
807    rcot = rcos/rsin
808    RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
809    varyList = covData['varyList']
810    covMatrix = covData['covMatrix']
811    vcov = G2mth.getVCov(RMnames,varyList,covMatrix)
812    Ax = np.array(A)
813    Ax[3:] /= 2.
814    drVdA = np.array([Ax[1]*Ax[2]-Ax[5]**2,Ax[0]*Ax[2]-Ax[4]**2,Ax[0]*Ax[1]-Ax[3]**2,
815        Ax[4]*Ax[5]-Ax[2]*Ax[3],Ax[3]*Ax[5]-Ax[1]*Ax[4],Ax[3]*Ax[4]-Ax[0]*Ax[5]])
816    srcvlsq = np.inner(drVdA,np.inner(vcov,drVdA.T))
817    Vol = 1/np.sqrt(rVsq)
818    sigVol = Vol**3*np.sqrt(srcvlsq)/2.
819    R123 = Ax[0]*Ax[1]*Ax[2]
820    dsasdg = np.zeros((3,6))
821    dadg = np.zeros((6,6))
822    for i0 in range(3):         #0  1   2
823        i1 = (i0+1)%3           #1  2   0
824        i2 = (i1+1)%3           #2  0   1
825        i3 = 5-i2               #3  5   4
826        i4 = 5-i1               #4  3   5
827        i5 = 5-i0               #5  4   3
828        dsasdg[i0][i1] = 0.5*scot[i0]*scos[i0]/Ax[i1]
829        dsasdg[i0][i2] = 0.5*scot[i0]*scos[i0]/Ax[i2]
830        dsasdg[i0][i5] = -scot[i0]/np.sqrt(Ax[i1]*Ax[i2])
831        denmsq = Ax[i0]*(R123-Ax[i1]*Ax[i4]**2-Ax[i2]*Ax[i3]**2+(Ax[i4]*Ax[i3])**2)
832        denom = np.sqrt(denmsq)
833        dadg[i5][i0] = -Ax[i5]/denom-rcos[i0]/denmsq*(R123-0.5*Ax[i1]*Ax[i4]**2-0.5*Ax[i2]*Ax[i3]**2)
834        dadg[i5][i1] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i2]-Ax[i0]*Ax[i4]**2)
835        dadg[i5][i2] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i1]-Ax[i0]*Ax[i3]**2)
836        dadg[i5][i3] = Ax[i4]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i2]*Ax[i3]-Ax[i3]*Ax[i4]**2)
837        dadg[i5][i4] = Ax[i3]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i1]*Ax[i4]-Ax[i3]**2*Ax[i4])
838        dadg[i5][i5] = -Ax[i0]/denom
839    for i0 in range(3):
840        i1 = (i0+1)%3
841        i2 = (i1+1)%3
842        i3 = 5-i2
843        for ij in range(6):
844            dadg[i0][ij] = cell[i0]*(rcot[i2]*dadg[i3][ij]/rsin[i2]-dsasdg[i1][ij]/ssin[i1])
845            if ij == i0:
846                dadg[i0][ij] = dadg[i0][ij]-0.5*cell[i0]/Ax[i0]
847            dadg[i3][ij] = -dadg[i3][ij]*rsin[2-i0]*dpr
848    sigMat = np.inner(dadg,np.inner(vcov,dadg.T))
849    var = np.diag(sigMat)
850    CS = np.where(var>0.,np.sqrt(var),0.)
851    cellSig = [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol]  #exchange sig(alp) & sig(gam) to get in right order
852    return cellSig           
853   
854def SetPhaseData(parmDict,sigDict,Phases,covData):
855   
856    def PrintAtomsAndSig(General,Atoms,atomsSig):
857        print '\n Atoms:'
858        line = '   name      x         y         z      frac   Uiso     U11     U22     U33     U12     U13     U23'
859        if General['Type'] == 'magnetic':
860            line += '   Mx     My     Mz'
861        elif General['Type'] == 'macromolecular':
862            line = ' res no  residue  chain '+line
863        print line
864        if General['Type'] == 'nuclear':
865            print 135*'-'
866            fmt = {0:'%7s',1:'%7s',3:'%10.5f',4:'%10.5f',5:'%10.5f',6:'%8.3f',10:'%8.5f',
867                11:'%8.5f',12:'%8.5f',13:'%8.5f',14:'%8.5f',15:'%8.5f',16:'%8.5f'}
868            noFXsig = {3:[10*' ','%10s'],4:[10*' ','%10s'],5:[10*' ','%10s'],6:[8*' ','%8s']}
869            for i,at in enumerate(Atoms):
870                name = fmt[0]%(at[0])+fmt[1]%(at[1])+':'
871                valstr = ' values:'
872                sigstr = ' sig   :'
873                for ind in [3,4,5,6]:
874                    sigind = str(i)+':'+str(ind)
875                    valstr += fmt[ind]%(at[ind])                   
876                    if sigind in atomsSig:
877                        sigstr += fmt[ind]%(atomsSig[sigind])
878                    else:
879                        sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
880                if at[9] == 'I':
881                    valstr += fmt[10]%(at[10])
882                    if str(i)+':10' in atomsSig:
883                        sigstr += fmt[10]%(atomsSig[str(i)+':10'])
884                    else:
885                        sigstr += 8*' '
886                else:
887                    valstr += 8*' '
888                    sigstr += 8*' '
889                    for ind in [11,12,13,14,15,16]:
890                        sigind = str(i)+':'+str(ind)
891                        valstr += fmt[ind]%(at[ind])
892                        if sigind in atomsSig:                       
893                            sigstr += fmt[ind]%(atomsSig[sigind])
894                        else:
895                            sigstr += 8*' '
896                print name
897                print valstr
898                print sigstr
899               
900    def PrintSHtextureAndSig(textureData,SHtextureSig):
901        print '\n Spherical harmonics texture: Order:' + str(textureData['Order'])
902        names = ['omega','chi','phi']
903        namstr = '  names :'
904        ptstr =  '  values:'
905        sigstr = '  esds  :'
906        for name in names:
907            namstr += '%12s'%(name)
908            ptstr += '%12.3f'%(textureData['Sample '+name][1])
909            if 'Sample '+name in SHtextureSig:
910                sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
911            else:
912                sigstr += 12*' '
913        print namstr
914        print ptstr
915        print sigstr
916        print '\n Texture coefficients:'
917        namstr = '  names :'
918        ptstr =  '  values:'
919        sigstr = '  esds  :'
920        SHcoeff = textureData['SH Coeff'][1]
921        for name in SHcoeff:
922            namstr += '%12s'%(name)
923            ptstr += '%12.3f'%(SHcoeff[name])
924            if name in SHtextureSig:
925                sigstr += '%12.3f'%(SHtextureSig[name])
926            else:
927                sigstr += 12*' '
928        print namstr
929        print ptstr
930        print sigstr
931       
932           
933    print '\n Phases:'
934    for phase in Phases:
935        print ' Result for phase: ',phase
936        Phase = Phases[phase]
937        General = Phase['General']
938        SGData = General['SGData']
939        Atoms = Phase['Atoms']
940        cell = General['Cell']
941        pId = Phase['pId']
942        pfx = str(pId)+'::'
943        if cell[0]:
944            A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
945            cellSig = getCellEsd(pfx,SGData,A,covData)  #includes sigVol
946            print ' Reciprocal metric tensor: '
947            ptfmt = "%15.9f"
948            names = ['A11','A22','A33','A12','A13','A23']
949            namstr = '  names :'
950            ptstr =  '  values:'
951            sigstr = '  esds  :'
952            for name,a,siga in zip(names,A,sigA):
953                namstr += '%15s'%(name)
954                ptstr += ptfmt%(a)
955                if siga:
956                    sigstr += ptfmt%(siga)
957                else:
958                    sigstr += 15*' '
959            print namstr
960            print ptstr
961            print sigstr
962            cell[1:7] = G2lat.A2cell(A)
963            cell[7] = G2lat.calc_V(A)
964            print ' New unit cell:'
965            ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
966            names = ['a','b','c','alpha','beta','gamma','Volume']
967            namstr = '  names :'
968            ptstr =  '  values:'
969            sigstr = '  esds  :'
970            for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
971                namstr += '%12s'%(name)
972                ptstr += fmt%(a)
973                if siga:
974                    sigstr += fmt%(siga)
975                else:
976                    sigstr += 12*' '
977            print namstr
978            print ptstr
979            print sigstr
980           
981        if Phase['General'].get('doPawley'):
982            pawleyRef = Phase['Pawley ref']
983            for i,refl in enumerate(pawleyRef):
984                key = pfx+'PWLref:'+str(i)
985                refl[6] = parmDict[key]
986                if key in sigDict:
987                    refl[7] = sigDict[key]
988                else:
989                    refl[7] = 0
990        else:
991            atomsSig = {}
992            if General['Type'] == 'nuclear':
993                for i,at in enumerate(Atoms):
994                    names = {3:pfx+'Ax:'+str(i),4:pfx+'Ay:'+str(i),5:pfx+'Az:'+str(i),6:pfx+'Afrac:'+str(i),
995                        10:pfx+'AUiso:'+str(i),11:pfx+'AU11:'+str(i),12:pfx+'AU22:'+str(i),13:pfx+'AU33:'+str(i),
996                        14:pfx+'AU12:'+str(i),15:pfx+'AU13:'+str(i),16:pfx+'AU23:'+str(i)}
997                    for ind in [3,4,5,6]:
998                        at[ind] = parmDict[names[ind]]
999                        if ind in [3,4,5]:
1000                            name = names[ind].replace('A','dA')
1001                        else:
1002                            name = names[ind]
1003                        if name in sigDict:
1004                            atomsSig[str(i)+':'+str(ind)] = sigDict[name]
1005                    if at[9] == 'I':
1006                        at[10] = parmDict[names[10]]
1007                        if names[10] in sigDict:
1008                            atomsSig[str(i)+':10'] = sigDict[names[10]]
1009                    else:
1010                        for ind in [11,12,13,14,15,16]:
1011                            at[ind] = parmDict[names[ind]]
1012                            if names[ind] in sigDict:
1013                                atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
1014            PrintAtomsAndSig(General,Atoms,atomsSig)
1015       
1016        textureData = General['SH Texture']   
1017        if textureData['Order']:
1018            SHtextureSig = {}
1019            for name in ['omega','chi','phi']:
1020                aname = pfx+'SH '+name
1021                textureData['Sample '+name][1] = parmDict[aname]
1022                if aname in sigDict:
1023                    SHtextureSig['Sample '+name] = sigDict[aname]
1024            for name in textureData['SH Coeff'][1]:
1025                aname = pfx+name
1026                textureData['SH Coeff'][1][name] = parmDict[aname]
1027                if aname in sigDict:
1028                    SHtextureSig[name] = sigDict[aname]
1029            PrintSHtextureAndSig(textureData,SHtextureSig)
1030
1031################################################################################
1032##### Histogram & Phase data
1033################################################################################       
1034                   
1035def GetHistogramPhaseData(Phases,Histograms,Print=True):
1036   
1037    def PrintSize(hapData):
1038        if hapData[0] in ['isotropic','uniaxial']:
1039            line = '\n Size model    : %9s'%(hapData[0])
1040            line += ' equatorial:'+'%12.3f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1041            if hapData[0] == 'uniaxial':
1042                line += ' axial:'+'%12.3f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1043            line += '\n\t LG mixing coeff.: %12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1044            print line
1045        else:
1046            print '\n Size model    : %s'%(hapData[0])+ \
1047                '\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1048            Snames = ['S11','S22','S33','S12','S13','S23']
1049            ptlbls = ' names :'
1050            ptstr =  ' values:'
1051            varstr = ' refine:'
1052            for i,name in enumerate(Snames):
1053                ptlbls += '%12s' % (name)
1054                ptstr += '%12.6f' % (hapData[4][i])
1055                varstr += '%12s' % (str(hapData[5][i]))
1056            print ptlbls
1057            print ptstr
1058            print varstr
1059       
1060    def PrintMuStrain(hapData,SGData):
1061        if hapData[0] in ['isotropic','uniaxial']:
1062            line = '\n Mustrain model: %9s'%(hapData[0])
1063            line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1064            if hapData[0] == 'uniaxial':
1065                line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1066            line +='\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1067            print line
1068        else:
1069            print '\n Mustrain model: %s'%(hapData[0])+ \
1070                '\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1071            Snames = G2spc.MustrainNames(SGData)
1072            ptlbls = ' names :'
1073            ptstr =  ' values:'
1074            varstr = ' refine:'
1075            for i,name in enumerate(Snames):
1076                ptlbls += '%12s' % (name)
1077                ptstr += '%12.6f' % (hapData[4][i])
1078                varstr += '%12s' % (str(hapData[5][i]))
1079            print ptlbls
1080            print ptstr
1081            print varstr
1082
1083    def PrintHStrain(hapData,SGData):
1084        print '\n Hydrostatic/elastic strain: '
1085        Hsnames = G2spc.HStrainNames(SGData)
1086        ptlbls = ' names :'
1087        ptstr =  ' values:'
1088        varstr = ' refine:'
1089        for i,name in enumerate(Hsnames):
1090            ptlbls += '%12s' % (name)
1091            ptstr += '%12.6f' % (hapData[0][i])
1092            varstr += '%12s' % (str(hapData[1][i]))
1093        print ptlbls
1094        print ptstr
1095        print varstr
1096
1097    def PrintSHPO(hapData):
1098        print '\n Spherical harmonics preferred orientation: Order:' + \
1099            str(hapData[4])+' Refine? '+str(hapData[2])
1100        ptlbls = ' names :'
1101        ptstr =  ' values:'
1102        for item in hapData[5]:
1103            ptlbls += '%12s'%(item)
1104            ptstr += '%12.3f'%(hapData[5][item]) 
1105        print ptlbls
1106        print ptstr
1107   
1108    hapDict = {}
1109    hapVary = []
1110    controlDict = {}
1111    poType = {}
1112    poAxes = {}
1113    spAxes = {}
1114    spType = {}
1115   
1116    for phase in Phases:
1117        HistoPhase = Phases[phase]['Histograms']
1118        SGData = Phases[phase]['General']['SGData']
1119        cell = Phases[phase]['General']['Cell'][1:7]
1120        A = G2lat.cell2A(cell)
1121        pId = Phases[phase]['pId']
1122        histoList = HistoPhase.keys()
1123        histoList.sort()
1124        for histogram in histoList:
1125            try:
1126                Histogram = Histograms[histogram]
1127            except KeyError:                       
1128                #skip if histogram not included e.g. in a sequential refinement
1129                continue
1130            hapData = HistoPhase[histogram]
1131            hId = Histogram['hId']
1132            if 'PWDR' in histogram:
1133                limits = Histogram['Limits'][1]
1134                inst = Histogram['Instrument Parameters']
1135                inst = dict(zip(inst[3],inst[1]))
1136                Zero = inst['Zero']
1137                if 'C' in inst['Type']:
1138                    try:
1139                        wave = inst['Lam']
1140                    except KeyError:
1141                        wave = inst['Lam1']
1142                    dmin = wave/(2.0*sind(limits[1]/2.0))
1143                pfx = str(pId)+':'+str(hId)+':'
1144                for item in ['Scale','Extinction']:
1145                    hapDict[pfx+item] = hapData[item][0]
1146                    if hapData[item][1]:
1147                        hapVary.append(pfx+item)
1148                names = G2spc.HStrainNames(SGData)
1149                for i,name in enumerate(names):
1150                    hapDict[pfx+name] = hapData['HStrain'][0][i]
1151                    if hapData['HStrain'][1][i]:
1152                        hapVary.append(pfx+name)
1153                controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
1154                if hapData['Pref.Ori.'][0] == 'MD':
1155                    hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
1156                    controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
1157                    if hapData['Pref.Ori.'][2]:
1158                        hapVary.append(pfx+'MD')
1159                else:                           #'SH' spherical harmonics
1160                    controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
1161                    controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
1162                    for item in hapData['Pref.Ori.'][5]:
1163                        hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
1164                        if hapData['Pref.Ori.'][2]:
1165                            hapVary.append(pfx+item)
1166                for item in ['Mustrain','Size']:
1167                    controlDict[pfx+item+'Type'] = hapData[item][0]
1168                    hapDict[pfx+item+':mx'] = hapData[item][1][2]
1169                    if hapData[item][2][2]:
1170                        hapVary.append(pfx+item+':mx')
1171                    if hapData[item][0] in ['isotropic','uniaxial']:
1172                        hapDict[pfx+item+':i'] = hapData[item][1][0]
1173                        if hapData[item][2][0]:
1174                            hapVary.append(pfx+item+':i')
1175                        if hapData[item][0] == 'uniaxial':
1176                            controlDict[pfx+item+'Axis'] = hapData[item][3]
1177                            hapDict[pfx+item+':a'] = hapData[item][1][1]
1178                            if hapData[item][2][1]:
1179                                hapVary.append(pfx+item+':a')
1180                    else:       #generalized for mustrain or ellipsoidal for size
1181                        Nterms = len(hapData[item][4])
1182                        if item == 'Mustrain':
1183                            names = G2spc.MustrainNames(SGData)
1184                            pwrs = []
1185                            for name in names:
1186                                h,k,l = name[1:]
1187                                pwrs.append([int(h),int(k),int(l)])
1188                            controlDict[pfx+'MuPwrs'] = pwrs
1189                        for i in range(Nterms):
1190                            sfx = ':'+str(i)
1191                            hapDict[pfx+item+sfx] = hapData[item][4][i]
1192                            if hapData[item][5][i]:
1193                                hapVary.append(pfx+item+sfx)
1194                               
1195                if Print: 
1196                    print '\n Phase: ',phase,' in histogram: ',histogram
1197                    print 135*'-'
1198                    print ' Phase fraction  : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1199                    print ' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
1200                    if hapData['Pref.Ori.'][0] == 'MD':
1201                        Ax = hapData['Pref.Ori.'][3]
1202                        print ' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
1203                            ' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
1204                    else: #'SH' for spherical harmonics
1205                        PrintSHPO(hapData['Pref.Ori.'])
1206                    PrintSize(hapData['Size'])
1207                    PrintMuStrain(hapData['Mustrain'],SGData)
1208                    PrintHStrain(hapData['HStrain'],SGData)
1209                HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
1210                refList = []
1211                for h,k,l,d in HKLd:
1212                    ext,mul,Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
1213                    mul *= 2      # for powder overlap of Friedel pairs
1214                    if ext:
1215                        continue
1216                    if 'C' in inst['Type']:
1217                        pos = 2.0*asind(wave/(2.0*d))+Zero
1218                        if limits[0] < pos < limits[1]:
1219                            refList.append([h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,0.0,Uniq,phi,0.0,{}])
1220                            #last item should contain form factor values by atom type
1221                    else:
1222                        raise ValueError 
1223                Histogram['Reflection Lists'][phase] = refList
1224            elif 'HKLF' in histogram:
1225                inst = Histogram['Instrument Parameters']
1226                inst = dict(zip(inst[2],inst[1]))
1227                hId = Histogram['hId']
1228                hfx = ':%d:'%(hId)
1229                for item in inst:
1230                    hapDict[hfx+item] = inst[item]
1231                pfx = str(pId)+':'+str(hId)+':'
1232                hapDict[pfx+'Scale'] = hapData['Scale'][0]
1233                if hapData['Scale'][1]:
1234                    hapVary.append(pfx+'Scale')
1235                extApprox,extType,extParms = hapData['Extinction']
1236                controlDict[pfx+'EType'] = extType
1237                controlDict[pfx+'EApprox'] = extApprox
1238                if 'Primary' in extType:
1239                    Ekey = ['Ep',]
1240                elif 'Secondary Type II' == extType:
1241                    Ekey = ['Es',]
1242                elif 'Secondary Type I' == extType:
1243                    Ekey = ['Eg',]
1244                else:
1245                    Ekey = ['Eg','Es']
1246                for eKey in Ekey:
1247                    hapDict[pfx+eKey] = extParms[eKey][0]
1248                    if extParms[eKey][1]:
1249                        hapVary.append(pfx+eKey)
1250                if Print: 
1251                    print '\n Phase: ',phase,' in histogram: ',histogram
1252                    print 135*'-'
1253                    print ' Scale factor     : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1254                    print ' Extinction approx: %10s'%(extApprox),' Type: %15s'%(extType),' tbar: %6.3f'%(extParms['Tbar'])
1255                    text = ' Parameters       :'
1256                    for eKey in Ekey:
1257                        text += ' %4s : %10.3g Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
1258                    print text
1259                Histogram['Reflection Lists'] = phase       
1260               
1261    return hapVary,hapDict,controlDict
1262   
1263def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,Print=True):
1264   
1265    def PrintSizeAndSig(hapData,sizeSig):
1266        line = '\n Size model:     %9s'%(hapData[0])
1267        refine = False
1268        if hapData[0] in ['isotropic','uniaxial']:
1269            line += ' equatorial:%12.3f'%(hapData[1][0])
1270            if sizeSig[0][0]:
1271                line += ', sig:%8.3f'%(sizeSig[0][0])
1272                refine = True
1273            if hapData[0] == 'uniaxial':
1274                line += ' axial:%12.3f'%(hapData[1][1])
1275                if sizeSig[0][1]:
1276                    refine = True
1277                    line += ', sig:%8.3f'%(sizeSig[0][1])
1278            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1279            if sizeSig[0][2]:
1280                refine = True
1281                line += ', sig:%8.3f'%(sizeSig[0][2])
1282            if refine:
1283                print line
1284        else:
1285            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1286            if sizeSig[0][2]:
1287                refine = True
1288                line += ', sig:%8.3f'%(sizeSig[0][2])
1289            Snames = ['S11','S22','S33','S12','S13','S23']
1290            ptlbls = ' name  :'
1291            ptstr =  ' value :'
1292            sigstr = ' sig   :'
1293            for i,name in enumerate(Snames):
1294                ptlbls += '%12s' % (name)
1295                ptstr += '%12.6f' % (hapData[4][i])
1296                if sizeSig[1][i]:
1297                    refine = True
1298                    sigstr += '%12.6f' % (sizeSig[1][i])
1299                else:
1300                    sigstr += 12*' '
1301            if refine:
1302                print line
1303                print ptlbls
1304                print ptstr
1305                print sigstr
1306       
1307    def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
1308        line = '\n Mustrain model: %9s'%(hapData[0])
1309        refine = False
1310        if hapData[0] in ['isotropic','uniaxial']:
1311            line += ' equatorial:%12.1f'%(hapData[1][0])
1312            if mustrainSig[0][0]:
1313                line += ', sig: %8.1f'%(mustrainSig[0][0])
1314                refine = True
1315            if hapData[0] == 'uniaxial':
1316                line += ' axial:%12.1f'%(hapData[1][1])
1317                if mustrainSig[0][1]:
1318                     line += ', sig:%8.1f'%(mustrainSig[0][1])
1319            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1320            if mustrainSig[0][2]:
1321                refine = True
1322                line += ', sig:%8.3f'%(mustrainSig[0][2])
1323            if refine:
1324                print line
1325        else:
1326            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1327            if mustrainSig[0][2]:
1328                refine = True
1329                line += ', sig:%8.3f'%(mustrainSig[0][2])
1330            Snames = G2spc.MustrainNames(SGData)
1331            ptlbls = ' name  :'
1332            ptstr =  ' value :'
1333            sigstr = ' sig   :'
1334            for i,name in enumerate(Snames):
1335                ptlbls += '%12s' % (name)
1336                ptstr += '%12.6f' % (hapData[4][i])
1337                if mustrainSig[1][i]:
1338                    refine = True
1339                    sigstr += '%12.6f' % (mustrainSig[1][i])
1340                else:
1341                    sigstr += 12*' '
1342            if refine:
1343                print line
1344                print ptlbls
1345                print ptstr
1346                print sigstr
1347           
1348    def PrintHStrainAndSig(hapData,strainSig,SGData):
1349        Hsnames = G2spc.HStrainNames(SGData)
1350        ptlbls = ' name  :'
1351        ptstr =  ' value :'
1352        sigstr = ' sig   :'
1353        refine = False
1354        for i,name in enumerate(Hsnames):
1355            ptlbls += '%12s' % (name)
1356            ptstr += '%12.6g' % (hapData[0][i])
1357            if name in strainSig:
1358                refine = True
1359                sigstr += '%12.6g' % (strainSig[name])
1360            else:
1361                sigstr += 12*' '
1362        if refine:
1363            print '\n Hydrostatic/elastic strain: '
1364            print ptlbls
1365            print ptstr
1366            print sigstr
1367       
1368    def PrintSHPOAndSig(hapData,POsig):
1369        print '\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
1370        ptlbls = ' names :'
1371        ptstr =  ' values:'
1372        sigstr = ' sig   :'
1373        for item in hapData[5]:
1374            ptlbls += '%12s'%(item)
1375            ptstr += '%12.3f'%(hapData[5][item])
1376            if item in POsig:
1377                sigstr += '%12.3f'%(POsig[item])
1378            else:
1379                sigstr += 12*' ' 
1380        print ptlbls
1381        print ptstr
1382        print sigstr
1383   
1384    for phase in Phases:
1385        HistoPhase = Phases[phase]['Histograms']
1386        SGData = Phases[phase]['General']['SGData']
1387        pId = Phases[phase]['pId']
1388        histoList = HistoPhase.keys()
1389        histoList.sort()
1390        for histogram in histoList:
1391            try:
1392                Histogram = Histograms[histogram]
1393            except KeyError:                       
1394                #skip if histogram not included e.g. in a sequential refinement
1395                continue
1396            print '\n Phase: ',phase,' in histogram: ',histogram
1397            print 130*'-'
1398            hapData = HistoPhase[histogram]
1399            hId = Histogram['hId']
1400            if 'PWDR' in histogram:
1401                pfx = str(pId)+':'+str(hId)+':'
1402                print ' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections'   \
1403                    %(Histogram[pfx+'Rf'],Histogram[pfx+'Rf^2'],Histogram[pfx+'Nref'])
1404               
1405                PhFrExtPOSig = {}
1406                for item in ['Scale','Extinction']:
1407                    hapData[item][0] = parmDict[pfx+item]
1408                    if pfx+item in sigDict:
1409                        PhFrExtPOSig[item] = sigDict[pfx+item]
1410                if hapData['Pref.Ori.'][0] == 'MD':
1411                    hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
1412                    if pfx+'MD' in sigDict:
1413                        PhFrExtPOSig['MD'] = sigDict[pfx+'MD']
1414                else:                           #'SH' spherical harmonics
1415                    for item in hapData['Pref.Ori.'][5]:
1416                        hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
1417                        if pfx+item in sigDict:
1418                            PhFrExtPOSig[item] = sigDict[pfx+item]
1419                if Print:
1420                    if 'Scale' in PhFrExtPOSig:
1421                        print ' Phase fraction  : %10.4f, sig %10.4f'%(hapData['Scale'][0],PhFrExtPOSig['Scale'])
1422                    if 'Extinction' in PhFrExtPOSig:
1423                        print ' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig['Extinction'])
1424                    if hapData['Pref.Ori.'][0] == 'MD':
1425                        if 'MD' in PhFrExtPOSig:
1426                            print ' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig['MD'])
1427                    else:
1428                        PrintSHPOAndSig(hapData['Pref.Ori.'],PhFrExtPOSig)
1429                SizeMuStrSig = {'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
1430                    'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
1431                    'HStrain':{}}                 
1432                for item in ['Mustrain','Size']:
1433                    hapData[item][1][2] = parmDict[pfx+item+':mx']
1434                    hapData[item][1][2] = min(1.,max(0.1,hapData[item][1][2]))
1435                    if pfx+item+':mx' in sigDict:
1436                        SizeMuStrSig[item][0][2] = sigDict[pfx+item+':mx']
1437                    if hapData[item][0] in ['isotropic','uniaxial']:                   
1438                        hapData[item][1][0] = parmDict[pfx+item+':i']
1439                        if item == 'Size':
1440                            hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
1441                        if pfx+item+':i' in sigDict: 
1442                            SizeMuStrSig[item][0][0] = sigDict[pfx+item+':i']
1443                        if hapData[item][0] == 'uniaxial':
1444                            hapData[item][1][1] = parmDict[pfx+item+':a']
1445                            if item == 'Size':
1446                                hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))                       
1447                            if pfx+item+':a' in sigDict:
1448                                SizeMuStrSig[item][0][1] = sigDict[pfx+item+':a']
1449                    else:       #generalized for mustrain or ellipsoidal for size
1450                        Nterms = len(hapData[item][4])
1451                        for i in range(Nterms):
1452                            sfx = ':'+str(i)
1453                            hapData[item][4][i] = parmDict[pfx+item+sfx]
1454                            if pfx+item+sfx in sigDict:
1455                                SizeMuStrSig[item][1][i] = sigDict[pfx+item+sfx]
1456                names = G2spc.HStrainNames(SGData)
1457                for i,name in enumerate(names):
1458                    hapData['HStrain'][0][i] = parmDict[pfx+name]
1459                    if pfx+name in sigDict:
1460                        SizeMuStrSig['HStrain'][name] = sigDict[pfx+name]
1461                if Print:
1462                    PrintSizeAndSig(hapData['Size'],SizeMuStrSig['Size'])
1463                    PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig['Mustrain'],SGData)
1464                    PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig['HStrain'],SGData)
1465               
1466            elif 'HKLF' in histogram:
1467                pfx = str(pId)+':'+str(hId)+':'
1468                print ' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections'   \
1469                    %(Histogram[pfx+'Rf'],Histogram[pfx+'Rf^2'],Histogram[pfx+'Nref'])
1470                ScalExtSig = {}
1471                for item in ['Scale','Ep','Eg','Es']:
1472                    if parmDict.get(pfx+item):
1473                        hapData[item][0] = parmDict[pfx+item]
1474                        if pfx+item in sigDict:
1475                            ScalExtSig[item] = sigDict[pfx+item]
1476                if Print: 
1477                    if 'Scale' in ScalExtSig:
1478                        print ' Scale factor : %10.4f, sig %10.4f'%(hapData['Scale'][0],ScalExtSig['Scale'])
1479# fix after it runs!               
1480#                    print '\n Phase: ',phase,' in histogram: ',histogram
1481#                    print 135*'-'
1482#                    print ' Scale factor     : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1483#                    print ' Extinction approx: %10s'%(extApprox),' Type: %15s'%(extType),' tbar: %6.3f'%(extParms['Tbar'])
1484#                    text = ' Parameters       :'
1485#                    for eKey in Ekey:
1486#                        text += ' %4s : %10.3g Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
1487#                    print text
1488   
1489################################################################################
1490##### Histogram data
1491################################################################################       
1492                   
1493def GetHistogramData(Histograms,Print=True):
1494   
1495    def GetBackgroundParms(hId,Background):
1496        Back = Background[0]
1497        DebyePeaks = Background[1]
1498        bakType,bakFlag = Back[:2]
1499        backVals = Back[3:]
1500        backNames = [':'+str(hId)+':Back:'+str(i) for i in range(len(backVals))]
1501        backDict = dict(zip(backNames,backVals))
1502        backVary = []
1503        if bakFlag:
1504            backVary = backNames
1505        backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
1506        backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
1507        debyeDict = {}
1508        debyeList = []
1509        for i in range(DebyePeaks['nDebye']):
1510            debyeNames = [':'+str(hId)+':DebyeA:'+str(i),':'+str(hId)+':DebyeR:'+str(i),':'+str(hId)+':DebyeU:'+str(i)]
1511            debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
1512            debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
1513        debyeVary = []
1514        for item in debyeList:
1515            if item[1]:
1516                debyeVary.append(item[0])
1517        backDict.update(debyeDict)
1518        backVary += debyeVary
1519        peakDict = {}
1520        peakList = []
1521        for i in range(DebyePeaks['nPeaks']):
1522            peakNames = [':'+str(hId)+':BkPkpos:'+str(i),':'+str(hId)+ \
1523                ':BkPkint:'+str(i),':'+str(hId)+':BkPksig:'+str(i),':'+str(hId)+':BkPkgam:'+str(i)]
1524            peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
1525            peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
1526        peakVary = []
1527        for item in peakList:
1528            if item[1]:
1529                peakVary.append(item[0])
1530        backDict.update(peakDict)
1531        backVary += peakVary
1532        return bakType,backDict,backVary           
1533       
1534    def GetInstParms(hId,Inst):
1535        insVals,insFlags,insNames = Inst[1:4]
1536        dataType = insVals[0]
1537        instDict = {}
1538        insVary = []
1539        pfx = ':'+str(hId)+':'
1540        for i,flag in enumerate(insFlags):
1541            insName = pfx+insNames[i]
1542            instDict[insName] = insVals[i]
1543            if flag:
1544                insVary.append(insName)
1545        instDict[pfx+'X'] = max(instDict[pfx+'X'],0.001)
1546        instDict[pfx+'Y'] = max(instDict[pfx+'Y'],0.001)
1547        instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
1548        return dataType,instDict,insVary
1549       
1550    def GetSampleParms(hId,Sample):
1551        sampVary = []
1552        hfx = ':'+str(hId)+':'       
1553        sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
1554            hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
1555        Type = Sample['Type']
1556        if 'Bragg' in Type:             #Bragg-Brentano
1557            for item in ['Scale','Shift','Transparency']:       #surface roughness?, diffuse scattering?
1558                sampDict[hfx+item] = Sample[item][0]
1559                if Sample[item][1]:
1560                    sampVary.append(hfx+item)
1561        elif 'Debye' in Type:        #Debye-Scherrer
1562            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
1563                sampDict[hfx+item] = Sample[item][0]
1564                if Sample[item][1]:
1565                    sampVary.append(hfx+item)
1566        return Type,sampDict,sampVary
1567       
1568    def PrintBackground(Background):
1569        Back = Background[0]
1570        DebyePeaks = Background[1]
1571        print '\n Background function: ',Back[0],' Refine?',bool(Back[1])
1572        line = ' Coefficients: '
1573        for i,back in enumerate(Back[3:]):
1574            line += '%10.3f'%(back)
1575            if i and not i%10:
1576                line += '\n'+15*' '
1577        print line
1578        if DebyePeaks['nDebye']:
1579            print '\n Debye diffuse scattering coefficients'
1580            parms = ['DebyeA','DebyeR','DebyeU']
1581            line = ' names :  '
1582            for parm in parms:
1583                line += '%8s refine?'%(parm)
1584            print line
1585            for j,term in enumerate(DebyePeaks['debyeTerms']):
1586                line = ' term'+'%2d'%(j)+':'
1587                for i in range(3):
1588                    line += '%10.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
1589                print line
1590        if DebyePeaks['nPeaks']:
1591            print '\n Single peak coefficients'
1592            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
1593            line = ' names :  '
1594            for parm in parms:
1595                line += '%8s refine?'%(parm)
1596            print line
1597            for j,term in enumerate(DebyePeaks['peaksList']):
1598                line = ' peak'+'%2d'%(j)+':'
1599                for i in range(4):
1600                    line += '%10.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
1601                print line
1602       
1603    def PrintInstParms(Inst):
1604        print '\n Instrument Parameters:'
1605        ptlbls = ' name  :'
1606        ptstr =  ' value :'
1607        varstr = ' refine:'
1608        instNames = Inst[3][1:]
1609        for i,name in enumerate(instNames):
1610            ptlbls += '%12s' % (name)
1611            ptstr += '%12.6f' % (Inst[1][i+1])
1612            if name in ['Lam1','Lam2','Azimuth']:
1613                varstr += 12*' '
1614            else:
1615                varstr += '%12s' % (str(bool(Inst[2][i+1])))
1616        print ptlbls
1617        print ptstr
1618        print varstr
1619       
1620    def PrintSampleParms(Sample):
1621        print '\n Sample Parameters:'
1622        print ' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \
1623            (Sample['Omega'],Sample['Chi'],Sample['Phi'])
1624        ptlbls = ' name  :'
1625        ptstr =  ' value :'
1626        varstr = ' refine:'
1627        if 'Bragg' in Sample['Type']:
1628            for item in ['Scale','Shift','Transparency']:
1629                ptlbls += '%14s'%(item)
1630                ptstr += '%14.4f'%(Sample[item][0])
1631                varstr += '%14s'%(str(bool(Sample[item][1])))
1632           
1633        elif 'Debye' in Type:        #Debye-Scherrer
1634            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
1635                ptlbls += '%14s'%(item)
1636                ptstr += '%14.4f'%(Sample[item][0])
1637                varstr += '%14s'%(str(bool(Sample[item][1])))
1638
1639        print ptlbls
1640        print ptstr
1641        print varstr
1642       
1643
1644    histDict = {}
1645    histVary = []
1646    controlDict = {}
1647    histoList = Histograms.keys()
1648    histoList.sort()
1649    for histogram in histoList:
1650        if 'PWDR' in histogram:
1651            Histogram = Histograms[histogram]
1652            hId = Histogram['hId']
1653            pfx = ':'+str(hId)+':'
1654            controlDict[pfx+'Limits'] = Histogram['Limits'][1]
1655           
1656            Background = Histogram['Background']
1657            Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
1658            controlDict[pfx+'bakType'] = Type
1659            histDict.update(bakDict)
1660            histVary += bakVary
1661           
1662            Inst = Histogram['Instrument Parameters']
1663            Type,instDict,insVary = GetInstParms(hId,Inst)
1664            controlDict[pfx+'histType'] = Type
1665            if pfx+'Lam1' in instDict:
1666                controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
1667            else:
1668                controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']           
1669            histDict.update(instDict)
1670            histVary += insVary
1671           
1672            Sample = Histogram['Sample Parameters']
1673            Type,sampDict,sampVary = GetSampleParms(hId,Sample)
1674            controlDict[pfx+'instType'] = Type
1675            histDict.update(sampDict)
1676            histVary += sampVary
1677   
1678            if Print: 
1679                print '\n Histogram: ',histogram,' histogram Id: ',hId
1680                print 135*'-'
1681                Units = {'C':' deg','T':' msec'}
1682                units = Units[controlDict[pfx+'histType'][2]]
1683                Limits = controlDict[pfx+'Limits']
1684                print ' Instrument type: ',Sample['Type']
1685                print ' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)     
1686                PrintSampleParms(Sample)
1687                PrintInstParms(Inst)
1688                PrintBackground(Background)
1689        elif 'HKLF' in histogram:
1690            Histogram = Histograms[histogram]
1691            hId = Histogram['hId']
1692            pfx = ':'+str(hId)+':'
1693            Inst = Histogram['Instrument Parameters']
1694            controlDict[pfx+'histType'] = Inst[1][0]
1695            histDict[pfx+'Lam'] = Inst[1][1]
1696            controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']                   
1697    return histVary,histDict,controlDict
1698   
1699def SetHistogramData(parmDict,sigDict,Histograms,Print=True):
1700   
1701    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
1702        Back = Background[0]
1703        DebyePeaks = Background[1]
1704        lenBack = len(Back[3:])
1705        backSig = [0 for i in range(lenBack+3*DebyePeaks['nDebye']+4*DebyePeaks['nPeaks'])]
1706        for i in range(lenBack):
1707            Back[3+i] = parmDict[pfx+'Back:'+str(i)]
1708            if pfx+'Back:'+str(i) in sigDict:
1709                backSig[i] = sigDict[pfx+'Back:'+str(i)]
1710        if DebyePeaks['nDebye']:
1711            for i in range(DebyePeaks['nDebye']):
1712                names = [pfx+'DebyeA:'+str(i),pfx+'DebyeR:'+str(i),pfx+'DebyeU:'+str(i)]
1713                for j,name in enumerate(names):
1714                    DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
1715                    if name in sigDict:
1716                        backSig[lenBack+3*i+j] = sigDict[name]           
1717        if DebyePeaks['nPeaks']:
1718            for i in range(DebyePeaks['nPeaks']):
1719                names = [pfx+'BkPkpos:'+str(i),pfx+'BkPkint:'+str(i),
1720                    pfx+'BkPksig:'+str(i),pfx+'BkPkgam:'+str(i)]
1721                for j,name in enumerate(names):
1722                    DebyePeaks['peaksList'][i][2*j] = parmDict[name]
1723                    if name in sigDict:
1724                        backSig[lenBack+3*DebyePeaks['nDebye']+4*i+j] = sigDict[name]
1725        return backSig
1726       
1727    def SetInstParms(pfx,Inst,parmDict,sigDict):
1728        insVals,insFlags,insNames = Inst[1:4]
1729        instSig = [0 for i in range(len(insVals))]
1730        for i,flag in enumerate(insFlags):
1731            insName = pfx+insNames[i]
1732            insVals[i] = parmDict[insName]
1733            if insName in sigDict:
1734                instSig[i] = sigDict[insName]
1735        return instSig
1736       
1737    def SetSampleParms(pfx,Sample,parmDict,sigDict):
1738        if 'Bragg' in Sample['Type']:             #Bragg-Brentano
1739            sampSig = [0 for i in range(3)]
1740            for i,item in enumerate(['Scale','Shift','Transparency']):       #surface roughness?, diffuse scattering?
1741                Sample[item][0] = parmDict[pfx+item]
1742                if pfx+item in sigDict:
1743                    sampSig[i] = sigDict[pfx+item]
1744        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
1745            sampSig = [0 for i in range(4)]
1746            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
1747                Sample[item][0] = parmDict[pfx+item]
1748                if pfx+item in sigDict:
1749                    sampSig[i] = sigDict[pfx+item]
1750        return sampSig
1751       
1752    def PrintBackgroundSig(Background,backSig):
1753        Back = Background[0]
1754        DebyePeaks = Background[1]
1755        lenBack = len(Back[3:])
1756        valstr = ' value : '
1757        sigstr = ' sig   : '
1758        refine = False
1759        for i,back in enumerate(Back[3:]):
1760            valstr += '%10.4g'%(back)
1761            if Back[1]:
1762                refine = True
1763                sigstr += '%10.4g'%(backSig[i])
1764            else:
1765                sigstr += 10*' '
1766        if refine:
1767            print '\n Background function: ',Back[0]
1768            print valstr
1769            print sigstr
1770        if DebyePeaks['nDebye']:
1771            ifAny = False
1772            ptfmt = "%12.3f"
1773            names =  ' names :'
1774            ptstr =  ' values:'
1775            sigstr = ' esds  :'
1776            for item in sigDict:
1777                if 'Debye' in item:
1778                    ifAny = True
1779                    names += '%12s'%(item)
1780                    ptstr += ptfmt%(parmDict[item])
1781                    sigstr += ptfmt%(sigDict[item])
1782            if ifAny:
1783                print '\n Debye diffuse scattering coefficients'
1784                print names
1785                print ptstr
1786                print sigstr
1787        if DebyePeaks['nPeaks']:
1788            ifAny = False
1789            ptfmt = "%14.3f"
1790            names =  ' names :'
1791            ptstr =  ' values:'
1792            sigstr = ' esds  :'
1793            for item in sigDict:
1794                if 'BkPk' in item:
1795                    ifAny = True
1796                    names += '%14s'%(item)
1797                    ptstr += ptfmt%(parmDict[item])
1798                    sigstr += ptfmt%(sigDict[item])
1799            if ifAny:
1800                print '\n Single peak coefficients'
1801                print names
1802                print ptstr
1803                print sigstr
1804       
1805    def PrintInstParmsSig(Inst,instSig):
1806        ptlbls = ' names :'
1807        ptstr =  ' value :'
1808        sigstr = ' sig   :'
1809        instNames = Inst[3][1:]
1810        refine = False
1811        for i,name in enumerate(instNames):
1812            ptlbls += '%12s' % (name)
1813            ptstr += '%12.6f' % (Inst[1][i+1])
1814            if instSig[i+1]:
1815                refine = True
1816                sigstr += '%12.6f' % (instSig[i+1])
1817            else:
1818                sigstr += 12*' '
1819        if refine:
1820            print '\n Instrument Parameters:'
1821            print ptlbls
1822            print ptstr
1823            print sigstr
1824       
1825    def PrintSampleParmsSig(Sample,sampleSig):
1826        ptlbls = ' names :'
1827        ptstr =  ' values:'
1828        sigstr = ' sig   :'
1829        refine = False
1830        if 'Bragg' in Sample['Type']:
1831            for i,item in enumerate(['Scale','Shift','Transparency']):
1832                ptlbls += '%14s'%(item)
1833                ptstr += '%14.4f'%(Sample[item][0])
1834                if sampleSig[i]:
1835                    refine = True
1836                    sigstr += '%14.4f'%(sampleSig[i])
1837                else:
1838                    sigstr += 14*' '
1839           
1840        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
1841            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
1842                ptlbls += '%14s'%(item)
1843                ptstr += '%14.4f'%(Sample[item][0])
1844                if sampleSig[i]:
1845                    refine = True
1846                    sigstr += '%14.4f'%(sampleSig[i])
1847                else:
1848                    sigstr += 14*' '
1849
1850        if refine:
1851            print '\n Sample Parameters:'
1852            print ptlbls
1853            print ptstr
1854            print sigstr
1855       
1856    histoList = Histograms.keys()
1857    histoList.sort()
1858    for histogram in histoList:
1859        if 'PWDR' in histogram:
1860            Histogram = Histograms[histogram]
1861            hId = Histogram['hId']
1862            pfx = ':'+str(hId)+':'
1863            Background = Histogram['Background']
1864            backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
1865           
1866            Inst = Histogram['Instrument Parameters']
1867            instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
1868       
1869            Sample = Histogram['Sample Parameters']
1870            sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
1871
1872            print '\n Histogram: ',histogram,' histogram Id: ',hId
1873            print 135*'-'
1874            print ' Final refinement wR = %.2f%% on %d observations in this histogram'%(Histogram['wR'],Histogram['Nobs'])
1875            if Print:
1876                print ' Instrument type: ',Sample['Type']
1877                PrintSampleParmsSig(Sample,sampSig)
1878                PrintInstParmsSig(Inst,instSig)
1879                PrintBackgroundSig(Background,backSig)
1880               
1881################################################################################
1882##### Penalty & restraint functions
1883################################################################################
1884
1885def penaltyFxn(parmDict,varyList):
1886    pFxn = np.zeros(len(varyList))
1887    for i,item in enumerate(varyList):
1888        if 'PWLref' in item and parmDict[item] < 0.:
1889            pFxn[i] = -parmDict[item]**2        #checked OK
1890    return pFxn
1891   
1892def penaltyDeriv(parmDict,varyList):
1893    pDerv = np.zeros(len(varyList))
1894    for i,item in enumerate(varyList):
1895        if 'PWLref' in item and parmDict[item] < 0.:
1896            pDerv[i] += 2.*parmDict[item]
1897    return pDerv/100.
1898
1899################################################################################
1900##### Function & derivative calculations
1901################################################################################       
1902                   
1903def GetAtomFXU(pfx,calcControls,parmDict):
1904    Natoms = calcControls['Natoms'][pfx]
1905    Tdata = Natoms*[' ',]
1906    Mdata = np.zeros(Natoms)
1907    IAdata = Natoms*[' ',]
1908    Fdata = np.zeros(Natoms)
1909    FFdata = []
1910    BLdata = []
1911    Xdata = np.zeros((3,Natoms))
1912    dXdata = np.zeros((3,Natoms))
1913    Uisodata = np.zeros(Natoms)
1914    Uijdata = np.zeros((6,Natoms))
1915    keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata,
1916        'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2],
1917        'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata,
1918        'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2],
1919        'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5]}
1920    for iatm in range(Natoms):
1921        for key in keys:
1922            parm = pfx+key+str(iatm)
1923            if parm in parmDict:
1924                keys[key][iatm] = parmDict[parm]
1925    return Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata
1926   
1927def getFFvalues(FFtables,SQ):
1928    FFvals = {}
1929    for El in FFtables:
1930        FFvals[El] = G2el.ScatFac(FFtables[El],SQ)[0]
1931    return FFvals
1932   
1933def getBLvalues(BLtables):
1934    BLvals = {}
1935    for El in BLtables:
1936        BLvals[El] = BLtables[El][1][1]
1937    return BLvals
1938       
1939def StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict):
1940    ''' Compute structure factors for all h,k,l for phase
1941    input:
1942        refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv]
1943        G:      reciprocal metric tensor
1944        pfx:    phase id string
1945        SGData: space group info. dictionary output from SpcGroup
1946        calcControls:
1947        ParmDict:
1948    puts result F^2 in each ref[8] in refList
1949    '''       
1950    twopi = 2.0*np.pi
1951    twopisq = 2.0*np.pi**2
1952    ast = np.sqrt(np.diag(G))
1953    Mast = twopisq*np.multiply.outer(ast,ast)
1954    FFtables = calcControls['FFtables']
1955    BLtables = calcControls['BLtables']
1956    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
1957    FF = np.zeros(len(Tdata))
1958    if 'N' in parmDict[hfx+'Type']:
1959        FP,FPP = G2el.BlenRes(Tdata,BLtables,parmDict[hfx+'Lam'])
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    for refl in refList:
1967        fbs = np.array([0,0])
1968        H = refl[:3]
1969        SQ = 1./(2.*refl[4])**2
1970        if not len(refl[-1]):                #no form factors
1971            if 'N' in parmDict[hfx+'Type']:
1972                refl[-1] = getBLvalues(BLtables)
1973            else:       #'X'
1974                refl[-1] = getFFvalues(FFtables,SQ)
1975        for i,El in enumerate(Tdata):
1976            FF[i] = refl[-1][El]           
1977        SQfactor = 4.0*SQ*twopisq
1978        Uniq = refl[11]
1979        phi = refl[12]
1980        phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+phi[:,np.newaxis])
1981        sinp = np.sin(phase)
1982        cosp = np.cos(phase)
1983        occ = Mdata*Fdata/len(Uniq)
1984        biso = -SQfactor*Uisodata
1985        Tiso = np.where(biso<1.,np.exp(biso),1.0)
1986        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
1987        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1988        Tcorr = Tiso*Tuij
1989        fa = np.array([(FF+FP)*occ*cosp*Tcorr,-FPP*occ*sinp*Tcorr])
1990        fas = np.sum(np.sum(fa,axis=1),axis=1)        #real
1991        if not SGData['SGInv']:
1992            fb = np.array([(FF+FP)*occ*sinp*Tcorr,FPP*occ*cosp*Tcorr])
1993            fbs = np.sum(np.sum(fb,axis=1),axis=1)
1994        fasq = fas**2
1995        fbsq = fbs**2        #imaginary
1996        refl[9] = np.sum(fasq)+np.sum(fbsq)
1997        refl[10] = atan2d(fbs[0],fas[0])
1998    return refList
1999   
2000def StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict):
2001    twopi = 2.0*np.pi
2002    twopisq = 2.0*np.pi**2
2003    ast = np.sqrt(np.diag(G))
2004    Mast = twopisq*np.multiply.outer(ast,ast)
2005    FFtables = calcControls['FFtables']
2006    BLtables = calcControls['BLtables']
2007    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
2008    FF = np.zeros(len(Tdata))
2009    if 'N' in parmDict[hfx+'Type']:
2010        FP = 0.
2011        FPP = 0.
2012    else:
2013        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
2014        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
2015    maxPos = len(SGData['SGOps'])       
2016    Uij = np.array(G2lat.U6toUij(Uijdata))
2017    bij = Mast*Uij.T
2018    dFdvDict = {}
2019    dFdfr = np.zeros((len(refList),len(Mdata)))
2020    dFdx = np.zeros((len(refList),len(Mdata),3))
2021    dFdui = np.zeros((len(refList),len(Mdata)))
2022    dFdua = np.zeros((len(refList),len(Mdata),6))
2023    for iref,refl in enumerate(refList):
2024        H = np.array(refl[:3])
2025        SQ = 1./(2.*refl[4])**2             # or (sin(theta)/lambda)**2
2026        for i,El in enumerate(Tdata):           
2027            FF[i] = refl[-1][El]           
2028        SQfactor = 8.0*SQ*np.pi**2
2029        Uniq = refl[11]
2030        phi = refl[12]
2031        phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+phi[np.newaxis,:])
2032        sinp = np.sin(phase)
2033        cosp = np.cos(phase)
2034        occ = Mdata*Fdata/len(Uniq)
2035        biso = -SQfactor*Uisodata
2036        Tiso = np.where(biso<1.,np.exp(biso),1.0)
2037        HbH = -np.inner(H,np.inner(bij,H))
2038        Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
2039        Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
2040        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
2041        Tcorr = Tiso*Tuij
2042        fot = (FF+FP)*occ*Tcorr
2043        fotp = FPP*occ*Tcorr
2044        fa = np.array([fot[:,np.newaxis]*cosp,fotp[:,np.newaxis]*cosp])       #non positions
2045        fb = np.array([fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])
2046       
2047        fas = np.sum(np.sum(fa,axis=1),axis=1)
2048        fbs = np.sum(np.sum(fb,axis=1),axis=1)
2049        fax = np.array([-fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])   #positions
2050        fbx = np.array([fot[:,np.newaxis]*cosp,-fot[:,np.newaxis]*cosp])
2051        #sum below is over Uniq
2052        dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)
2053        dfadx = np.sum(twopi*Uniq*fax[:,:,:,np.newaxis],axis=2)
2054        dfadui = np.sum(-SQfactor*fa,axis=2)
2055        dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
2056        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for     
2057        dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)
2058        dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])
2059        dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])
2060        dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])
2061        if not SGData['SGInv']:
2062            dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2)        #problem here if occ=0 for some atom
2063            dfbdx = np.sum(twopi*Uniq*fbx[:,:,:,np.newaxis],axis=2)         
2064            dfbdui = np.sum(-SQfactor*fb,axis=2)
2065            dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
2066            dFdfr[iref] += 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
2067            dFdx[iref] += 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
2068            dFdui[iref] += 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
2069            dFdua[iref] += 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
2070        #loop over atoms - each dict entry is list of derivatives for all the reflections
2071    for i in range(len(Mdata)):     
2072        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
2073        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
2074        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
2075        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
2076        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
2077        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
2078        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
2079        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
2080        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
2081        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
2082        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
2083    return dFdvDict
2084       
2085def Dict2Values(parmdict, varylist):
2086    '''Use before call to leastsq to setup list of values for the parameters
2087    in parmdict, as selected by key in varylist'''
2088    return [parmdict[key] for key in varylist] 
2089   
2090def Values2Dict(parmdict, varylist, values):
2091    ''' Use after call to leastsq to update the parameter dictionary with
2092    values corresponding to keys in varylist'''
2093    parmdict.update(zip(varylist,values))
2094   
2095def GetNewCellParms(parmDict,varyList):
2096    newCellDict = {}
2097    Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],['A'+str(i) for i in range(6)]))
2098    for item in varyList:
2099        keys = item.split(':')
2100        if keys[2] in Ddict:
2101            key = keys[0]+'::'+Ddict[keys[2]]
2102            parm = keys[0]+'::'+keys[2]
2103            newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
2104    return newCellDict
2105   
2106def ApplyXYZshifts(parmDict,varyList):
2107    ''' takes atom x,y,z shift and applies it to corresponding atom x,y,z value
2108        input:
2109            parmDict - parameter dictionary
2110            varyList - list of variables
2111        returns:
2112            newAtomDict - dictitemionary of new atomic coordinate names & values;
2113                key is parameter shift name
2114    '''
2115    newAtomDict = {}
2116    for item in parmDict:
2117        if 'dA' in item:
2118            parm = ''.join(item.split('d'))
2119            parmDict[parm] += parmDict[item]
2120            newAtomDict[item] = [parm,parmDict[parm]]
2121    return newAtomDict
2122   
2123def SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict):
2124    IFCoup = 'Bragg' in calcControls[hfx+'instType']
2125    odfCor = 1.0
2126    H = refl[:3]
2127    cell = G2lat.Gmat2cell(g)
2128    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
2129    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
2130    phi,beta = G2lat.CrsAng(H,cell,SGData)
2131    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
2132    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
2133    for item in SHnames:
2134        L,M,N = eval(item.strip('C'))
2135        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2136        Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
2137        Lnorm = G2lat.Lnorm(L)
2138        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
2139    return odfCor
2140   
2141def SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict):
2142    FORPI = 12.5663706143592
2143    IFCoup = 'Bragg' in calcControls[hfx+'instType']
2144    odfCor = 1.0
2145    dFdODF = {}
2146    dFdSA = [0,0,0]
2147    H = refl[:3]
2148    cell = G2lat.Gmat2cell(g)
2149    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
2150    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
2151    phi,beta = G2lat.CrsAng(H,cell,SGData)
2152    psi,gam,dPSdA,dGMdA = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup)
2153    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
2154    for item in SHnames:
2155        L,M,N = eval(item.strip('C'))
2156        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2157        Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
2158        Lnorm = G2lat.Lnorm(L)
2159        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
2160        dFdODF[pfx+item] = Lnorm*Kcl*Ksl
2161        for i in range(3):
2162            dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i])
2163    return odfCor,dFdODF,dFdSA
2164   
2165def SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict):
2166    odfCor = 1.0
2167    H = refl[:3]
2168    cell = G2lat.Gmat2cell(g)
2169    Sangl = [0.,0.,0.]
2170    if 'Bragg' in calcControls[hfx+'instType']:
2171        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
2172        IFCoup = True
2173    else:
2174        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
2175        IFCoup = False
2176    phi,beta = G2lat.CrsAng(H,cell,SGData)
2177    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
2178    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
2179    for item in SHnames:
2180        L,N = eval(item.strip('C'))
2181        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) 
2182        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
2183    return odfCor
2184   
2185def SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict):
2186    FORPI = 12.5663706143592
2187    odfCor = 1.0
2188    dFdODF = {}
2189    H = refl[:3]
2190    cell = G2lat.Gmat2cell(g)
2191    Sangl = [0.,0.,0.]
2192    if 'Bragg' in calcControls[hfx+'instType']:
2193        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
2194        IFCoup = True
2195    else:
2196        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
2197        IFCoup = False
2198    phi,beta = G2lat.CrsAng(H,cell,SGData)
2199    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
2200    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
2201    for item in SHnames:
2202        L,N = eval(item.strip('C'))
2203        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) 
2204        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
2205        dFdODF[phfx+item] = Kcsl*Lnorm
2206    return odfCor,dFdODF
2207   
2208def GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
2209    POcorr = 1.0
2210    if calcControls[phfx+'poType'] == 'MD':
2211        MD = parmDict[phfx+'MD']
2212        if MD != 1.0:
2213            MDAxis = calcControls[phfx+'MDAxis']
2214            sumMD = 0
2215            for H in refl[11]:           
2216                cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
2217                A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
2218                sumMD += A**3
2219            POcorr = sumMD/len(refl[11])
2220    else:   #spherical harmonics
2221        if calcControls[pfx+'SHord']:
2222            POcorr = SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict)
2223    return POcorr
2224   
2225def GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
2226    POcorr = 1.0
2227    POderv = {}
2228    if calcControls[phfx+'poType'] == 'MD':
2229        MD = parmDict[phfx+'MD']
2230        MDAxis = calcControls[phfx+'MDAxis']
2231        sumMD = 0
2232        sumdMD = 0
2233        for H in refl[11]:           
2234            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
2235            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
2236            sumMD += A**3
2237            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
2238        POcorr = sumMD/len(refl[11])
2239        POderv[phfx+'MD'] = sumdMD/len(refl[11])
2240    else:   #spherical harmonics
2241        if calcControls[pfx+'SHord']:
2242            POcorr,POderv = SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict)
2243    return POcorr,POderv
2244   
2245def GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
2246    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3]               #scale*multiplicity
2247    if 'X' in parmDict[hfx+'Type']:
2248        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
2249    Icorr *= GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
2250    if pfx+'SHorder' in parmDict:
2251        Icorr *= SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict)
2252    refl[13] = Icorr       
2253   
2254def GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
2255    dIdsh = 1./parmDict[hfx+'Scale']
2256    dIdsp = 1./parmDict[phfx+'Scale']
2257    if 'X' in parmDict[hfx+'Type']:
2258        pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
2259        dIdPola /= pola
2260    else:       #'N'
2261        dIdPola = 0.0
2262    POcorr,dIdPO = GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
2263    for iPO in dIdPO:
2264        dIdPO[iPO] /= POcorr
2265    dFdODF = {}
2266    dFdSA = [0,0,0]
2267    if pfx+'SHorder' in parmDict:
2268        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict)
2269        for iSH in dFdODF:
2270            dFdODF[iSH] /= odfCor
2271        for i in range(3):
2272            dFdSA[i] /= odfCor
2273    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA
2274       
2275def GetSampleSigGam(refl,wave,G,GB,phfx,calcControls,parmDict):
2276    costh = cosd(refl[5]/2.)
2277    #crystallite size
2278    if calcControls[phfx+'SizeType'] == 'isotropic':
2279        Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size:i']*costh)
2280    elif calcControls[phfx+'SizeType'] == 'uniaxial':
2281        H = np.array(refl[:3])
2282        P = np.array(calcControls[phfx+'SizeAxis'])
2283        cosP,sinP = G2lat.CosSinAngle(H,P,G)
2284        Sgam = (1.8*wave/np.pi)/(parmDict[phfx+'Size:i']*parmDict[phfx+'Size:a']*costh)
2285        Sgam *= np.sqrt((sinP*parmDict[phfx+'Size:a'])**2+(cosP*parmDict[phfx+'Size:i'])**2)
2286    else:           #ellipsoidal crystallites
2287        Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
2288        H = np.array(refl[:3])
2289        lenR = G2pwd.ellipseSize(H,Sij,GB)
2290        Sgam = 1.8*wave/(np.pi*costh*lenR)
2291    #microstrain               
2292    if calcControls[phfx+'MustrainType'] == 'isotropic':
2293        Mgam = 0.018*parmDict[phfx+'Mustrain:i']*tand(refl[5]/2.)/np.pi
2294    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2295        H = np.array(refl[:3])
2296        P = np.array(calcControls[phfx+'MustrainAxis'])
2297        cosP,sinP = G2lat.CosSinAngle(H,P,G)
2298        Si = parmDict[phfx+'Mustrain:i']
2299        Sa = parmDict[phfx+'Mustrain:a']
2300        Mgam = 0.018*Si*Sa*tand(refl[5]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
2301    else:       #generalized - P.W. Stephens model
2302        pwrs = calcControls[phfx+'MuPwrs']
2303        sum = 0
2304        for i,pwr in enumerate(pwrs):
2305            sum += parmDict[phfx+'Mustrain:'+str(i)]*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
2306        Mgam = 0.018*refl[4]**2*tand(refl[5]/2.)*sum
2307    gam = Sgam*parmDict[phfx+'Size:mx']+Mgam*parmDict[phfx+'Mustrain:mx']
2308    sig = (Sgam*(1.-parmDict[phfx+'Size:mx']))**2+(Mgam*(1.-parmDict[phfx+'Mustrain:mx']))**2
2309    sig /= ateln2
2310    return sig,gam
2311       
2312def GetSampleSigGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict):
2313    gamDict = {}
2314    sigDict = {}
2315    costh = cosd(refl[5]/2.)
2316    tanth = tand(refl[5]/2.)
2317    #crystallite size derivatives
2318    if calcControls[phfx+'SizeType'] == 'isotropic':
2319        Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size:i']*costh)
2320        gamDict[phfx+'Size:i'] = -1.80*wave*parmDict[phfx+'Size:mx']/(np.pi*costh)
2321        sigDict[phfx+'Size:i'] = -3.60*Sgam*wave*(1.-parmDict[phfx+'Size:mx'])**2/(np.pi*costh*ateln2)
2322    elif calcControls[phfx+'SizeType'] == 'uniaxial':
2323        H = np.array(refl[:3])
2324        P = np.array(calcControls[phfx+'SizeAxis'])
2325        cosP,sinP = G2lat.CosSinAngle(H,P,G)
2326        Si = parmDict[phfx+'Size:i']
2327        Sa = parmDict[phfx+'Size:a']
2328        gami = (1.8*wave/np.pi)/(Si*Sa)
2329        sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
2330        Sgam = gami*sqtrm
2331        gam = Sgam/costh
2332        dsi = (gami*Si*cosP**2/(sqtrm*costh)-gam/Si)
2333        dsa = (gami*Sa*sinP**2/(sqtrm*costh)-gam/Sa)
2334        gamDict[phfx+'Size:i'] = dsi*parmDict[phfx+'Size:mx']
2335        gamDict[phfx+'Size:a'] = dsa*parmDict[phfx+'Size:mx']
2336        sigDict[phfx+'Size:i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size:mx'])**2/ateln2
2337        sigDict[phfx+'Size:a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size:mx'])**2/ateln2
2338    else:           #ellipsoidal crystallites
2339        const = 1.8*wave/(np.pi*costh)
2340        Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
2341        H = np.array(refl[:3])
2342        lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
2343        Sgam = 1.8*wave/(np.pi*costh*lenR)
2344        for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
2345            gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size:mx']
2346            sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size:mx'])**2/ateln2
2347    gamDict[phfx+'Size:mx'] = Sgam
2348    sigDict[phfx+'Size:mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size:mx'])/ateln2
2349           
2350    #microstrain derivatives               
2351    if calcControls[phfx+'MustrainType'] == 'isotropic':
2352        Mgam = 0.018*parmDict[phfx+'Mustrain:i']*tand(refl[5]/2.)/np.pi
2353        gamDict[phfx+'Mustrain:i'] =  0.018*tanth*parmDict[phfx+'Mustrain:mx']/np.pi
2354        sigDict[phfx+'Mustrain:i'] =  0.036*Mgam*tanth*(1.-parmDict[phfx+'Mustrain:mx'])**2/(np.pi*ateln2)       
2355    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2356        H = np.array(refl[:3])
2357        P = np.array(calcControls[phfx+'MustrainAxis'])
2358        cosP,sinP = G2lat.CosSinAngle(H,P,G)
2359        Si = parmDict[phfx+'Mustrain:i']
2360        Sa = parmDict[phfx+'Mustrain:a']
2361        gami = 0.018*Si*Sa*tanth/np.pi
2362        sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
2363        Mgam = gami/sqtrm
2364        dsi = -gami*Si*cosP**2/sqtrm**3
2365        dsa = -gami*Sa*sinP**2/sqtrm**3
2366        gamDict[phfx+'Mustrain:i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain:mx']
2367        gamDict[phfx+'Mustrain:a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain:mx']
2368        sigDict[phfx+'Mustrain:i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain:mx'])**2/ateln2
2369        sigDict[phfx+'Mustrain:a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain:mx'])**2/ateln2       
2370    else:       #generalized - P.W. Stephens model
2371        pwrs = calcControls[phfx+'MuPwrs']
2372        const = 0.018*refl[4]**2*tanth
2373        sum = 0
2374        for i,pwr in enumerate(pwrs):
2375            term = refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
2376            sum += parmDict[phfx+'Mustrain:'+str(i)]*term
2377            gamDict[phfx+'Mustrain:'+str(i)] = const*term*parmDict[phfx+'Mustrain:mx']
2378            sigDict[phfx+'Mustrain:'+str(i)] = \
2379                2.*const*term*(1.-parmDict[phfx+'Mustrain:mx'])**2/ateln2
2380        Mgam = 0.018*refl[4]**2*tand(refl[5]/2.)*sum
2381        for i in range(len(pwrs)):
2382            sigDict[phfx+'Mustrain:'+str(i)] *= Mgam
2383    gamDict[phfx+'Mustrain:mx'] = Mgam
2384    sigDict[phfx+'Mustrain:mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain:mx'])/ateln2
2385    return sigDict,gamDict
2386       
2387def GetReflPos(refl,wave,G,hfx,calcControls,parmDict):
2388    h,k,l = refl[:3]
2389    dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G)
2390    d = np.sqrt(dsq)
2391
2392    refl[4] = d
2393    pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
2394    const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
2395    if 'Bragg' in calcControls[hfx+'instType']:
2396        pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
2397            parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
2398    else:               #Debye-Scherrer - simple but maybe not right
2399        pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
2400    return pos
2401
2402def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict):
2403    dpr = 180./np.pi
2404    h,k,l = refl[:3]
2405    dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
2406    dst = np.sqrt(dstsq)
2407    pos = refl[5]-parmDict[hfx+'Zero']
2408    const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0)
2409    dpdw = const*dst
2410    dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])
2411    dpdA *= const*wave/(2.0*dst)
2412    dpdZ = 1.0
2413    const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
2414    if 'Bragg' in calcControls[hfx+'instType']:
2415        dpdSh = -4.*const*cosd(pos/2.0)
2416        dpdTr = -const*sind(pos)*100.0
2417        return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.
2418    else:               #Debye-Scherrer - simple but maybe not right
2419        dpdXd = -const*cosd(pos)
2420        dpdYd = -const*sind(pos)
2421        return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd
2422           
2423def GetHStrainShift(refl,SGData,phfx,parmDict):
2424    laue = SGData['SGLaue']
2425    uniq = SGData['SGUniq']
2426    h,k,l = refl[:3]
2427    if laue in ['m3','m3m']:
2428        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
2429            refl[4]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
2430    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2431        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
2432    elif laue in ['3R','3mR']:
2433        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
2434    elif laue in ['4/m','4/mmm']:
2435        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
2436    elif laue in ['mmm']:
2437        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
2438    elif laue in ['2/m']:
2439        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
2440        if uniq == 'a':
2441            Dij += parmDict[phfx+'D23']*k*l
2442        elif uniq == 'b':
2443            Dij += parmDict[phfx+'D13']*h*l
2444        elif uniq == 'c':
2445            Dij += parmDict[phfx+'D12']*h*k
2446    else:
2447        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
2448            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
2449    return -Dij*refl[4]**2*tand(refl[5]/2.0)
2450           
2451def GetHStrainShiftDerv(refl,SGData,phfx):
2452    laue = SGData['SGLaue']
2453    uniq = SGData['SGUniq']
2454    h,k,l = refl[:3]
2455    if laue in ['m3','m3m']:
2456        dDijDict = {phfx+'D11':h**2+k**2+l**2,
2457            phfx+'eA':refl[4]**2*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
2458    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2459        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
2460    elif laue in ['3R','3mR']:
2461        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
2462    elif laue in ['4/m','4/mmm']:
2463        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
2464    elif laue in ['mmm']:
2465        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
2466    elif laue in ['2/m']:
2467        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
2468        if uniq == 'a':
2469            dDijDict[phfx+'D23'] = k*l
2470        elif uniq == 'b':
2471            dDijDict[phfx+'D13'] = h*l
2472        elif uniq == 'c':
2473            dDijDict[phfx+'D12'] = h*k
2474            names.append()
2475    else:
2476        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
2477            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
2478    for item in dDijDict:
2479        dDijDict[item] *= -refl[4]**2*tand(refl[5]/2.0)
2480    return dDijDict
2481   
2482def GetFprime(controlDict,Histograms):
2483    FFtables = controlDict['FFtables']
2484    if not FFtables:
2485        return
2486    histoList = Histograms.keys()
2487    histoList.sort()
2488    for histogram in histoList:
2489        if histogram[:4] in ['PWDR','HKLF']:
2490            Histogram = Histograms[histogram]
2491            hId = Histogram['hId']
2492            hfx = ':%d:'%(hId)
2493            keV = controlDict[hfx+'keV']
2494            for El in FFtables:
2495                Orbs = G2el.GetXsectionCoeff(El.split('+')[0].split('-')[0])
2496                FP,FPP,Mu = G2el.FPcalc(Orbs, keV)
2497                FFtables[El][hfx+'FP'] = FP
2498                FFtables[El][hfx+'FPP'] = FPP               
2499           
2500def GetFobsSq(Histograms,Phases,parmDict,calcControls):
2501    histoList = Histograms.keys()
2502    histoList.sort()
2503    for histogram in histoList:
2504        if 'PWDR' in histogram[:4]:
2505            Histogram = Histograms[histogram]
2506            hId = Histogram['hId']
2507            hfx = ':%d:'%(hId)
2508            Limits = calcControls[hfx+'Limits']
2509            shl = max(parmDict[hfx+'SH/L'],0.002)
2510            Ka2 = False
2511            kRatio = 0.0
2512            if hfx+'Lam1' in parmDict.keys():
2513                Ka2 = True
2514                lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2515                kRatio = parmDict[hfx+'I(L2)/I(L1)']
2516            x,y,w,yc,yb,yd = Histogram['Data']
2517            xB = np.searchsorted(x,Limits[0])
2518            xF = np.searchsorted(x,Limits[1])
2519            ymb = np.array(y-yb)
2520            ymb = np.where(ymb==0.,1.0,ymb)
2521            ycmb = np.array(yc-yb)
2522            ratio = ymb/ycmb           
2523            refLists = Histogram['Reflection Lists']
2524            for phase in refLists:
2525                Phase = Phases[phase]
2526                pId = Phase['pId']
2527                phfx = '%d:%d:'%(pId,hId)
2528                refList = refLists[phase]
2529                sumFo = 0.0
2530                sumdF = 0.0
2531                sumFosq = 0.0
2532                sumdFsq = 0.0
2533                for refl in refList:
2534                    if 'C' in calcControls[hfx+'histType']:
2535                        yp = np.zeros_like(yb)
2536                        Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
2537                        iBeg = max(xB,np.searchsorted(x,refl[5]-fmin))
2538                        iFin = max(xB,min(np.searchsorted(x,refl[5]+fmax),xF))
2539                        iFin2 = iFin
2540                        yp[iBeg:iFin] = refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
2541                        if Ka2:
2542                            pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
2543                            Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
2544                            iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
2545                            iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
2546                            if not iBeg2+iFin2:       #peak below low limit - skip peak
2547                                continue
2548                            elif not iBeg2-iFin2:     #peak above high limit - done
2549                                break
2550                            yp[iBeg2:iFin2] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])        #and here
2551                        refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[13]*(1.+kRatio)),0.0))
2552                    elif 'T' in calcControls[hfx+'histType']:
2553                        print 'TOF Undefined at present'
2554                        raise Exception    #no TOF yet
2555                    Fo = np.sqrt(np.abs(refl[8]))
2556                    Fc = np.sqrt(np.abs(refl[9]))
2557                    sumFo += Fo
2558                    sumFosq += refl[8]**2
2559                    sumdF += np.abs(Fo-Fc)
2560                    sumdFsq += (refl[8]-refl[9])**2
2561                Histogram[phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
2562                Histogram[phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
2563                Histogram[phfx+'Nref'] = len(refList)
2564               
2565def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
2566   
2567    def GetReflSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict):
2568        U = parmDict[hfx+'U']
2569        V = parmDict[hfx+'V']
2570        W = parmDict[hfx+'W']
2571        X = parmDict[hfx+'X']
2572        Y = parmDict[hfx+'Y']
2573        tanPos = tand(refl[5]/2.0)
2574        Ssig,Sgam = GetSampleSigGam(refl,wave,G,GB,phfx,calcControls,parmDict)
2575        sig = U*tanPos**2+V*tanPos+W+Ssig     #save peak sigma
2576        sig = max(0.001,sig)
2577        gam = X/cosd(refl[5]/2.0)+Y*tanPos+Sgam     #save peak gamma
2578        gam = max(0.001,gam)
2579        return sig,gam
2580               
2581    hId = Histogram['hId']
2582    hfx = ':%d:'%(hId)
2583    bakType = calcControls[hfx+'bakType']
2584    yb = G2pwd.getBackground(hfx,parmDict,bakType,x)
2585    yc = np.zeros_like(yb)
2586       
2587    if 'C' in calcControls[hfx+'histType']:   
2588        shl = max(parmDict[hfx+'SH/L'],0.002)
2589        Ka2 = False
2590        if hfx+'Lam1' in parmDict.keys():
2591            wave = parmDict[hfx+'Lam1']
2592            Ka2 = True
2593            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2594            kRatio = parmDict[hfx+'I(L2)/I(L1)']
2595        else:
2596            wave = parmDict[hfx+'Lam']
2597    else:
2598        print 'TOF Undefined at present'
2599        raise ValueError
2600    for phase in Histogram['Reflection Lists']:
2601        refList = Histogram['Reflection Lists'][phase]
2602        Phase = Phases[phase]
2603        pId = Phase['pId']
2604        pfx = '%d::'%(pId)
2605        phfx = '%d:%d:'%(pId,hId)
2606        hfx = ':%d:'%(hId)
2607        SGData = Phase['General']['SGData']
2608        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2609        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2610        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
2611        Vst = np.sqrt(nl.det(G))    #V*
2612        if not Phase['General'].get('doPawley'):
2613            refList = StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict)
2614        for refl in refList:
2615            if 'C' in calcControls[hfx+'histType']:
2616                h,k,l = refl[:3]
2617                refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict)         #corrected reflection position
2618                Lorenz = 1./(2.*sind(refl[5]/2.)**2*cosd(refl[5]/2.))           #Lorentz correction
2619                refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict)               #apply hydrostatic strain shift
2620                refl[6:8] = GetReflSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict)    #peak sig & gam
2621                GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)    #puts corrections in refl[13]
2622                refl[13] *= Vst*Lorenz
2623                if Phase['General'].get('doPawley'):
2624                    try:
2625                        pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
2626                        parmDict[pInd] = max(parmDict[pInd]/2.,parmDict[pInd])       
2627                        refl[9] = parmDict[pInd]
2628                    except KeyError:
2629#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
2630                        continue
2631                Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
2632                iBeg = np.searchsorted(x,refl[5]-fmin)
2633                iFin = np.searchsorted(x,refl[5]+fmax)
2634                if not iBeg+iFin:       #peak below low limit - skip peak
2635                    continue
2636                elif not iBeg-iFin:     #peak above high limit - done
2637                    break
2638                yc[iBeg:iFin] += refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
2639                if Ka2:
2640                    pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
2641                    Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
2642                    iBeg = np.searchsorted(x,pos2-fmin)
2643                    iFin = np.searchsorted(x,pos2+fmax)
2644                    if not iBeg+iFin:       #peak below low limit - skip peak
2645                        continue
2646                    elif not iBeg-iFin:     #peak above high limit - done
2647                        return yc,yb
2648                    yc[iBeg:iFin] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin])        #and here
2649            elif 'T' in calcControls[hfx+'histType']:
2650                print 'TOF Undefined at present'
2651                raise Exception    #no TOF yet
2652    return yc,yb
2653   
2654def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
2655   
2656    def cellVaryDerv(pfx,SGData,dpdA): 
2657        if SGData['SGLaue'] in ['-1',]:
2658            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
2659                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
2660        elif SGData['SGLaue'] in ['2/m',]:
2661            if SGData['SGUniq'] == 'a':
2662                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
2663            elif SGData['SGUniq'] == 'b':
2664                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
2665            else:
2666                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
2667        elif SGData['SGLaue'] in ['mmm',]:
2668            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
2669        elif SGData['SGLaue'] in ['4/m','4/mmm']:
2670#            return [[pfx+'A0',dpdA[0]+dpdA[1]],[pfx+'A2',dpdA[2]]]
2671            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2672        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
2673#            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[3]],[pfx+'A2',dpdA[2]]]
2674            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2675        elif SGData['SGLaue'] in ['3R', '3mR']:
2676            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
2677        elif SGData['SGLaue'] in ['m3m','m3']:
2678#            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]]]
2679            return [[pfx+'A0',dpdA[0]]]
2680    # create a list of dependent variables and set up a dictionary to hold their derivatives
2681    dependentVars = G2mv.GetDependentVars()
2682    depDerivDict = {}
2683    for j in dependentVars:
2684        depDerivDict[j] = np.zeros(shape=(len(x)))
2685    #print 'dependent vars',dependentVars
2686    lenX = len(x)               
2687    hId = Histogram['hId']
2688    hfx = ':%d:'%(hId)
2689    bakType = calcControls[hfx+'bakType']
2690    dMdv = np.zeros(shape=(len(varylist),len(x)))
2691    dMdb,dMddb,dMdpk = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x)
2692    if hfx+'Back:0' in varylist: # for now assume that Back:x vars to not appear in constraints
2693        bBpos =varylist.index(hfx+'Back:0')
2694        dMdv[bBpos:bBpos+len(dMdb)] = dMdb
2695    names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
2696    for name in varylist:
2697        if 'Debye' in name:
2698            id = int(name.split(':')[-1])
2699            parm = name[:int(name.rindex(':'))]
2700            ip = names.index(parm)
2701            dMdv[varylist.index(name)] = dMddb[3*id+ip]
2702    names = [hfx+'BkPkpos',hfx+'BkPkint',hfx+'BkPksig',hfx+'BkPkgam']
2703    for name in varylist:
2704        if 'BkPk' in name:
2705            id = int(name.split(':')[-1])
2706            parm = name[:int(name.rindex(':'))]
2707            ip = names.index(parm)
2708            dMdv[varylist.index(name)] = dMdpk[4*id+ip]
2709    if 'C' in calcControls[hfx+'histType']:   
2710        dx = x[1]-x[0]
2711        shl = max(parmDict[hfx+'SH/L'],0.002)
2712        Ka2 = False
2713        if hfx+'Lam1' in parmDict.keys():
2714            wave = parmDict[hfx+'Lam1']
2715            Ka2 = True
2716            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2717            kRatio = parmDict[hfx+'I(L2)/I(L1)']
2718        else:
2719            wave = parmDict[hfx+'Lam']
2720    else:
2721        print 'TOF Undefined at present'
2722        raise ValueError
2723    for phase in Histogram['Reflection Lists']:
2724        refList = Histogram['Reflection Lists'][phase]
2725        Phase = Phases[phase]
2726        SGData = Phase['General']['SGData']
2727        pId = Phase['pId']
2728        pfx = '%d::'%(pId)
2729        phfx = '%d:%d:'%(pId,hId)
2730        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2731        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2732        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
2733        if not Phase['General'].get('doPawley'):
2734            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)
2735        for iref,refl in enumerate(refList):
2736            if 'C' in calcControls[hfx+'histType']:        #CW powder
2737                h,k,l = refl[:3]
2738                dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA = GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2739                Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
2740                iBeg = np.searchsorted(x,refl[5]-fmin)
2741                iFin = np.searchsorted(x,refl[5]+fmax)
2742                if not iBeg+iFin:       #peak below low limit - skip peak
2743                    continue
2744                elif not iBeg-iFin:     #peak above high limit - done
2745                    break
2746                pos = refl[5]
2747                tanth = tand(pos/2.0)
2748                costh = cosd(pos/2.0)
2749                lenBF = iFin-iBeg
2750                dMdpk = np.zeros(shape=(6,lenBF))
2751                dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])
2752                for i in range(1,5):
2753                    dMdpk[i] += 100.*dx*refl[13]*refl[9]*dMdipk[i]
2754                dMdpk[0] += 100.*dx*refl[13]*refl[9]*dMdipk[0]
2755                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
2756                if Ka2:
2757                    pos2 = refl[5]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
2758                    kdelt = int((pos2-refl[5])/dx)               
2759                    iBeg2 = min(lenX,iBeg+kdelt)
2760                    iFin2 = min(lenX,iFin+kdelt)
2761                    if iBeg2-iFin2:
2762                        lenBF2 = iFin2-iBeg2
2763                        dMdpk2 = np.zeros(shape=(6,lenBF2))
2764                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])
2765                        for i in range(1,5):
2766                            dMdpk2[i] = 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[i]
2767                        dMdpk2[0] = 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[0]
2768                        dMdpk2[5] = 100.*dx*refl[13]*dMdipk2[0]
2769                        dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
2770                if Phase['General'].get('doPawley'):
2771                    dMdpw = np.zeros(len(x))
2772                    try:
2773                        pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
2774                        idx = varylist.index(pIdx)
2775                        dMdpw[iBeg:iFin] = dervDict['int']/refl[9]
2776                        if parmDict[pIdx] < 0.:
2777                            dMdpw[iBeg:iFin] = 2.*dervDict['int']/refl[9]
2778                        if Ka2:
2779                            dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9]
2780                            if parmDict[pIdx] < 0.:
2781                                dMdpw[iBeg2:iFin2] += 2.*dervDict['int']/refl[9]
2782                        dMdv[idx] = dMdpw
2783                    except ValueError:
2784                        pass
2785                dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
2786                names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
2787                    hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
2788                    hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
2789                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
2790                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
2791                    hfx+'DisplaceY':[dpdY,'pos'],}
2792                for name in names:
2793                    item = names[name]
2794                    if name in varylist:
2795                        dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
2796                        if Ka2:
2797                            dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
2798                    elif name in dependentVars:
2799                        if Ka2:
2800                            depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
2801                        depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
2802                for iPO in dIdPO:
2803                    if iPO in varylist:
2804                        dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
2805                        if Ka2:
2806                            dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
2807                    elif iPO in dependentVars:
2808                        depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
2809                        if Ka2:
2810                            depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
2811                for i,name in enumerate(['omega','chi','phi']):
2812                    aname = pfx+'SH '+name
2813                    if aname in varylist:
2814                        dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
2815                        if Ka2:
2816                            dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
2817                    elif aname in dependentVars:
2818                        depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
2819                        if Ka2:
2820                            depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
2821                for iSH in dFdODF:
2822                    if iSH in varylist:
2823                        dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
2824                        if Ka2:
2825                            dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
2826                    elif iSH in dependentVars:
2827                        depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
2828                        if Ka2:
2829                            depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
2830                cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
2831                for name,dpdA in cellDervNames:
2832                    if name in varylist:
2833                        dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
2834                        if Ka2:
2835                            dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
2836                    elif name in dependentVars:
2837                        depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
2838                        if Ka2:
2839                            depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
2840                dDijDict = GetHStrainShiftDerv(refl,SGData,phfx)
2841                for name in dDijDict:
2842                    if name in varylist:
2843                        dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
2844                        if Ka2:
2845                            dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
2846                    elif name in dependentVars:
2847                        depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
2848                        if Ka2:
2849                            depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
2850                sigDict,gamDict = GetSampleSigGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict)
2851                for name in gamDict:
2852                    if name in varylist:
2853                        dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
2854                        if Ka2:
2855                            dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
2856                    elif name in dependentVars:
2857                        depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
2858                        if Ka2:
2859                            depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
2860                for name in sigDict:
2861                    if name in varylist:
2862                        dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
2863                        if Ka2:
2864                            dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
2865                    elif name in dependentVars:
2866                        depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
2867                        if Ka2:
2868                            depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
2869                                               
2870            elif 'T' in calcControls[hfx+'histType']:
2871                print 'TOF Undefined at present'
2872                raise Exception    #no TOF yet
2873            #do atom derivatives -  for F,X & U so far             
2874            corr = dervDict['int']/refl[9]
2875            if Ka2:
2876                corr2 = dervDict2['int']/refl[9]
2877            for name in varylist+dependentVars:
2878                try:
2879                    aname = name.split(pfx)[1][:2]
2880                    if aname not in ['Af','dA','AU']: continue # skip anything not an atom param
2881                except IndexError:
2882                    continue
2883                if name in varylist:
2884                    dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
2885                    if Ka2:
2886                        dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
2887                elif name in dependentVars:
2888                    depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
2889                    if Ka2:
2890                        depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
2891    # now process derivatives in constraints
2892    G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
2893    return dMdv
2894
2895def dervRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
2896    parmdict.update(zip(varylist,values))
2897    G2mv.Dict2Map(parmdict,varylist)
2898    Histograms,Phases = HistoPhases
2899    nvar = len(varylist)
2900    dMdv = np.empty(0)
2901    histoList = Histograms.keys()
2902    histoList.sort()
2903    for histogram in histoList:
2904        if 'PWDR' in histogram[:4]:
2905            Histogram = Histograms[histogram]
2906            hId = Histogram['hId']
2907            hfx = ':%d:'%(hId)
2908            Limits = calcControls[hfx+'Limits']
2909            x,y,w,yc,yb,yd = Histogram['Data']
2910            xB = np.searchsorted(x,Limits[0])
2911            xF = np.searchsorted(x,Limits[1])
2912            dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF],
2913                varylist,Histogram,Phases,calcControls,pawleyLookup)
2914        elif 'HKLF' in histogram[:4]:
2915            Histogram = Histograms[histogram]
2916            nobs = Histogram['Nobs']
2917            phase = Histogram['Reflection Lists']
2918            Phase = Phases[phase]
2919            hId = Histogram['hId']
2920            hfx = ':%d:'%(hId)
2921            pfx = '%d::'%(Phase['pId'])
2922            phfx = '%d:%d:'%(Phase['pId'],hId)
2923            SGData = Phase['General']['SGData']
2924            A = [parmdict[pfx+'A%d'%(i)] for i in range(6)]
2925            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2926            refList = Histogram['Data']
2927            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmdict)
2928            dMdvh = np.zeros((len(varylist),len(refList)))
2929            for iref,ref in enumerate(refList):
2930                if ref[6] > 0:
2931                    if calcControls['F**2']:
2932                        if ref[5]/ref[6] >= calcControls['minF/sig']:
2933                            for j,var in enumerate(varylist):
2934                                if var in dFdvDict:
2935                                    dMdvh[j][iref] = dFdvDict[var][iref]/ref[6]
2936                            if phfx+'Scale' in varylist:
2937                                dMdvh[varylist.index(phfx+'Scale')][iref] = ref[9]/ref[6]
2938                    else:
2939                        Fo = np.sqrt(ref[5])
2940                        Fc = np.sqrt(ref[7])
2941                        sig = ref[6]/(2.0*Fo)
2942                        if Fo/sig >= calcControls['minF/sig']:
2943                            for j,var in enumerate(varylist):
2944                                if var in dFdvDict:
2945                                    dMdvh[j][iref] = dFdvDict[var][iref]/ref[6]
2946                            if phfx+'Scale' in varylist:
2947                                dMdvh[varylist.index(phfx+'Scale')][iref] = ref[9]/ref[6]                           
2948        else:
2949            continue        #skip non-histogram entries
2950        if len(dMdv):
2951            dMdv = np.concatenate((dMdv.T,dMdvh.T)).T
2952        else:
2953            dMdv = dMdvh
2954           
2955#    dpdv = penaltyDeriv(parmdict,varylist)
2956#    if np.any(dpdv):
2957#        dMdv = np.concatenate((dMdv.T,np.outer(dpdv,dpdv))).T
2958    return dMdv
2959
2960def HessRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
2961    parmdict.update(zip(varylist,values))
2962    G2mv.Dict2Map(parmdict,varylist)
2963    Histograms,Phases = HistoPhases
2964    nvar = len(varylist)
2965    Hess = np.empty(0)
2966    histoList = Histograms.keys()
2967    histoList.sort()
2968    for histogram in histoList:
2969        if 'PWDR' in histogram[:4]:
2970            Histogram = Histograms[histogram]
2971            hId = Histogram['hId']
2972            hfx = ':%d:'%(hId)
2973            Limits = calcControls[hfx+'Limits']
2974            x,y,w,yc,yb,yd = Histogram['Data']
2975            dy = y-yc
2976            xB = np.searchsorted(x,Limits[0])
2977            xF = np.searchsorted(x,Limits[1])
2978            dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF],
2979                varylist,Histogram,Phases,calcControls,pawleyLookup)
2980            if dlg:
2981                dlg.Update(Histogram['wR'],newmsg='Hessian for histogram %d Rw=%8.3f%s'%(hId,Histogram['wR'],'%'))[0]
2982            if len(Hess):
2983                Vec += np.sum(dMdvh*np.sqrt(w[xB:xF])*dy[xB:xF],axis=1)
2984                Hess += np.inner(dMdvh,dMdvh)
2985            else:
2986                Vec = np.sum(dMdvh*np.sqrt(w[xB:xF])*dy[xB:xF],axis=1)
2987                Hess = np.inner(dMdvh,dMdvh)
2988        elif 'HKLF' in histogram[:4]:
2989            Histogram = Histograms[histogram]
2990            nobs = Histogram['Nobs']
2991            phase = Histogram['Reflection Lists']
2992            Phase = Phases[phase]
2993            hId = Histogram['hId']
2994            hfx = ':%d:'%(hId)
2995            pfx = '%d::'%(Phase['pId'])
2996            phfx = '%d:%d:'%(Phase['pId'],hId)
2997            SGData = Phase['General']['SGData']
2998            A = [parmdict[pfx+'A%d'%(i)] for i in range(6)]
2999            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
3000            refList = Histogram['Data']
3001            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmdict)
3002            dMdvh = np.zeros((len(varylist),len(refList)))
3003            wdf = np.zeros(len(refList))
3004            for iref,ref in enumerate(refList):
3005                if ref[6] > 0:
3006                    if calcControls['F**2']:
3007                        if ref[5]/ref[6] >= calcControls['minF/sig']:
3008                            wdf[iref] = (ref[5]-ref[7])/ref[6]
3009                            for j,var in enumerate(varylist):
3010                                if var in dFdvDict:
3011                                    dMdvh[j][iref] = dFdvDict[var][iref]/ref[6]
3012                            if phfx+'Scale' in varylist:
3013                                dMdvh[varylist.index(phfx+'Scale')][iref] = ref[9]/ref[6]
3014                    else:
3015                        Fo = np.sqrt(ref[5])
3016                        Fc = np.sqrt(ref[7])
3017                        sig = ref[6]/(2.0*Fo)
3018                        wdf[iref] = (Fo-Fc)/sig
3019                        if Fo/sig >= calcControls['minF/sig']:
3020                            for j,var in enumerate(varylist):
3021                                if var in dFdvDict:
3022                                    dMdvh[j][iref] = dFdvDict[var][iref]/ref[6]
3023                            if phfx+'Scale' in varylist:
3024                                dMdvh[varylist.index(phfx+'Scale')][iref] = ref[9]/ref[6]                           
3025            if dlg:
3026                dlg.Update(Histogram['wR'],newmsg='Hessian for histogram %d Rw=%8.3f%s'%(hId,Histogram['wR'],'%'))[0]
3027            if len(Hess):
3028                Vec += np.sum(dMdvh*wdf,axis=1)
3029                Hess += np.inner(dMdvh,dMdvh)
3030            else:
3031                Vec = np.sum(dMdvh*wdf,axis=1)
3032                Hess = np.inner(dMdvh,dMdvh)
3033        else:
3034            continue        #skip non-histogram entries
3035#    dpdv = penaltyDeriv(parmdict,varylist)
3036#    if np.any(dpdv):
3037#        pVec = dpdv*penaltyFxn(parmdict,varylist)
3038#        Vec += pVec
3039#        pHess = np.zeros((len(varylist),len(varylist)))
3040#        for i,val in enumerate(dpdv):
3041#            pHess[i][i] = dpdv[i]**2
3042#        Hess += pHess
3043    return Vec,Hess
3044
3045def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):       
3046    parmdict.update(zip(varylist,values))
3047    Values2Dict(parmdict, varylist, values)
3048    G2mv.Dict2Map(parmdict,varylist)
3049    Histograms,Phases = HistoPhases
3050    M = np.empty(0)
3051    SumwYo = 0
3052    Nobs = 0
3053    histoList = Histograms.keys()
3054    histoList.sort()
3055    for histogram in histoList:
3056        if 'PWDR' in histogram[:4]:
3057            Histogram = Histograms[histogram]
3058            hId = Histogram['hId']
3059            hfx = ':%d:'%(hId)
3060            Limits = calcControls[hfx+'Limits']
3061            x,y,w,yc,yb,yd = Histogram['Data']
3062            yc *= 0.0                           #zero full calcd profiles
3063            yb *= 0.0
3064            yd *= 0.0
3065            xB = np.searchsorted(x,Limits[0])
3066            xF = np.searchsorted(x,Limits[1])
3067            Histogram['Nobs'] = xF-xB
3068            Nobs += Histogram['Nobs']
3069            Histogram['sumwYo'] = np.sum(w[xB:xF]*y[xB:xF]**2)
3070            SumwYo += Histogram['sumwYo']
3071            yc[xB:xF],yb[xB:xF] = getPowderProfile(parmdict,x[xB:xF],
3072                varylist,Histogram,Phases,calcControls,pawleyLookup)
3073            yc[xB:xF] += yb[xB:xF]
3074            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
3075            Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF]))
3076            wdy = -np.sqrt(w[xB:xF])*(yd[xB:xF])
3077            Histogram['wR'] = min(100.,np.sqrt(np.sum(wdy**2)/Histogram['sumwYo'])*100.)
3078            if dlg:
3079                dlg.Update(Histogram['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['wR'],'%'))[0]
3080            M = np.concatenate((M,wdy))
3081#end of PWDR processing
3082        elif 'HKLF' in histogram[:4]:
3083            Histogram = Histograms[histogram]
3084            phase = Histogram['Reflection Lists']
3085            Phase = Phases[phase]
3086            hId = Histogram['hId']
3087            hfx = ':%d:'%(hId)
3088            pfx = '%d::'%(Phase['pId'])
3089            phfx = '%d:%d:'%(Phase['pId'],hId)
3090            SGData = Phase['General']['SGData']
3091            A = [parmdict[pfx+'A%d'%(i)] for i in range(6)]
3092            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
3093            refList = Histogram['Data']
3094            refList = StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmdict)
3095            df = np.zeros(len(refList))
3096            sumwYo = 0
3097            sumFo = 0
3098            sumFo2 = 0
3099            sumdF = 0
3100            sumdF2 = 0
3101            nobs = 0
3102            for i,ref in enumerate(refList):
3103                if ref[6] > 0:
3104                    ref[7] = parmdict[phfx+'Scale']*ref[9]
3105                    ref[8] = ref[5]/parmdict[phfx+'Scale']
3106                    if calcControls['F**2']:
3107                        if ref[5]/ref[6] >= calcControls['minF/sig']:
3108                            sumFo2 += ref[5]
3109                            Fo = np.sqrt(ref[5])
3110                            sumFo += Fo
3111                            sumFo2 += ref[5]
3112                            sumdF += abs(Fo-np.sqrt(ref[7]))
3113                            sumdF2 += abs(ref[5]-ref[7])
3114                            nobs += 1
3115                            df[i] = -(ref[5]-ref[7])/ref[6]
3116                            sumwYo += (ref[5]/ref[6])**2
3117                    else:
3118                        Fo = np.sqrt(ref[5])
3119                        Fc = np.sqrt(ref[7])
3120                        sig = ref[6]/(2.0*Fo)
3121                        if Fo/sig >= calcControls['minF/sig']:
3122                            sumFo += Fo
3123                            sumFo2 += ref[5]
3124                            sumdF += abs(Fo-Fc)
3125                            sumdF2 += abs(ref[5]-ref[7])
3126                            nobs += 1
3127                            df[i] = -(Fo-Fc)/sig
3128                            sumwYo += (Fo/sig)**2
3129            Histogram['Nobs'] = nobs
3130            Histogram['sumwYo'] = sumwYo
3131            SumwYo += sumwYo
3132            Histogram['wR'] = min(100.,np.sqrt(np.sum(df**2)/Histogram['sumwYo'])*100.)
3133            Histogram[phfx+'Rf'] = 100.*sumdF/sumFo
3134            Histogram[phfx+'Rf^2'] = 100.*sumdF2/sumFo2
3135            Histogram[phfx+'Nref'] = nobs
3136            Nobs += nobs
3137            if dlg:
3138                dlg.Update(Histogram['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['wR'],'%'))[0]
3139            M = np.concatenate((M,df))
3140# end of HKLF processing
3141    Histograms['sumwYo'] = SumwYo
3142    Histograms['Nobs'] = Nobs
3143    Rw = min(100.,np.sqrt(np.sum(M**2)/SumwYo)*100.)
3144    if dlg:
3145        GoOn = dlg.Update(Rw,newmsg='%s%8.3f%s'%('All data Rw =',Rw,'%'))[0]
3146        if not GoOn:
3147            parmDict['saved values'] = values
3148            raise Exception         #Abort!!
3149#    pFunc = penaltyFxn(parmdict,varylist)
3150#    if np.any(pFunc):
3151#        M = np.concatenate((M,pFunc))
3152    return M
3153                       
3154def Refine(GPXfile,dlg):
3155    import pytexture as ptx
3156    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
3157   
3158    ShowBanner()
3159    varyList = []
3160    parmDict = {}
3161    G2mv.InitVars()   
3162    Controls = GetControls(GPXfile)
3163    ShowControls(Controls)
3164    calcControls = {}
3165    calcControls.update(Controls)           
3166    constrDict,fixedList = GetConstraints(GPXfile)
3167    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
3168    if not Phases:
3169        print ' *** ERROR - you have no phases! ***'
3170        print ' *** Refine aborted ***'
3171        raise Exception
3172    if not Histograms:
3173        print ' *** ERROR - you have no data to refine with! ***'
3174        print ' *** Refine aborted ***'
3175        raise Exception       
3176    Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases)
3177    calcControls['Natoms'] = Natoms
3178    calcControls['FFtables'] = FFtables
3179    calcControls['BLtables'] = BLtables
3180    hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms)
3181    calcControls.update(controlDict)
3182    histVary,histDict,controlDict = GetHistogramData(Histograms)
3183    calcControls.update(controlDict)
3184    varyList = phaseVary+hapVary+histVary
3185    parmDict.update(phaseDict)
3186    parmDict.update(hapDict)
3187    parmDict.update(histDict)
3188    GetFprime(calcControls,Histograms)
3189    # do constraint processing
3190    try:
3191        groups,parmlist = G2mv.GroupConstraints(constrDict)
3192        G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList)
3193    except:
3194        print ' *** ERROR - your constraints are internally inconsistent ***'
3195        # traceback for debug
3196        #print 'varyList',varyList
3197        #print 'constrDict',constrDict
3198        #print 'fixedList',fixedList
3199        #import traceback
3200        #print traceback.format_exc()
3201        raise Exception(' *** Refine aborted ***')
3202    # # check to see which generated parameters are fully varied
3203    # msg = G2mv.SetVaryFlags(varyList)
3204    # if msg:
3205    #     print ' *** ERROR - you have not set the refine flags for constraints consistently! ***'
3206    #     print msg
3207    #     raise Exception(' *** Refine aborted ***')
3208    #print G2mv.VarRemapShow(varyList)
3209    G2mv.Map2Dict(parmDict,varyList)
3210    Rvals = {}
3211    while True:
3212        begin = time.time()
3213        values =  np.array(Dict2Values(parmDict, varyList))
3214        Ftol = Controls['min dM/M']
3215        Factor = Controls['shift factor']
3216        maxCyc = Controls['max cyc']
3217        if 'Jacobian' in Controls['deriv type']:           
3218            result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
3219                ftol=Ftol,col_deriv=True,factor=Factor,
3220                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
3221            ncyc = int(result[2]['nfev']/2)
3222        elif 'Hessian' in Controls['deriv type']:
3223            result = G2mth.HessianLSQ(errRefine,values,Hess=HessRefine,ftol=Ftol,maxcyc=maxCyc,
3224                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
3225            ncyc = result[2]['num cyc']+1
3226            Rvals['lamMax'] = result[2]['lamMax']
3227        else:           #'numeric'
3228            result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
3229                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
3230            ncyc = int(result[2]['nfev']/len(varyList))
3231#        table = dict(zip(varyList,zip(values,result[0],(result[0]-values))))
3232#        for item in table: print item,table[item]               #useful debug - are things shifting?
3233        runtime = time.time()-begin
3234        Rvals['chisq'] = np.sum(result[2]['fvec']**2)
3235        Values2Dict(parmDict, varyList, result[0])
3236        G2mv.Dict2Map(parmDict,varyList)
3237       
3238        Rvals['Nobs'] = Histograms['Nobs']
3239        Rvals['Rwp'] = np.sqrt(Rvals['chisq']/Histograms['sumwYo'])*100.      #to %
3240        Rvals['GOF'] = Rvals['chisq']/(Histograms['Nobs']-len(varyList))
3241        print '\n Refinement results:'
3242        print 135*'-'
3243        print ' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histograms['Nobs'],' Number of parameters: ',len(varyList)
3244        print ' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles'%(runtime,runtime/ncyc,ncyc)
3245        print ' wR = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],Rvals['chisq'],Rvals['GOF'])
3246        try:
3247            covMatrix = result[1]*Rvals['GOF']
3248            sig = np.sqrt(np.diag(covMatrix))
3249            if np.any(np.isnan(sig)):
3250                print '*** Least squares aborted - some invalid esds possible ***'
3251#            table = dict(zip(varyList,zip(values,result[0],(result[0]-values)/sig)))
3252#            for item in table: print item,table[item]               #useful debug - are things shifting?
3253            break                   #refinement succeeded - finish up!
3254        except TypeError:          #result[1] is None on singular matrix
3255            print '**** Refinement failed - singular matrix ****'
3256            if 'Hessian' in Controls['deriv type']:
3257                num = len(varyList)-1
3258                for i,val in enumerate(np.flipud(result[2]['psing'])):
3259                    if val:
3260                        print 'Removing parameter: ',varyList[num-i]
3261                        del(varyList[num-i])                   
3262            else:
3263                Ipvt = result[2]['ipvt']
3264                for i,ipvt in enumerate(Ipvt):
3265                    if not np.sum(result[2]['fjac'],axis=1)[i]:
3266                        print 'Removing parameter: ',varyList[ipvt-1]
3267                        del(varyList[ipvt-1])
3268                        break
3269
3270#    print 'dependentParmList: ',G2mv.dependentParmList
3271#    print 'arrayList: ',G2mv.arrayList
3272#    print 'invarrayList: ',G2mv.invarrayList
3273#    print 'indParmList: ',G2mv.indParmList
3274#    print 'fixedDict: ',G2mv.fixedDict
3275#    print 'test1'
3276    GetFobsSq(Histograms,Phases,parmDict,calcControls)
3277#    print 'test2'
3278    sigDict = dict(zip(varyList,sig))
3279    newCellDict = GetNewCellParms(parmDict,varyList)
3280    newAtomDict = ApplyXYZshifts(parmDict,varyList)
3281    covData = {'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
3282        'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict,'newCellDict':newCellDict}
3283    # add the uncertainties into the esd dictionary (sigDict)
3284    sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
3285    SetPhaseData(parmDict,sigDict,Phases,covData)
3286    SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms)
3287    SetHistogramData(parmDict,sigDict,Histograms)
3288    G2mv.PrintIndependentVars(parmDict,varyList,sigDict)
3289    SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,covData)
3290   
3291#for testing purposes!!!
3292    if DEBUG:
3293        import cPickle
3294        fl = open('structTestdata.dat','wb')
3295        cPickle.dump(parmDict,fl,1)
3296        cPickle.dump(varyList,fl,1)
3297        for histogram in Histograms:
3298            if 'PWDR' in histogram[:4]:
3299                Histogram = Histograms[histogram]
3300        cPickle.dump(Histogram,fl,1)
3301        cPickle.dump(Phases,fl,1)
3302        cPickle.dump(calcControls,fl,1)
3303        cPickle.dump(pawleyLookup,fl,1)
3304        fl.close()
3305
3306    if dlg:
3307        return Rvals['Rwp']
3308
3309def SeqRefine(GPXfile,dlg):
3310    import pytexture as ptx
3311    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
3312   
3313    ShowBanner()
3314    print ' Sequential Refinement'
3315    G2mv.InitVars()   
3316    Controls = GetControls(GPXfile)
3317    ShowControls(Controls)           
3318    constrDict,fixedList = GetConstraints(GPXfile)
3319    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
3320    if not Phases:
3321        print ' *** ERROR - you have no histograms to refine! ***'
3322        print ' *** Refine aborted ***'
3323        raise Exception
3324    if not Histograms:
3325        print ' *** ERROR - you have no data to refine with! ***'
3326        print ' *** Refine aborted ***'
3327        raise Exception
3328    Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases,False)
3329    if 'Seq Data' in Controls:
3330        histNames = Controls['Seq Data']
3331    else:
3332        histNames = GetHistogramNames(GPXfile,['PWDR',])
3333    if 'Reverse Seq' in Controls:
3334        if Controls['Reverse Seq']:
3335            histNames.reverse()
3336    SeqResult = {'histNames':histNames}
3337    makeBack = True
3338    for ihst,histogram in enumerate(histNames):
3339        ifPrint = False
3340        if dlg:
3341            dlg.SetTitle('Residual for histogram '+str(ihst))
3342        calcControls = {}
3343        calcControls['Natoms'] = Natoms
3344        calcControls['FFtables'] = FFtables
3345        calcControls['BLtables'] = BLtables
3346        varyList = []
3347        parmDict = {}
3348        Histo = {histogram:Histograms[histogram],}
3349        hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histo,False)
3350        calcControls.update(controlDict)
3351        histVary,histDict,controlDict = GetHistogramData(Histo,False)
3352        calcControls.update(controlDict)
3353        varyList = phaseVary+hapVary+histVary
3354        if not ihst:
3355            saveVaryList = varyList[:]
3356            for i,item in enumerate(saveVaryList):
3357                items = item.split(':')
3358                if items[1]:
3359                    items[1] = ''
3360                item = ':'.join(items)
3361                saveVaryList[i] = item
3362            SeqResult['varyList'] = saveVaryList
3363        else:
3364            newVaryList = varyList[:]
3365            for i,item in enumerate(newVaryList):
3366                items = item.split(':')
3367                if items[1]:
3368                    items[1] = ''
3369                item = ':'.join(items)
3370                newVaryList[i] = item
3371            if newVaryList != SeqResult['varyList']:
3372                print newVaryList
3373                print SeqResult['varyList']
3374                print '**** ERROR - variable list for this histogram does not match previous'
3375                raise Exception
3376        parmDict.update(phaseDict)
3377        parmDict.update(hapDict)
3378        parmDict.update(histDict)
3379        GetFprime(calcControls,Histo)
3380        # do constraint processing
3381        try:
3382            groups,parmlist = G2mv.GroupConstraints(constrDict)
3383            G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList)
3384        except:
3385            print ' *** ERROR - your constraints are internally inconsistent ***'
3386            raise Exception(' *** Refine aborted ***')
3387        # check to see which generated parameters are fully varied
3388        # msg = G2mv.SetVaryFlags(varyList)
3389        # if msg:
3390        #     print ' *** ERROR - you have not set the refine flags for constraints consistently! ***'
3391        #     print msg
3392        #     raise Exception(' *** Refine aborted ***')
3393        #print G2mv.VarRemapShow(varyList)
3394        G2mv.Map2Dict(parmDict,varyList)
3395        Rvals = {}
3396        while True:
3397            begin = time.time()
3398            values =  np.array(Dict2Values(parmDict,varyList))
3399            Ftol = Controls['min dM/M']
3400            Factor = Controls['shift factor']
3401            maxCyc = Controls['max cyc']
3402
3403            if 'Jacobian' in Controls['deriv type']:           
3404                result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
3405                    ftol=Ftol,col_deriv=True,factor=Factor,
3406                    args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
3407                ncyc = int(result[2]['nfev']/2)
3408            elif 'Hessian' in Controls['deriv type']:
3409                result = G2mth.HessianLSQ(errRefine,values,Hess=HessRefine,ftol=Ftol,maxcyc=maxCyc,
3410                    args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
3411                ncyc = result[2]['num cyc']+1                           
3412            else:           #'numeric'
3413                result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
3414                    args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
3415                ncyc = int(result[2]['nfev']/len(varyList))
3416
3417
3418
3419            runtime = time.time()-begin
3420            Rvals['chisq'] = np.sum(result[2]['fvec']**2)
3421            Values2Dict(parmDict, varyList, result[0])
3422            G2mv.Dict2Map(parmDict,varyList)
3423           
3424            Rvals['Rwp'] = np.sqrt(Rvals['chisq']/Histo['sumwYo'])*100.      #to %
3425            Rvals['GOF'] = Rvals['Rwp']/(Histo['Nobs']-len(varyList))
3426            Rvals['Nobs'] = Histo['Nobs']
3427            print '\n Refinement results for histogram: v'+histogram
3428            print 135*'-'
3429            print ' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histo['Nobs'],' Number of parameters: ',len(varyList)
3430            print ' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles'%(runtime,runtime/ncyc,ncyc)
3431            print ' wRp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],Rvals['chisq'],Rvals['GOF'])
3432            try:
3433                covMatrix = result[1]*Rvals['GOF']
3434                sig = np.sqrt(np.diag(covMatrix))
3435                if np.any(np.isnan(sig)):
3436                    print '*** Least squares aborted - some invalid esds possible ***'
3437                    ifPrint = True
3438                break                   #refinement succeeded - finish up!
3439            except TypeError:          #result[1] is None on singular matrix
3440                print '**** Refinement failed - singular matrix ****'
3441                if 'Hessian' in Controls['deriv type']:
3442                    num = len(varyList)-1
3443                    for i,val in enumerate(np.flipud(result[2]['psing'])):
3444                        if val:
3445                            print 'Removing parameter: ',varyList[num-i]
3446                            del(varyList[num-i])                   
3447                else:
3448                    Ipvt = result[2]['ipvt']
3449                    for i,ipvt in enumerate(Ipvt):
3450                        if not np.sum(result[2]['fjac'],axis=1)[i]:
3451                            print 'Removing parameter: ',varyList[ipvt-1]
3452                            del(varyList[ipvt-1])
3453                            break
3454   
3455        GetFobsSq(Histo,Phases,parmDict,calcControls)
3456        sigDict = dict(zip(varyList,sig))
3457        newCellDict = GetNewCellParms(parmDict,varyList)
3458        newAtomDict = ApplyXYZshifts(parmDict,varyList)
3459        covData = {'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
3460            'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,'newCellDict':newCellDict}
3461        SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,ifPrint)
3462        SetHistogramData(parmDict,sigDict,Histo,ifPrint)
3463        SeqResult[histogram] = covData
3464        SetUsedHistogramsAndPhases(GPXfile,Histo,Phases,covData,makeBack)
3465        makeBack = False
3466    SetSeqResult(GPXfile,Histograms,SeqResult)
3467
3468def DistAngle(DisAglCtls,DisAglData):
3469    import numpy.ma as ma
3470   
3471    def ShowBanner(name):
3472        print 80*'*'
3473        print '   Interatomic Distances and Angles for phase '+name
3474        print 80*'*','\n'
3475
3476    ShowBanner(DisAglCtls['Name'])
3477    SGData = DisAglData['SGData']
3478    SGtext = G2spc.SGPrint(SGData)
3479    for line in SGtext: print line
3480    Cell = DisAglData['Cell']
3481   
3482    Amat,Bmat = G2lat.cell2AB(Cell[:6])
3483    covData = {}
3484    if 'covData' in DisAglData:   
3485        covData = DisAglData['covData']
3486        covMatrix = covData['covMatrix']
3487        varyList = covData['varyList']
3488        pfx = str(DisAglData['pId'])+'::'
3489        A = G2lat.cell2A(Cell[:6])
3490        cellSig = getCellEsd(pfx,SGData,A,covData)
3491        names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']
3492        valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]
3493        line = '\n Unit cell:'
3494        for name,vals in zip(names,valEsd):
3495            line += name+vals 
3496        print line
3497    else: 
3498        print '\n Unit cell: a = ','%.5f'%(Cell[0]),' b = ','%.5f'%(Cell[1]),' c = ','%.5f'%(Cell[2]), \
3499            ' alpha = ','%.3f'%(Cell[3]),' beta = ','%.3f'%(Cell[4]),' gamma = ', \
3500            '%.3f'%(Cell[5]),' volume = ','%.3f'%(Cell[6])
3501    Factor = DisAglCtls['Factors']
3502    Radii = dict(zip(DisAglCtls['AtomTypes'],zip(DisAglCtls['BondRadii'],DisAglCtls['AngleRadii'])))
3503    indices = (-1,0,1)
3504    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
3505    origAtoms = DisAglData['OrigAtoms']
3506    targAtoms = DisAglData['TargAtoms']
3507    for Oatom in origAtoms:
3508        OxyzNames = ''
3509        IndBlist = []
3510        Dist = []
3511        Vect = []
3512        VectA = []
3513        angles = []
3514        for Tatom in targAtoms:
3515            Xvcov = []
3516            TxyzNames = ''
3517            if 'covData' in DisAglData:
3518                OxyzNames = [pfx+'dAx:%d'%(Oatom[0]),pfx+'dAy:%d'%(Oatom[0]),pfx+'dAz:%d'%(Oatom[0])]
3519                TxyzNames = [pfx+'dAx:%d'%(Tatom[0]),pfx+'dAy:%d'%(Tatom[0]),pfx+'dAz:%d'%(Tatom[0])]
3520                Xvcov = G2mth.getVCov(OxyzNames+TxyzNames,varyList,covMatrix)
3521            result = G2spc.GenAtom(Tatom[3:6],SGData,False,Move=False)
3522            BsumR = (Radii[Oatom[2]][0]+Radii[Tatom[2]][0])*Factor[0]
3523            AsumR = (Radii[Oatom[2]][1]+Radii[Tatom[2]][1])*Factor[1]
3524            for Txyz,Top,Tunit in result:
3525                Dx = (Txyz-np.array(Oatom[3:6]))+Units
3526                dx = np.inner(Amat,Dx)
3527                dist = ma.masked_less(np.sqrt(np.sum(dx**2,axis=0)),0.5)
3528                IndB = ma.nonzero(ma.masked_greater(dist-BsumR,0.))
3529                if np.any(IndB):
3530                    for indb in IndB:
3531                        for i in range(len(indb)):
3532                            if str(dx.T[indb][i]) not in IndBlist:
3533                                IndBlist.append(str(dx.T[indb][i]))
3534                                unit = Units[indb][i]
3535                                tunit = '[%2d%2d%2d]'%(unit[0]+Tunit[0],unit[1]+Tunit[1],unit[2]+Tunit[2])
3536                                pdpx = G2mth.getDistDerv(Oatom[3:6],Tatom[3:6],Amat,unit,Top,SGData)
3537                                sig = 0.0
3538                                if len(Xvcov):
3539                                    sig = np.sqrt(np.inner(pdpx,np.inner(Xvcov,pdpx)))
3540                                Dist.append([Oatom[1],Tatom[1],tunit,Top,ma.getdata(dist[indb])[i],sig])
3541                                if (Dist[-1][-1]-AsumR) <= 0.:
3542                                    Vect.append(dx.T[indb][i]/Dist[-1][-2])
3543                                    VectA.append([OxyzNames,np.array(Oatom[3:6]),TxyzNames,np.array(Tatom[3:6]),unit,Top])
3544                                else:
3545                                    Vect.append([0.,0.,0.])
3546                                    VectA.append([])
3547        Vect = np.array(Vect)
3548        angles = np.zeros((len(Vect),len(Vect)))
3549        angsig = np.zeros((len(Vect),len(Vect)))
3550        for i,veca in enumerate(Vect):
3551            if np.any(veca):
3552                for j,vecb in enumerate(Vect):
3553                    if np.any(vecb):
3554                        angles[i][j],angsig[i][j] = G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)
3555        line = ''
3556        for i,x in enumerate(Oatom[3:6]):
3557            if len(Xvcov):
3558                line += '%12s'%(G2mth.ValEsd(x,np.sqrt(Xvcov[i][i]),True))
3559            else:
3560                line += '%12.5f'%(x)
3561        print '\n Distances & angles for ',Oatom[1],' at ',line
3562        print 80*'*'
3563        line = ''
3564        for dist in Dist[:-1]:
3565            line += '%12s'%(dist[1].center(12))
3566        print '  To       cell +(sym. op.)      dist.  ',line
3567        for i,dist in enumerate(Dist):
3568            line = ''
3569            for j,angle in enumerate(angles[i][0:i]):
3570                sig = angsig[i][j]
3571                if angle:
3572                    if sig:
3573                        line += '%12s'%(G2mth.ValEsd(angle,sig,True).center(12))
3574                    else:
3575                        val = '%.3f'%(angle)
3576                        line += '%12s'%(val.center(12))
3577                else:
3578                    line += 12*' '
3579            if dist[5]:            #sig exists!
3580                val = G2mth.ValEsd(dist[4],dist[5])
3581            else:
3582                val = '%8.4f'%(dist[4])
3583            print %8s%10s+(%4d) %12s'%(dist[1].ljust(8),dist[2].ljust(10),dist[3],val.center(12)),line
3584
3585def DisAglTor(DATData):
3586    SGData = DATData['SGData']
3587    Cell = DATData['Cell']
3588   
3589    Amat,Bmat = G2lat.cell2AB(Cell[:6])
3590    covData = {}
3591    pfx = ''
3592    if 'covData' in DATData:   
3593        covData = DATData['covData']
3594        covMatrix = covData['covMatrix']
3595        varyList = covData['varyList']
3596        pfx = str(DATData['pId'])+'::'
3597    Datoms = []
3598    Oatoms = []
3599    for i,atom in enumerate(DATData['Datoms']):
3600        symop = atom[-1].split('+')
3601        if len(symop) == 1:
3602            symop.append('0,0,0')       
3603        symop[0] = int(symop[0])
3604        symop[1] = eval(symop[1])
3605        atom.append(symop)
3606        Datoms.append(atom)
3607        oatom = DATData['Oatoms'][i]
3608        names = ['','','']
3609        if pfx:
3610            names = [pfx+'dAx:'+str(oatom[0]),pfx+'dAy:'+str(oatom[0]),pfx+'dAz:'+str(oatom[0])]
3611        oatom += [names,]
3612        Oatoms.append(oatom)
3613    atmSeq = [atom[1]+'('+atom[-2]+')' for atom in Datoms]
3614    if DATData['Natoms'] == 4:  #torsion
3615        Tors,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
3616        print ' Torsion angle for '+DATData['Name']+' atom sequence: ',atmSeq,'=',G2mth.ValEsd(Tors,sig)
3617        print ' NB: Atom sequence determined by selection order'
3618        return      # done with torsion
3619    elif DATData['Natoms'] == 3:  #angle
3620        Ang,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
3621        print ' Angle in '+DATData['Name']+' for atom sequence: ',atmSeq,'=',G2mth.ValEsd(Ang,sig)
3622        print ' NB: Atom sequence determined by selection order'
3623    else:   #2 atoms - distance
3624        Dist,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
3625        print ' Distance in '+DATData['Name']+' for atom sequence: ',atmSeq,'=',G2mth.ValEsd(Dist,sig)
3626               
3627def BestPlane(PlaneData):
3628
3629    def ShowBanner(name):
3630        print 80*'*'
3631        print '   Best plane result for phase '+name
3632        print 80*'*','\n'
3633
3634    ShowBanner(PlaneData['Name'])
3635
3636    Cell = PlaneData['Cell']   
3637    Amat,Bmat = G2lat.cell2AB(Cell[:6])       
3638    Atoms = PlaneData['Atoms']
3639    sumXYZ = np.zeros(3)
3640    XYZ = []
3641    Natoms = len(Atoms)
3642    for atom in Atoms:
3643        xyz = np.array(atom[3:6])
3644        XYZ.append(xyz)
3645        sumXYZ += xyz
3646    sumXYZ /= Natoms
3647    XYZ = np.array(XYZ)-sumXYZ
3648    XYZ = np.inner(Amat,XYZ).T
3649    Zmat = np.zeros((3,3))
3650    for i,xyz in enumerate(XYZ):
3651        Zmat += np.outer(xyz.T,xyz)
3652    print ' Selected atoms centered at %10.5f %10.5f %10.5f'%(sumXYZ[0],sumXYZ[1],sumXYZ[2])
3653    Evec,Emat = nl.eig(Zmat)
3654    Evec = np.sqrt(Evec)/(Natoms-3)
3655    Order = np.argsort(Evec)
3656    XYZ = np.inner(XYZ,Emat.T).T
3657    XYZ = np.array([XYZ[Order[2]],XYZ[Order[1]],XYZ[Order[0]]]).T
3658    print ' Atoms in Cartesian best plane coordinates:'
3659    print ' Name         X         Y         Z'
3660    for i,xyz in enumerate(XYZ):
3661        print ' %6s%10.3f%10.3f%10.3f'%(Atoms[i][1].ljust(6),xyz[0],xyz[1],xyz[2])
3662    print '\n Best plane RMS X =%8.3f, Y =%8.3f, Z =%8.3f'%(Evec[Order[2]],Evec[Order[1]],Evec[Order[0]])   
3663
3664           
3665def main():
3666    arg = sys.argv
3667    if len(arg) > 1:
3668        GPXfile = arg[1]
3669        if not ospath.exists(GPXfile):
3670            print 'ERROR - ',GPXfile," doesn't exist!"
3671            exit()
3672        GPXpath = ospath.dirname(arg[1])
3673        Refine(GPXfile,None)
3674    else:
3675        print 'ERROR - missing filename'
3676        exit()
3677    print "Done"
3678         
3679if __name__ == '__main__':
3680    main()
Note: See TracBrowser for help on using the repository browser.