source: trunk/GSASIIstruct.py @ 694

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

remove error check on use of variable in constraint & equivalence - this is now allowed; code commented out in 3 places in mapvars

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