source: trunk/GSASIIstruct.py @ 874

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

add citation to banners

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