source: trunk/GSASIIstruct.py @ 888

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

fix deriv for vector RB mags.

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