source: trunk/GSASIIstruct.py @ 850

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

fix problem with atoms with An+ style of valence
more places for size results
more rigid body work
allow for skinnier images

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