source: trunk/GSASIIstruct.py @ 860

Last change on this file since 860 was 860, checked in by vondreele, 9 years ago

more fixes to rigid body stuff

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