source: trunk/GSASIIstruct.py @ 660

Last change on this file since 660 was 660, checked in by toby, 11 years ago

try revision number tracking

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