source: trunk/GSASIIstruct.py @ 869

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

a fix to restraint refinement

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 188.5 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASIIstructure - structure computation routines
3########### SVN repository information ###################
4# $Date: 2013-03-21 15:54:41 +0000 (Thu, 21 Mar 2013) $
5# $Author: vondreele $
6# $Revision: 869 $
7# $URL: trunk/GSASIIstruct.py $
8# $Id: GSASIIstruct.py 869 2013-03-21 15:54:41Z 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: 869 $")
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                    if hapData['Babinet']['BabA'][0]:
1319                        PrintBabinet(hapData['Babinet'])
1320                HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
1321                refList = []
1322                for h,k,l,d in HKLd:
1323                    ext,mul,Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
1324                    mul *= 2      # for powder overlap of Friedel pairs
1325                    if ext:
1326                        continue
1327                    if 'C' in inst['Type'][0]:
1328                        pos = 2.0*asind(wave/(2.0*d))+Zero
1329                        if limits[0] < pos < limits[1]:
1330                            refList.append([h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,0.0,Uniq,phi,0.0,{}])
1331                            #last item should contain form factor values by atom type
1332                    else:
1333                        raise ValueError 
1334                Histogram['Reflection Lists'][phase] = refList
1335            elif 'HKLF' in histogram:
1336                inst = Histogram['Instrument Parameters'][0]
1337                hId = Histogram['hId']
1338                hfx = ':%d:'%(hId)
1339                for item in inst:
1340                    hapDict[hfx+item] = inst[item][1]
1341                pfx = str(pId)+':'+str(hId)+':'
1342                hapDict[pfx+'Scale'] = hapData['Scale'][0]
1343                if hapData['Scale'][1]:
1344                    hapVary.append(pfx+'Scale')
1345                               
1346                extApprox,extType,extParms = hapData['Extinction']
1347                controlDict[pfx+'EType'] = extType
1348                controlDict[pfx+'EApprox'] = extApprox
1349                controlDict[pfx+'Tbar'] = extParms['Tbar']
1350                controlDict[pfx+'Cos2TM'] = extParms['Cos2TM']
1351                if 'Primary' in extType:
1352                    Ekey = ['Ep',]
1353                elif 'I & II' in extType:
1354                    Ekey = ['Eg','Es']
1355                elif 'Secondary Type II' == extType:
1356                    Ekey = ['Es',]
1357                elif 'Secondary Type I' == extType:
1358                    Ekey = ['Eg',]
1359                else:   #'None'
1360                    Ekey = []
1361                for eKey in Ekey:
1362                    hapDict[pfx+eKey] = extParms[eKey][0]
1363                    if extParms[eKey][1]:
1364                        hapVary.append(pfx+eKey)
1365                for bab in ['BabA','BabU']:
1366                    hapDict[pfx+bab] = hapData['Babinet'][bab][0]
1367                    if hapData['Babinet'][bab][1]:
1368                        hapVary.append(pfx+bab)
1369                if Print: 
1370                    print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
1371                    print >>pFile,135*'-'
1372                    print >>pFile,' Scale factor     : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1373                    if extType != 'None':
1374                        print >>pFile,' Extinction  Type: %15s'%(extType),' approx: %10s'%(extApprox),' tbar: %6.3f'%(extParms['Tbar'])
1375                        text = ' Parameters       :'
1376                        for eKey in Ekey:
1377                            text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
1378                        print >>pFile,text
1379                    if hapData['Babinet']['BabA'][0]:
1380                        PrintBabinet(hapData['Babinet'])
1381                Histogram['Reflection Lists'] = phase       
1382               
1383    return hapVary,hapDict,controlDict
1384   
1385def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,Print=True,pFile=None):
1386   
1387    def PrintSizeAndSig(hapData,sizeSig):
1388        line = '\n Size model:     %9s'%(hapData[0])
1389        refine = False
1390        if hapData[0] in ['isotropic','uniaxial']:
1391            line += ' equatorial:%12.4f'%(hapData[1][0])
1392            if sizeSig[0][0]:
1393                line += ', sig:%8.4f'%(sizeSig[0][0])
1394                refine = True
1395            if hapData[0] == 'uniaxial':
1396                line += ' axial:%12.4f'%(hapData[1][1])
1397                if sizeSig[0][1]:
1398                    refine = True
1399                    line += ', sig:%8.4f'%(sizeSig[0][1])
1400            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1401            if sizeSig[0][2]:
1402                refine = True
1403                line += ', sig:%8.4f'%(sizeSig[0][2])
1404            if refine:
1405                print >>pFile,line
1406        else:
1407            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1408            if sizeSig[0][2]:
1409                refine = True
1410                line += ', sig:%8.4f'%(sizeSig[0][2])
1411            Snames = ['S11','S22','S33','S12','S13','S23']
1412            ptlbls = ' name  :'
1413            ptstr =  ' value :'
1414            sigstr = ' sig   :'
1415            for i,name in enumerate(Snames):
1416                ptlbls += '%12s' % (name)
1417                ptstr += '%12.6f' % (hapData[4][i])
1418                if sizeSig[1][i]:
1419                    refine = True
1420                    sigstr += '%12.6f' % (sizeSig[1][i])
1421                else:
1422                    sigstr += 12*' '
1423            if refine:
1424                print >>pFile,line
1425                print >>pFile,ptlbls
1426                print >>pFile,ptstr
1427                print >>pFile,sigstr
1428       
1429    def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
1430        line = '\n Mustrain model: %9s'%(hapData[0])
1431        refine = False
1432        if hapData[0] in ['isotropic','uniaxial']:
1433            line += ' equatorial:%12.1f'%(hapData[1][0])
1434            if mustrainSig[0][0]:
1435                line += ', sig:%8.1f'%(mustrainSig[0][0])
1436                refine = True
1437            if hapData[0] == 'uniaxial':
1438                line += ' axial:%12.1f'%(hapData[1][1])
1439                if mustrainSig[0][1]:
1440                     line += ', sig:%8.1f'%(mustrainSig[0][1])
1441            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1442            if mustrainSig[0][2]:
1443                refine = True
1444                line += ', sig:%8.3f'%(mustrainSig[0][2])
1445            if refine:
1446                print >>pFile,line
1447        else:
1448            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1449            if mustrainSig[0][2]:
1450                refine = True
1451                line += ', sig:%8.3f'%(mustrainSig[0][2])
1452            Snames = G2spc.MustrainNames(SGData)
1453            ptlbls = ' name  :'
1454            ptstr =  ' value :'
1455            sigstr = ' sig   :'
1456            for i,name in enumerate(Snames):
1457                ptlbls += '%12s' % (name)
1458                ptstr += '%12.6f' % (hapData[4][i])
1459                if mustrainSig[1][i]:
1460                    refine = True
1461                    sigstr += '%12.6f' % (mustrainSig[1][i])
1462                else:
1463                    sigstr += 12*' '
1464            if refine:
1465                print >>pFile,line
1466                print >>pFile,ptlbls
1467                print >>pFile,ptstr
1468                print >>pFile,sigstr
1469           
1470    def PrintHStrainAndSig(hapData,strainSig,SGData):
1471        Hsnames = G2spc.HStrainNames(SGData)
1472        ptlbls = ' name  :'
1473        ptstr =  ' value :'
1474        sigstr = ' sig   :'
1475        refine = False
1476        for i,name in enumerate(Hsnames):
1477            ptlbls += '%12s' % (name)
1478            ptstr += '%12.6g' % (hapData[0][i])
1479            if name in strainSig:
1480                refine = True
1481                sigstr += '%12.6g' % (strainSig[name])
1482            else:
1483                sigstr += 12*' '
1484        if refine:
1485            print >>pFile,'\n Hydrostatic/elastic strain: '
1486            print >>pFile,ptlbls
1487            print >>pFile,ptstr
1488            print >>pFile,sigstr
1489       
1490    def PrintSHPOAndSig(pfx,hapData,POsig):
1491        print >>pFile,'\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
1492        ptlbls = ' names :'
1493        ptstr =  ' values:'
1494        sigstr = ' sig   :'
1495        for item in hapData[5]:
1496            ptlbls += '%12s'%(item)
1497            ptstr += '%12.3f'%(hapData[5][item])
1498            if pfx+item in POsig:
1499                sigstr += '%12.3f'%(POsig[pfx+item])
1500            else:
1501                sigstr += 12*' ' 
1502        print >>pFile,ptlbls
1503        print >>pFile,ptstr
1504        print >>pFile,sigstr
1505       
1506    def PrintExtAndSig(pfx,hapData,ScalExtSig):
1507        print >>pFile,'\n Single crystal extinction: Type: ',hapData[0],' Approx: ',hapData[1]
1508        text = ''
1509        for item in hapData[2]:
1510            if pfx+item in ScalExtSig:
1511                text += '       %s: '%(item)
1512                text += '%12.2e'%(hapData[2][item][0])
1513                if pfx+item in ScalExtSig:
1514                    text += ' sig: %12.2e'%(ScalExtSig[pfx+item])
1515        print >>pFile,text       
1516       
1517    def PrintBabinetAndSig(pfx,hapData,BabSig):
1518        print >>pFile,'\n Babinet form factor modification: '
1519        ptlbls = ' names :'
1520        ptstr =  ' values:'
1521        sigstr = ' sig   :'
1522        for item in hapData:
1523            ptlbls += '%12s'%(item)
1524            ptstr += '%12.3f'%(hapData[item][0])
1525            if pfx+item in BabSig:
1526                sigstr += '%12.3f'%(BabSig[pfx+item])
1527            else:
1528                sigstr += 12*' ' 
1529        print >>pFile,ptlbls
1530        print >>pFile,ptstr
1531        print >>pFile,sigstr
1532   
1533    PhFrExtPOSig = {}
1534    SizeMuStrSig = {}
1535    ScalExtSig = {}
1536    BabSig = {}
1537    wtFrSum = {}
1538    for phase in Phases:
1539        HistoPhase = Phases[phase]['Histograms']
1540        General = Phases[phase]['General']
1541        SGData = General['SGData']
1542        pId = Phases[phase]['pId']
1543        histoList = HistoPhase.keys()
1544        histoList.sort()
1545        for histogram in histoList:
1546            try:
1547                Histogram = Histograms[histogram]
1548            except KeyError:                       
1549                #skip if histogram not included e.g. in a sequential refinement
1550                continue
1551            hapData = HistoPhase[histogram]
1552            hId = Histogram['hId']
1553            pfx = str(pId)+':'+str(hId)+':'
1554            if hId not in wtFrSum:
1555                wtFrSum[hId] = 0.
1556            if 'PWDR' in histogram:
1557                for item in ['Scale','Extinction']:
1558                    hapData[item][0] = parmDict[pfx+item]
1559                    if pfx+item in sigDict:
1560                        PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
1561                wtFrSum[hId] += hapData['Scale'][0]*General['Mass']
1562                if hapData['Pref.Ori.'][0] == 'MD':
1563                    hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
1564                    if pfx+'MD' in sigDict:
1565                        PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],})
1566                else:                           #'SH' spherical harmonics
1567                    for item in hapData['Pref.Ori.'][5]:
1568                        hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
1569                        if pfx+item in sigDict:
1570                            PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
1571                SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
1572                    pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
1573                    pfx+'HStrain':{}})                 
1574                for item in ['Mustrain','Size']:
1575                    hapData[item][1][2] = parmDict[pfx+item+';mx']
1576                    hapData[item][1][2] = min(1.,max(0.1,hapData[item][1][2]))
1577                    if pfx+item+';mx' in sigDict:
1578                        SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
1579                    if hapData[item][0] in ['isotropic','uniaxial']:                   
1580                        hapData[item][1][0] = parmDict[pfx+item+';i']
1581                        if item == 'Size':
1582                            hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
1583                        if pfx+item+';i' in sigDict: 
1584                            SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
1585                        if hapData[item][0] == 'uniaxial':
1586                            hapData[item][1][1] = parmDict[pfx+item+';a']
1587                            if item == 'Size':
1588                                hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))                       
1589                            if pfx+item+';a' in sigDict:
1590                                SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
1591                    else:       #generalized for mustrain or ellipsoidal for size
1592                        Nterms = len(hapData[item][4])
1593                        for i in range(Nterms):
1594                            sfx = ':'+str(i)
1595                            hapData[item][4][i] = parmDict[pfx+item+sfx]
1596                            if pfx+item+sfx in sigDict:
1597                                SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx]
1598                names = G2spc.HStrainNames(SGData)
1599                for i,name in enumerate(names):
1600                    hapData['HStrain'][0][i] = parmDict[pfx+name]
1601                    if pfx+name in sigDict:
1602                        SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
1603                for name in ['BabA','BabU']:
1604                    hapData['Babinet'][name][0] = parmDict[pfx+name]
1605                    if pfx+name in sigDict:
1606                        BabSig[pfx+name] = sigDict[pfx+name]               
1607               
1608            elif 'HKLF' in histogram:
1609                for item in ['Scale',]:
1610                    if parmDict.get(pfx+item):
1611                        hapData[item][0] = parmDict[pfx+item]
1612                        if pfx+item in sigDict:
1613                            ScalExtSig[pfx+item] = sigDict[pfx+item]
1614                for item in ['Ep','Eg','Es']:
1615                    if parmDict.get(pfx+item):
1616                        hapData['Extinction'][2][item][0] = parmDict[pfx+item]
1617                        if pfx+item in sigDict:
1618                            ScalExtSig[pfx+item] = sigDict[pfx+item]
1619                for name in ['BabA','BabU']:
1620                    hapData['Babinet'][name][0] = parmDict[pfx+name]
1621                    if pfx+name in sigDict:
1622                        BabSig[pfx+name] = sigDict[pfx+name]               
1623
1624    if Print:
1625        for phase in Phases:
1626            HistoPhase = Phases[phase]['Histograms']
1627            General = Phases[phase]['General']
1628            SGData = General['SGData']
1629            pId = Phases[phase]['pId']
1630            histoList = HistoPhase.keys()
1631            histoList.sort()
1632            for histogram in histoList:
1633                try:
1634                    Histogram = Histograms[histogram]
1635                except KeyError:                       
1636                    #skip if histogram not included e.g. in a sequential refinement
1637                    continue
1638                print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
1639                print >>pFile,130*'-'
1640                hapData = HistoPhase[histogram]
1641                hId = Histogram['hId']
1642                pfx = str(pId)+':'+str(hId)+':'
1643                if 'PWDR' in histogram:
1644                    print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections'   \
1645                        %(Histogram[pfx+'Rf'],Histogram[pfx+'Rf^2'],Histogram[pfx+'Nref'])
1646               
1647                    if pfx+'Scale' in PhFrExtPOSig:
1648                        wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum[hId]
1649                        sigwtFr = PhFrExtPOSig[pfx+'Scale']*wtFr/hapData['Scale'][0]
1650                        print >>pFile,' Phase fraction  : %10.5f, sig %10.5f Weight fraction  : %8.5f, sig %10.5f' \
1651                            %(hapData['Scale'][0],PhFrExtPOSig[pfx+'Scale'],wtFr,sigwtFr)
1652                    if pfx+'Extinction' in PhFrExtPOSig:
1653                        print >>pFile,' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction'])
1654                    if hapData['Pref.Ori.'][0] == 'MD':
1655                        if pfx+'MD' in PhFrExtPOSig:
1656                            print >>pFile,' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD'])
1657                    else:
1658                        PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig)
1659                    PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size'])
1660                    PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
1661                    PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
1662                    PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
1663                   
1664                elif 'HKLF' in histogram:
1665                    print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections'   \
1666                        %(Histogram[pfx+'Rf'],Histogram[pfx+'Rf^2'],Histogram[pfx+'Nref'])
1667                    print >>pFile,' HKLF histogram weight factor = ','%.3f'%(Histogram['wtFactor'])
1668                    if pfx+'Scale' in ScalExtSig:
1669                        print >>pFile,' Scale factor : %10.4f, sig %10.4f'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale'])
1670                    if hapData['Extinction'][0] != 'None':
1671                        PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig)
1672                    if len(BabSig):
1673                        PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
1674
1675################################################################################
1676##### Histogram data
1677################################################################################       
1678                   
1679def GetHistogramData(Histograms,Print=True,pFile=None):
1680   
1681    def GetBackgroundParms(hId,Background):
1682        Back = Background[0]
1683        DebyePeaks = Background[1]
1684        bakType,bakFlag = Back[:2]
1685        backVals = Back[3:]
1686        backNames = [':'+str(hId)+':Back:'+str(i) for i in range(len(backVals))]
1687        backDict = dict(zip(backNames,backVals))
1688        backVary = []
1689        if bakFlag:
1690            backVary = backNames
1691        backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
1692        backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
1693        debyeDict = {}
1694        debyeList = []
1695        for i in range(DebyePeaks['nDebye']):
1696            debyeNames = [':'+str(hId)+':DebyeA:'+str(i),':'+str(hId)+':DebyeR:'+str(i),':'+str(hId)+':DebyeU:'+str(i)]
1697            debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
1698            debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
1699        debyeVary = []
1700        for item in debyeList:
1701            if item[1]:
1702                debyeVary.append(item[0])
1703        backDict.update(debyeDict)
1704        backVary += debyeVary
1705        peakDict = {}
1706        peakList = []
1707        for i in range(DebyePeaks['nPeaks']):
1708            peakNames = [':'+str(hId)+':BkPkpos:'+str(i),':'+str(hId)+ \
1709                ':BkPkint:'+str(i),':'+str(hId)+':BkPksig:'+str(i),':'+str(hId)+':BkPkgam:'+str(i)]
1710            peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
1711            peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
1712        peakVary = []
1713        for item in peakList:
1714            if item[1]:
1715                peakVary.append(item[0])
1716        backDict.update(peakDict)
1717        backVary += peakVary
1718        return bakType,backDict,backVary           
1719       
1720    def GetInstParms(hId,Inst):     
1721        dataType = Inst['Type'][0]
1722        instDict = {}
1723        insVary = []
1724        pfx = ':'+str(hId)+':'
1725        insKeys = Inst.keys()
1726        insKeys.sort()
1727        for item in insKeys:
1728            insName = pfx+item
1729            instDict[insName] = Inst[item][1]
1730            if Inst[item][2]:
1731                insVary.append(insName)
1732#        instDict[pfx+'X'] = max(instDict[pfx+'X'],0.001)
1733#        instDict[pfx+'Y'] = max(instDict[pfx+'Y'],0.001)
1734        instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
1735        return dataType,instDict,insVary
1736       
1737    def GetSampleParms(hId,Sample):
1738        sampVary = []
1739        hfx = ':'+str(hId)+':'       
1740        sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
1741            hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
1742        Type = Sample['Type']
1743        if 'Bragg' in Type:             #Bragg-Brentano
1744            for item in ['Scale','Shift','Transparency']:       #surface roughness?, diffuse scattering?
1745                sampDict[hfx+item] = Sample[item][0]
1746                if Sample[item][1]:
1747                    sampVary.append(hfx+item)
1748        elif 'Debye' in Type:        #Debye-Scherrer
1749            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
1750                sampDict[hfx+item] = Sample[item][0]
1751                if Sample[item][1]:
1752                    sampVary.append(hfx+item)
1753        return Type,sampDict,sampVary
1754       
1755    def PrintBackground(Background):
1756        Back = Background[0]
1757        DebyePeaks = Background[1]
1758        print >>pFile,'\n Background function: ',Back[0],' Refine?',bool(Back[1])
1759        line = ' Coefficients: '
1760        for i,back in enumerate(Back[3:]):
1761            line += '%10.3f'%(back)
1762            if i and not i%10:
1763                line += '\n'+15*' '
1764        print >>pFile,line
1765        if DebyePeaks['nDebye']:
1766            print >>pFile,'\n Debye diffuse scattering coefficients'
1767            parms = ['DebyeA','DebyeR','DebyeU']
1768            line = ' names :  '
1769            for parm in parms:
1770                line += '%8s refine?'%(parm)
1771            print >>pFile,line
1772            for j,term in enumerate(DebyePeaks['debyeTerms']):
1773                line = ' term'+'%2d'%(j)+':'
1774                for i in range(3):
1775                    line += '%10.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
1776                print >>pFile,line
1777        if DebyePeaks['nPeaks']:
1778            print >>pFile,'\n Single peak coefficients'
1779            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
1780            line = ' names :  '
1781            for parm in parms:
1782                line += '%8s refine?'%(parm)
1783            print >>pFile,line
1784            for j,term in enumerate(DebyePeaks['peaksList']):
1785                line = ' peak'+'%2d'%(j)+':'
1786                for i in range(4):
1787                    line += '%10.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
1788                print >>pFile,line
1789       
1790    def PrintInstParms(Inst):
1791        print >>pFile,'\n Instrument Parameters:'
1792        ptlbls = ' name  :'
1793        ptstr =  ' value :'
1794        varstr = ' refine:'
1795        insKeys = Inst.keys()
1796        insKeys.sort()
1797        for item in insKeys:
1798            if item != 'Type':
1799                ptlbls += '%12s' % (item)
1800                ptstr += '%12.6f' % (Inst[item][1])
1801                if item in ['Lam1','Lam2','Azimuth']:
1802                    varstr += 12*' '
1803                else:
1804                    varstr += '%12s' % (str(bool(Inst[item][2])))
1805        print >>pFile,ptlbls
1806        print >>pFile,ptstr
1807        print >>pFile,varstr
1808       
1809    def PrintSampleParms(Sample):
1810        print >>pFile,'\n Sample Parameters:'
1811        print >>pFile,' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \
1812            (Sample['Omega'],Sample['Chi'],Sample['Phi'])
1813        ptlbls = ' name  :'
1814        ptstr =  ' value :'
1815        varstr = ' refine:'
1816        if 'Bragg' in Sample['Type']:
1817            for item in ['Scale','Shift','Transparency']:
1818                ptlbls += '%14s'%(item)
1819                ptstr += '%14.4f'%(Sample[item][0])
1820                varstr += '%14s'%(str(bool(Sample[item][1])))
1821           
1822        elif 'Debye' in Type:        #Debye-Scherrer
1823            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
1824                ptlbls += '%14s'%(item)
1825                ptstr += '%14.4f'%(Sample[item][0])
1826                varstr += '%14s'%(str(bool(Sample[item][1])))
1827
1828        print >>pFile,ptlbls
1829        print >>pFile,ptstr
1830        print >>pFile,varstr
1831       
1832    histDict = {}
1833    histVary = []
1834    controlDict = {}
1835    histoList = Histograms.keys()
1836    histoList.sort()
1837    for histogram in histoList:
1838        if 'PWDR' in histogram:
1839            Histogram = Histograms[histogram]
1840            hId = Histogram['hId']
1841            pfx = ':'+str(hId)+':'
1842            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
1843            controlDict[pfx+'Limits'] = Histogram['Limits'][1]
1844           
1845            Background = Histogram['Background']
1846            Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
1847            controlDict[pfx+'bakType'] = Type
1848            histDict.update(bakDict)
1849            histVary += bakVary
1850           
1851            Inst = Histogram['Instrument Parameters'][0]
1852            Type,instDict,insVary = GetInstParms(hId,Inst)
1853            controlDict[pfx+'histType'] = Type
1854            if pfx+'Lam1' in instDict:
1855                controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
1856            else:
1857                controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']           
1858            histDict.update(instDict)
1859            histVary += insVary
1860           
1861            Sample = Histogram['Sample Parameters']
1862            Type,sampDict,sampVary = GetSampleParms(hId,Sample)
1863            controlDict[pfx+'instType'] = Type
1864            histDict.update(sampDict)
1865            histVary += sampVary
1866           
1867   
1868            if Print: 
1869                print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
1870                print >>pFile,135*'-'
1871                Units = {'C':' deg','T':' msec'}
1872                units = Units[controlDict[pfx+'histType'][2]]
1873                Limits = controlDict[pfx+'Limits']
1874                print >>pFile,' Instrument type: ',Sample['Type']
1875                print >>pFile,' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)     
1876                PrintSampleParms(Sample)
1877                PrintInstParms(Inst)
1878                PrintBackground(Background)
1879        elif 'HKLF' in histogram:
1880            Histogram = Histograms[histogram]
1881            hId = Histogram['hId']
1882            pfx = ':'+str(hId)+':'
1883            controlDict[pfx+'wtFactor'] =Histogram['wtFactor']
1884            Inst = Histogram['Instrument Parameters'][0]
1885            controlDict[pfx+'histType'] = Inst['Type'][0]
1886            histDict[pfx+'Lam'] = Inst['Lam'][1]
1887            controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']                   
1888    return histVary,histDict,controlDict
1889   
1890def SetHistogramData(parmDict,sigDict,Histograms,Print=True,pFile=None):
1891   
1892    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
1893        Back = Background[0]
1894        DebyePeaks = Background[1]
1895        lenBack = len(Back[3:])
1896        backSig = [0 for i in range(lenBack+3*DebyePeaks['nDebye']+4*DebyePeaks['nPeaks'])]
1897        for i in range(lenBack):
1898            Back[3+i] = parmDict[pfx+'Back:'+str(i)]
1899            if pfx+'Back:'+str(i) in sigDict:
1900                backSig[i] = sigDict[pfx+'Back:'+str(i)]
1901        if DebyePeaks['nDebye']:
1902            for i in range(DebyePeaks['nDebye']):
1903                names = [pfx+'DebyeA:'+str(i),pfx+'DebyeR:'+str(i),pfx+'DebyeU:'+str(i)]
1904                for j,name in enumerate(names):
1905                    DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
1906                    if name in sigDict:
1907                        backSig[lenBack+3*i+j] = sigDict[name]           
1908        if DebyePeaks['nPeaks']:
1909            for i in range(DebyePeaks['nPeaks']):
1910                names = [pfx+'BkPkpos:'+str(i),pfx+'BkPkint:'+str(i),
1911                    pfx+'BkPksig:'+str(i),pfx+'BkPkgam:'+str(i)]
1912                for j,name in enumerate(names):
1913                    DebyePeaks['peaksList'][i][2*j] = parmDict[name]
1914                    if name in sigDict:
1915                        backSig[lenBack+3*DebyePeaks['nDebye']+4*i+j] = sigDict[name]
1916        return backSig
1917       
1918    def SetInstParms(pfx,Inst,parmDict,sigDict):
1919        instSig = {}
1920        insKeys = Inst.keys()
1921        insKeys.sort()
1922        for item in insKeys:
1923            insName = pfx+item
1924            Inst[item][1] = parmDict[insName]
1925            if insName in sigDict:
1926                instSig[item] = sigDict[insName]
1927            else:
1928                instSig[item] = 0
1929        return instSig
1930       
1931    def SetSampleParms(pfx,Sample,parmDict,sigDict):
1932        if 'Bragg' in Sample['Type']:             #Bragg-Brentano
1933            sampSig = [0 for i in range(3)]
1934            for i,item in enumerate(['Scale','Shift','Transparency']):       #surface roughness?, diffuse scattering?
1935                Sample[item][0] = parmDict[pfx+item]
1936                if pfx+item in sigDict:
1937                    sampSig[i] = sigDict[pfx+item]
1938        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
1939            sampSig = [0 for i in range(4)]
1940            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
1941                Sample[item][0] = parmDict[pfx+item]
1942                if pfx+item in sigDict:
1943                    sampSig[i] = sigDict[pfx+item]
1944        return sampSig
1945       
1946    def PrintBackgroundSig(Background,backSig):
1947        Back = Background[0]
1948        DebyePeaks = Background[1]
1949        lenBack = len(Back[3:])
1950        valstr = ' value : '
1951        sigstr = ' sig   : '
1952        refine = False
1953        for i,back in enumerate(Back[3:]):
1954            valstr += '%10.4g'%(back)
1955            if Back[1]:
1956                refine = True
1957                sigstr += '%10.4g'%(backSig[i])
1958            else:
1959                sigstr += 10*' '
1960        if refine:
1961            print >>pFile,'\n Background function: ',Back[0]
1962            print >>pFile,valstr
1963            print >>pFile,sigstr
1964        if DebyePeaks['nDebye']:
1965            ifAny = False
1966            ptfmt = "%12.3f"
1967            names =  ' names :'
1968            ptstr =  ' values:'
1969            sigstr = ' esds  :'
1970            for item in sigDict:
1971                if 'Debye' in item:
1972                    ifAny = True
1973                    names += '%12s'%(item)
1974                    ptstr += ptfmt%(parmDict[item])
1975                    sigstr += ptfmt%(sigDict[item])
1976            if ifAny:
1977                print >>pFile,'\n Debye diffuse scattering coefficients'
1978                print >>pFile,names
1979                print >>pFile,ptstr
1980                print >>pFile,sigstr
1981        if DebyePeaks['nPeaks']:
1982            ifAny = False
1983            ptfmt = "%14.3f"
1984            names =  ' names :'
1985            ptstr =  ' values:'
1986            sigstr = ' esds  :'
1987            for item in sigDict:
1988                if 'BkPk' in item:
1989                    ifAny = True
1990                    names += '%14s'%(item)
1991                    ptstr += ptfmt%(parmDict[item])
1992                    sigstr += ptfmt%(sigDict[item])
1993            if ifAny:
1994                print >>pFile,'\n Single peak coefficients'
1995                print >>pFile,names
1996                print >>pFile,ptstr
1997                print >>pFile,sigstr
1998       
1999    def PrintInstParmsSig(Inst,instSig):
2000        ptlbls = ' names :'
2001        ptstr =  ' value :'
2002        sigstr = ' sig   :'
2003        refine = False
2004        insKeys = instSig.keys()
2005        insKeys.sort()
2006        for name in insKeys:
2007            if name not in  ['Type','Lam1','Lam2','Azimuth']:
2008                ptlbls += '%12s' % (name)
2009                ptstr += '%12.6f' % (Inst[name][1])
2010                if instSig[name]:
2011                    refine = True
2012                    sigstr += '%12.6f' % (instSig[name])
2013                else:
2014                    sigstr += 12*' '
2015        if refine:
2016            print >>pFile,'\n Instrument Parameters:'
2017            print >>pFile,ptlbls
2018            print >>pFile,ptstr
2019            print >>pFile,sigstr
2020       
2021    def PrintSampleParmsSig(Sample,sampleSig):
2022        ptlbls = ' names :'
2023        ptstr =  ' values:'
2024        sigstr = ' sig   :'
2025        refine = False
2026        if 'Bragg' in Sample['Type']:
2027            for i,item in enumerate(['Scale','Shift','Transparency']):
2028                ptlbls += '%14s'%(item)
2029                ptstr += '%14.4f'%(Sample[item][0])
2030                if sampleSig[i]:
2031                    refine = True
2032                    sigstr += '%14.4f'%(sampleSig[i])
2033                else:
2034                    sigstr += 14*' '
2035           
2036        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
2037            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
2038                ptlbls += '%14s'%(item)
2039                ptstr += '%14.4f'%(Sample[item][0])
2040                if sampleSig[i]:
2041                    refine = True
2042                    sigstr += '%14.4f'%(sampleSig[i])
2043                else:
2044                    sigstr += 14*' '
2045
2046        if refine:
2047            print >>pFile,'\n Sample Parameters:'
2048            print >>pFile,ptlbls
2049            print >>pFile,ptstr
2050            print >>pFile,sigstr
2051       
2052    histoList = Histograms.keys()
2053    histoList.sort()
2054    for histogram in histoList:
2055        if 'PWDR' in histogram:
2056            Histogram = Histograms[histogram]
2057            hId = Histogram['hId']
2058            pfx = ':'+str(hId)+':'
2059            Background = Histogram['Background']
2060            backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
2061           
2062            Inst = Histogram['Instrument Parameters'][0]
2063            instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
2064       
2065            Sample = Histogram['Sample Parameters']
2066            sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
2067
2068            print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
2069            print >>pFile,135*'-'
2070            print >>pFile,' Final refinement wR = %.2f%% on %d observations in this histogram'%(Histogram['wR'],Histogram['Nobs'])
2071            print >>pFile,' PWDR histogram weight factor = '+'%.3f'%(Histogram['wtFactor'])
2072            if Print:
2073                print >>pFile,' Instrument type: ',Sample['Type']
2074                PrintSampleParmsSig(Sample,sampSig)
2075                PrintInstParmsSig(Inst,instSig)
2076                PrintBackgroundSig(Background,backSig)
2077               
2078################################################################################
2079##### Penalty & restraint functions
2080################################################################################
2081
2082def penaltyFxn(HistoPhases,parmDict,varyList):
2083    Histograms,Phases,restraintDict = HistoPhases
2084    pNames = []
2085    pVals = []
2086    pWt = []
2087    negWt = {}
2088    for phase in Phases:
2089        pId = Phases[phase]['pId']
2090        negWt[pId] = Phases[phase]['General']['Pawley neg wt']
2091        General = Phases[phase]['General']
2092        textureData = General['SH Texture']
2093        SGData = General['SGData']
2094        AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'])
2095        cell = General['Cell'][1:7]
2096        Amat,Bmat = G2lat.cell2AB(cell)
2097        if phase not in restraintDict:
2098            continue
2099        phaseRest = restraintDict[phase]
2100        names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
2101            ['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
2102            ['ChemComp','Sites'],['Texture','HKLs']]
2103        for name,rest in names:
2104            itemRest = phaseRest[name]
2105            if itemRest[rest] and itemRest['Use']:
2106                wt = itemRest['wtFactor']
2107                if name in ['Bond','Angle','Plane','Chiral']:
2108                    for i,[indx,ops,obs,esd] in enumerate(itemRest[rest]):
2109                        pNames.append(str(pId)+':'+name+':'+str(i))
2110                        XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
2111                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
2112                        if name == 'Bond':
2113                            calc = G2mth.getRestDist(XYZ,Amat)
2114                        elif name == 'Angle':
2115                            calc = G2mth.getRestAngle(XYZ,Amat)
2116                        elif name == 'Plane':
2117                            calc = G2mth.getRestPlane(XYZ,Amat)
2118                        elif name == 'Chiral':
2119                            calc = G2mth.getRestChiral(XYZ,Amat)
2120                        pVals.append(obs-calc)
2121                        pWt.append(wt/esd**2)
2122                elif name in ['Torsion','Rama']:
2123                    coeffDict = itemRest['Coeff']
2124                    for i,[indx,ops,cofName,esd] in enumerate(torsionList):
2125                        pNames.append(str(pId)+':'+name+':'+str(i))
2126                        XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
2127                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
2128                        if name == 'Torsion':
2129                            tor = G2mth.getRestTorsion(XYZ,Amat)
2130                            restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
2131                        else:
2132                            phi,psi = G2mth.getRestRama(XYZ,Amat)
2133                            restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])                               
2134                        pVals.append(obs-calc)
2135                        pWt.append(wt/esd**2)
2136                elif name == 'ChemComp':
2137                    for i,[indx,factors,obs,esd] in enumerate(itemRest[rest]):
2138                        pNames.append(str(pId)+':'+name+':'+str(i))
2139                        mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs+1))
2140                        frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs-1))
2141                        calc = np.sum(mul*frac*factors)
2142                        pVals.append(obs-calc)
2143                        pWt.append(wt/esd**2)                   
2144                elif name == 'Texture':
2145                    SHkeys = textureData['SH Coeff'][1].keys()
2146                    SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
2147                    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
2148                    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
2149                    for i,[hkl,grid,esd1,ifesd2,esd2] in enumerate(itemRest[rest]):
2150                        PH = np.array(hkl)
2151                        phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
2152                        ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
2153                        R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
2154                        Z1 = -ma.masked_greater(Z,0.0)
2155                        IndZ1 = np.array(ma.nonzero(Z1))
2156                        for ind in IndZ1.T:
2157                            pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name,i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
2158                            pVals.append(Z1[ind[0]][ind[1]])
2159                            pWt.append(wt/esd1**2)
2160                        if ifesd2:
2161                            Z2 = 1.-Z
2162                            for ind in np.ndindex(grid,grid):
2163                                pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name+'-unit',i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
2164                                pVals.append(Z1[ind[0]][ind[1]])
2165                                pWt.append(wt/esd2**2)
2166         
2167    for item in varyList:
2168        if 'PWLref' in item and parmDict[item] < 0.:
2169            pId = int(item.split(':')[0])
2170            if negWt[pId]:
2171                pNames.append(item)
2172                pVals.append(-parmDict[item])
2173                pWt.append(negWt[pId])
2174    pVals = np.array(pVals)
2175    pWt = np.array(pWt)         #should this be np.sqrt?
2176    return pNames,pVals,pWt
2177   
2178def penaltyDeriv(pNames,pVal,HistoPhases,parmDict,varyList):
2179    Histograms,Phases,restraintDict = HistoPhases
2180    pDerv = np.zeros((len(varyList),len(pVal)))
2181    for phase in Phases:
2182        if phase not in restraintDict:
2183            continue
2184        pId = Phases[phase]['pId']
2185        General = Phases[phase]['General']
2186        SGData = General['SGData']
2187        AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'])
2188        cell = General['Cell'][1:7]
2189        Amat,Bmat = G2lat.cell2AB(cell)
2190        textureData = General['SH Texture']
2191
2192        SHkeys = textureData['SH Coeff'][1].keys()
2193        SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
2194        shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
2195        SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
2196        sam = SamSym[textureData['Model']]
2197        phaseRest = restraintDict[phase]
2198        names = {'Bond':'Bonds','Angle':'Angles','Plane':'Planes',
2199            'Chiral':'Volumes','Torsion':'Torsions','Rama':'Ramas',
2200            'ChemComp':'Sites','Texture':'HKLs'}
2201        lasthkl = np.array([0,0,0])
2202        for ip,pName in enumerate(pNames):
2203            pnames = pName.split(':')
2204            if pId == int(pnames[0]):
2205                name = pnames[1]
2206                if not name:        #empty for Pawley restraints; pName has '::' in it
2207                    pDerv[varyList.index(pName)][ip] += 1.
2208                    continue
2209                id = int(pnames[2]) 
2210                itemRest = phaseRest[name]
2211                if name in ['Bond','Angle','Plane','Chiral']:
2212                    indx,ops,obs,esd = itemRest[names[name]][id]
2213                    dNames = []
2214                    for ind in indx:
2215                        dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
2216                    XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
2217                    if name == 'Bond':
2218                        deriv = G2mth.getRestDeriv(G2mth.getRestDist,XYZ,Amat,ops,SGData)
2219                    elif name == 'Angle':
2220                        deriv = G2mth.getRestDeriv(G2mth.getRestAngle,XYZ,Amat,ops,SGData)
2221                    elif name == 'Plane':
2222                        deriv = G2mth.getRestDeriv(G2mth.getRestPlane,XYZ,Amat,ops,SGData)
2223                    elif name == 'Chiral':
2224                        deriv = G2mth.getRestDeriv(G2mth.getRestChiral,XYZ,Amat,ops,SGData)
2225                elif name in ['Torsion','Rama']:
2226                    coffDict = itemRest['Coeff']
2227                    indx,ops,cofName,esd = itemRest[names[name]][id]
2228                    dNames = []
2229                    for ind in indx:
2230                        dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
2231                    XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
2232                    if name == 'Torsion':
2233                        deriv = G2mth.getTorsionDeriv(XYZ,Amat,coffDict[cofName])
2234                    else:
2235                        deriv = G2mth.getRamaDeriv(XYZ,Amat,coffDict[cofName])
2236                elif name == 'ChemComp':
2237                    indx,factors,obs,esd = itemRest[names[name]][id]
2238                    dNames = []
2239                    for ind in indx:
2240                        dNames += [str(pId)+'::Afrac:'+str(AtLookup[ind])]
2241                        mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs+1))
2242                        deriv = mul*factors
2243                elif 'Texture' in name:
2244                    deriv = []
2245                    dNames = []
2246                    hkl,grid,esd1,ifesd2,esd2 = itemRest[names[name]][id]
2247                    hkl = np.array(hkl)
2248                    if np.any(lasthkl-hkl):
2249                        PH = np.array(hkl)
2250                        phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
2251                        ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
2252                        lasthkl = copy.copy(hkl)                       
2253                    if 'unit' in name:
2254                        pass
2255                    else:
2256                        gam = float(pnames[3])
2257                        psi = float(pnames[4])
2258                        for SHname in ODFln:
2259                            l,m,n = eval(SHname[1:])
2260                            Ksl = G2lat.GetKsl(l,m,sam,psi,gam)[0]
2261                            dNames += [str(pId)+'::'+SHname]
2262                            deriv.append(-ODFln[SHname][0]*Ksl/SHCoef[SHname])
2263                for dName,drv in zip(dNames,deriv):
2264                    try:
2265                        ind = varyList.index(dName)
2266                        pDerv[ind][ip] += drv
2267                    except ValueError:
2268                        pass
2269    return pDerv
2270
2271################################################################################
2272##### Function & derivative calculations
2273################################################################################       
2274                   
2275def GetAtomFXU(pfx,calcControls,parmDict):
2276    Natoms = calcControls['Natoms'][pfx]
2277    Tdata = Natoms*[' ',]
2278    Mdata = np.zeros(Natoms)
2279    IAdata = Natoms*[' ',]
2280    Fdata = np.zeros(Natoms)
2281    FFdata = []
2282    BLdata = []
2283    Xdata = np.zeros((3,Natoms))
2284    dXdata = np.zeros((3,Natoms))
2285    Uisodata = np.zeros(Natoms)
2286    Uijdata = np.zeros((6,Natoms))
2287    keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata,
2288        'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2],
2289        'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata,
2290        'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2],
2291        'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5]}
2292    for iatm in range(Natoms):
2293        for key in keys:
2294            parm = pfx+key+str(iatm)
2295            if parm in parmDict:
2296                keys[key][iatm] = parmDict[parm]
2297    return Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata
2298   
2299def getFFvalues(FFtables,SQ):
2300    FFvals = {}
2301    for El in FFtables:
2302        FFvals[El] = G2el.ScatFac(FFtables[El],SQ)[0]
2303    return FFvals
2304   
2305def getBLvalues(BLtables):
2306    BLvals = {}
2307    for El in BLtables:
2308        BLvals[El] = BLtables[El][1][1]
2309    return BLvals
2310       
2311def StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict):
2312    ''' Compute structure factors for all h,k,l for phase
2313    input:
2314        refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv]
2315        G:      reciprocal metric tensor
2316        pfx:    phase id string
2317        SGData: space group info. dictionary output from SpcGroup
2318        calcControls:
2319        ParmDict:
2320    puts result F^2 in each ref[8] in refList
2321    '''       
2322    twopi = 2.0*np.pi
2323    twopisq = 2.0*np.pi**2
2324    phfx = pfx.split(':')[0]+hfx
2325    ast = np.sqrt(np.diag(G))
2326    Mast = twopisq*np.multiply.outer(ast,ast)
2327    FFtables = calcControls['FFtables']
2328    BLtables = calcControls['BLtables']
2329    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
2330    FF = np.zeros(len(Tdata))
2331    if 'N' in parmDict[hfx+'Type']:
2332        FP,FPP = G2el.BlenRes(Tdata,BLtables,parmDict[hfx+'Lam'])
2333    else:
2334        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
2335        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
2336    maxPos = len(SGData['SGOps'])
2337    Uij = np.array(G2lat.U6toUij(Uijdata))
2338    bij = Mast*Uij.T
2339    for refl in refList:
2340        fbs = np.array([0,0])
2341        H = refl[:3]
2342        SQ = 1./(2.*refl[4])**2
2343        SQfactor = 4.0*SQ*twopisq
2344        Bab = parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor)
2345        if not len(refl[-1]):                #no form factors
2346            if 'N' in parmDict[hfx+'Type']:
2347                refl[-1] = getBLvalues(BLtables)
2348            else:       #'X'
2349                refl[-1] = getFFvalues(FFtables,SQ)
2350        for i,El in enumerate(Tdata):
2351            FF[i] = refl[-1][El]           
2352        Uniq = refl[11]
2353        phi = refl[12]
2354        phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+phi[:,np.newaxis])
2355        sinp = np.sin(phase)
2356        cosp = np.cos(phase)
2357        occ = Mdata*Fdata/len(Uniq)
2358        biso = -SQfactor*Uisodata
2359        Tiso = np.where(biso<1.,np.exp(biso),1.0)
2360        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
2361        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
2362        Tcorr = Tiso*Tuij
2363        fa = np.array([(FF+FP-Bab)*occ*cosp*Tcorr,-FPP*occ*sinp*Tcorr])
2364        fas = np.sum(np.sum(fa,axis=1),axis=1)        #real
2365        if not SGData['SGInv']:
2366            fb = np.array([(FF+FP-Bab)*occ*sinp*Tcorr,FPP*occ*cosp*Tcorr])
2367            fbs = np.sum(np.sum(fb,axis=1),axis=1)
2368        fasq = fas**2
2369        fbsq = fbs**2        #imaginary
2370        refl[9] = np.sum(fasq)+np.sum(fbsq)
2371        refl[10] = atan2d(fbs[0],fas[0])
2372    return refList
2373   
2374def StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict):
2375    twopi = 2.0*np.pi
2376    twopisq = 2.0*np.pi**2
2377    phfx = pfx.split(':')[0]+hfx
2378    ast = np.sqrt(np.diag(G))
2379    Mast = twopisq*np.multiply.outer(ast,ast)
2380    FFtables = calcControls['FFtables']
2381    BLtables = calcControls['BLtables']
2382    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
2383    FF = np.zeros(len(Tdata))
2384    if 'N' in parmDict[hfx+'Type']:
2385        FP = 0.
2386        FPP = 0.
2387    else:
2388        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
2389        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
2390    maxPos = len(SGData['SGOps'])       
2391    Uij = np.array(G2lat.U6toUij(Uijdata))
2392    bij = Mast*Uij.T
2393    dFdvDict = {}
2394    dFdfr = np.zeros((len(refList),len(Mdata)))
2395    dFdx = np.zeros((len(refList),len(Mdata),3))
2396    dFdui = np.zeros((len(refList),len(Mdata)))
2397    dFdua = np.zeros((len(refList),len(Mdata),6))
2398    dFdbab = np.zeros((len(refList),2))
2399    for iref,refl in enumerate(refList):
2400        H = np.array(refl[:3])
2401        SQ = 1./(2.*refl[4])**2             # or (sin(theta)/lambda)**2
2402        SQfactor = 8.0*SQ*np.pi**2
2403        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
2404        Bab = parmDict[phfx+'BabA']*dBabdA
2405        for i,El in enumerate(Tdata):           
2406            FF[i] = refl[-1][El]           
2407        Uniq = refl[11]
2408        phi = refl[12]
2409        phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+phi[np.newaxis,:])
2410        sinp = np.sin(phase)
2411        cosp = np.cos(phase)
2412        occ = Mdata*Fdata/len(Uniq)
2413        biso = -SQfactor*Uisodata
2414        Tiso = np.where(biso<1.,np.exp(biso),1.0)
2415        HbH = -np.inner(H,np.inner(bij,H))
2416        Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
2417        Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
2418        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
2419        Tcorr = Tiso*Tuij
2420        fot = (FF+FP-Bab)*occ*Tcorr
2421        fotp = FPP*occ*Tcorr
2422        fa = np.array([fot[:,np.newaxis]*cosp,fotp[:,np.newaxis]*cosp])       #non positions
2423        fb = np.array([fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])
2424       
2425        fas = np.sum(np.sum(fa,axis=1),axis=1)
2426        fbs = np.sum(np.sum(fb,axis=1),axis=1)
2427        fax = np.array([-fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])   #positions
2428        fbx = np.array([fot[:,np.newaxis]*cosp,-fot[:,np.newaxis]*cosp])
2429        #sum below is over Uniq
2430        dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)
2431        dfadx = np.sum(twopi*Uniq*fax[:,:,:,np.newaxis],axis=2)
2432        dfadui = np.sum(-SQfactor*fa,axis=2)
2433        dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
2434        dfadba = np.sum(-cosp*(occ*Tcorr)[:,np.newaxis],axis=1)
2435        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for     
2436        dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)
2437        dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])
2438        dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])
2439        dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])
2440        dFdbab[iref] = np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
2441        if not SGData['SGInv']:
2442            dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2)        #problem here if occ=0 for some atom
2443            dfbdx = np.sum(twopi*Uniq*fbx[:,:,:,np.newaxis],axis=2)         
2444            dfbdui = np.sum(-SQfactor*fb,axis=2)
2445            dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
2446            dfbdba = np.sum(-sinp*(occ*Tcorr)[:,np.newaxis],axis=1)
2447            dFdfr[iref] += 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
2448            dFdx[iref] += 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
2449            dFdui[iref] += 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
2450            dFdua[iref] += 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
2451            dFdbab[iref] += np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
2452        #loop over atoms - each dict entry is list of derivatives for all the reflections
2453    for i in range(len(Mdata)):     
2454        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
2455        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
2456        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
2457        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
2458        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
2459        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
2460        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
2461        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
2462        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
2463        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
2464        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
2465        dFdvDict[pfx+'BabA'] = dFdbab.T[0]
2466        dFdvDict[pfx+'BabU'] = dFdbab.T[1]
2467    return dFdvDict
2468   
2469def SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varyList):
2470    ''' Single crystal extinction function; puts correction in ref[13] and returns
2471    corrections needed for derivatives
2472    '''
2473    ref[13] = 1.0
2474    dervCor = 1.0
2475    dervDict = {}
2476    if calcControls[phfx+'EType'] != 'None':
2477        cos2T = 1.0-0.5*(parmDict[hfx+'Lam']/ref[4])**2         #cos(2theta)
2478        if 'SXC' in parmDict[hfx+'Type']:
2479            AV = 7.9406e5/parmDict[pfx+'Vol']**2
2480            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
2481            P12 = (calcControls[phfx+'Cos2TM']+cos2T**4)/(calcControls[phfx+'Cos2TM']+cos2T**2)
2482        elif 'SNT' in parmDict[hfx+'Type']:
2483            AV = 1.e7/parmDict[pfx+'Vol']**2
2484            PL = 1./(4.*refl[4]**2)
2485            P12 = 1.0
2486        elif 'SNC' in parmDict[hfx+'Type']:
2487            AV = 1.e7/parmDict[pfx+'Vol']**2
2488            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
2489            P12 = 1.0
2490           
2491        PLZ = AV*P12*parmDict[hfx+'Lam']**2*ref[7]
2492        if 'Primary' in calcControls[phfx+'EType']:
2493            PLZ *= 1.5
2494        else:
2495            PLZ *= calcControls[phfx+'Tbar']
2496                       
2497        if 'Primary' in calcControls[phfx+'EType']:
2498            PSIG = parmDict[phfx+'Ep']
2499        elif 'I & II' in calcControls[phfx+'EType']:
2500            PSIG = parmDict[phfx+'Eg']/np.sqrt(1.+(parmDict[phfx+'Es']*PL/parmDict[phfx+'Eg'])**2)
2501        elif 'Type II' in calcControls[phfx+'EType']:
2502            PSIG = parmDict[phfx+'Es']
2503        else:       # 'Secondary Type I'
2504            PSIG = parmDict[phfx+'Eg']/PL
2505           
2506        AG = 0.58+0.48*cos2T+0.24*cos2T**2
2507        AL = 0.025+0.285*cos2T
2508        BG = 0.02-0.025*cos2T
2509        BL = 0.15-0.2*(0.75-cos2T)**2
2510        if cos2T < 0.:
2511            BL = -0.45*cos2T
2512        CG = 2.
2513        CL = 2.
2514        PF = PLZ*PSIG
2515       
2516        if 'Gaussian' in calcControls[phfx+'EApprox']:
2517            PF4 = 1.+CG*PF+AG*PF**2/(1.+BG*PF)
2518            extCor = np.sqrt(PF4)
2519            PF3 = 0.5*(CG+2.*AG*PF/(1.+BG*PF)-AG*PF**2*BG/(1.+BG*PF)**2)/(PF4*extCor)
2520        else:
2521            PF4 = 1.+CL*PF+AL*PF**2/(1.+BL*PF)
2522            extCor = np.sqrt(PF4)
2523            PF3 = 0.5*(CL+2.*AL*PF/(1.+BL*PF)-AL*PF**2*BL/(1.+BL*PF)**2)/(PF4*extCor)
2524
2525        dervCor = (1.+PF)*PF3
2526        if 'Primary' in calcControls[phfx+'EType'] and phfx+'Ep' in varyList:
2527            dervDict[phfx+'Ep'] = -ref[7]*PLZ*PF3
2528        if 'II' in calcControls[phfx+'EType'] and phfx+'Es' in varyList:
2529            dervDict[phfx+'Es'] = -ref[7]*PLZ*PF3*(PSIG/parmDict[phfx+'Es'])**3
2530        if 'I' in calcControls[phfx+'EType'] and phfx+'Eg' in varyList:
2531            dervDict[phfx+'Eg'] = -ref[7]*PLZ*PF3*(PSIG/parmDict[phfx+'Eg'])**3*PL**2
2532               
2533        ref[13] = 1./extCor
2534    return dervCor,dervDict
2535       
2536   
2537def Dict2Values(parmdict, varylist):
2538    '''Use before call to leastsq to setup list of values for the parameters
2539    in parmdict, as selected by key in varylist'''
2540    return [parmdict[key] for key in varylist] 
2541   
2542def Values2Dict(parmdict, varylist, values):
2543    ''' Use after call to leastsq to update the parameter dictionary with
2544    values corresponding to keys in varylist'''
2545    parmdict.update(zip(varylist,values))
2546   
2547def GetNewCellParms(parmDict,varyList):
2548    newCellDict = {}
2549    Anames = ['A'+str(i) for i in range(6)]
2550    Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],Anames))
2551    for item in varyList:
2552        keys = item.split(':')
2553        if keys[2] in Ddict:
2554            key = keys[0]+'::'+Ddict[keys[2]]       #key is e.g. '0::A0'
2555            parm = keys[0]+'::'+keys[2]             #parm is e.g. '0::D11'
2556            newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
2557    return newCellDict          # is e.g. {'0::D11':A0+D11}
2558   
2559def ApplyXYZshifts(parmDict,varyList):
2560    ''' takes atom x,y,z shift and applies it to corresponding atom x,y,z value
2561        input:
2562            parmDict - parameter dictionary
2563            varyList - list of variables
2564        returns:
2565            newAtomDict - dictionary of new atomic coordinate names & values;
2566                key is parameter shift name
2567    '''
2568    newAtomDict = {}
2569    for item in parmDict:
2570        if 'dA' in item:
2571            parm = ''.join(item.split('d'))
2572            parmDict[parm] += parmDict[item]
2573            newAtomDict[item] = [parm,parmDict[parm]]
2574    return newAtomDict
2575   
2576def SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict):
2577    #Spherical harmonics texture
2578    IFCoup = 'Bragg' in calcControls[hfx+'instType']
2579    odfCor = 1.0
2580    H = refl[:3]
2581    cell = G2lat.Gmat2cell(g)
2582    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
2583    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
2584    phi,beta = G2lat.CrsAng(H,cell,SGData)
2585    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
2586    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
2587    for item in SHnames:
2588        L,M,N = eval(item.strip('C'))
2589        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2590        Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
2591        Lnorm = G2lat.Lnorm(L)
2592        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
2593    return odfCor
2594   
2595def SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict):
2596    #Spherical harmonics texture derivatives
2597    FORPI = 4.0*np.pi
2598    IFCoup = 'Bragg' in calcControls[hfx+'instType']
2599    odfCor = 1.0
2600    dFdODF = {}
2601    dFdSA = [0,0,0]
2602    H = refl[:3]
2603    cell = G2lat.Gmat2cell(g)
2604    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
2605    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
2606    phi,beta = G2lat.CrsAng(H,cell,SGData)
2607    psi,gam,dPSdA,dGMdA = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup)
2608    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
2609    for item in SHnames:
2610        L,M,N = eval(item.strip('C'))
2611        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2612        Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
2613        Lnorm = G2lat.Lnorm(L)
2614        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
2615        dFdODF[pfx+item] = Lnorm*Kcl*Ksl
2616        for i in range(3):
2617            dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i])
2618    return odfCor,dFdODF,dFdSA
2619   
2620def SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict):
2621    #sphericaql harmonics preferred orientation (cylindrical symmetry only)
2622    odfCor = 1.0
2623    H = refl[:3]
2624    cell = G2lat.Gmat2cell(g)
2625    Sangl = [0.,0.,0.]
2626    if 'Bragg' in calcControls[hfx+'instType']:
2627        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
2628        IFCoup = True
2629    else:
2630        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
2631        IFCoup = False
2632    phi,beta = G2lat.CrsAng(H,cell,SGData)
2633    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
2634    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
2635    for item in SHnames:
2636        L,N = eval(item.strip('C'))
2637        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) 
2638        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
2639    return odfCor
2640   
2641def SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict):
2642    #spherical harmonics preferred orientation derivatives (cylindrical symmetry only)
2643    FORPI = 12.5663706143592
2644    odfCor = 1.0
2645    dFdODF = {}
2646    H = refl[:3]
2647    cell = G2lat.Gmat2cell(g)
2648    Sangl = [0.,0.,0.]
2649    if 'Bragg' in calcControls[hfx+'instType']:
2650        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
2651        IFCoup = True
2652    else:
2653        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
2654        IFCoup = False
2655    phi,beta = G2lat.CrsAng(H,cell,SGData)
2656    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
2657    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
2658    for item in SHnames:
2659        L,N = eval(item.strip('C'))
2660        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) 
2661        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
2662        dFdODF[phfx+item] = Kcsl*Lnorm
2663    return odfCor,dFdODF
2664   
2665def GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
2666    POcorr = 1.0
2667    if calcControls[phfx+'poType'] == 'MD':
2668        MD = parmDict[phfx+'MD']
2669        if MD != 1.0:
2670            MDAxis = calcControls[phfx+'MDAxis']
2671            sumMD = 0
2672            for H in refl[11]:           
2673                cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
2674                A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
2675                sumMD += A**3
2676            POcorr = sumMD/len(refl[11])
2677    else:   #spherical harmonics
2678        if calcControls[phfx+'SHord']:
2679            POcorr = SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict)
2680    return POcorr
2681   
2682def GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
2683    POcorr = 1.0
2684    POderv = {}
2685    if calcControls[phfx+'poType'] == 'MD':
2686        MD = parmDict[phfx+'MD']
2687        MDAxis = calcControls[phfx+'MDAxis']
2688        sumMD = 0
2689        sumdMD = 0
2690        for H in refl[11]:           
2691            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
2692            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
2693            sumMD += A**3
2694            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
2695        POcorr = sumMD/len(refl[11])
2696        POderv[phfx+'MD'] = sumdMD/len(refl[11])
2697    else:   #spherical harmonics
2698        if calcControls[phfx+'SHord']:
2699            POcorr,POderv = SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict)
2700    return POcorr,POderv
2701   
2702def GetAbsorb(refl,hfx,calcControls,parmDict):
2703    if 'Debye' in calcControls[hfx+'instType']:
2704        return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0)
2705    else:
2706        return 1.0
2707   
2708def GetAbsorbDerv(refl,hfx,calcControls,parmDict):
2709    if 'Debye' in calcControls[hfx+'instType']:
2710        return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0)
2711    else:
2712        return 0.0
2713   
2714def GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
2715    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3]               #scale*multiplicity
2716    if 'X' in parmDict[hfx+'Type']:
2717        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
2718    Icorr *= GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
2719    if pfx+'SHorder' in parmDict:
2720        Icorr *= SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict)
2721    Icorr *= GetAbsorb(refl,hfx,calcControls,parmDict)
2722    refl[13] = Icorr       
2723   
2724def GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
2725    dIdsh = 1./parmDict[hfx+'Scale']
2726    dIdsp = 1./parmDict[phfx+'Scale']
2727    if 'X' in parmDict[hfx+'Type']:
2728        pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
2729        dIdPola /= pola
2730    else:       #'N'
2731        dIdPola = 0.0
2732    POcorr,dIdPO = GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
2733    for iPO in dIdPO:
2734        dIdPO[iPO] /= POcorr
2735    dFdODF = {}
2736    dFdSA = [0,0,0]
2737    if pfx+'SHorder' in parmDict:
2738        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict)
2739        for iSH in dFdODF:
2740            dFdODF[iSH] /= odfCor
2741        for i in range(3):
2742            dFdSA[i] /= odfCor
2743    dFdAb = GetAbsorbDerv(refl,hfx,calcControls,parmDict)
2744    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb
2745       
2746def GetSampleSigGam(refl,wave,G,GB,phfx,calcControls,parmDict):
2747    costh = cosd(refl[5]/2.)
2748    #crystallite size
2749    if calcControls[phfx+'SizeType'] == 'isotropic':
2750        Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
2751    elif calcControls[phfx+'SizeType'] == 'uniaxial':
2752        H = np.array(refl[:3])
2753        P = np.array(calcControls[phfx+'SizeAxis'])
2754        cosP,sinP = G2lat.CosSinAngle(H,P,G)
2755        Sgam = (1.8*wave/np.pi)/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a']*costh)
2756        Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
2757    else:           #ellipsoidal crystallites
2758        Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
2759        H = np.array(refl[:3])
2760        lenR = G2pwd.ellipseSize(H,Sij,GB)
2761        Sgam = 1.8*wave/(np.pi*costh*lenR)
2762    #microstrain               
2763    if calcControls[phfx+'MustrainType'] == 'isotropic':
2764        Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5]/2.)/np.pi
2765    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2766        H = np.array(refl[:3])
2767        P = np.array(calcControls[phfx+'MustrainAxis'])
2768        cosP,sinP = G2lat.CosSinAngle(H,P,G)
2769        Si = parmDict[phfx+'Mustrain;i']
2770        Sa = parmDict[phfx+'Mustrain;a']
2771        Mgam = 0.018*Si*Sa*tand(refl[5]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
2772    else:       #generalized - P.W. Stephens model
2773        pwrs = calcControls[phfx+'MuPwrs']
2774        sum = 0
2775        for i,pwr in enumerate(pwrs):
2776            sum += parmDict[phfx+'Mustrain:'+str(i)]*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
2777        Mgam = 0.018*refl[4]**2*tand(refl[5]/2.)*sum
2778    gam = Sgam*parmDict[phfx+'Size;mx']+Mgam*parmDict[phfx+'Mustrain;mx']
2779    sig = (Sgam*(1.-parmDict[phfx+'Size;mx']))**2+(Mgam*(1.-parmDict[phfx+'Mustrain;mx']))**2
2780    sig /= ateln2
2781    return sig,gam
2782       
2783def GetSampleSigGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict):
2784    gamDict = {}
2785    sigDict = {}
2786    costh = cosd(refl[5]/2.)
2787    tanth = tand(refl[5]/2.)
2788    #crystallite size derivatives
2789    if calcControls[phfx+'SizeType'] == 'isotropic':
2790        Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
2791        gamDict[phfx+'Size;i'] = -1.8*wave*parmDict[phfx+'Size;mx']/(np.pi*costh)
2792        sigDict[phfx+'Size;i'] = -3.6*Sgam*wave*(1.-parmDict[phfx+'Size;mx'])**2/(np.pi*costh*ateln2)
2793    elif calcControls[phfx+'SizeType'] == 'uniaxial':
2794        H = np.array(refl[:3])
2795        P = np.array(calcControls[phfx+'SizeAxis'])
2796        cosP,sinP = G2lat.CosSinAngle(H,P,G)
2797        Si = parmDict[phfx+'Size;i']
2798        Sa = parmDict[phfx+'Size;a']
2799        gami = (1.8*wave/np.pi)/(Si*Sa)
2800        sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
2801        Sgam = gami*sqtrm
2802        gam = Sgam/costh
2803        dsi = (gami*Si*cosP**2/(sqtrm*costh)-gam/Si)
2804        dsa = (gami*Sa*sinP**2/(sqtrm*costh)-gam/Sa)
2805        gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
2806        gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
2807        sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2808        sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2809    else:           #ellipsoidal crystallites
2810        const = 1.8*wave/(np.pi*costh)
2811        Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
2812        H = np.array(refl[:3])
2813        lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
2814        Sgam = 1.8*wave/(np.pi*costh*lenR)
2815        for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
2816            gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
2817            sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2818    gamDict[phfx+'Size;mx'] = Sgam
2819    sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2
2820           
2821    #microstrain derivatives               
2822    if calcControls[phfx+'MustrainType'] == 'isotropic':
2823        Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5]/2.)/np.pi
2824        gamDict[phfx+'Mustrain;i'] =  0.018*tanth*parmDict[phfx+'Mustrain;mx']/np.pi
2825        sigDict[phfx+'Mustrain;i'] =  0.036*Mgam*tanth*(1.-parmDict[phfx+'Mustrain;mx'])**2/(np.pi*ateln2)       
2826    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2827        H = np.array(refl[:3])
2828        P = np.array(calcControls[phfx+'MustrainAxis'])
2829        cosP,sinP = G2lat.CosSinAngle(H,P,G)
2830        Si = parmDict[phfx+'Mustrain;i']
2831        Sa = parmDict[phfx+'Mustrain;a']
2832        gami = 0.018*Si*Sa*tanth/np.pi
2833        sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
2834        Mgam = gami/sqtrm
2835        dsi = -gami*Si*cosP**2/sqtrm**3
2836        dsa = -gami*Sa*sinP**2/sqtrm**3
2837        gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
2838        gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
2839        sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
2840        sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
2841    else:       #generalized - P.W. Stephens model
2842        pwrs = calcControls[phfx+'MuPwrs']
2843        const = 0.018*refl[4]**2*tanth
2844        sum = 0
2845        for i,pwr in enumerate(pwrs):
2846            term = refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
2847            sum += parmDict[phfx+'Mustrain:'+str(i)]*term
2848            gamDict[phfx+'Mustrain:'+str(i)] = const*term*parmDict[phfx+'Mustrain;mx']
2849            sigDict[phfx+'Mustrain:'+str(i)] = \
2850                2.*const*term*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
2851        Mgam = 0.018*refl[4]**2*tand(refl[5]/2.)*sum
2852        for i in range(len(pwrs)):
2853            sigDict[phfx+'Mustrain:'+str(i)] *= Mgam
2854    gamDict[phfx+'Mustrain;mx'] = Mgam
2855    sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
2856    return sigDict,gamDict
2857       
2858def GetReflPos(refl,wave,G,hfx,calcControls,parmDict):
2859    h,k,l = refl[:3]
2860    dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G)
2861    d = np.sqrt(dsq)
2862
2863    refl[4] = d
2864    pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
2865    const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
2866    if 'Bragg' in calcControls[hfx+'instType']:
2867        pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
2868            parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
2869    else:               #Debye-Scherrer - simple but maybe not right
2870        pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
2871    return pos
2872
2873def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict):
2874    dpr = 180./np.pi
2875    h,k,l = refl[:3]
2876    dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
2877    dst = np.sqrt(dstsq)
2878    pos = refl[5]-parmDict[hfx+'Zero']
2879    const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0)
2880    dpdw = const*dst
2881    dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])
2882    dpdA *= const*wave/(2.0*dst)
2883    dpdZ = 1.0
2884    const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
2885    if 'Bragg' in calcControls[hfx+'instType']:
2886        dpdSh = -4.*const*cosd(pos/2.0)
2887        dpdTr = -const*sind(pos)*100.0
2888        return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.
2889    else:               #Debye-Scherrer - simple but maybe not right
2890        dpdXd = -const*cosd(pos)
2891        dpdYd = -const*sind(pos)
2892        return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd
2893           
2894def GetHStrainShift(refl,SGData,phfx,parmDict):
2895    laue = SGData['SGLaue']
2896    uniq = SGData['SGUniq']
2897    h,k,l = refl[:3]
2898    if laue in ['m3','m3m']:
2899        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
2900            refl[4]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
2901    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2902        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
2903    elif laue in ['3R','3mR']:
2904        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
2905    elif laue in ['4/m','4/mmm']:
2906        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
2907    elif laue in ['mmm']:
2908        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
2909    elif laue in ['2/m']:
2910        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
2911        if uniq == 'a':
2912            Dij += parmDict[phfx+'D23']*k*l
2913        elif uniq == 'b':
2914            Dij += parmDict[phfx+'D13']*h*l
2915        elif uniq == 'c':
2916            Dij += parmDict[phfx+'D12']*h*k
2917    else:
2918        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
2919            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
2920    return -Dij*refl[4]**2*tand(refl[5]/2.0)
2921           
2922def GetHStrainShiftDerv(refl,SGData,phfx):
2923    laue = SGData['SGLaue']
2924    uniq = SGData['SGUniq']
2925    h,k,l = refl[:3]
2926    if laue in ['m3','m3m']:
2927        dDijDict = {phfx+'D11':h**2+k**2+l**2,
2928            phfx+'eA':refl[4]**2*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
2929    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2930        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
2931    elif laue in ['3R','3mR']:
2932        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
2933    elif laue in ['4/m','4/mmm']:
2934        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
2935    elif laue in ['mmm']:
2936        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
2937    elif laue in ['2/m']:
2938        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
2939        if uniq == 'a':
2940            dDijDict[phfx+'D23'] = k*l
2941        elif uniq == 'b':
2942            dDijDict[phfx+'D13'] = h*l
2943        elif uniq == 'c':
2944            dDijDict[phfx+'D12'] = h*k
2945            names.append()
2946    else:
2947        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
2948            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
2949    for item in dDijDict:
2950        dDijDict[item] *= -refl[4]**2*tand(refl[5]/2.0)
2951    return dDijDict
2952   
2953def GetFprime(controlDict,Histograms):
2954    FFtables = controlDict['FFtables']
2955    if not FFtables:
2956        return
2957    histoList = Histograms.keys()
2958    histoList.sort()
2959    for histogram in histoList:
2960        if histogram[:4] in ['PWDR','HKLF']:
2961            Histogram = Histograms[histogram]
2962            hId = Histogram['hId']
2963            hfx = ':%d:'%(hId)
2964            keV = controlDict[hfx+'keV']
2965            for El in FFtables:
2966                Orbs = G2el.GetXsectionCoeff(El.split('+')[0].split('-')[0])
2967                FP,FPP,Mu = G2el.FPcalc(Orbs, keV)
2968                FFtables[El][hfx+'FP'] = FP
2969                FFtables[El][hfx+'FPP'] = FPP               
2970           
2971def GetFobsSq(Histograms,Phases,parmDict,calcControls):
2972    histoList = Histograms.keys()
2973    histoList.sort()
2974    for histogram in histoList:
2975        if 'PWDR' in histogram[:4]:
2976            Histogram = Histograms[histogram]
2977            hId = Histogram['hId']
2978            hfx = ':%d:'%(hId)
2979            Limits = calcControls[hfx+'Limits']
2980            shl = max(parmDict[hfx+'SH/L'],0.0005)
2981            Ka2 = False
2982            kRatio = 0.0
2983            if hfx+'Lam1' in parmDict.keys():
2984                Ka2 = True
2985                lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2986                kRatio = parmDict[hfx+'I(L2)/I(L1)']
2987            x,y,w,yc,yb,yd = Histogram['Data']
2988            xB = np.searchsorted(x,Limits[0])
2989            xF = np.searchsorted(x,Limits[1])
2990            ymb = np.array(y-yb)
2991            ymb = np.where(ymb,ymb,1.0)
2992            ycmb = np.array(yc-yb)
2993            ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)         
2994            refLists = Histogram['Reflection Lists']
2995            for phase in refLists:
2996                Phase = Phases[phase]
2997                pId = Phase['pId']
2998                phfx = '%d:%d:'%(pId,hId)
2999                refList = refLists[phase]
3000                sumFo = 0.0
3001                sumdF = 0.0
3002                sumFosq = 0.0
3003                sumdFsq = 0.0
3004                for refl in refList:
3005                    if 'C' in calcControls[hfx+'histType']:
3006                        yp = np.zeros_like(yb)
3007                        Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
3008                        iBeg = max(xB,np.searchsorted(x,refl[5]-fmin))
3009                        iFin = max(xB,min(np.searchsorted(x,refl[5]+fmax),xF))
3010                        iFin2 = iFin
3011                        yp[iBeg:iFin] = refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
3012                        if Ka2:
3013                            pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
3014                            Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
3015                            iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
3016                            iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
3017                            if not iBeg2+iFin2:       #peak below low limit - skip peak
3018                                continue
3019                            elif not iBeg2-iFin2:     #peak above high limit - done
3020                                break
3021                            yp[iBeg2:iFin2] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])        #and here
3022                        refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[13]*(1.+kRatio)),0.0))
3023                    elif 'T' in calcControls[hfx+'histType']:
3024                        print 'TOF Undefined at present'
3025                        raise Exception    #no TOF yet
3026                    Fo = np.sqrt(np.abs(refl[8]))
3027                    Fc = np.sqrt(np.abs(refl[9]))
3028                    sumFo += Fo
3029                    sumFosq += refl[8]**2
3030                    sumdF += np.abs(Fo-Fc)
3031                    sumdFsq += (refl[8]-refl[9])**2
3032                Histogram[phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
3033                Histogram[phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
3034                Histogram[phfx+'Nref'] = len(refList)
3035               
3036def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
3037   
3038    def GetReflSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict):
3039        U = parmDict[hfx+'U']
3040        V = parmDict[hfx+'V']
3041        W = parmDict[hfx+'W']
3042        X = parmDict[hfx+'X']
3043        Y = parmDict[hfx+'Y']
3044        tanPos = tand(refl[5]/2.0)
3045        Ssig,Sgam = GetSampleSigGam(refl,wave,G,GB,phfx,calcControls,parmDict)
3046        sig = U*tanPos**2+V*tanPos+W+Ssig     #save peak sigma
3047        sig = max(0.001,sig)
3048        gam = X/cosd(refl[5]/2.0)+Y*tanPos+Sgam     #save peak gamma
3049        gam = max(0.001,gam)
3050        return sig,gam
3051               
3052    hId = Histogram['hId']
3053    hfx = ':%d:'%(hId)
3054    bakType = calcControls[hfx+'bakType']
3055    yb = G2pwd.getBackground(hfx,parmDict,bakType,x)
3056    yc = np.zeros_like(yb)
3057       
3058    if 'C' in calcControls[hfx+'histType']:   
3059        shl = max(parmDict[hfx+'SH/L'],0.002)
3060        Ka2 = False
3061        if hfx+'Lam1' in parmDict.keys():
3062            wave = parmDict[hfx+'Lam1']
3063            Ka2 = True
3064            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
3065            kRatio = parmDict[hfx+'I(L2)/I(L1)']
3066        else:
3067            wave = parmDict[hfx+'Lam']
3068    else:
3069        print 'TOF Undefined at present'
3070        raise ValueError
3071    for phase in Histogram['Reflection Lists']:
3072        refList = Histogram['Reflection Lists'][phase]
3073        Phase = Phases[phase]
3074        pId = Phase['pId']
3075        pfx = '%d::'%(pId)
3076        phfx = '%d:%d:'%(pId,hId)
3077        hfx = ':%d:'%(hId)
3078        SGData = Phase['General']['SGData']
3079        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
3080        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
3081        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
3082        Vst = np.sqrt(nl.det(G))    #V*
3083        if not Phase['General'].get('doPawley'):
3084            refList = StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict)
3085        for refl in refList:
3086            if 'C' in calcControls[hfx+'histType']:
3087                h,k,l = refl[:3]
3088                refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict)         #corrected reflection position
3089                Lorenz = 1./(2.*sind(refl[5]/2.)**2*cosd(refl[5]/2.))           #Lorentz correction
3090                refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict)               #apply hydrostatic strain shift
3091                refl[6:8] = GetReflSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict)    #peak sig & gam
3092                GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)    #puts corrections in refl[13]
3093                refl[13] *= Vst*Lorenz
3094                if Phase['General'].get('doPawley'):
3095                    try:
3096                        pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
3097                        refl[9] = parmDict[pInd]
3098                    except KeyError:
3099#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
3100                        continue
3101                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
3102                iBeg = np.searchsorted(x,refl[5]-fmin)
3103                iFin = np.searchsorted(x,refl[5]+fmax)
3104                if not iBeg+iFin:       #peak below low limit - skip peak
3105                    continue
3106                elif not iBeg-iFin:     #peak above high limit - done
3107                    break
3108                yc[iBeg:iFin] += refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
3109                if Ka2:
3110                    pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
3111                    Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
3112                    iBeg = np.searchsorted(x,pos2-fmin)
3113                    iFin = np.searchsorted(x,pos2+fmax)
3114                    if not iBeg+iFin:       #peak below low limit - skip peak
3115                        continue
3116                    elif not iBeg-iFin:     #peak above high limit - done
3117                        return yc,yb
3118                    yc[iBeg:iFin] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin])        #and here
3119            elif 'T' in calcControls[hfx+'histType']:
3120                print 'TOF Undefined at present'
3121                raise Exception    #no TOF yet
3122    return yc,yb
3123   
3124def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
3125   
3126    def cellVaryDerv(pfx,SGData,dpdA): 
3127        if SGData['SGLaue'] in ['-1',]:
3128            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
3129                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
3130        elif SGData['SGLaue'] in ['2/m',]:
3131            if SGData['SGUniq'] == 'a':
3132                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
3133            elif SGData['SGUniq'] == 'b':
3134                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
3135            else:
3136                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
3137        elif SGData['SGLaue'] in ['mmm',]:
3138            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
3139        elif SGData['SGLaue'] in ['4/m','4/mmm']:
3140            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
3141        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
3142            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
3143        elif SGData['SGLaue'] in ['3R', '3mR']:
3144            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
3145        elif SGData['SGLaue'] in ['m3m','m3']:
3146            return [[pfx+'A0',dpdA[0]]]
3147    # create a list of dependent variables and set up a dictionary to hold their derivatives
3148    dependentVars = G2mv.GetDependentVars()
3149    depDerivDict = {}
3150    for j in dependentVars:
3151        depDerivDict[j] = np.zeros(shape=(len(x)))
3152    #print 'dependent vars',dependentVars
3153    lenX = len(x)               
3154    hId = Histogram['hId']
3155    hfx = ':%d:'%(hId)
3156    bakType = calcControls[hfx+'bakType']
3157    dMdv = np.zeros(shape=(len(varylist),len(x)))
3158    dMdb,dMddb,dMdpk = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x)
3159    if hfx+'Back:0' in varylist: # for now assume that Back:x vars to not appear in constraints
3160        bBpos =varylist.index(hfx+'Back:0')
3161        dMdv[bBpos:bBpos+len(dMdb)] = dMdb
3162    names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
3163    for name in varylist:
3164        if 'Debye' in name:
3165            id = int(name.split(':')[-1])
3166            parm = name[:int(name.rindex(':'))]
3167            ip = names.index(parm)
3168            dMdv[varylist.index(name)] = dMddb[3*id+ip]
3169    names = [hfx+'BkPkpos',hfx+'BkPkint',hfx+'BkPksig',hfx+'BkPkgam']
3170    for name in varylist:
3171        if 'BkPk' in name:
3172            id = int(name.split(':')[-1])
3173            parm = name[:int(name.rindex(':'))]
3174            ip = names.index(parm)
3175            dMdv[varylist.index(name)] = dMdpk[4*id+ip]
3176    cw = np.diff(x)
3177    cw = np.append(cw,cw[-1])
3178    if 'C' in calcControls[hfx+'histType']:   
3179        shl = max(parmDict[hfx+'SH/L'],0.002)
3180        Ka2 = False
3181        if hfx+'Lam1' in parmDict.keys():
3182            wave = parmDict[hfx+'Lam1']
3183            Ka2 = True
3184            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
3185            kRatio = parmDict[hfx+'I(L2)/I(L1)']
3186        else:
3187            wave = parmDict[hfx+'Lam']
3188    else:
3189        print 'TOF Undefined at present'
3190        raise ValueError
3191    for phase in Histogram['Reflection Lists']:
3192        refList = Histogram['Reflection Lists'][phase]
3193        Phase = Phases[phase]
3194        SGData = Phase['General']['SGData']
3195        pId = Phase['pId']
3196        pfx = '%d::'%(pId)
3197        phfx = '%d:%d:'%(pId,hId)
3198        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
3199        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
3200        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
3201        if not Phase['General'].get('doPawley'):
3202            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)
3203        for iref,refl in enumerate(refList):
3204            if 'C' in calcControls[hfx+'histType']:        #CW powder
3205                h,k,l = refl[:3]
3206                dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb = GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
3207                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
3208                iBeg = np.searchsorted(x,refl[5]-fmin)
3209                iFin = np.searchsorted(x,refl[5]+fmax)
3210                if not iBeg+iFin:       #peak below low limit - skip peak
3211                    continue
3212                elif not iBeg-iFin:     #peak above high limit - done
3213                    break
3214                pos = refl[5]
3215                tanth = tand(pos/2.0)
3216                costh = cosd(pos/2.0)
3217                lenBF = iFin-iBeg
3218                dMdpk = np.zeros(shape=(6,lenBF))
3219                dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])
3220                for i in range(5):
3221                    dMdpk[i] += 100.*cw[iBeg:iFin]*refl[13]*refl[9]*dMdipk[i]
3222                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
3223                if Ka2:
3224                    pos2 = refl[5]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
3225                    iBeg2 = np.searchsorted(x,pos2-fmin)
3226                    iFin2 = np.searchsorted(x,pos2+fmax)
3227                    if iBeg2-iFin2:
3228                        lenBF2 = iFin2-iBeg2
3229                        dMdpk2 = np.zeros(shape=(6,lenBF2))
3230                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])
3231                        for i in range(5):
3232                            dMdpk2[i] = 100.*cw[iBeg2:iFin2]*refl[13]*refl[9]*kRatio*dMdipk2[i]
3233                        dMdpk2[5] = 100.*cw[iBeg2:iFin2]*refl[13]*dMdipk2[0]
3234                        dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
3235                if Phase['General'].get('doPawley'):
3236                    dMdpw = np.zeros(len(x))
3237                    try:
3238                        pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
3239                        idx = varylist.index(pIdx)
3240                        dMdpw[iBeg:iFin] = dervDict['int']/refl[9]
3241                        if Ka2:
3242                            dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9]
3243                        dMdv[idx] = dMdpw
3244                    except: # ValueError:
3245                        pass
3246                dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
3247                names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
3248                    hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
3249                    hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
3250                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
3251                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
3252                    hfx+'DisplaceY':[dpdY,'pos'],hfx+'Absorption':[dFdAb,'int'],}
3253                for name in names:
3254                    item = names[name]
3255                    if name in varylist:
3256                        dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
3257                        if Ka2:
3258                            dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
3259                    elif name in dependentVars:
3260                        if Ka2:
3261                            depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
3262                        depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
3263                for iPO in dIdPO:
3264                    if iPO in varylist:
3265                        dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
3266                        if Ka2:
3267                            dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
3268                    elif iPO in dependentVars:
3269                        depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
3270                        if Ka2:
3271                            depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
3272                for i,name in enumerate(['omega','chi','phi']):
3273                    aname = pfx+'SH '+name
3274                    if aname in varylist:
3275                        dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
3276                        if Ka2:
3277                            dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
3278                    elif aname in dependentVars:
3279                        depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
3280                        if Ka2:
3281                            depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
3282                for iSH in dFdODF:
3283                    if iSH in varylist:
3284                        dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
3285                        if Ka2:
3286                            dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
3287                    elif iSH in dependentVars:
3288                        depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
3289                        if Ka2:
3290                            depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
3291                cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
3292                for name,dpdA in cellDervNames:
3293                    if name in varylist:
3294                        dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
3295                        if Ka2:
3296                            dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
3297                    elif name in dependentVars:
3298                        depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
3299                        if Ka2:
3300                            depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
3301                dDijDict = GetHStrainShiftDerv(refl,SGData,phfx)
3302                for name in dDijDict:
3303                    if name in varylist:
3304                        dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
3305                        if Ka2:
3306                            dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
3307                    elif name in dependentVars:
3308                        depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
3309                        if Ka2:
3310                            depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
3311                sigDict,gamDict = GetSampleSigGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict)
3312                for name in gamDict:
3313                    if name in varylist:
3314                        dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
3315                        if Ka2:
3316                            dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
3317                    elif name in dependentVars:
3318                        depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
3319                        if Ka2:
3320                            depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
3321                for name in sigDict:
3322                    if name in varylist:
3323                        dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
3324                        if Ka2:
3325                            dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
3326                    elif name in dependentVars:
3327                        depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
3328                        if Ka2:
3329                            depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
3330                for name in ['BabA','BabU']:
3331                    if phfx+name in varylist:
3332                        dMdv[varylist.index(phfx+name)][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']*cw[iBeg:iFin]
3333                        if Ka2:
3334                            dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']*cw[iBeg2:iFin2]
3335                    elif phfx+name in dependentVars:                   
3336                        depDerivDict[phfx+name][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']*cw[iBeg:iFin]
3337                        if Ka2:
3338                            depDerivDict[phfx+name][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']*cw[iBeg2:iFin2]                 
3339            elif 'T' in calcControls[hfx+'histType']:
3340                print 'TOF Undefined at present'
3341                raise Exception    #no TOF yet
3342            #do atom derivatives -  for F,X & U so far             
3343            corr = dervDict['int']/refl[9]
3344            if Ka2:
3345                corr2 = dervDict2['int']/refl[9]
3346            for name in varylist+dependentVars:
3347                try:
3348                    aname = name.split(pfx)[1][:2]
3349                    if aname not in ['Af','dA','AU']: continue # skip anything not an atom param
3350                except IndexError:
3351                    continue
3352                if name in varylist:
3353                    dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
3354                    if Ka2:
3355                        dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
3356                elif name in dependentVars:
3357                    depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
3358                    if Ka2:
3359                        depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
3360    # now process derivatives in constraints
3361    G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
3362    return dMdv
3363
3364def dervRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
3365    parmdict.update(zip(varylist,values))
3366    G2mv.Dict2Map(parmdict,varylist)
3367    Histograms,Phases,restraintDict = HistoPhases
3368    nvar = len(varylist)
3369    dMdv = np.empty(0)
3370    histoList = Histograms.keys()
3371    histoList.sort()
3372    for histogram in histoList:
3373        if 'PWDR' in histogram[:4]:
3374            Histogram = Histograms[histogram]
3375            hId = Histogram['hId']
3376            hfx = ':%d:'%(hId)
3377            wtFactor = calcControls[hfx+'wtFactor']
3378            Limits = calcControls[hfx+'Limits']
3379            x,y,w,yc,yb,yd = Histogram['Data']
3380            W = wtFactor*w
3381            xB = np.searchsorted(x,Limits[0])
3382            xF = np.searchsorted(x,Limits[1])
3383            dMdvh = np.sqrt(W[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF],
3384                varylist,Histogram,Phases,calcControls,pawleyLookup)
3385        elif 'HKLF' in histogram[:4]:
3386            Histogram = Histograms[histogram]
3387            nobs = Histogram['Nobs']
3388            phase = Histogram['Reflection Lists']
3389            Phase = Phases[phase]
3390            hId = Histogram['hId']
3391            hfx = ':%d:'%(hId)
3392            wtFactor = calcControls[hfx+'wtFactor']
3393            pfx = '%d::'%(Phase['pId'])
3394            phfx = '%d:%d:'%(Phase['pId'],hId)
3395            SGData = Phase['General']['SGData']
3396            A = [parmdict[pfx+'A%d'%(i)] for i in range(6)]
3397            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
3398            refList = Histogram['Data']
3399            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmdict)
3400            dMdvh = np.zeros((len(varylist),len(refList)))
3401            for iref,ref in enumerate(refList):
3402                if ref[6] > 0:
3403                    dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmdict,varylist) #puts correction in refl[13]
3404                    if calcControls['F**2']:
3405                        if ref[5]/ref[6] >= calcControls['minF/sig']:
3406                            w = wtFactor/ref[6]
3407                            for j,var in enumerate(varylist):
3408                                if var in dFdvDict:
3409                                    dMdvh[j][iref] = w*dFdvDict[var][iref]*dervCor
3410                            if phfx+'Scale' in varylist:
3411                                dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*dervCor
3412                    else:
3413                        Fo = np.sqrt(ref[5])
3414                        Fc = np.sqrt(ref[7])
3415                        sig = ref[6]/(2.0*Fo)
3416                        if Fo/sig >= calcControls['minF/sig']:
3417                            w = wtFactor/sig
3418                            for j,var in enumerate(varylist):
3419                                if var in dFdvDict:
3420                                    dMdvh[j][iref] = w*dFdvDict[var][iref]*np.sqrt(dervCor)
3421                            if phfx+'Scale' in varylist:
3422                                dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*np.sqrt(dervCor)                           
3423                    for item in ['Ep','Es','Eg']:
3424                        if phfx+item in varylist:
3425                            dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]
3426        else:
3427            continue        #skip non-histogram entries
3428        if len(dMdv):
3429            dMdv = np.concatenate((dMdv.T,dMdvh.T)).T
3430        else:
3431            dMdv = dMdvh
3432           
3433    pNames,pVals,pWt = penaltyFxn(HistoPhases,parmdict,varylist)
3434    if np.any(pVals):
3435        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmdict,varylist)
3436        dMdv = np.concatenate((dMdv.T,dpdv.T)).T
3437       
3438    return dMdv
3439
3440def HessRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
3441    parmdict.update(zip(varylist,values))
3442    G2mv.Dict2Map(parmdict,varylist)
3443    Histograms,Phases,restraintDict = HistoPhases
3444    nvar = len(varylist)
3445    Hess = np.empty(0)
3446    histoList = Histograms.keys()
3447    histoList.sort()
3448    for histogram in histoList:
3449        if 'PWDR' in histogram[:4]:
3450            Histogram = Histograms[histogram]
3451            hId = Histogram['hId']
3452            hfx = ':%d:'%(hId)
3453            wtFactor = calcControls[hfx+'wtFactor']
3454            Limits = calcControls[hfx+'Limits']
3455            x,y,w,yc,yb,yd = Histogram['Data']
3456            W = wtFactor*w
3457            dy = y-yc
3458            xB = np.searchsorted(x,Limits[0])
3459            xF = np.searchsorted(x,Limits[1])
3460            dMdvh = getPowderProfileDerv(parmdict,x[xB:xF],
3461                varylist,Histogram,Phases,calcControls,pawleyLookup)
3462            Wt = np.sqrt(W[xB:xF])[np.newaxis,:]
3463            Dy = dy[xB:xF][np.newaxis,:]
3464            dMdvh *= Wt
3465            if dlg:
3466                dlg.Update(Histogram['wR'],newmsg='Hessian for histogram %d\nAll data Rw=%8.3f%s'%(hId,Histogram['wR'],'%'))[0]
3467            if len(Hess):
3468                Hess += np.inner(dMdvh,dMdvh)
3469                dMdvh *= Wt*Dy
3470                Vec += np.sum(dMdvh,axis=1)
3471            else:
3472                Hess = np.inner(dMdvh,dMdvh)
3473                dMdvh *= Wt*Dy
3474                Vec = np.sum(dMdvh,axis=1)
3475        elif 'HKLF' in histogram[:4]:
3476            Histogram = Histograms[histogram]
3477            nobs = Histogram['Nobs']
3478            phase = Histogram['Reflection Lists']
3479            Phase = Phases[phase]
3480            hId = Histogram['hId']
3481            hfx = ':%d:'%(hId)
3482            wtFactor = calcControls[hfx+'wtFactor']
3483            pfx = '%d::'%(Phase['pId'])
3484            phfx = '%d:%d:'%(Phase['pId'],hId)
3485            SGData = Phase['General']['SGData']
3486            A = [parmdict[pfx+'A%d'%(i)] for i in range(6)]
3487            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
3488            refList = Histogram['Data']
3489            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmdict)
3490            dMdvh = np.zeros((len(varylist),len(refList)))
3491            wdf = np.zeros(len(refList))
3492            for iref,ref in enumerate(refList):
3493                if ref[6] > 0:
3494                    dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmdict,varylist) #puts correction in refl[13]
3495                    if calcControls['F**2']:
3496                        if ref[5]/ref[6] >= calcControls['minF/sig']:
3497                            w =  wtFactor/ref[6]
3498                            wdf[iref] = w*(ref[5]-ref[7])
3499                            for j,var in enumerate(varylist):
3500                                if var in dFdvDict:
3501                                    dMdvh[j][iref] = w*dFdvDict[var][iref]*dervCor
3502                            if phfx+'Scale' in varylist:
3503                                dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*dervCor
3504                    else:
3505                        if ref[5] > 0.:
3506                            Fo = np.sqrt(ref[5])
3507                            Fc = np.sqrt(ref[7])
3508                            sig = ref[6]/(2.0*Fo)
3509                            w = wtFactor/sig
3510                            wdf[iref] = w*(Fo-Fc)
3511                            if Fo/sig >= calcControls['minF/sig']:
3512                                for j,var in enumerate(varylist):
3513                                    if var in dFdvDict:
3514                                        dMdvh[j][iref] = w*dFdvDict[var][iref]*np.sqrt(dervCor)
3515                                if phfx+'Scale' in varylist:
3516                                    dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*np.sqrt(dervCor)                           
3517                    for item in ['Ep','Es','Eg']:
3518                        if phfx+item in varylist:
3519                            dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]
3520            if dlg:
3521                dlg.Update(Histogram['wR'],newmsg='Hessian for histogram %d Rw=%8.3f%s'%(hId,Histogram['wR'],'%'))[0]
3522            if len(Hess):
3523                Vec += np.sum(dMdvh*wdf,axis=1)
3524                Hess += np.inner(dMdvh,dMdvh)
3525            else:
3526                Vec = np.sum(dMdvh*wdf,axis=1)
3527                Hess = np.inner(dMdvh,dMdvh)
3528        else:
3529            continue        #skip non-histogram entries
3530    pNames,pVals,pWt = penaltyFxn(HistoPhases,parmdict,varylist)
3531    if np.any(pVals):
3532        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmdict,varylist)
3533        Vec += np.sum(dpdv*pWt*pVals,axis=1)
3534        Hess += np.inner(dpdv*pWt,dpdv)
3535    return Vec,Hess
3536
3537def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):       
3538    parmdict.update(zip(varylist,values))
3539    Values2Dict(parmdict, varylist, values)
3540    G2mv.Dict2Map(parmdict,varylist)
3541    Histograms,Phases,restraintDict = HistoPhases
3542    M = np.empty(0)
3543    SumwYo = 0
3544    Nobs = 0
3545    histoList = Histograms.keys()
3546    histoList.sort()
3547    for histogram in histoList:
3548        if 'PWDR' in histogram[:4]:
3549            Histogram = Histograms[histogram]
3550            hId = Histogram['hId']
3551            hfx = ':%d:'%(hId)
3552            wtFactor = calcControls[hfx+'wtFactor']
3553            Limits = calcControls[hfx+'Limits']
3554            x,y,w,yc,yb,yd = Histogram['Data']
3555            W = wtFactor*w
3556            yc *= 0.0                           #zero full calcd profiles
3557            yb *= 0.0
3558            yd *= 0.0
3559            xB = np.searchsorted(x,Limits[0])
3560            xF = np.searchsorted(x,Limits[1])
3561            Histogram['Nobs'] = xF-xB
3562            Nobs += Histogram['Nobs']
3563            Histogram['sumwYo'] = np.sum(W[xB:xF]*y[xB:xF]**2)
3564            SumwYo += Histogram['sumwYo']
3565            yc[xB:xF],yb[xB:xF] = getPowderProfile(parmdict,x[xB:xF],
3566                varylist,Histogram,Phases,calcControls,pawleyLookup)
3567            yc[xB:xF] += yb[xB:xF]
3568            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
3569            Histogram['sumwYd'] = np.sum(np.sqrt(W[xB:xF])*(yd[xB:xF]))
3570            wdy = -np.sqrt(W[xB:xF])*(yd[xB:xF])
3571            Histogram['wR'] = min(100.,np.sqrt(np.sum(wdy**2)/Histogram['sumwYo'])*100.)
3572            if dlg:
3573                dlg.Update(Histogram['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['wR'],'%'))[0]
3574            M = np.concatenate((M,wdy))
3575#end of PWDR processing
3576        elif 'HKLF' in histogram[:4]:
3577            Histogram = Histograms[histogram]
3578            phase = Histogram['Reflection Lists']
3579            Phase = Phases[phase]
3580            hId = Histogram['hId']
3581            hfx = ':%d:'%(hId)
3582            wtFactor = calcControls[hfx+'wtFactor']
3583            pfx = '%d::'%(Phase['pId'])
3584            phfx = '%d:%d:'%(Phase['pId'],hId)
3585            SGData = Phase['General']['SGData']
3586            A = [parmdict[pfx+'A%d'%(i)] for i in range(6)]
3587            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
3588            refList = Histogram['Data']
3589            refList = StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmdict)
3590            df = np.zeros(len(refList))
3591            sumwYo = 0
3592            sumFo = 0
3593            sumFo2 = 0
3594            sumdF = 0
3595            sumdF2 = 0
3596            nobs = 0
3597            for i,ref in enumerate(refList):
3598                if ref[6] > 0:
3599                    SCExtinction(ref,phfx,hfx,pfx,calcControls,parmdict,varylist) #puts correction in refl[13]
3600                    ref[7] = parmdict[phfx+'Scale']*ref[9]
3601                    ref[7] *= ref[13]
3602                    ref[8] = ref[5]/parmdict[phfx+'Scale']
3603                    if calcControls['F**2']:
3604                        if ref[5]/ref[6] >= calcControls['minF/sig']:
3605                            sumFo2 += ref[5]
3606                            Fo = np.sqrt(ref[5])
3607                            sumFo += Fo
3608                            sumFo2 += ref[5]
3609                            sumdF += abs(Fo-np.sqrt(ref[7]))
3610                            sumdF2 += abs(ref[5]-ref[7])
3611                            nobs += 1
3612                            df[i] = -np.sqrt(wtFactor)*(ref[5]-ref[7])/ref[6]
3613                            sumwYo += wtFactor*(ref[5]/ref[6])**2
3614                    else:
3615                        Fo = np.sqrt(ref[5])
3616                        Fc = np.sqrt(ref[7])
3617                        sig = ref[6]/(2.0*Fo)
3618                        if Fo/sig >= calcControls['minF/sig']:
3619                            sumFo += Fo
3620                            sumFo2 += ref[5]
3621                            sumdF += abs(Fo-Fc)
3622                            sumdF2 += abs(ref[5]-ref[7])
3623                            nobs += 1
3624                            df[i] = -np.sqrt(wtFactor)*(Fo-Fc)/sig
3625                            sumwYo += wtFactor*(Fo/sig)**2
3626            Histogram['Nobs'] = nobs
3627            Histogram['sumwYo'] = sumwYo
3628            SumwYo += sumwYo
3629            Histogram['wR'] = min(100.,np.sqrt(np.sum(df**2)/Histogram['sumwYo'])*100.)
3630            Histogram[phfx+'Rf'] = 100.*sumdF/sumFo
3631            Histogram[phfx+'Rf^2'] = 100.*sumdF2/sumFo2
3632            Histogram[phfx+'Nref'] = nobs
3633            Nobs += nobs
3634            if dlg:
3635                dlg.Update(Histogram['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['wR'],'%'))[0]
3636            M = np.concatenate((M,df))
3637# end of HKLF processing
3638    Histograms['sumwYo'] = SumwYo
3639    Histograms['Nobs'] = Nobs
3640    Rw = min(100.,np.sqrt(np.sum(M**2)/SumwYo)*100.)
3641    if dlg:
3642        GoOn = dlg.Update(Rw,newmsg='%s%8.3f%s'%('All data Rw =',Rw,'%'))[0]
3643        if not GoOn:
3644            parmdict['saved values'] = values
3645            dlg.Destroy()
3646            raise Exception         #Abort!!
3647    pDict,pVals,pWt = penaltyFxn(HistoPhases,parmdict,varylist)
3648    if np.any(pVals):
3649        pSum = np.sum(pWt*pVals**2)
3650        print 'Penalty function: %.3f on %d terms'%(pSum,len(pVals))
3651        Nobs += len(pVals)
3652        M = np.concatenate((M,np.sqrt(pWt)*pVals))
3653    return M
3654                       
3655def Refine(GPXfile,dlg):
3656    import pytexture as ptx
3657    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
3658   
3659    printFile = open(ospath.splitext(GPXfile)[0]+'.lst','w')
3660    ShowBanner(printFile)
3661    varyList = []
3662    parmDict = {}
3663    G2mv.InitVars()   
3664    Controls = GetControls(GPXfile)
3665    ShowControls(Controls,printFile)
3666    calcControls = {}
3667    calcControls.update(Controls)           
3668    constrDict,fixedList = GetConstraints(GPXfile)
3669    restraintDict = GetRestraints(GPXfile)
3670    rigidbodyDict = GetRigidBodies(GPXfile)
3671    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
3672    if not Phases:
3673        print ' *** ERROR - you have no phases! ***'
3674        print ' *** Refine aborted ***'
3675        raise Exception
3676    if not Histograms:
3677        print ' *** ERROR - you have no data to refine with! ***'
3678        print ' *** Refine aborted ***'
3679        raise Exception       
3680    Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases,restraintDict,pFile=printFile)
3681    calcControls['atomIndx'] = atomIndx
3682    calcControls['Natoms'] = Natoms
3683    calcControls['FFtables'] = FFtables
3684    calcControls['BLtables'] = BLtables
3685    hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms,pFile=printFile)
3686    calcControls.update(controlDict)
3687    histVary,histDict,controlDict = GetHistogramData(Histograms,pFile=printFile)
3688    calcControls.update(controlDict)
3689    varyList = phaseVary+hapVary+histVary
3690    parmDict.update(phaseDict)
3691    parmDict.update(hapDict)
3692    parmDict.update(histDict)
3693    GetFprime(calcControls,Histograms)
3694    # do constraint processing
3695    try:
3696        groups,parmlist = G2mv.GroupConstraints(constrDict)
3697        G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList)
3698    except:
3699        print ' *** ERROR - your constraints are internally inconsistent ***'
3700        # traceback for debug
3701        #print 'varyList',varyList
3702        #print 'constrDict',constrDict
3703        #print 'fixedList',fixedList
3704        #import traceback
3705        #print traceback.format_exc()
3706        raise Exception(' *** Refine aborted ***')
3707    # # check to see which generated parameters are fully varied
3708    # msg = G2mv.SetVaryFlags(varyList)
3709    # if msg:
3710    #     print ' *** ERROR - you have not set the refine flags for constraints consistently! ***'
3711    #     print msg
3712    #     raise Exception(' *** Refine aborted ***')
3713    #print G2mv.VarRemapShow(varyList)
3714    G2mv.Map2Dict(parmDict,varyList)
3715    Rvals = {}
3716    while True:
3717        begin = time.time()
3718        values =  np.array(Dict2Values(parmDict, varyList))
3719        Ftol = Controls['min dM/M']
3720        Factor = Controls['shift factor']
3721        maxCyc = Controls['max cyc']
3722        if 'Jacobian' in Controls['deriv type']:           
3723            result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
3724                ftol=Ftol,col_deriv=True,factor=Factor,
3725                args=([Histograms,Phases,restraintDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
3726            ncyc = int(result[2]['nfev']/2)
3727        elif 'Hessian' in Controls['deriv type']:
3728            result = G2mth.HessianLSQ(errRefine,values,Hess=HessRefine,ftol=Ftol,maxcyc=maxCyc,
3729                args=([Histograms,Phases,restraintDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
3730            ncyc = result[2]['num cyc']+1
3731            Rvals['lamMax'] = result[2]['lamMax']
3732        else:           #'numeric'
3733            result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
3734                args=([Histograms,Phases,restraintDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
3735            ncyc = int(result[2]['nfev']/len(varyList))
3736#        table = dict(zip(varyList,zip(values,result[0],(result[0]-values))))
3737#        for item in table: print item,table[item]               #useful debug - are things shifting?
3738        runtime = time.time()-begin
3739        Rvals['chisq'] = np.sum(result[2]['fvec']**2)
3740        Values2Dict(parmDict, varyList, result[0])
3741        G2mv.Dict2Map(parmDict,varyList)
3742       
3743        Rvals['Nobs'] = Histograms['Nobs']
3744        Rvals['Rwp'] = np.sqrt(Rvals['chisq']/Histograms['sumwYo'])*100.      #to %
3745        Rvals['GOF'] = Rvals['chisq']/(Histograms['Nobs']-len(varyList))
3746        print >>printFile,'\n Refinement results:'
3747        print >>printFile,135*'-'
3748        print >>printFile,' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histograms['Nobs'],' Number of parameters: ',len(varyList)
3749        print >>printFile,' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles'%(runtime,runtime/ncyc,ncyc)
3750        print >>printFile,' wR = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],Rvals['chisq'],Rvals['GOF'])
3751        try:
3752            covMatrix = result[1]*Rvals['GOF']
3753            sig = np.sqrt(np.diag(covMatrix))
3754            if np.any(np.isnan(sig)):
3755                print '*** Least squares aborted - some invalid esds possible ***'
3756#            table = dict(zip(varyList,zip(values,result[0],(result[0]-values)/sig)))
3757#            for item in table: print item,table[item]               #useful debug - are things shifting?
3758            break                   #refinement succeeded - finish up!
3759        except TypeError:          #result[1] is None on singular matrix
3760            print '**** Refinement failed - singular matrix ****'
3761            if 'Hessian' in Controls['deriv type']:
3762                num = len(varyList)-1
3763                for i,val in enumerate(np.flipud(result[2]['psing'])):
3764                    if val:
3765                        print 'Removing parameter: ',varyList[num-i]
3766                        del(varyList[num-i])                   
3767            else:
3768                Ipvt = result[2]['ipvt']
3769                for i,ipvt in enumerate(Ipvt):
3770                    if not np.sum(result[2]['fjac'],axis=1)[i]:
3771                        print 'Removing parameter: ',varyList[ipvt-1]
3772                        del(varyList[ipvt-1])
3773                        break
3774
3775#    print 'dependentParmList: ',G2mv.dependentParmList
3776#    print 'arrayList: ',G2mv.arrayList
3777#    print 'invarrayList: ',G2mv.invarrayList
3778#    print 'indParmList: ',G2mv.indParmList
3779#    print 'fixedDict: ',G2mv.fixedDict
3780#    print 'test1'
3781    GetFobsSq(Histograms,Phases,parmDict,calcControls)
3782#    print 'test2'
3783    sigDict = dict(zip(varyList,sig))
3784    newCellDict = GetNewCellParms(parmDict,varyList)
3785    newAtomDict = ApplyXYZshifts(parmDict,varyList)
3786    covData = {'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
3787        'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict,'newCellDict':newCellDict}
3788    # add the uncertainties into the esd dictionary (sigDict)
3789    sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
3790    G2mv.PrintIndependentVars(parmDict,varyList,sigDict,pFile=printFile)
3791    SetPhaseData(parmDict,sigDict,Phases,covData,restraintDict,printFile)
3792    SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,pFile=printFile)
3793    SetHistogramData(parmDict,sigDict,Histograms,pFile=printFile)
3794    SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,covData)
3795    printFile.close()
3796    print ' Refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst'
3797    print ' ***** Refinement successful *****'
3798   
3799#for testing purposes!!!
3800    if DEBUG:
3801        import cPickle
3802        fl = open('structTestdata.dat','wb')
3803        cPickle.dump(parmDict,fl,1)
3804        cPickle.dump(varyList,fl,1)
3805        for histogram in Histograms:
3806            if 'PWDR' in histogram[:4]:
3807                Histogram = Histograms[histogram]
3808        cPickle.dump(Histogram,fl,1)
3809        cPickle.dump(Phases,fl,1)
3810        cPickle.dump(calcControls,fl,1)
3811        cPickle.dump(pawleyLookup,fl,1)
3812        fl.close()
3813
3814    if dlg:
3815        return Rvals['Rwp']
3816
3817def SeqRefine(GPXfile,dlg):
3818    import pytexture as ptx
3819    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
3820   
3821    printFile = open(ospath.splitext(GPXfile)[0]+'.lst','w')
3822    print ' Sequential Refinement'
3823    ShowBanner(printFile)
3824    G2mv.InitVars()   
3825    Controls = GetControls(GPXfile)
3826    ShowControls(Controls,printFile)           
3827    constrDict,fixedList = GetConstraints(GPXfile)
3828    restraintDict = GetRestraints(GPXfile)
3829    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
3830    if not Phases:
3831        print ' *** ERROR - you have no histograms to refine! ***'
3832        print ' *** Refine aborted ***'
3833        raise Exception
3834    if not Histograms:
3835        print ' *** ERROR - you have no data to refine with! ***'
3836        print ' *** Refine aborted ***'
3837        raise Exception
3838    Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases,restraintDict,False,printFile)
3839    for item in phaseVary:
3840        if '::A0' in item:
3841            print '**** WARNING - lattice parameters should not be refined in a sequential refinement ****'
3842            print '****           instead use the Dij parameters for each powder histogram            ****'
3843    if 'Seq Data' in Controls:
3844        histNames = Controls['Seq Data']
3845    else:
3846        histNames = GetHistogramNames(GPXfile,['PWDR',])
3847    if 'Reverse Seq' in Controls:
3848        if Controls['Reverse Seq']:
3849            histNames.reverse()
3850    SeqResult = {'histNames':histNames}
3851    makeBack = True
3852    for ihst,histogram in enumerate(histNames):
3853        ifPrint = False
3854        if dlg:
3855            dlg.SetTitle('Residual for histogram '+str(ihst))
3856        calcControls = {}
3857        calcControls['atomIndx'] = atomIndx
3858        calcControls['Natoms'] = Natoms
3859        calcControls['FFtables'] = FFtables
3860        calcControls['BLtables'] = BLtables
3861        varyList = []
3862        parmDict = {}
3863        Histo = {histogram:Histograms[histogram],}
3864        hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histo,False)
3865        calcControls.update(controlDict)
3866        histVary,histDict,controlDict = GetHistogramData(Histo,False)
3867        calcControls.update(controlDict)
3868        varyList = phaseVary+hapVary+histVary
3869        if not ihst:
3870            saveVaryList = varyList[:]
3871            for i,item in enumerate(saveVaryList):
3872                items = item.split(':')
3873                if items[1]:
3874                    items[1] = ''
3875                item = ':'.join(items)
3876                saveVaryList[i] = item
3877            SeqResult['varyList'] = saveVaryList
3878        else:
3879            newVaryList = varyList[:]
3880            for i,item in enumerate(newVaryList):
3881                items = item.split(':')
3882                if items[1]:
3883                    items[1] = ''
3884                item = ':'.join(items)
3885                newVaryList[i] = item
3886            if newVaryList != SeqResult['varyList']:
3887                print newVaryList
3888                print SeqResult['varyList']
3889                print '**** ERROR - variable list for this histogram does not match previous'
3890                raise Exception
3891        parmDict.update(phaseDict)
3892        parmDict.update(hapDict)
3893        parmDict.update(histDict)
3894        GetFprime(calcControls,Histo)
3895        # do constraint processing
3896        try:
3897            groups,parmlist = G2mv.GroupConstraints(constrDict)
3898            G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList)
3899        except:
3900            print ' *** ERROR - your constraints are internally inconsistent ***'
3901            raise Exception(' *** Refine aborted ***')
3902        # check to see which generated parameters are fully varied
3903        # msg = G2mv.SetVaryFlags(varyList)
3904        # if msg:
3905        #     print ' *** ERROR - you have not set the refine flags for constraints consistently! ***'
3906        #     print msg
3907        #     raise Exception(' *** Refine aborted ***')
3908        #print G2mv.VarRemapShow(varyList)
3909        G2mv.Map2Dict(parmDict,varyList)
3910        Rvals = {}
3911        while True:
3912            begin = time.time()
3913            values =  np.array(Dict2Values(parmDict,varyList))
3914            Ftol = Controls['min dM/M']
3915            Factor = Controls['shift factor']
3916            maxCyc = Controls['max cyc']
3917
3918            if 'Jacobian' in Controls['deriv type']:           
3919                result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
3920                    ftol=Ftol,col_deriv=True,factor=Factor,
3921                    args=([Histo,Phases,restraintDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
3922                ncyc = int(result[2]['nfev']/2)
3923            elif 'Hessian' in Controls['deriv type']:
3924                result = G2mth.HessianLSQ(errRefine,values,Hess=HessRefine,ftol=Ftol,maxcyc=maxCyc,
3925                    args=([Histo,Phases,restraintDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
3926                ncyc = result[2]['num cyc']+1                           
3927            else:           #'numeric'
3928                result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
3929                    args=([Histo,Phases,restraintDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
3930                ncyc = int(result[2]['nfev']/len(varyList))
3931
3932            runtime = time.time()-begin
3933            Rvals['chisq'] = np.sum(result[2]['fvec']**2)
3934            Values2Dict(parmDict, varyList, result[0])
3935            G2mv.Dict2Map(parmDict,varyList)
3936           
3937            Rvals['Rwp'] = np.sqrt(Rvals['chisq']/Histo['sumwYo'])*100.      #to %
3938            Rvals['GOF'] = Rvals['Rwp']/(Histo['Nobs']-len(varyList))
3939            Rvals['Nobs'] = Histo['Nobs']
3940            print >>printFile,'\n Refinement results for histogram: v'+histogram
3941            print >>printFile,135*'-'
3942            print >>printFile,' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histo['Nobs'],' Number of parameters: ',len(varyList)
3943            print >>printFile,' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles'%(runtime,runtime/ncyc,ncyc)
3944            print >>printFile,' wRp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],Rvals['chisq'],Rvals['GOF'])
3945            try:
3946                covMatrix = result[1]*Rvals['GOF']
3947                sig = np.sqrt(np.diag(covMatrix))
3948                if np.any(np.isnan(sig)):
3949                    print '*** Least squares aborted - some invalid esds possible ***'
3950                    ifPrint = True
3951                break                   #refinement succeeded - finish up!
3952            except TypeError:          #result[1] is None on singular matrix
3953                print '**** Refinement failed - singular matrix ****'
3954                if 'Hessian' in Controls['deriv type']:
3955                    num = len(varyList)-1
3956                    for i,val in enumerate(np.flipud(result[2]['psing'])):
3957                        if val:
3958