source: trunk/GSASIIstrIO.py @ 2483

Last change on this file since 2483 was 2483, checked in by vondreele, 5 years ago

fix site symmetry issues for magnetic moments on special positions

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 144.9 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2016-10-03 18:24:31 +0000 (Mon, 03 Oct 2016) $
4# $Author: vondreele $
5# $Revision: 2483 $
6# $URL: trunk/GSASIIstrIO.py $
7# $Id: GSASIIstrIO.py 2483 2016-10-03 18:24:31Z vondreele $
8########### SVN repository information ###################
9'''
10*GSASIIstrIO: structure I/O routines*
11-------------------------------------
12
13'''
14import sys
15import os
16import os.path as ospath
17import time
18import math
19import copy
20import random
21import cPickle
22import numpy as np
23import numpy.ma as ma
24import numpy.linalg as nl
25import scipy.optimize as so
26import GSASIIpath
27GSASIIpath.SetVersionNumber("$Revision: 2483 $")
28import GSASIIElem as G2el
29import GSASIIgrid as G2gd
30import GSASIIlattice as G2lat
31import GSASIIspc as G2spc
32import GSASIIobj as G2obj
33import GSASIImapvars as G2mv
34import GSASIImath as G2mth
35
36sind = lambda x: np.sin(x*np.pi/180.)
37cosd = lambda x: np.cos(x*np.pi/180.)
38tand = lambda x: np.tan(x*np.pi/180.)
39asind = lambda x: 180.*np.arcsin(x)/np.pi
40acosd = lambda x: 180.*np.arccos(x)/np.pi
41atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
42   
43ateln2 = 8.0*math.log(2.0)
44
45def GetControls(GPXfile):
46    ''' Returns dictionary of control items found in GSASII gpx file
47
48    :param str GPXfile: full .gpx file name
49    :return: dictionary of control items
50    '''
51    Controls = copy.copy(G2obj.DefaultControls)
52    fl = open(GPXfile,'rb')
53    while True:
54        try:
55            data = cPickle.load(fl)
56        except EOFError:
57            break
58        datum = data[0]
59        if datum[0] == 'Controls':
60            Controls.update(datum[1])
61    fl.close()
62    return Controls
63   
64def GetConstraints(GPXfile):
65    '''Read the constraints from the GPX file and interpret them
66
67    called in :func:`ReadCheckConstraints`, :func:`GSASIIstrMain.Refine`
68    and :func:`GSASIIstrMain.SeqRefine`.
69    '''
70    fl = open(GPXfile,'rb')
71    while True:
72        try:
73            data = cPickle.load(fl)
74        except EOFError:
75            break
76        datum = data[0]
77        if datum[0] == 'Constraints':
78            constList = []
79            for item in datum[1]:
80                if item.startswith('_'): continue
81                constList += datum[1][item]
82            fl.close()
83            constDict,fixedList,ignored = ProcessConstraints(constList)
84            if ignored:
85                print ignored,'Constraints were rejected. Was a constrained phase, histogram or atom deleted?'
86            return constDict,fixedList
87    fl.close()
88    raise Exception,"No constraints in GPX file"
89   
90def ProcessConstraints(constList):
91    """Interpret the constraints in the constList input into a dictionary, etc.
92    All :class:`GSASIIobj.G2VarObj` objects are mapped to the appropriate
93    phase/hist/atoms based on the object internals (random Ids). If this can't be
94    done (if a phase has been deleted, etc.), the variable is ignored.
95    If the constraint cannot be used due to too many dropped variables,
96    it is counted as ignored.
97   
98    :param list constList: a list of lists where each item in the outer list
99      specifies a constraint of some form, as described in the :mod:`GSASIIobj`
100      :ref:`Constraint definition <Constraint_definitions_table>`.
101
102    :returns:  a tuple of (constDict,fixedList,ignored) where:
103     
104      * constDict (list of dicts) contains the constraint relationships
105      * fixedList (list) contains the fixed values for each type
106        of constraint.
107      * ignored (int) counts the number of invalid constraint items
108        (should always be zero!)
109    """
110    constDict = []
111    fixedList = []
112    ignored = 0
113    for item in constList:
114        if item[-1] == 'h':
115            # process a hold
116            fixedList.append('0')
117            var = str(item[0][1])
118            if '?' not in var:
119                constDict.append({var:0.0})
120            else:
121                ignored += 1
122        elif item[-1] == 'f':
123            # process a new variable
124            fixedList.append(None)
125            D = {}
126            varyFlag = item[-2]
127            varname = item[-3]
128            for term in item[:-3]:
129                var = str(term[1])
130                if '?' not in var:
131                    D[var] = term[0]
132            if len(D) > 1:
133                # add extra dict terms for input variable name and vary flag
134                if varname is not None:                   
135                    if varname.startswith('::'):
136                        varname = varname[2:].replace(':',';')
137                    else:
138                        varname = varname.replace(':',';')
139                    D['_name'] = '::' + varname
140                D['_vary'] = varyFlag == True # force to bool
141                constDict.append(D)
142            else:
143                ignored += 1
144            #constFlag[-1] = ['Vary']
145        elif item[-1] == 'c': 
146            # process a contraint relationship
147            D = {}
148            for term in item[:-3]:
149                var = str(term[1])
150                if '?' not in var:
151                    D[var] = term[0]
152            if len(D) >= 1:
153                fixedList.append(str(item[-3]))
154                constDict.append(D)
155            else:
156                ignored += 1
157        elif item[-1] == 'e':
158            # process an equivalence
159            firstmult = None
160            eqlist = []
161            for term in item[:-3]:
162                if term[0] == 0: term[0] = 1.0
163                var = str(term[1])
164                if '?' in var: continue
165                if firstmult is None:
166                    firstmult = term[0]
167                    firstvar = var
168                else:
169                    eqlist.append([var,firstmult/term[0]])
170            if len(eqlist) > 0:
171                G2mv.StoreEquivalence(firstvar,eqlist)
172            else:
173                ignored += 1
174        else:
175            ignored += 1
176    return constDict,fixedList,ignored
177
178def ReadCheckConstraints(GPXfile):
179    '''Load constraints and related info and return any error or warning messages'''
180    # init constraints
181    G2mv.InitVars()   
182    # get variables
183    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
184    if not Phases:
185        return 'Error: No phases or no histograms in phases!',''
186    if not Histograms:
187        return 'Error: no diffraction data',''
188    constrDict,fixedList = GetConstraints(GPXfile) # load user constraints before internally generated ones
189    rigidbodyDict = GetRigidBodies(GPXfile)
190    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
191    rbVary,rbDict = GetRigidBodyModels(rigidbodyDict,Print=False)
192    Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,maxSSwave = \
193        GetPhaseData(Phases,RestraintDict=None,rbIds=rbIds,Print=False) # generates atom symmetry constraints
194    hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms,Print=False)
195    histVary,histDict,controlDict = GetHistogramData(Histograms,Print=False)
196    varyList = rbVary+phaseVary+hapVary+histVary
197    errmsg, warnmsg = G2mv.CheckConstraints(varyList,constrDict,fixedList)
198    if errmsg:
199        # print some diagnostic info on the constraints
200        print('Error in constraints:\n'+errmsg+
201              '\nRefinement not possible due to conflict in constraints, see below:\n')
202        print G2mv.VarRemapShow(varyList,True)
203    return errmsg, warnmsg
204   
205def makeTwinFrConstr(Phases,Histograms,hapVary):
206    TwConstr = []
207    TwFixed = []
208    for Phase in Phases:
209        pId = Phases[Phase]['pId']
210        for Histogram in Phases[Phase]['Histograms']:
211            try:
212                hId = Histograms[Histogram]['hId']
213                phfx = '%d:%d:'%(pId,hId)
214                if phfx+'TwinFr:0' in hapVary:
215                    TwFixed.append('1.0')     #constraint value
216                    nTwin = len(Phases[Phase]['Histograms'][Histogram]['Twins'])
217                    TwConstr.append({phfx+'TwinFr:'+str(i):'1.0' for i in range(nTwin)})
218            except KeyError:    #unused histograms?
219                pass
220    return TwConstr,TwFixed   
221   
222def GetRestraints(GPXfile):
223    '''Read the restraints from the GPX file.
224    Throws an exception if not found in the .GPX file
225    '''
226    fl = open(GPXfile,'rb')
227    while True:
228        try:
229            data = cPickle.load(fl)
230        except EOFError:
231            break
232        datum = data[0]
233        if datum[0] == 'Restraints':
234            restraintDict = datum[1]
235    fl.close()
236    return restraintDict
237   
238def GetRigidBodies(GPXfile):
239    '''Read the rigid body models from the GPX file
240    '''
241    fl = open(GPXfile,'rb')
242    while True:
243        try:
244            data = cPickle.load(fl)
245        except EOFError:
246            break
247        datum = data[0]
248        if datum[0] == 'Rigid bodies':
249            rigidbodyDict = datum[1]
250    fl.close()
251    return rigidbodyDict
252       
253def GetFprime(controlDict,Histograms):
254    'Needs a doc string'
255    FFtables = controlDict['FFtables']
256    if not FFtables:
257        return
258    histoList = Histograms.keys()
259    histoList.sort()
260    for histogram in histoList:
261        if histogram[:4] in ['PWDR','HKLF']:
262            Histogram = Histograms[histogram]
263            hId = Histogram['hId']
264            hfx = ':%d:'%(hId)
265            if 'X' in controlDict[hfx+'histType']:
266                keV = controlDict[hfx+'keV']
267                for El in FFtables:
268                    Orbs = G2el.GetXsectionCoeff(El.split('+')[0].split('-')[0])
269                    FP,FPP,Mu = G2el.FPcalc(Orbs, keV)
270                    FFtables[El][hfx+'FP'] = FP
271                    FFtables[El][hfx+'FPP'] = FPP
272                   
273def PrintFprime(FFtables,pfx,pFile):
274    print >>pFile,'\n Resonant form factors:'
275    Elstr = ' Element:'
276    FPstr = " f'     :"
277    FPPstr = ' f"     :'
278    for El in FFtables:
279        Elstr += ' %8s'%(El)
280        FPstr += ' %8.3f'%(FFtables[El][pfx+'FP'])
281        FPPstr += ' %8.3f'%(FFtables[El][pfx+'FPP'])
282    print >>pFile,Elstr
283    print >>pFile,FPstr
284    print >>pFile,FPPstr
285           
286def GetPhaseNames(GPXfile):
287    ''' Returns a list of phase names found under 'Phases' in GSASII gpx file
288
289    :param str GPXfile: full .gpx file name
290    :return: list of phase names
291    '''
292    fl = open(GPXfile,'rb')
293    PhaseNames = []
294    while True:
295        try:
296            data = cPickle.load(fl)
297        except EOFError:
298            break
299        datum = data[0]
300        if 'Phases' == datum[0]:
301            for datus in data[1:]:
302                PhaseNames.append(datus[0])
303    fl.close()
304    return PhaseNames
305
306def GetAllPhaseData(GPXfile,PhaseName):
307    ''' Returns the entire dictionary for PhaseName from GSASII gpx file
308
309    :param str GPXfile: full .gpx file name
310    :param str PhaseName: phase name
311    :return: phase dictionary
312    '''       
313    fl = open(GPXfile,'rb')
314    General = {}
315    Atoms = []
316    while True:
317        try:
318            data = cPickle.load(fl)
319        except EOFError:
320            break
321        datum = data[0]
322        if 'Phases' == datum[0]:
323            for datus in data[1:]:
324                if datus[0] == PhaseName:
325                    break
326    fl.close()
327    return datus[1]
328   
329def GetHistograms(GPXfile,hNames):
330    """ Returns a dictionary of histograms found in GSASII gpx file
331
332    :param str GPXfile: full .gpx file name
333    :param str hNames: list of histogram names
334    :return: dictionary of histograms (types = PWDR & HKLF)
335
336    """
337    fl = open(GPXfile,'rb')
338    Histograms = {}
339    while True:
340        try:
341            data = cPickle.load(fl)
342        except EOFError:
343            break
344        datum = data[0]
345        hist = datum[0]
346        if hist in hNames:
347            if 'PWDR' in hist[:4]:
348                PWDRdata = {}
349                PWDRdata.update(datum[1][0])        #weight factor
350                PWDRdata['Data'] = ma.array(ma.getdata(datum[1][1]))          #masked powder data arrays/clear previous masks
351                PWDRdata[data[2][0]] = data[2][1]       #Limits & excluded regions (if any)
352                PWDRdata[data[3][0]] = data[3][1]       #Background
353                PWDRdata[data[4][0]] = data[4][1]       #Instrument parameters
354                PWDRdata[data[5][0]] = data[5][1]       #Sample parameters
355                try:
356                    PWDRdata[data[9][0]] = data[9][1]       #Reflection lists might be missing
357                except IndexError:
358                    PWDRdata['Reflection Lists'] = {}
359                PWDRdata['Residuals'] = {}
360   
361                Histograms[hist] = PWDRdata
362            elif 'HKLF' in hist[:4]:
363                HKLFdata = {}
364                HKLFdata.update(datum[1][0])        #weight factor
365#patch
366                if 'list' in str(type(datum[1][1])):
367                #if isinstance(datum[1][1],list):
368                    RefData = {'RefList':[],'FF':{}}
369                    for ref in datum[1][1]:
370                        RefData['RefList'].append(ref[:11]+[ref[13],])
371                    RefData['RefList'] = np.array(RefData['RefList'])
372                    datum[1][1] = RefData
373#end patch
374                datum[1][1]['FF'] = {}
375                HKLFdata['Data'] = datum[1][1]
376                HKLFdata[data[1][0]] = data[1][1]       #Instrument parameters
377                HKLFdata['Reflection Lists'] = None
378                HKLFdata['Residuals'] = {}
379                Histograms[hist] = HKLFdata           
380    fl.close()
381    return Histograms
382   
383def GetHistogramNames(GPXfile,hType):
384    """ Returns a list of histogram names found in GSASII gpx file
385
386    :param str GPXfile: full .gpx file name
387    :param str hType: list of histogram types
388    :return: list of histogram names (types = PWDR & HKLF)
389
390    """
391    fl = open(GPXfile,'rb')
392    HistogramNames = []
393    while True:
394        try:
395            data = cPickle.load(fl)
396        except EOFError:
397            break
398        datum = data[0]
399        if datum[0][:4] in hType:
400            HistogramNames.append(datum[0])
401    fl.close()
402    return HistogramNames
403   
404def GetUsedHistogramsAndPhases(GPXfile):
405    ''' Returns all histograms that are found in any phase
406    and any phase that uses a histogram. This also
407    assigns numbers to used phases and histograms by the
408    order they appear in the file.
409
410    :param str GPXfile: full .gpx file name
411    :returns: (Histograms,Phases)
412
413     * Histograms = dictionary of histograms as {name:data,...}
414     * Phases = dictionary of phases that use histograms
415
416    '''
417    phaseNames = GetPhaseNames(GPXfile)
418    histoList = GetHistogramNames(GPXfile,['PWDR','HKLF'])
419    allHistograms = GetHistograms(GPXfile,histoList)
420    phaseData = {}
421    for name in phaseNames: 
422        phaseData[name] =  GetAllPhaseData(GPXfile,name)
423    Histograms = {}
424    Phases = {}
425    for phase in phaseData:
426        Phase = phaseData[phase]
427        if Phase['General']['Type'] == 'faulted': continue      #don't use faulted phases!
428        if Phase['Histograms']:
429            if phase not in Phases:
430                pId = phaseNames.index(phase)
431                Phase['pId'] = pId
432                Phases[phase] = Phase
433            for hist in Phase['Histograms']:
434                if 'Use' not in Phase['Histograms'][hist]:      #patch
435                    Phase['Histograms'][hist]['Use'] = True         
436                if hist not in Histograms and Phase['Histograms'][hist]['Use']:
437                    try:
438                        Histograms[hist] = allHistograms[hist]
439                        hId = histoList.index(hist)
440                        Histograms[hist]['hId'] = hId
441                    except KeyError: # would happen if a referenced histogram were
442                        # renamed or deleted
443                        print('For phase "'+str(phase)+
444                              '" unresolved reference to histogram "'+str(hist)+'"')
445    G2obj.IndexAllIds(Histograms=Histograms,Phases=Phases)
446    return Histograms,Phases
447   
448def getBackupName(GPXfile,makeBack):
449    '''
450    Get the name for the backup .gpx file name
451   
452    :param str GPXfile: full .gpx file name
453    :param bool makeBack: if True the name of a new file is returned, if
454      False the name of the last file that exists is returned
455    :returns: the name of a backup file
456   
457    '''
458    GPXpath,GPXname = ospath.split(GPXfile)
459    if GPXpath == '': GPXpath = '.'
460    Name = ospath.splitext(GPXname)[0]
461    files = os.listdir(GPXpath)
462    last = 0
463    for name in files:
464        name = name.split('.')
465        if len(name) == 3 and name[0] == Name and 'bak' in name[1]:
466            if makeBack:
467                last = max(last,int(name[1].strip('bak'))+1)
468            else:
469                last = max(last,int(name[1].strip('bak')))
470    GPXback = ospath.join(GPXpath,ospath.splitext(GPXname)[0]+'.bak'+str(last)+'.gpx')
471    return GPXback   
472       
473def GPXBackup(GPXfile,makeBack=True):
474    '''
475    makes a backup of the current .gpx file (?)
476   
477    :param str GPXfile: full .gpx file name
478    :param bool makeBack: if True (default), the backup is written to
479      a new file; if False, the last backup is overwritten
480    :returns: the name of the backup file that was written
481    '''
482    import distutils.file_util as dfu
483    GPXback = getBackupName(GPXfile,makeBack)
484    tries = 0
485    while True:
486        try:
487            dfu.copy_file(GPXfile,GPXback)
488            break
489        except:
490            tries += 1
491            if tries > 10:
492                return GPXfile  #failed!
493            time.sleep(1)           #just wait a second!         
494    return GPXback
495
496def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,RigidBodies,CovData,makeBack=True):
497    ''' Updates gpxfile from all histograms that are found in any phase
498    and any phase that used a histogram. Also updates rigid body definitions.
499
500
501    :param str GPXfile: full .gpx file name
502    :param dict Histograms: dictionary of histograms as {name:data,...}
503    :param dict Phases: dictionary of phases that use histograms
504    :param dict RigidBodies: dictionary of rigid bodies
505    :param dict CovData: dictionary of refined variables, varyList, & covariance matrix
506    :param bool makeBack: True if new backup of .gpx file is to be made; else use the last one made
507
508    '''
509                       
510    import distutils.file_util as dfu
511    GPXback = GPXBackup(GPXfile,makeBack)
512    print 'Read from file:',GPXback
513    print 'Save to file  :',GPXfile
514    infile = open(GPXback,'rb')
515    outfile = open(GPXfile,'wb')
516    while True:
517        try:
518            data = cPickle.load(infile)
519        except EOFError:
520            break
521        datum = data[0]
522#        print 'read: ',datum[0]
523        if datum[0] == 'Phases':
524            for iphase in range(len(data)):
525                if data[iphase][0] in Phases:
526                    phaseName = data[iphase][0]
527                    data[iphase][1].update(Phases[phaseName])
528        elif datum[0] == 'Covariance':
529            data[0][1] = CovData
530        elif datum[0] == 'Rigid bodies':
531            data[0][1] = RigidBodies
532        try:
533            histogram = Histograms[datum[0]]
534#            print 'found ',datum[0]
535            data[0][1][1] = histogram['Data']
536            data[0][1][0].update(histogram['Residuals'])
537            for datus in data[1:]:
538#                print '    read: ',datus[0]
539                if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
540                    datus[1] = histogram[datus[0]]
541        except KeyError:
542            pass
543        try:                       
544            cPickle.dump(data,outfile,1)
545        except AttributeError:
546            print 'ERROR - bad data in least squares result'
547            infile.close()
548            outfile.close()
549            dfu.copy_file(GPXback,GPXfile)
550            print 'GPX file save failed - old version retained'
551            return
552           
553    print 'GPX file save successful'
554   
555def SetSeqResult(GPXfile,Histograms,SeqResult):
556    '''
557    Needs doc string
558   
559    :param str GPXfile: full .gpx file name
560    '''
561    GPXback = GPXBackup(GPXfile)
562    print 'Read from file:',GPXback
563    print 'Save to file  :',GPXfile
564    infile = open(GPXback,'rb')
565    outfile = open(GPXfile,'wb')
566    while True:
567        try:
568            data = cPickle.load(infile)
569        except EOFError:
570            break
571        datum = data[0]
572        if datum[0] == 'Sequential results':
573            data[0][1] = SeqResult
574        # reset the Copy Next flag, since it should not be needed twice in a row
575        if datum[0] == 'Controls':
576            data[0][1]['Copy2Next'] = False
577        try:
578            histogram = Histograms[datum[0]]
579            data[0][1][1] = list(histogram['Data'])
580            for datus in data[1:]:
581                if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
582                    datus[1] = histogram[datus[0]]
583        except KeyError:
584            pass
585                               
586        cPickle.dump(data,outfile,1)
587    infile.close()
588    outfile.close()
589    print 'GPX file save successful'
590                       
591def ShowBanner(pFile=None):
592    'Print authorship, copyright and citation notice'
593    print >>pFile,80*'*'
594    print >>pFile,'   General Structure Analysis System-II Crystal Structure Refinement'
595    print >>pFile,'              by Robert B. Von Dreele & Brian H. Toby'
596    print >>pFile,'                Argonne National Laboratory(C), 2010'
597    print >>pFile,' This product includes software developed by the UChicago Argonne, LLC,' 
598    print >>pFile,'            as Operator of Argonne National Laboratory.'
599    print >>pFile,'                          Please cite:'
600    print >>pFile,'   B.H. Toby & R.B. Von Dreele, J. Appl. Cryst. 46, 544-549 (2013)'
601
602    print >>pFile,80*'*','\n'
603
604def ShowControls(Controls,pFile=None,SeqRef=False):
605    'Print controls information'
606    print >>pFile,' Least squares controls:'
607    print >>pFile,' Refinement type: ',Controls['deriv type']
608    if 'Hessian' in Controls['deriv type']:
609        print >>pFile,' Maximum number of cycles:',Controls['max cyc']
610    else:
611        print >>pFile,' Minimum delta-M/M for convergence: ','%.2g'%(Controls['min dM/M'])
612    print >>pFile,' Regularize hydrogens (if any):',Controls.get('HatomFix',False)
613    print >>pFile,' Initial shift factor: ','%.3f'%(Controls['shift factor'])
614    if SeqRef:
615        print >>pFile,' Sequential refinement controls:'
616        print >>pFile,' Copy of histogram results to next: ',Controls['Copy2Next']
617        print >>pFile,' Process histograms in reverse order: ',Controls['Reverse Seq']
618   
619def GetPawleyConstr(SGLaue,PawleyRef,im,pawleyVary):
620    'needs a doc string'
621#    if SGLaue in ['-1','2/m','mmm']:
622#        return                      #no Pawley symmetry required constraints
623    eqvDict = {}
624    for i,varyI in enumerate(pawleyVary):
625        eqvDict[varyI] = []
626        refI = int(varyI.split(':')[-1])
627        ih,ik,il = PawleyRef[refI][:3]
628        dspI = PawleyRef[refI][4+im]
629        for varyJ in pawleyVary[i+1:]:
630            refJ = int(varyJ.split(':')[-1])
631            jh,jk,jl = PawleyRef[refJ][:3]
632            dspJ = PawleyRef[refJ][4+im]
633            if SGLaue in ['4/m','4/mmm']:
634                isum = ih**2+ik**2
635                jsum = jh**2+jk**2
636                if abs(il) == abs(jl) and isum == jsum:
637                    eqvDict[varyI].append(varyJ) 
638            elif SGLaue in ['3R','3mR']:
639                isum = ih**2+ik**2+il**2
640                jsum = jh**2+jk**2+jl**2
641                isum2 = ih*ik+ih*il+ik*il
642                jsum2 = jh*jk+jh*jl+jk*jl
643                if isum == jsum and isum2 == jsum2:
644                    eqvDict[varyI].append(varyJ) 
645            elif SGLaue in ['3','3m1','31m','6/m','6/mmm']:
646                isum = ih**2+ik**2+ih*ik
647                jsum = jh**2+jk**2+jh*jk
648                if abs(il) == abs(jl) and isum == jsum:
649                    eqvDict[varyI].append(varyJ) 
650            elif SGLaue in ['m3','m3m']:
651                isum = ih**2+ik**2+il**2
652                jsum = jh**2+jk**2+jl**2
653                if isum == jsum:
654                    eqvDict[varyI].append(varyJ)
655            elif abs(dspI-dspJ)/dspI < 1.e-4:
656                eqvDict[varyI].append(varyJ) 
657    for item in pawleyVary:
658        if eqvDict[item]:
659            for item2 in pawleyVary:
660                if item2 in eqvDict[item]:
661                    eqvDict[item2] = []
662            G2mv.StoreEquivalence(item,eqvDict[item])
663                   
664def cellVary(pfx,SGData): 
665    'needs a doc string'
666    if SGData['SGLaue'] in ['-1',]:
667        return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
668    elif SGData['SGLaue'] in ['2/m',]:
669        if SGData['SGUniq'] == 'a':
670            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A5']
671        elif SGData['SGUniq'] == 'b':
672            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A4']
673        else:
674            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3']
675    elif SGData['SGLaue'] in ['mmm',]:
676        return [pfx+'A0',pfx+'A1',pfx+'A2']
677    elif SGData['SGLaue'] in ['4/m','4/mmm']:
678        G2mv.StoreEquivalence(pfx+'A0',(pfx+'A1',))
679        return [pfx+'A0',pfx+'A1',pfx+'A2']
680    elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
681        G2mv.StoreEquivalence(pfx+'A0',(pfx+'A1',pfx+'A3',))
682        return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3']
683    elif SGData['SGLaue'] in ['3R', '3mR']:
684        G2mv.StoreEquivalence(pfx+'A0',(pfx+'A1',pfx+'A2',))
685        G2mv.StoreEquivalence(pfx+'A3',(pfx+'A4',pfx+'A5',))
686        return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']                       
687    elif SGData['SGLaue'] in ['m3m','m3']:
688        G2mv.StoreEquivalence(pfx+'A0',(pfx+'A1',pfx+'A2',))
689        return [pfx+'A0',pfx+'A1',pfx+'A2']
690       
691def modVary(pfx,SSGData):
692    vary = []
693    for i,item in enumerate(SSGData['modSymb']):
694        if item in ['a','b','g']:
695            vary.append(pfx+'mV%d'%(i))
696    return vary
697       
698################################################################################
699##### Rigid Body Models and not General.get('doPawley')
700################################################################################
701       
702def GetRigidBodyModels(rigidbodyDict,Print=True,pFile=None):
703    'needs a doc string'
704   
705    def PrintResRBModel(RBModel):
706        atNames = RBModel['atNames']
707        rbRef = RBModel['rbRef']
708        rbSeq = RBModel['rbSeq']
709        print >>pFile,'Residue RB name: ',RBModel['RBname'],' no.atoms: ',len(RBModel['rbTypes']), \
710            'No. times used: ',RBModel['useCount']
711        print >>pFile,'    At name       x          y          z'
712        for name,xyz in zip(atNames,RBModel['rbXYZ']):
713            print >>pFile,%8s %10.4f %10.4f %10.4f'%(name,xyz[0],xyz[1],xyz[2])
714        print >>pFile,'Orientation defined by:',atNames[rbRef[0]],' -> ',atNames[rbRef[1]], \
715            ' & ',atNames[rbRef[0]],' -> ',atNames[rbRef[2]]
716        if rbSeq:
717            for i,rbseq in enumerate(rbSeq):
718                print >>pFile,'Torsion sequence ',i,' Bond: '+atNames[rbseq[0]],' - ', \
719                    atNames[rbseq[1]],' riding: ',[atNames[i] for i in rbseq[3]]
720       
721    def PrintVecRBModel(RBModel):
722        rbRef = RBModel['rbRef']
723        atTypes = RBModel['rbTypes']
724        print >>pFile,'Vector RB name: ',RBModel['RBname'],' no.atoms: ',len(RBModel['rbTypes']), \
725            'No. times used: ',RBModel['useCount']
726        for i in range(len(RBModel['VectMag'])):
727            print >>pFile,'Vector no.: ',i,' Magnitude: ', \
728                '%8.4f'%(RBModel['VectMag'][i]),' Refine? ',RBModel['VectRef'][i]
729            print >>pFile,'  No. Type     vx         vy         vz'
730            for j,[name,xyz] in enumerate(zip(atTypes,RBModel['rbVect'][i])):
731                print >>pFile,%d   %2s %10.4f %10.4f %10.4f'%(j,name,xyz[0],xyz[1],xyz[2])
732        print >>pFile,'  No. Type      x          y          z'
733        for i,[name,xyz] in enumerate(zip(atTypes,RBModel['rbXYZ'])):
734            print >>pFile,%d   %2s %10.4f %10.4f %10.4f'%(i,name,xyz[0],xyz[1],xyz[2])
735        print >>pFile,'Orientation defined by: atom ',rbRef[0],' -> atom ',rbRef[1], \
736            ' & atom ',rbRef[0],' -> atom ',rbRef[2]
737    rbVary = []
738    rbDict = {}
739    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
740    if len(rbIds['Vector']):
741        for irb,item in enumerate(rbIds['Vector']):
742            if rigidbodyDict['Vector'][item]['useCount']:
743                RBmags = rigidbodyDict['Vector'][item]['VectMag']
744                RBrefs = rigidbodyDict['Vector'][item]['VectRef']
745                for i,[mag,ref] in enumerate(zip(RBmags,RBrefs)):
746                    pid = '::RBV;'+str(i)+':'+str(irb)
747                    rbDict[pid] = mag
748                    if ref:
749                        rbVary.append(pid)
750                if Print:
751                    print >>pFile,'\nVector rigid body model:'
752                    PrintVecRBModel(rigidbodyDict['Vector'][item])
753    if len(rbIds['Residue']):
754        for item in rbIds['Residue']:
755            if rigidbodyDict['Residue'][item]['useCount']:
756                if Print:
757                    print >>pFile,'\nResidue rigid body model:'
758                    PrintResRBModel(rigidbodyDict['Residue'][item])
759    return rbVary,rbDict
760   
761def SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,pFile=None):
762    'needs a doc string'
763   
764    def PrintRBVectandSig(VectRB,VectSig):
765        print >>pFile,'\n Rigid body vector magnitudes for '+VectRB['RBname']+':'
766        namstr = '  names :'
767        valstr = '  values:'
768        sigstr = '  esds  :'
769        for i,[val,sig] in enumerate(zip(VectRB['VectMag'],VectSig)):
770            namstr += '%12s'%('Vect '+str(i))
771            valstr += '%12.4f'%(val)
772            if sig:
773                sigstr += '%12.4f'%(sig)
774            else:
775                sigstr += 12*' '
776        print >>pFile,namstr
777        print >>pFile,valstr
778        print >>pFile,sigstr       
779       
780    RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})  #these are lists of rbIds
781    if not RBIds['Vector']:
782        return
783    for irb,item in enumerate(RBIds['Vector']):
784        if rigidbodyDict['Vector'][item]['useCount']:
785            VectSig = []
786            RBmags = rigidbodyDict['Vector'][item]['VectMag']
787            for i,mag in enumerate(RBmags):
788                name = '::RBV;'+str(i)+':'+str(irb)
789                mag = parmDict[name]
790                if name in sigDict:
791                    VectSig.append(sigDict[name])
792            PrintRBVectandSig(rigidbodyDict['Vector'][item],VectSig)   
793       
794################################################################################
795##### Phase data
796################################################################################       
797                   
798def GetPhaseData(PhaseData,RestraintDict={},rbIds={},Print=True,pFile=None,seqRef=False):
799    'needs a doc string'
800           
801    def PrintFFtable(FFtable):
802        print >>pFile,'\n X-ray scattering factors:'
803        print >>pFile,'   Symbol     fa                                      fb                                      fc'
804        print >>pFile,99*'-'
805        for Ename in FFtable:
806            ffdata = FFtable[Ename]
807            fa = ffdata['fa']
808            fb = ffdata['fb']
809            print >>pFile,' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f' %  \
810                (Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],ffdata['fc'])
811               
812    def PrintMFtable(MFtable):
813        print >>pFile,'\n <j0> Magnetic scattering factors:'
814        print >>pFile,'   Symbol     mfa                                    mfb                                     mfc'
815        print >>pFile,99*'-'
816        for Ename in MFtable:
817            mfdata = MFtable[Ename]
818            fa = mfdata['mfa']
819            fb = mfdata['mfb']
820            print >>pFile,' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f' %  \
821                (Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],mfdata['mfc'])
822        print >>pFile,'\n <j2> Magnetic scattering factors:'
823        print >>pFile,'   Symbol     nfa                                    nfb                                     nfc'
824        print >>pFile,99*'-'
825        for Ename in MFtable:
826            mfdata = MFtable[Ename]
827            fa = mfdata['nfa']
828            fb = mfdata['nfb']
829            print >>pFile,' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f' %  \
830                (Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],mfdata['nfc'])
831               
832    def PrintBLtable(BLtable):
833        print >>pFile,'\n Neutron scattering factors:'
834        print >>pFile,'   Symbol   isotope       mass       b       resonant terms'
835        print >>pFile,99*'-'
836        for Ename in BLtable:
837            bldata = BLtable[Ename]
838            isotope = bldata[0]
839            mass = bldata[1]['Mass']
840            if 'SL' in bldata[1]:
841                blen = bldata[1]['SL'][0]
842                bres = []
843            else:
844                blen = 0
845                bres = bldata[1]['BW-LS']
846            line = ' %8s%11s %10.3f %8.3f'%(Ename.ljust(8),isotope.center(11),mass,blen)
847            for item in bres:
848                line += '%10.5g'%(item)
849            print >>pFile,line
850           
851    def PrintRBObjects(resRBData,vecRBData):
852       
853        def PrintRBThermals():
854            tlstr = ['11','22','33','12','13','23']
855            sstr = ['12','13','21','23','31','32','AA','BB']
856            TLS = RB['ThermalMotion'][1]
857            TLSvar = RB['ThermalMotion'][2]
858            if 'T' in RB['ThermalMotion'][0]:
859                print >>pFile,'TLS data'
860                text = ''
861                for i in range(6):
862                    text += 'T'+tlstr[i]+' %8.4f %s '%(TLS[i],str(TLSvar[i])[0])
863                print >>pFile,text
864                if 'L' in RB['ThermalMotion'][0]: 
865                    text = ''
866                    for i in range(6,12):
867                        text += 'L'+tlstr[i-6]+' %8.2f %s '%(TLS[i],str(TLSvar[i])[0])
868                    print >>pFile,text
869                if 'S' in RB['ThermalMotion'][0]:
870                    text = ''
871                    for i in range(12,20):
872                        text += 'S'+sstr[i-12]+' %8.3f %s '%(TLS[i],str(TLSvar[i])[0])
873                    print >>pFile,text
874            if 'U' in RB['ThermalMotion'][0]:
875                print >>pFile,'Uiso data'
876                text = 'Uiso'+' %10.3f %s'%(TLS[0],str(TLSvar[0])[0])           
877           
878        if len(resRBData):
879            for RB in resRBData:
880                Oxyz = RB['Orig'][0]
881                Qrijk = RB['Orient'][0]
882                Angle = 2.0*acosd(Qrijk[0])
883                print >>pFile,'\nRBObject ',RB['RBname'],' at ',      \
884                    '%10.4f %10.4f %10.4f'%(Oxyz[0],Oxyz[1],Oxyz[2]),' Refine?',RB['Orig'][1]
885                print >>pFile,'Orientation angle,vector:',      \
886                    '%10.3f %10.4f %10.4f %10.4f'%(Angle,Qrijk[1],Qrijk[2],Qrijk[3]),' Refine? ',RB['Orient'][1]
887                Torsions = RB['Torsions']
888                if len(Torsions):
889                    text = 'Torsions: '
890                    for torsion in Torsions:
891                        text += '%10.4f Refine? %s'%(torsion[0],torsion[1])
892                    print >>pFile,text
893                PrintRBThermals()
894        if len(vecRBData):
895            for RB in vecRBData:
896                Oxyz = RB['Orig'][0]
897                Qrijk = RB['Orient'][0]
898                Angle = 2.0*acosd(Qrijk[0])
899                print >>pFile,'\nRBObject ',RB['RBname'],' at ',      \
900                    '%10.4f %10.4f %10.4f'%(Oxyz[0],Oxyz[1],Oxyz[2]),' Refine?',RB['Orig'][1]           
901                print >>pFile,'Orientation angle,vector:',      \
902                    '%10.3f %10.4f %10.4f %10.4f'%(Angle,Qrijk[1],Qrijk[2],Qrijk[3]),' Refine? ',RB['Orient'][1]
903                PrintRBThermals()
904               
905    def PrintAtoms(General,Atoms):
906        cx,ct,cs,cia = General['AtomPtrs']
907        print >>pFile,'\n Atoms:'
908        line = '   name    type  refine?   x         y         z    '+ \
909            '  frac site sym  mult I/A   Uiso     U11     U22     U33     U12     U13     U23'
910        if General['Type'] == 'macromolecular':
911            line = ' res no residue chain'+line
912        print >>pFile,line
913        if General['Type'] in ['nuclear','magnetic','faulted',]:
914            print >>pFile,135*'-'
915            for i,at in enumerate(Atoms):
916                line = '%7s'%(at[ct-1])+'%7s'%(at[ct])+'%7s'%(at[ct+1])+'%10.5f'%(at[cx])+'%10.5f'%(at[cx+1])+ \
917                    '%10.5f'%(at[cx+2])+'%8.3f'%(at[cx+3])+'%7s'%(at[cs])+'%5d'%(at[cs+1])+'%5s'%(at[cia])
918                if at[cia] == 'I':
919                    line += '%8.5f'%(at[cia+1])+48*' '
920                else:
921                    line += 8*' '
922                    for j in range(6):
923                        line += '%8.5f'%(at[cia+2+j])
924                print >>pFile,line
925        elif General['Type'] == 'macromolecular':
926            print >>pFile,135*'-'           
927            for i,at in enumerate(Atoms):
928                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])+ \
929                    '%10.5f'%(at[cx+2])+'%8.3f'%(at[cx+3])+'%7s'%(at[cs])+'%5d'%(at[cs+1])+'%5s'%(at[cia])
930                if at[cia] == 'I':
931                    line += '%8.4f'%(at[cia+1])+48*' '
932                else:
933                    line += 8*' '
934                    for j in range(6):
935                        line += '%8.4f'%(at[cia+2+j])
936                print >>pFile,line
937               
938    def PrintMoments(General,Atoms):
939        cx,ct,cs,cia = General['AtomPtrs']
940        cmx = cx+4
941        AtInfo = dict(zip(General['AtomTypes'],General['Lande g']))
942        print >>pFile,'\n Magnetic moments:'
943        line = '   name    type  refine?  Mx        My        Mz    '
944        print >>pFile,line
945        print >>pFile,135*'-'
946        for i,at in enumerate(Atoms):
947            if AtInfo[at[ct]]:
948                line = '%7s'%(at[ct-1])+'%7s'%(at[ct])+'%7s'%(at[ct+1])+'%10.5f'%(at[cmx])+'%10.5f'%(at[cmx+1])+ \
949                    '%10.5f'%(at[cmx+2])
950                print >>pFile,line
951       
952               
953    def PrintWaves(General,Atoms):
954        cx,ct,cs,cia = General['AtomPtrs']
955        print >>pFile,'\n Modulation waves'
956        names = {'Sfrac':['Fsin','Fcos','Fzero','Fwid'],'Spos':['Xsin','Ysin','Zsin','Xcos','Ycos','Zcos','Tmin','Tmax','Xmax','Ymax','Zmax'],
957            'Sadp':['U11sin','U22sin','U33sin','U12sin','U13sin','U23sin','U11cos','U22cos',
958            'U33cos','U12cos','U13cos','U23cos'],'Smag':['MXsin','MYsin','MZsin','MXcos','MYcos','MZcos']}
959        print >>pFile,135*'-'
960        for i,at in enumerate(Atoms):
961            AtomSS = at[-1]['SS1']
962            for Stype in ['Sfrac','Spos','Sadp','Smag']:
963                Waves = AtomSS[Stype]
964                if len(Waves):
965                    print >>pFile,' atom: %s, site sym: %s, type: %s wave type: %s:'    \
966                        %(at[ct-1],at[cs],Stype,AtomSS['waveType'])
967                for iw,wave in enumerate(Waves):                   
968                    line = ''
969                    if AtomSS['waveType'] in ['Block','ZigZag'] and Stype == 'Spos' and not iw:
970                        for item in names[Stype][6:]:
971                            line += '%8s '%(item)                       
972                    else:
973                        if Stype == 'Spos':
974                            for item in names[Stype][:6]:
975                                line += '%8s '%(item)
976                        else:
977                            for item in names[Stype]:
978                                line += '%8s '%(item)
979                    print >>pFile,line
980                    line = ''
981                    for item in wave[0]:
982                        line += '%8.4f '%(item)
983                    line += ' Refine? '+str(wave[1])
984                    print >>pFile,line
985       
986    def PrintTexture(textureData):
987        topstr = '\n Spherical harmonics texture: Order:' + \
988            str(textureData['Order'])
989        if textureData['Order']:
990            print >>pFile,topstr+' Refine? '+str(textureData['SH Coeff'][0])
991        else:
992            print >>pFile,topstr
993            return
994        names = ['omega','chi','phi']
995        line = '\n'
996        for name in names:
997            line += ' SH '+name+':'+'%12.4f'%(textureData['Sample '+name][1])+' Refine? '+str(textureData['Sample '+name][0])
998        print >>pFile,line
999        print >>pFile,'\n Texture coefficients:'
1000        SHcoeff = textureData['SH Coeff'][1]
1001        SHkeys = SHcoeff.keys()
1002        nCoeff = len(SHcoeff)
1003        nBlock = nCoeff/10+1
1004        iBeg = 0
1005        iFin = min(iBeg+10,nCoeff)
1006        for block in range(nBlock):
1007            ptlbls = ' names :'
1008            ptstr =  ' values:'
1009            for item in SHkeys[iBeg:iFin]:
1010                ptlbls += '%12s'%(item)
1011                ptstr += '%12.4f'%(SHcoeff[item]) 
1012            print >>pFile,ptlbls
1013            print >>pFile,ptstr
1014            iBeg += 10
1015            iFin = min(iBeg+10,nCoeff)
1016       
1017    def MakeRBParms(rbKey,phaseVary,phaseDict):
1018        rbid = str(rbids.index(RB['RBId']))
1019        pfxRB = pfx+'RB'+rbKey+'P'
1020        pstr = ['x','y','z']
1021        ostr = ['a','i','j','k']
1022        for i in range(3):
1023            name = pfxRB+pstr[i]+':'+str(iRB)+':'+rbid
1024            phaseDict[name] = RB['Orig'][0][i]
1025            if RB['Orig'][1]:
1026                phaseVary += [name,]
1027        pfxRB = pfx+'RB'+rbKey+'O'
1028        for i in range(4):
1029            name = pfxRB+ostr[i]+':'+str(iRB)+':'+rbid
1030            phaseDict[name] = RB['Orient'][0][i]
1031            if RB['Orient'][1] == 'AV' and i:
1032                phaseVary += [name,]
1033            elif RB['Orient'][1] == 'A' and not i:
1034                phaseVary += [name,]
1035           
1036    def MakeRBThermals(rbKey,phaseVary,phaseDict):
1037        rbid = str(rbids.index(RB['RBId']))
1038        tlstr = ['11','22','33','12','13','23']
1039        sstr = ['12','13','21','23','31','32','AA','BB']
1040        if 'T' in RB['ThermalMotion'][0]:
1041            pfxRB = pfx+'RB'+rbKey+'T'
1042            for i in range(6):
1043                name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid
1044                phaseDict[name] = RB['ThermalMotion'][1][i]
1045                if RB['ThermalMotion'][2][i]:
1046                    phaseVary += [name,]
1047        if 'L' in RB['ThermalMotion'][0]:
1048            pfxRB = pfx+'RB'+rbKey+'L'
1049            for i in range(6):
1050                name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid
1051                phaseDict[name] = RB['ThermalMotion'][1][i+6]
1052                if RB['ThermalMotion'][2][i+6]:
1053                    phaseVary += [name,]
1054        if 'S' in RB['ThermalMotion'][0]:
1055            pfxRB = pfx+'RB'+rbKey+'S'
1056            for i in range(8):
1057                name = pfxRB+sstr[i]+':'+str(iRB)+':'+rbid
1058                phaseDict[name] = RB['ThermalMotion'][1][i+12]
1059                if RB['ThermalMotion'][2][i+12]:
1060                    phaseVary += [name,]
1061        if 'U' in RB['ThermalMotion'][0]:
1062            name = pfx+'RB'+rbKey+'U:'+str(iRB)+':'+rbid
1063            phaseDict[name] = RB['ThermalMotion'][1][0]
1064            if RB['ThermalMotion'][2][0]:
1065                phaseVary += [name,]
1066               
1067    def MakeRBTorsions(rbKey,phaseVary,phaseDict):
1068        rbid = str(rbids.index(RB['RBId']))
1069        pfxRB = pfx+'RB'+rbKey+'Tr;'
1070        for i,torsion in enumerate(RB['Torsions']):
1071            name = pfxRB+str(i)+':'+str(iRB)+':'+rbid
1072            phaseDict[name] = torsion[0]
1073            if torsion[1]:
1074                phaseVary += [name,]
1075                   
1076    if Print:
1077        print  >>pFile,'\n Phases:'
1078    phaseVary = []
1079    phaseDict = {}
1080    phaseConstr = {}
1081    pawleyLookup = {}
1082    FFtables = {}                   #scattering factors - xrays
1083    MFtables = {}                   #Mag. form factors
1084    BLtables = {}                   # neutrons
1085    Natoms = {}
1086    AtMults = {}
1087    AtIA = {}
1088    maxSSwave = {}
1089    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1090    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
1091    atomIndx = {}
1092    for name in PhaseData:
1093        General = PhaseData[name]['General']
1094        pId = PhaseData[name]['pId']
1095        pfx = str(pId)+'::'
1096        FFtable = G2el.GetFFtable(General['AtomTypes'])
1097        BLtable = G2el.GetBLtable(General)
1098        FFtables.update(FFtable)
1099        BLtables.update(BLtable)
1100        phaseDict[pfx+'isMag'] = False
1101        SGData = General['SGData']
1102        SGtext,SGtable = G2spc.SGPrint(SGData)
1103        if General['Type'] == 'magnetic':
1104            MFtable = G2el.GetMFtable(General['AtomTypes'],General['Lande g'])
1105            MFtables.update(MFtable)
1106            phaseDict[pfx+'isMag'] = True
1107            SpnFlp = SGData['SpnFlp']
1108        Atoms = PhaseData[name]['Atoms']
1109        if Atoms and not General.get('doPawley'):
1110            cx,ct,cs,cia = General['AtomPtrs']
1111            AtLookup = G2mth.FillAtomLookUp(Atoms,cia+8)
1112        PawleyRef = PhaseData[name].get('Pawley ref',[])
1113        cell = General['Cell']
1114        A = G2lat.cell2A(cell[1:7])
1115        phaseDict.update({pfx+'A0':A[0],pfx+'A1':A[1],pfx+'A2':A[2],
1116            pfx+'A3':A[3],pfx+'A4':A[4],pfx+'A5':A[5],pfx+'Vol':G2lat.calc_V(A)})
1117        if cell[0]:
1118            phaseVary += cellVary(pfx,SGData)       #also fills in symmetry required constraints
1119        SSGtext = []    #no superstructure
1120        im = 0
1121        if General.get('Modulated',False):
1122            im = 1      #refl offset
1123            Vec,vRef,maxH = General['SuperVec']
1124            phaseDict.update({pfx+'mV0':Vec[0],pfx+'mV1':Vec[1],pfx+'mV2':Vec[2]})
1125            SSGData = General['SSGData']
1126            SSGtext,SSGtable = G2spc.SSGPrint(SGData,SSGData)
1127            if vRef:
1128                phaseVary += modVary(pfx,SSGData)       
1129        resRBData = PhaseData[name]['RBModels'].get('Residue',[])
1130        if resRBData:
1131            rbids = rbIds['Residue']    #NB: used in the MakeRB routines
1132            for iRB,RB in enumerate(resRBData):
1133                MakeRBParms('R',phaseVary,phaseDict)
1134                MakeRBThermals('R',phaseVary,phaseDict)
1135                MakeRBTorsions('R',phaseVary,phaseDict)
1136       
1137        vecRBData = PhaseData[name]['RBModels'].get('Vector',[])
1138        if vecRBData:
1139            rbids = rbIds['Vector']    #NB: used in the MakeRB routines
1140            for iRB,RB in enumerate(vecRBData):
1141                MakeRBParms('V',phaseVary,phaseDict)
1142                MakeRBThermals('V',phaseVary,phaseDict)
1143                   
1144        Natoms[pfx] = 0
1145        maxSSwave[pfx] = {'Sfrac':0,'Spos':0,'Sadp':0,'Smag':0}
1146        if Atoms and not General.get('doPawley'):
1147            cx,ct,cs,cia = General['AtomPtrs']
1148            Natoms[pfx] = len(Atoms)
1149            for i,at in enumerate(Atoms):
1150                atomIndx[at[cia+8]] = [pfx,i]      #lookup table for restraints
1151                phaseDict.update({pfx+'Atype:'+str(i):at[ct],pfx+'Afrac:'+str(i):at[cx+3],pfx+'Amul:'+str(i):at[cs+1],
1152                    pfx+'Ax:'+str(i):at[cx],pfx+'Ay:'+str(i):at[cx+1],pfx+'Az:'+str(i):at[cx+2],
1153                    pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx+'dAz:'+str(i):0.,         #refined shifts for x,y,z
1154                    pfx+'AI/A:'+str(i):at[cia],})
1155                if at[cia] == 'I':
1156                    phaseDict[pfx+'AUiso:'+str(i)] = at[cia+1]
1157                else:
1158                    phaseDict.update({pfx+'AU11:'+str(i):at[cia+2],pfx+'AU22:'+str(i):at[cia+3],pfx+'AU33:'+str(i):at[cia+4],
1159                        pfx+'AU12:'+str(i):at[cia+5],pfx+'AU13:'+str(i):at[cia+6],pfx+'AU23:'+str(i):at[cia+7]})
1160                if General['Type'] == 'magnetic':
1161                    phaseDict.update({pfx+'AMx:'+str(i):at[cx+4],pfx+'AMy:'+str(i):at[cx+5],pfx+'AMz:'+str(i):at[cx+6]})
1162                if 'F' in at[ct+1]:
1163                    phaseVary.append(pfx+'Afrac:'+str(i))
1164                if 'X' in at[ct+1]:
1165                    try:    #patch for sytsym name changes
1166                        xId,xCoef = G2spc.GetCSxinel(at[cs])
1167                    except KeyError:
1168                        Sytsym = G2spc.SytSym(at[cx:cx+3],SGData)[0]
1169                        at[cs] = Sytsym
1170                        xId,xCoef = G2spc.GetCSxinel(at[cs])
1171                    xId,xCoef = G2spc.GetCSxinel(at[cs])
1172                    names = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)]
1173                    equivs = [[],[],[]]
1174                    for j in range(3):
1175                        if xId[j] > 0:                               
1176                            phaseVary.append(names[j])
1177                            equivs[xId[j]-1].append([names[j],xCoef[j]])
1178                    for equiv in equivs:
1179                        if len(equiv) > 1:
1180                            name = equiv[0][0]
1181                            coef = equiv[0][1]
1182                            for eqv in equiv[1:]:
1183                                eqv[1] /= coef
1184                                G2mv.StoreEquivalence(name,(eqv,))
1185                if 'U' in at[ct+1]:
1186                    if at[cia] == 'I':
1187                        phaseVary.append(pfx+'AUiso:'+str(i))
1188                    else:
1189                        try:    #patch for sytsym name changes
1190                            uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
1191                        except KeyError:
1192                            Sytsym = G2spc.SytSym(at[cx:cx+3],SGData)[0]
1193                            at[cs] = Sytsym
1194                            uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
1195                        uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
1196                        names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i),
1197                            pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)]
1198                        equivs = [[],[],[],[],[],[]]
1199                        for j in range(6):
1200                            if uId[j] > 0:                               
1201                                phaseVary.append(names[j])
1202                                equivs[uId[j]-1].append([names[j],uCoef[j]])
1203                        for equiv in equivs:
1204                            if len(equiv) > 1:
1205                                name = equiv[0][0]
1206                                coef = equiv[0][1]
1207                                for eqv in equiv[1:]:
1208                                    eqv[1] /= coef
1209                                G2mv.StoreEquivalence(name,equiv[1:])
1210                if 'M' in at[ct+1]:
1211                    SytSym,Mul,Nop,dupDir = G2spc.SytSym(at[cx:cx+3],SGData)
1212                    mId,mCoef = G2spc.GetCSpqinel(SytSym,SpnFlp,dupDir)
1213                    names = [pfx+'AMx:'+str(i),pfx+'AMy:'+str(i),pfx+'AMz:'+str(i)]
1214                    equivs = [[],[],[]]
1215                    for j in range(3):
1216                        if mId[j] > 0:
1217                            phaseVary.append(names[j])
1218                            equivs[mId[j]-1].append([names[j],mCoef[j]])
1219                    for equiv in equivs:
1220                        if len(equiv) > 1:
1221                            name = equiv[0][0]
1222                            coef = equiv[0][1]
1223                            for eqv in equiv[1:]:
1224                                eqv[1] /= coef
1225                                G2mv.StoreEquivalence(name,(eqv,))
1226                if General.get('Modulated',False):
1227                    AtomSS = at[-1]['SS1']
1228                    waveType = AtomSS['waveType']
1229                    phaseDict[pfx+'waveType:'+str(i)] = waveType
1230                    for Stype in ['Sfrac','Spos','Sadp','Smag']:
1231                        Waves = AtomSS[Stype]
1232                        nx = 0
1233                        for iw,wave in enumerate(Waves):
1234                            if not iw:
1235                                if waveType in ['ZigZag','Block']:
1236                                    nx = 1
1237                                CSI = G2spc.GetSSfxuinel(waveType,1,at[cx:cx+3],SGData,SSGData)
1238                            else:
1239                                CSI = G2spc.GetSSfxuinel('Fourier',iw+1-nx,at[cx:cx+3],SGData,SSGData)
1240                            uId,uCoef = CSI[Stype]
1241                            stiw = str(i)+':'+str(iw)
1242                            if Stype == 'Spos':
1243                                if waveType in ['ZigZag','Block',] and not iw:
1244                                    names = [pfx+'Tmin:'+stiw,pfx+'Tmax:'+stiw,pfx+'Xmax:'+stiw,pfx+'Ymax:'+stiw,pfx+'Zmax:'+stiw]
1245                                    equivs = [[],[], [],[],[]]
1246                                else:
1247                                    names = [pfx+'Xsin:'+stiw,pfx+'Ysin:'+stiw,pfx+'Zsin:'+stiw,
1248                                        pfx+'Xcos:'+stiw,pfx+'Ycos:'+stiw,pfx+'Zcos:'+stiw]
1249                                    equivs = [[],[],[], [],[],[]]
1250                            elif Stype == 'Sadp':
1251                                names = [pfx+'U11sin:'+stiw,pfx+'U22sin:'+stiw,pfx+'U33sin:'+stiw,
1252                                    pfx+'U12sin:'+stiw,pfx+'U13sin:'+stiw,pfx+'U23sin:'+stiw,
1253                                    pfx+'U11cos:'+stiw,pfx+'U22cos:'+stiw,pfx+'U33cos:'+stiw,
1254                                    pfx+'U12cos:'+stiw,pfx+'U13cos:'+stiw,pfx+'U23cos:'+stiw]
1255                                equivs = [[],[],[],[],[],[], [],[],[],[],[],[]]
1256                            elif Stype == 'Sfrac':
1257                                equivs = [[],[]]
1258                                if 'Crenel' in waveType and not iw:
1259                                    names = [pfx+'Fzero:'+stiw,pfx+'Fwid:'+stiw]
1260                                else:
1261                                    names = [pfx+'Fsin:'+stiw,pfx+'Fcos:'+stiw]
1262                            elif Stype == 'Smag':
1263                                equivs = [[],[],[], [],[],[]]
1264                                names = [pfx+'MXsin:'+stiw,pfx+'MYsin:'+stiw,pfx+'MZsin:'+stiw,
1265                                    pfx+'MXcos:'+stiw,pfx+'MYcos:'+stiw,pfx+'MZcos:'+stiw]
1266                            phaseDict.update(dict(zip(names,wave[0])))
1267                            if wave[1]: #what do we do here for multiple terms in modulation constraints?
1268                                for j in range(len(equivs)):
1269                                    if uId[j][0] > 0:                               
1270                                        phaseVary.append(names[j])
1271                                        equivs[uId[j][0]-1].append([names[j],uCoef[j][0]])
1272                                for equiv in equivs:
1273                                    if len(equiv) > 1:
1274                                        name = equiv[0][0]
1275                                        coef = equiv[0][1]
1276                                        for eqv in equiv[1:]:
1277                                            eqv[1] /= coef
1278                                        G2mv.StoreEquivalence(name,equiv[1:])
1279                            maxSSwave[pfx][Stype] = max(maxSSwave[pfx][Stype],iw+1)
1280            textureData = General['SH Texture']
1281            if textureData['Order'] and not seqRef:
1282                phaseDict[pfx+'SHorder'] = textureData['Order']
1283                phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']]
1284                for item in ['omega','chi','phi']:
1285                    phaseDict[pfx+'SH '+item] = textureData['Sample '+item][1]
1286                    if textureData['Sample '+item][0]:
1287                        phaseVary.append(pfx+'SH '+item)
1288                for item in textureData['SH Coeff'][1]:
1289                    phaseDict[pfx+item] = textureData['SH Coeff'][1][item]
1290                    if textureData['SH Coeff'][0]:
1291                        phaseVary.append(pfx+item)
1292               
1293            if Print:
1294                print >>pFile,'\n Phase name: ',General['Name']
1295                print >>pFile,135*'-'
1296                PrintFFtable(FFtable)
1297                PrintBLtable(BLtable)
1298                if General['Type'] == 'magnetic':
1299                    PrintMFtable(MFtable)
1300                print >>pFile,''
1301                #how do we print magnetic symmetry table? TBD
1302                if len(SSGtext):    #if superstructure
1303                    for line in SSGtext: print >>pFile,line
1304                    if len(SSGtable):
1305                        for item in SSGtable:
1306                            line = ' %s '%(item)
1307                            print >>pFile,line   
1308                    else:
1309                        print >>pFile,' ( 1)    %s'%(SSGtable[0])
1310                else:
1311                    for line in SGtext: print >>pFile,line
1312                    if len(SGtable):
1313                        for item in SGtable:
1314                            line = ' %s '%(item)
1315                            print >>pFile,line   
1316                    else:
1317                        print >>pFile,' ( 1)    %s'%(SGtable[0])
1318                PrintRBObjects(resRBData,vecRBData)
1319                PrintAtoms(General,Atoms)
1320                if General['Type'] == 'magnetic':
1321                    PrintMoments(General,Atoms)
1322                if General.get('Modulated',False):
1323                    PrintWaves(General,Atoms)
1324                print >>pFile,'\n Unit cell: a = %.5f'%(cell[1]),' b = %.5f'%(cell[2]),' c = %.5f'%(cell[3]), \
1325                    ' alpha = %.3f'%(cell[4]),' beta = %.3f'%(cell[5]),' gamma = %.3f'%(cell[6]), \
1326                    ' volume = %.3f'%(cell[7]),' Refine?',cell[0]
1327                if len(SSGtext):    #if superstructure
1328                    print >>pFile,'\n Modulation vector: mV0 = %.4f'%(Vec[0]),' mV1 = %.4f'%(Vec[1]),   \
1329                        ' mV2 = %.4f'%(Vec[2]),' max mod. index = %d'%(maxH),' Refine?',vRef
1330                if not seqRef:
1331                    PrintTexture(textureData)
1332                if name in RestraintDict:
1333                    PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
1334                        textureData,RestraintDict[name],pFile)
1335                   
1336        elif PawleyRef:
1337            if Print:
1338                print >>pFile,'\n Phase name: ',General['Name']
1339                print >>pFile,135*'-'
1340                print >>pFile,''
1341                if len(SSGtext):    #if superstructure
1342                    for line in SSGtext: print >>pFile,line
1343                    if len(SSGtable):
1344                        for item in SSGtable:
1345                            line = ' %s '%(item)
1346                            print >>pFile,line   
1347                    else:
1348                        print >>pFile,' ( 1)    %s'%(SSGtable[0])
1349                else:
1350                    for line in SGtext: print >>pFile,line
1351                    if len(SGtable):
1352                        for item in SGtable:
1353                            line = ' %s '%(item)
1354                            print >>pFile,line   
1355                    else:
1356                        print >>pFile,' ( 1)    %s'%(SGtable[0])
1357                print >>pFile,'\n Unit cell: a = %.5f'%(cell[1]),' b = %.5f'%(cell[2]),' c = %.5f'%(cell[3]), \
1358                    ' alpha = %.3f'%(cell[4]),' beta = %.3f'%(cell[5]),' gamma = %.3f'%(cell[6]), \
1359                    ' volume = %.3f'%(cell[7]),' Refine?',cell[0]
1360                if len(SSGtext):    #if superstructure
1361                    print >>pFile,'\n Modulation vector: mV0 = %.4f'%(Vec[0]),' mV1 = %.4f'%(Vec[1]),   \
1362                        ' mV2 = %.4f'%(Vec[2]),' max mod. index = %d'%(maxH),' Refine?',vRef
1363            pawleyVary = []
1364            for i,refl in enumerate(PawleyRef):
1365                phaseDict[pfx+'PWLref:'+str(i)] = refl[6+im]
1366                if im:
1367                    pawleyLookup[pfx+'%d,%d,%d,%d'%(refl[0],refl[1],refl[2],refl[3])] = i
1368                else:
1369                    pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
1370                if refl[5+im]:
1371                    pawleyVary.append(pfx+'PWLref:'+str(i))
1372            GetPawleyConstr(SGData['SGLaue'],PawleyRef,im,pawleyVary)      #does G2mv.StoreEquivalence
1373            phaseVary += pawleyVary
1374               
1375    return Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,maxSSwave
1376   
1377def cellFill(pfx,SGData,parmDict,sigDict): 
1378    '''Returns the filled-out reciprocal cell (A) terms and their uncertainties
1379    from the parameter and sig dictionaries.
1380
1381    :param str pfx: parameter prefix ("n::", where n is a phase number)
1382    :param dict SGdata: a symmetry object
1383    :param dict parmDict: a dictionary of parameters
1384    :param dict sigDict:  a dictionary of uncertainties on parameters
1385
1386    :returns: A,sigA where each is a list of six terms with the A terms
1387    '''
1388    if SGData['SGLaue'] in ['-1',]:
1389        A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1390            parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
1391    elif SGData['SGLaue'] in ['2/m',]:
1392        if SGData['SGUniq'] == 'a':
1393            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1394                0,0,parmDict[pfx+'A5']]
1395        elif SGData['SGUniq'] == 'b':
1396            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1397                0,parmDict[pfx+'A4'],0]
1398        else:
1399            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1400                parmDict[pfx+'A3'],0,0]
1401    elif SGData['SGLaue'] in ['mmm',]:
1402        A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
1403    elif SGData['SGLaue'] in ['4/m','4/mmm']:
1404        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
1405    elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1406        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
1407            parmDict[pfx+'A0'],0,0]
1408    elif SGData['SGLaue'] in ['3R', '3mR']:
1409        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
1410            parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
1411    elif SGData['SGLaue'] in ['m3m','m3']:
1412        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
1413
1414    try:
1415        if SGData['SGLaue'] in ['-1',]:
1416            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1417                sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
1418        elif SGData['SGLaue'] in ['2/m',]:
1419            if SGData['SGUniq'] == 'a':
1420                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1421                    0,0,sigDict[pfx+'A5']]
1422            elif SGData['SGUniq'] == 'b':
1423                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1424                    0,sigDict[pfx+'A4'],0]
1425            else:
1426                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1427                    sigDict[pfx+'A3'],0,0]
1428        elif SGData['SGLaue'] in ['mmm',]:
1429            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
1430        elif SGData['SGLaue'] in ['4/m','4/mmm']:
1431            sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1432        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1433            sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1434        elif SGData['SGLaue'] in ['3R', '3mR']:
1435            sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
1436        elif SGData['SGLaue'] in ['m3m','m3']:
1437            sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
1438    except KeyError:
1439        sigA = [0,0,0,0,0,0]
1440    return A,sigA
1441       
1442def PrintRestraints(cell,SGData,AtPtrs,Atoms,AtLookup,textureData,phaseRest,pFile):
1443    'needs a doc string'
1444    if phaseRest:
1445        Amat = G2lat.cell2AB(cell)[0]
1446        cx,ct,cs = AtPtrs[:3]
1447        names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
1448            ['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
1449            ['ChemComp','Sites'],['Texture','HKLs']]
1450        for name,rest in names:
1451            itemRest = phaseRest[name]
1452            if itemRest[rest] and itemRest['Use']:
1453                print >>pFile,'\n  %s %10.3f Use: %s'%(name+' restraint weight factor',itemRest['wtFactor'],str(itemRest['Use']))
1454                if name in ['Bond','Angle','Plane','Chiral']:
1455                    print >>pFile,'     calc       obs      sig   delt/sig  atoms(symOp): '
1456                    for indx,ops,obs,esd in itemRest[rest]:
1457                        try:
1458                            AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1459                            AtName = ''
1460                            for i,Aname in enumerate(AtNames):
1461                                AtName += Aname
1462                                if ops[i] == '1':
1463                                    AtName += '-'
1464                                else:
1465                                    AtName += '+('+ops[i]+')-'
1466                            XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1467                            XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1468                            if name == 'Bond':
1469                                calc = G2mth.getRestDist(XYZ,Amat)
1470                            elif name == 'Angle':
1471                                calc = G2mth.getRestAngle(XYZ,Amat)
1472                            elif name == 'Plane':
1473                                calc = G2mth.getRestPlane(XYZ,Amat)
1474                            elif name == 'Chiral':
1475                                calc = G2mth.getRestChiral(XYZ,Amat)
1476                            print >>pFile,' %9.3f %9.3f %8.3f %8.3f   %s'%(calc,obs,esd,(obs-calc)/esd,AtName[:-1])
1477                        except KeyError:
1478                            del itemRest[rest]
1479                elif name in ['Torsion','Rama']:
1480                    print >>pFile,'  atoms(symOp)  calc  obs  sig  delt/sig  torsions: '
1481                    coeffDict = itemRest['Coeff']
1482                    for indx,ops,cofName,esd in itemRest[rest]:
1483                        AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1484                        AtName = ''
1485                        for i,Aname in enumerate(AtNames):
1486                            AtName += Aname+'+('+ops[i]+')-'
1487                        XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1488                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1489                        if name == 'Torsion':
1490                            tor = G2mth.getRestTorsion(XYZ,Amat)
1491                            restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
1492                            print >>pFile,' %8.3f %8.3f %.3f %8.3f %8.3f %s'%(calc,obs,esd,(obs-calc)/esd,tor,AtName[:-1])
1493                        else:
1494                            phi,psi = G2mth.getRestRama(XYZ,Amat)
1495                            restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])                               
1496                            print >>pFile,' %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %s'%(calc,obs,esd,(obs-calc)/esd,phi,psi,AtName[:-1])
1497                elif name == 'ChemComp':
1498                    print >>pFile,'     atoms   mul*frac  factor     prod'
1499                    for indx,factors,obs,esd in itemRest[rest]:
1500                        try:
1501                            atoms = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1502                            mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs+1))
1503                            frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs-1))
1504                            mulfrac = mul*frac
1505                            calcs = mul*frac*factors
1506                            for iatm,[atom,mf,fr,clc] in enumerate(zip(atoms,mulfrac,factors,calcs)):
1507                                print >>pFile,' %10s %8.3f %8.3f %8.3f'%(atom,mf,fr,clc)
1508                            print >>pFile,' Sum:                   calc: %8.3f obs: %8.3f esd: %8.3f'%(np.sum(calcs),obs,esd)
1509                        except KeyError:
1510                            del itemRest[rest]
1511                elif name == 'Texture' and textureData['Order']:
1512                    Start = False
1513                    SHCoef = textureData['SH Coeff'][1]
1514                    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1515                    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
1516                    print '    HKL  grid  neg esd  sum neg  num neg use unit?  unit esd  '
1517                    for hkl,grid,esd1,ifesd2,esd2 in itemRest[rest]:
1518                        PH = np.array(hkl)
1519                        phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
1520                        ODFln = G2lat.Flnh(Start,SHCoef,phi,beta,SGData)
1521                        R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
1522                        Z = ma.masked_greater(Z,0.0)
1523                        num = ma.count(Z)
1524                        sum = 0
1525                        if num:
1526                            sum = np.sum(Z)
1527                        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)
1528       
1529def getCellEsd(pfx,SGData,A,covData):
1530    'needs a doc string'
1531    dpr = 180./np.pi
1532    rVsq = G2lat.calc_rVsq(A)
1533    G,g = G2lat.A2Gmat(A)       #get recip. & real metric tensors
1534    cell = np.array(G2lat.Gmat2cell(g))   #real cell
1535    cellst = np.array(G2lat.Gmat2cell(G)) #recip. cell
1536    scos = cosd(cellst[3:6])
1537    ssin = sind(cellst[3:6])
1538    scot = scos/ssin
1539    rcos = cosd(cell[3:6])
1540    rsin = sind(cell[3:6])
1541    rcot = rcos/rsin
1542    RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
1543    varyList = covData['varyList']
1544    covMatrix = covData['covMatrix']
1545    vcov = G2mth.getVCov(RMnames,varyList,covMatrix)
1546    Ax = np.array(A)
1547    Ax[3:] /= 2.
1548    drVdA = np.array([
1549        Ax[1]*Ax[2]-Ax[5]**2,
1550        Ax[0]*Ax[2]-Ax[4]**2,
1551        Ax[0]*Ax[1]-Ax[3]**2,
1552        Ax[4]*Ax[5]-Ax[2]*Ax[3],
1553        Ax[3]*Ax[5]-Ax[1]*Ax[4],
1554        Ax[3]*Ax[4]-Ax[0]*Ax[5]])
1555    srcvlsq = np.inner(drVdA,np.inner(drVdA,vcov))
1556    Vol = 1/np.sqrt(rVsq)
1557    sigVol = Vol**3*np.sqrt(srcvlsq)/2.         #ok - checks with GSAS
1558    R123 = Ax[0]*Ax[1]*Ax[2]
1559    dsasdg = np.zeros((3,6))
1560    dadg = np.zeros((6,6))
1561    for i0 in range(3):         #0  1   2
1562        i1 = (i0+1)%3           #1  2   0
1563        i2 = (i1+1)%3           #2  0   1
1564        i3 = 5-i2               #3  5   4
1565        i4 = 5-i1               #4  3   5
1566        i5 = 5-i0               #5  4   3
1567        dsasdg[i0][i1] = 0.5*scot[i0]*scos[i0]/Ax[i1]
1568        dsasdg[i0][i2] = 0.5*scot[i0]*scos[i0]/Ax[i2]
1569        dsasdg[i0][i5] = -scot[i0]/np.sqrt(Ax[i1]*Ax[i2])
1570        denmsq = Ax[i0]*(R123-Ax[i1]*Ax[i4]**2-Ax[i2]*Ax[i3]**2+(Ax[i4]*Ax[i3])**2)
1571        denom = np.sqrt(denmsq)
1572        dadg[i5][i0] = -Ax[i5]/denom-rcos[i0]/denmsq*(R123-0.5*Ax[i1]*Ax[i4]**2-0.5*Ax[i2]*Ax[i3]**2)
1573        dadg[i5][i1] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i2]-Ax[i0]*Ax[i4]**2)
1574        dadg[i5][i2] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i1]-Ax[i0]*Ax[i3]**2)
1575        dadg[i5][i3] = Ax[i4]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i2]*Ax[i3]-Ax[i3]*Ax[i4]**2)
1576        dadg[i5][i4] = Ax[i3]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i1]*Ax[i4]-Ax[i3]**2*Ax[i4])
1577        dadg[i5][i5] = -Ax[i0]/denom
1578    for i0 in range(3):
1579        i1 = (i0+1)%3
1580        i2 = (i1+1)%3
1581        i3 = 5-i2
1582        for ij in range(6):
1583            dadg[i0][ij] = cell[i0]*(rcot[i2]*dadg[i3][ij]/rsin[i2]-dsasdg[i1][ij]/ssin[i1])
1584            if ij == i0:
1585                dadg[i0][ij] = dadg[i0][ij]-0.5*cell[i0]/Ax[i0]
1586            dadg[i3][ij] = -dadg[i3][ij]*rsin[2-i0]*dpr
1587    sigMat = np.inner(dadg,np.inner(dadg,vcov))
1588    var = np.diag(sigMat)
1589    CS = np.where(var>0.,np.sqrt(var),0.)
1590    return [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol]
1591   
1592def SetPhaseData(parmDict,sigDict,Phases,RBIds,covData,RestraintDict=None,pFile=None):
1593    'needs a doc string'
1594   
1595    def PrintAtomsAndSig(General,Atoms,atomsSig):
1596        print >>pFile,'\n Atoms:'
1597        line = '   name      x         y         z      frac   Uiso     U11     U22     U33     U12     U13     U23'
1598        if General['Type'] == 'macromolecular':
1599            line = ' res no residue chain '+line
1600        cx,ct,cs,cia = General['AtomPtrs']
1601        print >>pFile,line
1602        print >>pFile,135*'-'
1603        fmt = {0:'%7s',ct:'%7s',cx:'%10.5f',cx+1:'%10.5f',cx+2:'%10.5f',cx+3:'%8.3f',cia+1:'%8.5f',
1604            cia+2:'%8.5f',cia+3:'%8.5f',cia+4:'%8.5f',cia+5:'%8.5f',cia+6:'%8.5f',cia+7:'%8.5f'}
1605        noFXsig = {cx:[10*' ','%10s'],cx+1:[10*' ','%10s'],cx+2:[10*' ','%10s'],cx+3:[8*' ','%8s']}
1606        for atyp in General['AtomTypes']:       #zero composition
1607            General['NoAtoms'][atyp] = 0.
1608        for i,at in enumerate(Atoms):
1609            General['NoAtoms'][at[ct]] += at[cx+3]*float(at[cx+5])     #new composition
1610            if General['Type'] == 'macromolecular':
1611                name = ' %s %s %s %s:'%(at[0],at[1],at[2],at[3])
1612                valstr = ' values:          '
1613                sigstr = ' sig   :          '
1614            else:
1615                name = fmt[0]%(at[ct-1])+fmt[1]%(at[ct])+':'
1616                valstr = ' values:'
1617                sigstr = ' sig   :'
1618            for ind in range(cx,cx+4):
1619                sigind = str(i)+':'+str(ind)
1620                valstr += fmt[ind]%(at[ind])                   
1621                if sigind in atomsSig:
1622                    sigstr += fmt[ind]%(atomsSig[sigind])
1623                else:
1624                    sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
1625            if at[cia] == 'I':
1626                valstr += fmt[cia+1]%(at[cia+1])
1627                if '%d:%d'%(i,cia+1) in atomsSig:
1628                    sigstr += fmt[cia+1]%(atomsSig['%d:%d'%(i,cia+1)])
1629                else:
1630                    sigstr += 8*' '
1631            else:
1632                valstr += 8*' '
1633                sigstr += 8*' '
1634                for ind in range(cia+2,cia+8):
1635                    sigind = str(i)+':'+str(ind)
1636                    valstr += fmt[ind]%(at[ind])
1637                    if sigind in atomsSig:                       
1638                        sigstr += fmt[ind]%(atomsSig[sigind])
1639                    else:
1640                        sigstr += 8*' '
1641            print >>pFile,name
1642            print >>pFile,valstr
1643            print >>pFile,sigstr
1644           
1645    def PrintMomentsAndSig(General,Atoms,atomsSig):
1646        print >>pFile,'\n Magnetic Moments:'    #add magnitude & angle, etc.? TBD
1647        line = '   name      Mx        My        Mz'
1648        cx,ct,cs,cia = General['AtomPtrs']
1649        cmx = cx+4
1650        AtInfo = dict(zip(General['AtomTypes'],General['Lande g']))
1651        print >>pFile,line
1652        print >>pFile,135*'-'
1653        fmt = {0:'%7s',ct:'%7s',cmx:'%10.5f',cmx+1:'%10.5f',cmx+2:'%10.5f'}
1654        noFXsig = {cmx:[10*' ','%10s'],cmx+1:[10*' ','%10s'],cmx+2:[10*' ','%10s']}
1655        for i,at in enumerate(Atoms):
1656            if AtInfo[at[ct]]:
1657                name = fmt[0]%(at[ct-1])+fmt[1]%(at[ct])+':'
1658                valstr = ' values:'
1659                sigstr = ' sig   :'
1660                for ind in range(cmx,cmx+3):
1661                    sigind = str(i)+':'+str(ind)
1662                    valstr += fmt[ind]%(at[ind])                   
1663                    if sigind in atomsSig:
1664                        sigstr += fmt[ind]%(atomsSig[sigind])
1665                    else:
1666                        sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
1667                print >>pFile,name
1668                print >>pFile,valstr
1669                print >>pFile,sigstr
1670           
1671    def PrintWavesAndSig(General,Atoms,wavesSig):
1672        cx,ct,cs,cia = General['AtomPtrs']
1673        print >>pFile,'\n Modulation waves'
1674        names = {'Sfrac':['Fsin','Fcos','Fzero','Fwid'],'Spos':['Xsin','Ysin','Zsin','Xcos','Ycos','Zcos','Tmin','Tmax','Xmax','Ymax','Zmax'],
1675            'Sadp':['U11sin','U22sin','U33sin','U12sin','U13sin','U23sin','U11cos','U22cos',
1676            'U33cos','U12cos','U13cos','U23cos'],'Smag':['MXsin','MYsin','MZsin','MXcos','MYcos','MZcos']}
1677        print >>pFile,135*'-'
1678        for i,at in enumerate(Atoms):
1679            AtomSS = at[-1]['SS1']
1680            waveType = AtomSS['waveType']
1681            for Stype in ['Sfrac','Spos','Sadp','Smag']:
1682                Waves = AtomSS[Stype]
1683                if len(Waves):
1684                    print >>pFile,' atom: %s, site sym: %s, type: %s wave type: %s:'    \
1685                        %(at[ct-1],at[cs],Stype,waveType)
1686                    line = ''
1687                    for iw,wave in enumerate(Waves):
1688                        stiw = ':'+str(i)+':'+str(iw)
1689                        namstr = '  names :'
1690                        valstr = '  values:'
1691                        sigstr = '  esds  :'
1692                        if Stype == 'Spos':
1693                            nt = 6
1694                            ot = 0
1695                            if waveType in ['ZigZag','Block',] and not iw:
1696                                nt = 5
1697                                ot = 6
1698                            for j in range(nt):
1699                                name = names['Spos'][j+ot]
1700                                namstr += '%12s'%(name)
1701                                valstr += '%12.4f'%(wave[0][j])
1702                                if name+stiw in wavesSig:
1703                                    sigstr += '%12.4f'%(wavesSig[name+stiw])
1704                                else:
1705                                    sigstr += 12*' '
1706                        elif Stype == 'Sfrac':
1707                            ot = 0
1708                            if 'Crenel' in waveType and not iw:
1709                                ot = 2
1710                            for j in range(2):
1711                                name = names['Sfrac'][j+ot]
1712                                namstr += '%12s'%(names['Sfrac'][j+ot])
1713                                valstr += '%12.4f'%(wave[0][j])
1714                                if name+stiw in wavesSig:
1715                                    sigstr += '%12.4f'%(wavesSig[name+stiw])
1716                                else:
1717                                    sigstr += 12*' '
1718                        elif Stype == 'Sadp':
1719                            for j in range(12):
1720                                name = names['Sadp'][j]
1721                                namstr += '%10s'%(names['Sadp'][j])
1722                                valstr += '%10.6f'%(wave[0][j])
1723                                if name+stiw in wavesSig:
1724                                    sigstr += '%10.6f'%(wavesSig[name+stiw])
1725                                else:
1726                                    sigstr += 10*' '
1727                        elif Stype == 'Smag':
1728                            for j in range(6):
1729                                name = names['Smag'][j]
1730                                namstr += '%12s'%(names['Smag'][j])
1731                                valstr += '%12.4f'%(wave[0][j])
1732                                if name+stiw in wavesSig:
1733                                    sigstr += '%12.4f'%(wavesSig[name+stiw])
1734                                else:
1735                                    sigstr += 12*' '
1736                               
1737                    print >>pFile,namstr
1738                    print >>pFile,valstr
1739                    print >>pFile,sigstr
1740       
1741               
1742    def PrintRBObjPOAndSig(rbfx,rbsx):
1743        namstr = '  names :'
1744        valstr = '  values:'
1745        sigstr = '  esds  :'
1746        for i,px in enumerate(['Px:','Py:','Pz:']):
1747            name = pfx+rbfx+px+rbsx
1748            namstr += '%12s'%('Pos '+px[1])
1749            valstr += '%12.5f'%(parmDict[name])
1750            if name in sigDict:
1751                sigstr += '%12.5f'%(sigDict[name])
1752            else:
1753                sigstr += 12*' '
1754        for i,po in enumerate(['Oa:','Oi:','Oj:','Ok:']):
1755            name = pfx+rbfx+po+rbsx
1756            namstr += '%12s'%('Ori '+po[1])
1757            valstr += '%12.5f'%(parmDict[name])
1758            if name in sigDict:
1759                sigstr += '%12.5f'%(sigDict[name])
1760            else:
1761                sigstr += 12*' '
1762        print >>pFile,namstr
1763        print >>pFile,valstr
1764        print >>pFile,sigstr
1765       
1766    def PrintRBObjTLSAndSig(rbfx,rbsx,TLS):
1767        namstr = '  names :'
1768        valstr = '  values:'
1769        sigstr = '  esds  :'
1770        if 'N' not in TLS:
1771            print >>pFile,' Thermal motion:'
1772        if 'T' in TLS:
1773            for i,pt in enumerate(['T11:','T22:','T33:','T12:','T13:','T23:']):
1774                name = pfx+rbfx+pt+rbsx
1775                namstr += '%12s'%(pt[:3])
1776                valstr += '%12.5f'%(parmDict[name])
1777                if name in sigDict:
1778                    sigstr += '%12.5f'%(sigDict[name])
1779                else:
1780                    sigstr += 12*' '
1781            print >>pFile,namstr
1782            print >>pFile,valstr
1783            print >>pFile,sigstr
1784        if 'L' in TLS:
1785            namstr = '  names :'
1786            valstr = '  values:'
1787            sigstr = '  esds  :'
1788            for i,pt in enumerate(['L11:','L22:','L33:','L12:','L13:','L23:']):
1789                name = pfx+rbfx+pt+rbsx
1790                namstr += '%12s'%(pt[:3])
1791                valstr += '%12.3f'%(parmDict[name])
1792                if name in sigDict:
1793                    sigstr += '%12.3f'%(sigDict[name])
1794                else:
1795                    sigstr += 12*' '
1796            print >>pFile,namstr
1797            print >>pFile,valstr
1798            print >>pFile,sigstr
1799        if 'S' in TLS:
1800            namstr = '  names :'
1801            valstr = '  values:'
1802            sigstr = '  esds  :'
1803            for i,pt in enumerate(['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']):
1804                name = pfx+rbfx+pt+rbsx
1805                namstr += '%12s'%(pt[:3])
1806                valstr += '%12.4f'%(parmDict[name])
1807                if name in sigDict:
1808                    sigstr += '%12.4f'%(sigDict[name])
1809                else:
1810                    sigstr += 12*' '
1811            print >>pFile,namstr
1812            print >>pFile,valstr
1813            print >>pFile,sigstr
1814        if 'U' in TLS:
1815            name = pfx+rbfx+'U:'+rbsx
1816            namstr = '  names :'+'%12s'%('Uiso')
1817            valstr = '  values:'+'%12.5f'%(parmDict[name])
1818            if name in sigDict:
1819                sigstr = '  esds  :'+'%12.5f'%(sigDict[name])
1820            else:
1821                sigstr = '  esds  :'+12*' '
1822            print >>pFile,namstr
1823            print >>pFile,valstr
1824            print >>pFile,sigstr
1825       
1826    def PrintRBObjTorAndSig(rbsx):
1827        namstr = '  names :'
1828        valstr = '  values:'
1829        sigstr = '  esds  :'
1830        nTors = len(RBObj['Torsions'])
1831        if nTors:
1832            print >>pFile,' Torsions:'
1833            for it in range(nTors):
1834                name = pfx+'RBRTr;'+str(it)+':'+rbsx
1835                namstr += '%12s'%('Tor'+str(it))
1836                valstr += '%12.4f'%(parmDict[name])
1837                if name in sigDict:
1838                    sigstr += '%12.4f'%(sigDict[name])
1839            print >>pFile,namstr
1840            print >>pFile,valstr
1841            print >>pFile,sigstr
1842               
1843    def PrintSHtextureAndSig(textureData,SHtextureSig):
1844        print >>pFile,'\n Spherical harmonics texture: Order:' + str(textureData['Order'])
1845        names = ['omega','chi','phi']
1846        namstr = '  names :'
1847        ptstr =  '  values:'
1848        sigstr = '  esds  :'
1849        for name in names:
1850            namstr += '%12s'%(name)
1851            ptstr += '%12.3f'%(textureData['Sample '+name][1])
1852            if 'Sample '+name in SHtextureSig:
1853                sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
1854            else:
1855                sigstr += 12*' '
1856        print >>pFile,namstr
1857        print >>pFile,ptstr
1858        print >>pFile,sigstr
1859        print >>pFile,'\n Texture coefficients:'
1860        SHcoeff = textureData['SH Coeff'][1]
1861        SHkeys = SHcoeff.keys()
1862        nCoeff = len(SHcoeff)
1863        nBlock = nCoeff/10+1
1864        iBeg = 0
1865        iFin = min(iBeg+10,nCoeff)
1866        for block in range(nBlock):
1867            namstr = '  names :'
1868            ptstr =  '  values:'
1869            sigstr = '  esds  :'
1870            for name in SHkeys[iBeg:iFin]:
1871                namstr += '%12s'%(name)
1872                ptstr += '%12.3f'%(SHcoeff[name])
1873                if name in SHtextureSig:
1874                    sigstr += '%12.3f'%(SHtextureSig[name])
1875                else:
1876                    sigstr += 12*' '
1877            print >>pFile,namstr
1878            print >>pFile,ptstr
1879            print >>pFile,sigstr
1880            iBeg += 10
1881            iFin = min(iBeg+10,nCoeff)
1882           
1883    print >>pFile,'\n Phases:'
1884    for phase in Phases:
1885        print >>pFile,' Result for phase: ',phase
1886        Phase = Phases[phase]
1887        General = Phase['General']
1888        SGData = General['SGData']
1889        Atoms = Phase['Atoms']
1890        if Atoms and not General.get('doPawley'):
1891            cx,ct,cs,cia = General['AtomPtrs']
1892            AtLookup = G2mth.FillAtomLookUp(Atoms,cia+8)
1893        cell = General['Cell']
1894        pId = Phase['pId']
1895        pfx = str(pId)+'::'
1896        if cell[0]:
1897            A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
1898            cellSig = getCellEsd(pfx,SGData,A,covData)  #includes sigVol
1899            print >>pFile,' Reciprocal metric tensor: '
1900            ptfmt = "%15.9f"
1901            names = ['A11','A22','A33','A12','A13','A23']
1902            namstr = '  names :'
1903            ptstr =  '  values:'
1904            sigstr = '  esds  :'
1905            for name,a,siga in zip(names,A,sigA):
1906                namstr += '%15s'%(name)
1907                ptstr += ptfmt%(a)
1908                if siga:
1909                    sigstr += ptfmt%(siga)
1910                else:
1911                    sigstr += 15*' '
1912            print >>pFile,namstr
1913            print >>pFile,ptstr
1914            print >>pFile,sigstr
1915            cell[1:7] = G2lat.A2cell(A)
1916            cell[7] = G2lat.calc_V(A)
1917            print >>pFile,' New unit cell:'
1918            ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
1919            names = ['a','b','c','alpha','beta','gamma','Volume']
1920            namstr = '  names :'
1921            ptstr =  '  values:'
1922            sigstr = '  esds  :'
1923            for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
1924                namstr += '%12s'%(name)
1925                ptstr += fmt%(a)
1926                if siga:
1927                    sigstr += fmt%(siga)
1928                else:
1929                    sigstr += 12*' '
1930            print >>pFile,namstr
1931            print >>pFile,ptstr
1932            print >>pFile,sigstr
1933        ik = 6  #for Pawley stuff below
1934        if General.get('Modulated',False):
1935            ik = 7
1936            Vec,vRef,maxH = General['SuperVec']
1937            if vRef:
1938                print >>pFile,' New modulation vector:'
1939                namstr = '  names :'
1940                ptstr =  '  values:'
1941                sigstr = '  esds  :'
1942                for var in ['mV0','mV1','mV2']:
1943                    namstr += '%12s'%(pfx+var)
1944                    ptstr += '%12.6f'%(parmDict[pfx+var])
1945                    if pfx+var in sigDict:
1946                        sigstr += '%12.6f'%(sigDict[pfx+var])
1947                    else:
1948                        sigstr += 12*' '
1949                print >>pFile,namstr
1950                print >>pFile,ptstr
1951                print >>pFile,sigstr
1952           
1953        General['Mass'] = 0.
1954        if Phase['General'].get('doPawley'):
1955            pawleyRef = Phase['Pawley ref']
1956            for i,refl in enumerate(pawleyRef):
1957                key = pfx+'PWLref:'+str(i)
1958                refl[ik] = parmDict[key]
1959                if key in sigDict:
1960                    refl[ik+1] = sigDict[key]
1961                else:
1962                    refl[ik+1] = 0
1963        else:
1964            VRBIds = RBIds['Vector']
1965            RRBIds = RBIds['Residue']
1966            RBModels = Phase['RBModels']
1967            for irb,RBObj in enumerate(RBModels.get('Vector',[])):
1968                jrb = VRBIds.index(RBObj['RBId'])
1969                rbsx = str(irb)+':'+str(jrb)
1970                print >>pFile,' Vector rigid body parameters:'
1971                PrintRBObjPOAndSig('RBV',rbsx)
1972                PrintRBObjTLSAndSig('RBV',rbsx,RBObj['ThermalMotion'][0])
1973            for irb,RBObj in enumerate(RBModels.get('Residue',[])):
1974                jrb = RRBIds.index(RBObj['RBId'])
1975                rbsx = str(irb)+':'+str(jrb)
1976                print >>pFile,' Residue rigid body parameters:'
1977                PrintRBObjPOAndSig('RBR',rbsx)
1978                PrintRBObjTLSAndSig('RBR',rbsx,RBObj['ThermalMotion'][0])
1979                PrintRBObjTorAndSig(rbsx)
1980            atomsSig = {}
1981            wavesSig = {}
1982            cx,ct,cs,cia = General['AtomPtrs']
1983            for i,at in enumerate(Atoms):
1984                names = {cx:pfx+'Ax:'+str(i),cx+1:pfx+'Ay:'+str(i),cx+2:pfx+'Az:'+str(i),cx+3:pfx+'Afrac:'+str(i),
1985                    cia+1:pfx+'AUiso:'+str(i),cia+2:pfx+'AU11:'+str(i),cia+3:pfx+'AU22:'+str(i),cia+4:pfx+'AU33:'+str(i),
1986                    cia+5:pfx+'AU12:'+str(i),cia+6:pfx+'AU13:'+str(i),cia+7:pfx+'AU23:'+str(i),
1987                    cx+4:pfx+'AMx:'+str(i),cx+5:pfx+'AMy:'+str(i),cx+6:pfx+'AMz:'+str(i)}
1988                for ind in range(cx,cx+4):
1989                    at[ind] = parmDict[names[ind]]
1990                    if ind in range(cx,cx+3):
1991                        name = names[ind].replace('A','dA')
1992                    else:
1993                        name = names[ind]
1994                    if name in sigDict:
1995                        atomsSig[str(i)+':'+str(ind)] = sigDict[name]
1996                if at[cia] == 'I':
1997                    at[cia+1] = parmDict[names[cia+1]]
1998                    if names[cia+1] in sigDict:
1999                        atomsSig['%d:%d'%(i,cia+1)] = sigDict[names[cia+1]]
2000                else:
2001                    for ind in range(cia+2,cia+8):
2002                        at[ind] = parmDict[names[ind]]
2003                        if names[ind] in sigDict:
2004                            atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
2005                if General['Type'] == 'magnetic':
2006                    for ind in range(cx+4,cx+7):
2007                        at[ind] = parmDict[names[ind]]
2008                        if names[ind] in sigDict:
2009                            atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
2010                ind = General['AtomTypes'].index(at[ct])
2011                General['Mass'] += General['AtomMass'][ind]*at[cx+3]*at[cx+5]
2012                if General.get('Modulated',False):
2013                    AtomSS = at[-1]['SS1']
2014                    waveType = AtomSS['waveType']
2015                    for Stype in ['Sfrac','Spos','Sadp','Smag']:
2016                        Waves = AtomSS[Stype]
2017                        for iw,wave in enumerate(Waves):
2018                            stiw = str(i)+':'+str(iw)
2019                            if Stype == 'Spos':
2020                                if waveType in ['ZigZag','Block',] and not iw:
2021                                    names = ['Tmin:'+stiw,'Tmax:'+stiw,'Xmax:'+stiw,'Ymax:'+stiw,'Zmax:'+stiw]
2022                                else:
2023                                    names = ['Xsin:'+stiw,'Ysin:'+stiw,'Zsin:'+stiw,
2024                                        'Xcos:'+stiw,'Ycos:'+stiw,'Zcos:'+stiw]
2025                            elif Stype == 'Sadp':
2026                                names = ['U11sin:'+stiw,'U22sin:'+stiw,'U33sin:'+stiw,
2027                                    'U12sin:'+stiw,'U13sin:'+stiw,'U23sin:'+stiw,
2028                                    'U11cos:'+stiw,'U22cos:'+stiw,'U33cos:'+stiw,
2029                                    'U12cos:'+stiw,'U13cos:'+stiw,'U23cos:'+stiw]
2030                            elif Stype == 'Sfrac':
2031                                if 'Crenel' in waveType and not iw:
2032                                    names = ['Fzero:'+stiw,'Fwid:'+stiw]
2033                                else:
2034                                    names = ['Fsin:'+stiw,'Fcos:'+stiw]
2035                            elif Stype == 'Smag':
2036                                names = ['MXsin:'+stiw,'MYsin:'+stiw,'MZsin:'+stiw,
2037                                    'MXcos:'+stiw,'MYcos:'+stiw,'MZcos:'+stiw]
2038                            for iname,name in enumerate(names):
2039                                AtomSS[Stype][iw][0][iname] = parmDict[pfx+name]
2040                                if pfx+name in sigDict:
2041                                    wavesSig[name] = sigDict[pfx+name]
2042                   
2043            PrintAtomsAndSig(General,Atoms,atomsSig)
2044            if General['Type'] == 'magnetic':
2045                PrintMomentsAndSig(General,Atoms,atomsSig)
2046            if General.get('Modulated',False):
2047                PrintWavesAndSig(General,Atoms,wavesSig)
2048           
2049       
2050        textureData = General['SH Texture']   
2051        if textureData['Order']:
2052            SHtextureSig = {}
2053            for name in ['omega','chi','phi']:
2054                aname = pfx+'SH '+name
2055                textureData['Sample '+name][1] = parmDict[aname]
2056                if aname in sigDict:
2057                    SHtextureSig['Sample '+name] = sigDict[aname]
2058            for name in textureData['SH Coeff'][1]:
2059                aname = pfx+name
2060                textureData['SH Coeff'][1][name] = parmDict[aname]
2061                if aname in sigDict:
2062                    SHtextureSig[name] = sigDict[aname]
2063            PrintSHtextureAndSig(textureData,SHtextureSig)
2064        if phase in RestraintDict and not Phase['General'].get('doPawley'):
2065            PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
2066                textureData,RestraintDict[phase],pFile)
2067                   
2068################################################################################
2069##### Histogram & Phase data
2070################################################################################       
2071                   
2072def GetHistogramPhaseData(Phases,Histograms,Print=True,pFile=None,resetRefList=True):
2073    '''Loads the HAP histogram/phase information into dicts
2074
2075    :param dict Phases: phase information
2076    :param dict Histograms: Histogram information
2077    :param bool Print: prints information as it is read
2078    :param file pFile: file object to print to (the default, None causes printing to the console)
2079    :param bool resetRefList: Should the contents of the Reflection List be initialized
2080      on loading. The default, True, initializes the Reflection List as it is loaded.
2081
2082    :returns: (hapVary,hapDict,controlDict)
2083      * hapVary: list of refined variables
2084      * hapDict: dict with refined variables and their values
2085      * controlDict: dict with fixed parameters
2086    '''
2087   
2088    def PrintSize(hapData):
2089        if hapData[0] in ['isotropic','uniaxial']:
2090            line = '\n Size model    : %9s'%(hapData[0])
2091            line += ' equatorial:'+'%12.5f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
2092            if hapData[0] == 'uniaxial':
2093                line += ' axial:'+'%12.5f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
2094            line += '\n\t LG mixing coeff.: %12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2095            print >>pFile,line
2096        else:
2097            print >>pFile,'\n Size model    : %s'%(hapData[0])+ \
2098                '\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2099            Snames = ['S11','S22','S33','S12','S13','S23']
2100            ptlbls = ' names :'
2101            ptstr =  ' values:'
2102            varstr = ' refine:'
2103            for i,name in enumerate(Snames):
2104                ptlbls += '%12s' % (name)
2105                ptstr += '%12.6f' % (hapData[4][i])
2106                varstr += '%12s' % (str(hapData[5][i]))
2107            print >>pFile,ptlbls
2108            print >>pFile,ptstr
2109            print >>pFile,varstr
2110       
2111    def PrintMuStrain(hapData,SGData):
2112        if hapData[0] in ['isotropic','uniaxial']:
2113            line = '\n Mustrain model: %9s'%(hapData[0])
2114            line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
2115            if hapData[0] == 'uniaxial':
2116                line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
2117            line +='\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2118            print >>pFile,line
2119        else:
2120            print >>pFile,'\n Mustrain model: %s'%(hapData[0])+ \
2121                '\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2122            Snames = G2spc.MustrainNames(SGData)
2123            ptlbls = ' names :'
2124            ptstr =  ' values:'
2125            varstr = ' refine:'
2126            for i,name in enumerate(Snames):
2127                ptlbls += '%12s' % (name)
2128                ptstr += '%12.6f' % (hapData[4][i])
2129                varstr += '%12s' % (str(hapData[5][i]))
2130            print >>pFile,ptlbls
2131            print >>pFile,ptstr
2132            print >>pFile,varstr
2133
2134    def PrintHStrain(hapData,SGData):
2135        print >>pFile,'\n Hydrostatic/elastic strain: '
2136        Hsnames = G2spc.HStrainNames(SGData)
2137        ptlbls = ' names :'
2138        ptstr =  ' values:'
2139        varstr = ' refine:'
2140        for i,name in enumerate(Hsnames):
2141            ptlbls += '%12s' % (name)
2142            ptstr += '%12.4g' % (hapData[0][i])
2143            varstr += '%12s' % (str(hapData[1][i]))
2144        print >>pFile,ptlbls
2145        print >>pFile,ptstr
2146        print >>pFile,varstr
2147
2148    def PrintSHPO(hapData):
2149        print >>pFile,'\n Spherical harmonics preferred orientation: Order:' + \
2150            str(hapData[4])+' Refine? '+str(hapData[2])
2151        ptlbls = ' names :'
2152        ptstr =  ' values:'
2153        for item in hapData[5]:
2154            ptlbls += '%12s'%(item)
2155            ptstr += '%12.3f'%(hapData[5][item]) 
2156        print >>pFile,ptlbls
2157        print >>pFile,ptstr
2158   
2159    def PrintBabinet(hapData):
2160        print >>pFile,'\n Babinet form factor modification: '
2161        ptlbls = ' names :'
2162        ptstr =  ' values:'
2163        varstr = ' refine:'
2164        for name in ['BabA','BabU']:
2165            ptlbls += '%12s' % (name)
2166            ptstr += '%12.6f' % (hapData[name][0])
2167            varstr += '%12s' % (str(hapData[name][1]))
2168        print >>pFile,ptlbls
2169        print >>pFile,ptstr
2170        print >>pFile,varstr
2171       
2172    hapDict = {}
2173    hapVary = []
2174    controlDict = {}
2175    poType = {}
2176    poAxes = {}
2177    spAxes = {}
2178    spType = {}
2179   
2180    for phase in Phases:
2181        HistoPhase = Phases[phase]['Histograms']
2182        SGData = Phases[phase]['General']['SGData']
2183        cell = Phases[phase]['General']['Cell'][1:7]
2184        A = G2lat.cell2A(cell)
2185        if Phases[phase]['General'].get('Modulated',False):
2186            SSGData = Phases[phase]['General']['SSGData']
2187            Vec,x,maxH = Phases[phase]['General']['SuperVec']
2188        pId = Phases[phase]['pId']
2189        histoList = HistoPhase.keys()
2190        histoList.sort()
2191        for histogram in histoList:
2192            try:
2193                Histogram = Histograms[histogram]
2194            except KeyError:                       
2195                #skip if histogram not included e.g. in a sequential refinement
2196                continue
2197            hapData = HistoPhase[histogram]
2198            hId = Histogram['hId']
2199            if 'PWDR' in histogram:
2200                limits = Histogram['Limits'][1]
2201                inst = Histogram['Instrument Parameters'][0]
2202                Zero = inst['Zero'][0]
2203                if 'C' in inst['Type'][1]:
2204                    try:
2205                        wave = inst['Lam'][1]
2206                    except KeyError:
2207                        wave = inst['Lam1'][1]
2208                    dmin = wave/(2.0*sind(limits[1]/2.0))
2209                elif 'T' in inst['Type'][0]:
2210                    dmin = limits[0]/inst['difC'][1]
2211                if Phases[phase]['General']['Type'] == 'magnetic':
2212                    dmin = max(dmin,Phases[phase]['General']['MagDmin'])
2213                pfx = str(pId)+':'+str(hId)+':'
2214                for item in ['Scale','Extinction']:
2215                    hapDict[pfx+item] = hapData[item][0]
2216                    if hapData[item][1]:
2217                        hapVary.append(pfx+item)
2218                names = G2spc.HStrainNames(SGData)
2219                HSvals = []
2220                for i,name in enumerate(names):
2221                    hapDict[pfx+name] = hapData['HStrain'][0][i]
2222                    HSvals.append(hapDict[pfx+name])
2223                    if hapData['HStrain'][1][i]:
2224                        hapVary.append(pfx+name)
2225                DIJS = G2spc.HStrainVals(HSvals,SGData)
2226                controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
2227                if hapData['Pref.Ori.'][0] == 'MD':
2228                    hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
2229                    controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
2230                    if hapData['Pref.Ori.'][2]:
2231                        hapVary.append(pfx+'MD')
2232                else:                           #'SH' spherical harmonics
2233                    controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
2234                    controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
2235                    controlDict[pfx+'SHnames'] = G2lat.GenSHCoeff(SGData['SGLaue'],'0',controlDict[pfx+'SHord'],False)
2236                    controlDict[pfx+'SHhkl'] = []
2237                    try: #patch for old Pref.Ori. items
2238                        controlDict[pfx+'SHtoler'] = 0.1
2239                        if hapData['Pref.Ori.'][6][0] != '':
2240                            controlDict[pfx+'SHhkl'] = [eval(a.replace(' ',',')) for a in hapData['Pref.Ori.'][6]]
2241                        controlDict[pfx+'SHtoler'] = hapData['Pref.Ori.'][7]
2242                    except IndexError:
2243                        pass
2244                    for item in hapData['Pref.Ori.'][5]:
2245                        hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
2246                        if hapData['Pref.Ori.'][2]:
2247                            hapVary.append(pfx+item)
2248                for item in ['Mustrain','Size']:
2249                    controlDict[pfx+item+'Type'] = hapData[item][0]
2250                    hapDict[pfx+item+';mx'] = hapData[item][1][2]
2251                    if hapData[item][2][2]:
2252                        hapVary.append(pfx+item+';mx')
2253                    if hapData[item][0] in ['isotropic','uniaxial']:
2254                        hapDict[pfx+item+';i'] = hapData[item][1][0]
2255                        if hapData[item][2][0]:
2256                            hapVary.append(pfx+item+';i')
2257                        if hapData[item][0] == 'uniaxial':
2258                            controlDict[pfx+item+'Axis'] = hapData[item][3]
2259                            hapDict[pfx+item+';a'] = hapData[item][1][1]
2260                            if hapData[item][2][1]:
2261                                hapVary.append(pfx+item+';a')
2262                    else:       #generalized for mustrain or ellipsoidal for size
2263                        Nterms = len(hapData[item][4])
2264                        if item == 'Mustrain':
2265                            names = G2spc.MustrainNames(SGData)
2266                            pwrs = []
2267                            for name in names:
2268                                h,k,l = name[1:]
2269                                pwrs.append([int(h),int(k),int(l)])
2270                            controlDict[pfx+'MuPwrs'] = pwrs
2271                        for i in range(Nterms):
2272                            sfx = ';'+str(i)
2273                            hapDict[pfx+item+sfx] = hapData[item][4][i]
2274                            if hapData[item][5][i]:
2275                                hapVary.append(pfx+item+sfx)
2276                for bab in ['BabA','BabU']:
2277                    hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2278                    if hapData['Babinet'][bab][1]:
2279                        hapVary.append(pfx+bab)
2280                               
2281                if Print: 
2282                    print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
2283                    print >>pFile,135*'-'
2284                    print >>pFile,' Phase fraction  : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
2285                    print >>pFile,' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
2286                    if hapData['Pref.Ori.'][0] == 'MD':
2287                        Ax = hapData['Pref.Ori.'][3]
2288                        print >>pFile,' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
2289                            ' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
2290                    else: #'SH' for spherical harmonics
2291                        PrintSHPO(hapData['Pref.Ori.'])
2292                        print >>pFile,' Penalty hkl list: '+str(controlDict[pfx+'SHhkl'])+' tolerance: %.2f'%(controlDict[pfx+'SHtoler'])
2293                    PrintSize(hapData['Size'])
2294                    PrintMuStrain(hapData['Mustrain'],SGData)
2295                    PrintHStrain(hapData['HStrain'],SGData)
2296                    if hapData['Babinet']['BabA'][0]:
2297                        PrintBabinet(hapData['Babinet'])                       
2298                if resetRefList:
2299                    refList = []
2300                    Uniq = []
2301                    Phi = []
2302                    useExt = 'magnetic' in Phases[phase]['General']['Type'] and 'N' in inst['Type'][0]
2303                    if Phases[phase]['General'].get('Modulated',False):
2304                        ifSuper = True
2305                        HKLd = np.array(G2lat.GenSSHLaue(dmin,SGData,SSGData,Vec,maxH,A))
2306                        HKLd = G2mth.sortArray(HKLd,4,reverse=True)
2307                        for h,k,l,m,d in HKLd:
2308                            ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2309                            mul *= 2      # for powder overlap of Friedel pairs
2310                            if m or not ext or useExt:
2311                                if 'C' in inst['Type'][0]:
2312                                    pos = G2lat.Dsp2pos(inst,d)
2313                                    if limits[0] < pos < limits[1]:
2314                                        refList.append([h,k,l,m,mul,d, pos,0.0,0.0,0.0,0.0, 0.0,0.0,1.0,1.0,1.0])
2315                                        #... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2316                                        Uniq.append(uniq)
2317                                        Phi.append(phi)
2318                                elif 'T' in inst['Type'][0]:
2319                                    pos = G2lat.Dsp2pos(inst,d)
2320                                    if limits[0] < pos < limits[1]:
2321                                        wave = inst['difC'][1]*d/(252.816*inst['fltPath'][0])
2322                                        refList.append([h,k,l,m,mul,d, pos,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,wave, 1.0,1.0,1.0])
2323                                        # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2324                                        Uniq.append(uniq)
2325                                        Phi.append(phi)
2326                    else:
2327                        ifSuper = False
2328                        HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
2329                        HKLd = G2mth.sortArray(HKLd,3,reverse=True)
2330                        for h,k,l,d in HKLd:
2331                            ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2332                            mul *= 2      # for powder overlap of Friedel pairs
2333                            if ext and not useExt:
2334                                continue
2335                            if 'C' in inst['Type'][0]:
2336                                pos = G2lat.Dsp2pos(inst,d)
2337                                if limits[0] < pos < limits[1]:
2338                                    refList.append([h,k,l,mul,d, pos,0.0,0.0,0.0,0.0, 0.0,0.0,1.0,1.0,1.0])
2339                                    #... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2340                                    Uniq.append(uniq)
2341                                    Phi.append(phi)
2342                            elif 'T' in inst['Type'][0]:
2343                                pos = G2lat.Dsp2pos(inst,d)
2344                                if limits[0] < pos < limits[1]:
2345                                    wave = inst['difC'][1]*d/(252.816*inst['fltPath'][0])
2346                                    refList.append([h,k,l,mul,d, pos,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,wave, 1.0,1.0,1.0])
2347                                    # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2348                                    Uniq.append(uniq)
2349                                    Phi.append(phi)
2350                    Histogram['Reflection Lists'][phase] = {'RefList':np.array(refList),'FF':{},'Type':inst['Type'][0],'Super':ifSuper}
2351            elif 'HKLF' in histogram:
2352                inst = Histogram['Instrument Parameters'][0]
2353                hId = Histogram['hId']
2354                hfx = ':%d:'%(hId)
2355                for item in inst:
2356                    if type(inst) is not list and item != 'Type': continue # skip over non-refined items (such as InstName)
2357                    hapDict[hfx+item] = inst[item][1]
2358                pfx = str(pId)+':'+str(hId)+':'
2359                hapDict[pfx+'Scale'] = hapData['Scale'][0]
2360                if hapData['Scale'][1]:
2361                    hapVary.append(pfx+'Scale')
2362                               
2363                extApprox,extType,extParms = hapData['Extinction']
2364                controlDict[pfx+'EType'] = extType
2365                controlDict[pfx+'EApprox'] = extApprox
2366                if 'C' in inst['Type'][0]:
2367                    controlDict[pfx+'Tbar'] = extParms['Tbar']
2368                    controlDict[pfx+'Cos2TM'] = extParms['Cos2TM']
2369                if 'Primary' in extType:
2370                    Ekey = ['Ep',]
2371                elif 'I & II' in extType:
2372                    Ekey = ['Eg','Es']
2373                elif 'Secondary Type II' == extType:
2374                    Ekey = ['Es',]
2375                elif 'Secondary Type I' == extType:
2376                    Ekey = ['Eg',]
2377                else:   #'None'
2378                    Ekey = []
2379                for eKey in Ekey:
2380                    hapDict[pfx+eKey] = extParms[eKey][0]
2381                    if extParms[eKey][1]:
2382                        hapVary.append(pfx+eKey)
2383                for bab in ['BabA','BabU']:
2384                    hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2385                    if hapData['Babinet'][bab][1]:
2386                        hapVary.append(pfx+bab)
2387                Twins = hapData.get('Twins',[[np.array([[1,0,0],[0,1,0],[0,0,1]]),[1.0,False,0]],])
2388                if len(Twins) == 1:
2389                    hapDict[pfx+'Flack'] = hapData.get('Flack',[0.,False])[0]
2390                    if hapData.get('Flack',[0,False])[1]:
2391                        hapVary.append(pfx+'Flack')
2392                sumTwFr = 0.
2393                controlDict[pfx+'TwinLaw'] = []
2394                controlDict[pfx+'TwinInv'] = []
2395                NTL = 0           
2396                for it,twin in enumerate(Twins):
2397                    if 'bool' in str(type(twin[0])):
2398                        controlDict[pfx+'TwinInv'].append(twin[0])
2399                        controlDict[pfx+'TwinLaw'].append(np.zeros((3,3)))
2400                    else:
2401                        NTL += 1
2402                        controlDict[pfx+'TwinInv'].append(False)
2403                        controlDict[pfx+'TwinLaw'].append(twin[0])
2404                    if it:
2405                        hapDict[pfx+'TwinFr:'+str(it)] = twin[1]
2406                        sumTwFr += twin[1]
2407                    else:
2408                        hapDict[pfx+'TwinFr:0'] = twin[1][0]
2409                        controlDict[pfx+'TwinNMN'] = twin[1][1]
2410                    if Twins[0][1][1]:
2411                        hapVary.append(pfx+'TwinFr:'+str(it))
2412                controlDict[pfx+'NTL'] = NTL
2413                #need to make constraint on TwinFr
2414                controlDict[pfx+'TwinLaw'] = np.array(controlDict[pfx+'TwinLaw'])
2415                if len(Twins) > 1:    #force sum to unity
2416                    hapDict[pfx+'TwinFr:0'] = 1.-sumTwFr
2417                if Print: 
2418                    print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
2419                    print >>pFile,135*'-'
2420                    print >>pFile,' Scale factor     : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
2421                    if extType != 'None':
2422                        print >>pFile,' Extinction  Type: %15s'%(extType),' approx: %10s'%(extApprox)
2423                        text = ' Parameters       :'
2424                        for eKey in Ekey:
2425                            text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
2426                        print >>pFile,text
2427                    if hapData['Babinet']['BabA'][0]:
2428                        PrintBabinet(hapData['Babinet'])
2429                    if not SGData['SGInv'] and len(Twins) == 1:
2430                        print >>pFile,' Flack parameter: %10.3f'%(hapData['Flack'][0]),' Refine?',hapData['Flack'][1]
2431                    if len(Twins) > 1:
2432                        for it,twin in enumerate(Twins):
2433                            if 'bool' in str(type(twin[0])):
2434                                print >>pFile,' Nonmerohedral twin fr.: %5.3f Inv? %s Refine? '%(hapDict[pfx+'TwinFr:'+str(it)],str(controlDict[pfx+'TwinInv'][it])),Twins[0][1][1] 
2435                            else:
2436                                print >>pFile,' Twin law: %s'%(str(twin[0]).replace('\n',',')),' Twin fr.: %5.3f Refine? '%(hapDict[pfx+'TwinFr:'+str(it)]),Twins[0][1][1] 
2437                       
2438                Histogram['Reflection Lists'] = phase       
2439               
2440    return hapVary,hapDict,controlDict
2441   
2442def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,FFtables,Print=True,pFile=None):
2443    'needs a doc string'
2444   
2445    def PrintSizeAndSig(hapData,sizeSig):
2446        line = '\n Size model:     %9s'%(hapData[0])
2447        refine = False
2448        if hapData[0] in ['isotropic','uniaxial']:
2449            line += ' equatorial:%12.5f'%(hapData[1][0])
2450            if sizeSig[0][0]:
2451                line += ', sig:%8.4f'%(sizeSig[0][0])
2452                refine = True
2453            if hapData[0] == 'uniaxial':
2454                line += ' axial:%12.4f'%(hapData[1][1])
2455                if sizeSig[0][1]:
2456                    refine = True
2457                    line += ', sig:%8.4f'%(sizeSig[0][1])
2458            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2459            if sizeSig[0][2]:
2460                refine = True
2461                line += ', sig:%8.4f'%(sizeSig[0][2])
2462            if refine:
2463                print >>pFile,line
2464        else:
2465            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2466            if sizeSig[0][2]:
2467                refine = True
2468                line += ', sig:%8.4f'%(sizeSig[0][2])
2469            Snames = ['S11','S22','S33','S12','S13','S23']
2470            ptlbls = ' name  :'
2471            ptstr =  ' value :'
2472            sigstr = ' sig   :'
2473            for i,name in enumerate(Snames):
2474                ptlbls += '%12s' % (name)
2475                ptstr += '%12.6f' % (hapData[4][i])
2476                if sizeSig[1][i]:
2477                    refine = True
2478                    sigstr += '%12.6f' % (sizeSig[1][i])
2479                else:
2480                    sigstr += 12*' '
2481            if refine:
2482                print >>pFile,line
2483                print >>pFile,ptlbls
2484                print >>pFile,ptstr
2485                print >>pFile,sigstr
2486       
2487    def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
2488        line = '\n Mustrain model: %9s'%(hapData[0])
2489        refine = False
2490        if hapData[0] in ['isotropic','uniaxial']:
2491            line += ' equatorial:%12.1f'%(hapData[1][0])
2492            if mustrainSig[0][0]:
2493                line += ', sig:%8.1f'%(mustrainSig[0][0])
2494                refine = True
2495            if hapData[0] == 'uniaxial':
2496                line += ' axial:%12.1f'%(hapData[1][1])
2497                if mustrainSig[0][1]:
2498                     line += ', sig:%8.1f'%(mustrainSig[0][1])
2499            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2500            if mustrainSig[0][2]:
2501                refine = True
2502                line += ', sig:%8.3f'%(mustrainSig[0][2])
2503            if refine:
2504                print >>pFile,line
2505        else:
2506            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2507            if mustrainSig[0][2]:
2508                refine = True
2509                line += ', sig:%8.3f'%(mustrainSig[0][2])
2510            Snames = G2spc.MustrainNames(SGData)
2511            ptlbls = ' name  :'
2512            ptstr =  ' value :'
2513            sigstr = ' sig   :'
2514            for i,name in enumerate(Snames):
2515                ptlbls += '%12s' % (name)
2516                ptstr += '%12.1f' % (hapData[4][i])
2517                if mustrainSig[1][i]:
2518                    refine = True
2519                    sigstr += '%12.1f' % (mustrainSig[1][i])
2520                else:
2521                    sigstr += 12*' '
2522            if refine:
2523                print >>pFile,line
2524                print >>pFile,ptlbls
2525                print >>pFile,ptstr
2526                print >>pFile,sigstr
2527           
2528    def PrintHStrainAndSig(hapData,strainSig,SGData):
2529        Hsnames = G2spc.HStrainNames(SGData)
2530        ptlbls = ' name  :'
2531        ptstr =  ' value :'
2532        sigstr = ' sig   :'
2533        refine = False
2534        for i,name in enumerate(Hsnames):
2535            ptlbls += '%12s' % (name)
2536            ptstr += '%12.4g' % (hapData[0][i])
2537            if name in strainSig:
2538                refine = True
2539                sigstr += '%12.4g' % (strainSig[name])
2540            else:
2541                sigstr += 12*' '
2542        if refine:
2543            print >>pFile,'\n Hydrostatic/elastic strain: '
2544            print >>pFile,ptlbls
2545            print >>pFile,ptstr
2546            print >>pFile,sigstr
2547       
2548    def PrintSHPOAndSig(pfx,hapData,POsig):
2549        print >>pFile,'\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
2550        ptlbls = ' names :'
2551        ptstr =  ' values:'
2552        sigstr = ' sig   :'
2553        for item in hapData[5]:
2554            ptlbls += '%12s'%(item)
2555            ptstr += '%12.3f'%(hapData[5][item])
2556            if pfx+item in POsig:
2557                sigstr += '%12.3f'%(POsig[pfx+item])
2558            else:
2559                sigstr += 12*' ' 
2560        print >>pFile,ptlbls
2561        print >>pFile,ptstr
2562        print >>pFile,sigstr
2563       
2564    def PrintExtAndSig(pfx,hapData,ScalExtSig):
2565        print >>pFile,'\n Single crystal extinction: Type: ',hapData[0],' Approx: ',hapData[1]
2566        text = ''
2567        for item in hapData[2]:
2568            if pfx+item in ScalExtSig:
2569                text += '       %s: '%(item)
2570                text += '%12.2e'%(hapData[2][item][0])
2571                if pfx+item in ScalExtSig:
2572                    text += ' sig: %12.2e'%(ScalExtSig[pfx+item])
2573        print >>pFile,text       
2574       
2575    def PrintBabinetAndSig(pfx,hapData,BabSig):
2576        print >>pFile,'\n Babinet form factor modification: '
2577        ptlbls = ' names :'
2578        ptstr =  ' values:'
2579        sigstr = ' sig   :'
2580        for item in hapData:
2581            ptlbls += '%12s'%(item)
2582            ptstr += '%12.3f'%(hapData[item][0])
2583            if pfx+item in BabSig:
2584                sigstr += '%12.3f'%(BabSig[pfx+item])
2585            else:
2586                sigstr += 12*' ' 
2587        print >>pFile,ptlbls
2588        print >>pFile,ptstr
2589        print >>pFile,sigstr
2590       
2591    def PrintTwinsAndSig(pfx,twinData,TwinSig):
2592        print >>pFile,'\n Twin Law fractions : '
2593        ptlbls = ' names :'
2594        ptstr =  ' values:'
2595        sigstr = ' sig   :'
2596        for it,item in enumerate(twinData):
2597            ptlbls += '     twin #%d'%(it)
2598            if it:
2599                ptstr += '%12.3f'%(item[1])
2600            else:
2601                ptstr += '%12.3f'%(item[1][0])
2602            if pfx+'TwinFr:'+str(it) in TwinSig:
2603                sigstr += '%12.3f'%(TwinSig[pfx+'TwinFr:'+str(it)])
2604            else:
2605                sigstr += 12*' ' 
2606        print >>pFile,ptlbls
2607        print >>pFile,ptstr
2608        print >>pFile,sigstr
2609       
2610   
2611    PhFrExtPOSig = {}
2612    SizeMuStrSig = {}
2613    ScalExtSig = {}
2614    BabSig = {}
2615    TwinFrSig = {}
2616    wtFrSum = {}
2617    for phase in Phases:
2618        HistoPhase = Phases[phase]['Histograms']
2619        General = Phases[phase]['General']
2620        SGData = General['SGData']
2621        pId = Phases[phase]['pId']
2622        histoList = HistoPhase.keys()
2623        histoList.sort()
2624        for histogram in histoList:
2625            try:
2626                Histogram = Histograms[histogram]
2627            except KeyError:                       
2628                #skip if histogram not included e.g. in a sequential refinement
2629                continue
2630            hapData = HistoPhase[histogram]
2631            hId = Histogram['hId']
2632            pfx = str(pId)+':'+str(hId)+':'
2633            if hId not in wtFrSum:
2634                wtFrSum[hId] = 0.
2635            if 'PWDR' in histogram:
2636                for item in ['Scale','Extinction']:
2637                    hapData[item][0] = parmDict[pfx+item]
2638                    if pfx+item in sigDict:
2639                        PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2640                wtFrSum[hId] += hapData['Scale'][0]*General['Mass']
2641                if hapData['Pref.Ori.'][0] == 'MD':
2642                    hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
2643                    if pfx+'MD' in sigDict:
2644                        PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],})
2645                else:                           #'SH' spherical harmonics
2646                    for item in hapData['Pref.Ori.'][5]:
2647                        hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
2648                        if pfx+item in sigDict:
2649                            PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2650                SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
2651                    pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
2652                    pfx+'HStrain':{}})                 
2653                for item in ['Mustrain','Size']:
2654                    hapData[item][1][2] = parmDict[pfx+item+';mx']
2655#                    hapData[item][1][2] = min(1.,max(0.,hapData[item][1][2]))
2656                    if pfx+item+';mx' in sigDict:
2657                        SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
2658                    if hapData[item][0] in ['isotropic','uniaxial']:                   
2659                        hapData[item][1][0] = parmDict[pfx+item+';i']
2660                        if item == 'Size':
2661                            hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
2662                        if pfx+item+';i' in sigDict: 
2663                            SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
2664                        if hapData[item][0] == 'uniaxial':
2665                            hapData[item][1][1] = parmDict[pfx+item+';a']
2666                            if item == 'Size':
2667                                hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))                       
2668                            if pfx+item+';a' in sigDict:
2669                                SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
2670                    else:       #generalized for mustrain or ellipsoidal for size
2671                        Nterms = len(hapData[item][4])
2672                        for i in range(Nterms):
2673                            sfx = ';'+str(i)
2674                            hapData[item][4][i] = parmDict[pfx+item+sfx]
2675                            if pfx+item+sfx in sigDict:
2676                                SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx]
2677                names = G2spc.HStrainNames(SGData)
2678                for i,name in enumerate(names):
2679                    hapData['HStrain'][0][i] = parmDict[pfx+name]
2680                    if pfx+name in sigDict:
2681                        SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
2682                for name in ['BabA','BabU']:
2683                    hapData['Babinet'][name][0] = parmDict[pfx+name]
2684                    if pfx+name in sigDict:
2685                        BabSig[pfx+name] = sigDict[pfx+name]               
2686               
2687            elif 'HKLF' in histogram:
2688                for item in ['Scale','Flack']:
2689                    if parmDict.get(pfx+item):
2690                        hapData[item][0] = parmDict[pfx+item]
2691                        if pfx+item in sigDict:
2692                            ScalExtSig[pfx+item] = sigDict[pfx+item]
2693                for item in ['Ep','Eg','Es']:
2694                    if parmDict.get(pfx+item):
2695                        hapData['Extinction'][2][item][0] = parmDict[pfx+item]
2696                        if pfx+item in sigDict:
2697                            ScalExtSig[pfx+item] = sigDict[pfx+item]
2698                for item in ['BabA','BabU']:
2699                    hapData['Babinet'][item][0] = parmDict[pfx+item]
2700                    if pfx+item in sigDict:
2701                        BabSig[pfx+item] = sigDict[pfx+item]
2702                if 'Twins' in hapData:
2703                    it = 1
2704                    sumTwFr = 0.
2705                    while True:
2706                        try:
2707                            hapData['Twins'][it][1] = parmDict[pfx+'TwinFr:'+str(it)]
2708                            if pfx+'TwinFr:'+str(it) in sigDict:
2709                                TwinFrSig[pfx+'TwinFr:'+str(it)] = sigDict[pfx+'TwinFr:'+str(it)]
2710                            if it:
2711                                sumTwFr += hapData['Twins'][it][1]
2712                            it += 1
2713                        except KeyError:
2714                            break
2715                    hapData['Twins'][0][1][0] = 1.-sumTwFr
2716
2717    if Print:
2718        for phase in Phases:
2719            HistoPhase = Phases[phase]['Histograms']
2720            General = Phases[phase]['General']
2721            SGData = General['SGData']
2722            pId = Phases[phase]['pId']
2723            histoList = HistoPhase.keys()
2724            histoList.sort()
2725            for histogram in histoList:
2726                try:
2727                    Histogram = Histograms[histogram]
2728                except KeyError:                       
2729                    #skip if histogram not included e.g. in a sequential refinement
2730                    continue
2731                print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
2732                print >>pFile,130*'-'
2733                hapData = HistoPhase[histogram]
2734                hId = Histogram['hId']
2735                Histogram['Residuals'][str(pId)+'::Name'] = phase
2736                pfx = str(pId)+':'+str(hId)+':'
2737                hfx = ':%s:'%(hId)
2738                if 'PWDR' in histogram:
2739                    print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections'   \
2740                        %(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'])
2741                    print >>pFile,' Bragg intensity sum = %.3g'%(Histogram['Residuals'][pfx+'sumInt'])
2742               
2743                    if pfx+'Scale' in PhFrExtPOSig:
2744                        wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum[hId]
2745                        sigwtFr = PhFrExtPOSig[pfx+'Scale']*wtFr/hapData['Scale'][0]
2746                        print >>pFile,' Phase fraction  : %10.5f, sig %10.5f Weight fraction  : %8.5f, sig %10.5f' \
2747                            %(hapData['Scale'][0],PhFrExtPOSig[pfx+'Scale'],wtFr,sigwtFr)
2748                    if pfx+'Extinction' in PhFrExtPOSig:
2749                        print >>pFile,' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction'])
2750                    if hapData['Pref.Ori.'][0] == 'MD':
2751                        if pfx+'MD' in PhFrExtPOSig:
2752                            print >>pFile,' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD'])
2753                    else:
2754                        PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig)
2755                    PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size'])
2756                    PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
2757                    PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
2758                    if len(BabSig):
2759                        PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2760                   
2761                elif 'HKLF' in histogram:
2762                    Inst = Histogram['Instrument Parameters'][0]
2763                    print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections (%d user rejected, %d sp.gp.extinct)'   \
2764                        %(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'],
2765                        Histogram['Residuals'][pfx+'Nrej'],Histogram['Residuals'][pfx+'Next'])
2766                    if FFtables != None and 'N' not in Inst['Type'][0]:
2767                        PrintFprime(FFtables,hfx,pFile)
2768                    print >>pFile,' HKLF histogram weight factor = ','%.3f'%(Histogram['wtFactor'])
2769                    if pfx+'Scale' in ScalExtSig:
2770                        print >>pFile,' Scale factor : %10.4f, sig %10.4f'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale'])
2771                    if hapData['Extinction'][0] != 'None':
2772                        PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig)
2773                    if len(BabSig):
2774                        PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2775                    if pfx+'Flack' in ScalExtSig:
2776                        print >>pFile,' Flack parameter : %10.3f, sig %10.3f'%(hapData['Flack'][0],ScalExtSig[pfx+'Flack'])
2777                    if len(TwinFrSig):
2778                        PrintTwinsAndSig(pfx,hapData['Twins'],TwinFrSig)
2779
2780################################################################################
2781##### Histogram data
2782################################################################################       
2783                   
2784def GetHistogramData(Histograms,Print=True,pFile=None):
2785    'needs a doc string'
2786   
2787    def GetBackgroundParms(hId,Background):
2788        Back = Background[0]
2789        DebyePeaks = Background[1]
2790        bakType,bakFlag = Back[:2]
2791        backVals = Back[3:]
2792        backNames = [':'+str(hId)+':Back;'+str(i) for i in range(len(backVals))]
2793        backDict = dict(zip(backNames,backVals))
2794        backVary = []
2795        if bakFlag:
2796            backVary = backNames
2797        backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
2798        backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
2799        debyeDict = {}
2800        debyeList = []
2801        for i in range(DebyePeaks['nDebye']):
2802            debyeNames = [':'+str(hId)+':DebyeA;'+str(i),':'+str(hId)+':DebyeR;'+str(i),':'+str(hId)+':DebyeU;'+str(i)]
2803            debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
2804            debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
2805        debyeVary = []
2806        for item in debyeList:
2807            if item[1]:
2808                debyeVary.append(item[0])
2809        backDict.update(debyeDict)
2810        backVary += debyeVary
2811        peakDict = {}
2812        peakList = []
2813        for i in range(DebyePeaks['nPeaks']):
2814            peakNames = [':'+str(hId)+':BkPkpos;'+str(i),':'+str(hId)+ \
2815                ':BkPkint;'+str(i),':'+str(hId)+':BkPksig;'+str(i),':'+str(hId)+':BkPkgam;'+str(i)]
2816            peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
2817            peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
2818        peakVary = []
2819        for item in peakList:
2820            if item[1]:
2821                peakVary.append(item[0])
2822        backDict.update(peakDict)
2823        backVary += peakVary
2824        return bakType,backDict,backVary           
2825       
2826    def GetInstParms(hId,Inst):     
2827        dataType = Inst['Type'][0]
2828        instDict = {}
2829        insVary = []
2830        pfx = ':'+str(hId)+':'
2831        insKeys = Inst.keys()
2832        insKeys.sort()
2833        for item in insKeys:
2834            insName = pfx+item
2835            instDict[insName] = Inst[item][1]
2836            if len(Inst[item]) > 2 and Inst[item][2]:
2837                insVary.append(insName)
2838        if 'C' in dataType:
2839            instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
2840        return dataType,instDict,insVary
2841       
2842    def GetSampleParms(hId,Sample):
2843        sampVary = []
2844        hfx = ':'+str(hId)+':'       
2845        sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
2846            hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
2847        for key in ('Temperature','Pressure','FreePrm1','FreePrm2','FreePrm3'):
2848            if key in Sample:
2849                sampDict[hfx+key] = Sample[key]
2850        Type = Sample['Type']
2851        if 'Bragg' in Type:             #Bragg-Brentano
2852            for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2853                sampDict[hfx+item] = Sample[item][0]
2854                if Sample[item][1]:
2855                    sampVary.append(hfx+item)
2856        elif 'Debye' in Type:        #Debye-Scherrer
2857            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2858                sampDict[hfx+item] = Sample[item][0]
2859                if Sample[item][1]:
2860                    sampVary.append(hfx+item)
2861        return Type,sampDict,sampVary
2862       
2863    def PrintBackground(Background):
2864        Back = Background[0]
2865        DebyePeaks = Background[1]
2866        print >>pFile,'\n Background function: ',Back[0],' Refine?',bool(Back[1])
2867        line = ' Coefficients: '
2868        for i,back in enumerate(Back[3:]):
2869            line += '%10.3f'%(back)
2870            if i and not i%10:
2871                line += '\n'+15*' '
2872        print >>pFile,line
2873        if DebyePeaks['nDebye']:
2874            print >>pFile,'\n Debye diffuse scattering coefficients'
2875            parms = ['DebyeA','DebyeR','DebyeU']
2876            line = ' names :  '
2877            for parm in parms:
2878                line += '%8s refine?'%(parm)
2879            print >>pFile,line
2880            for j,term in enumerate(DebyePeaks['debyeTerms']):
2881                line = ' term'+'%2d'%(j)+':'
2882                for i in range(3):
2883                    line += '%10.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
2884                print >>pFile,line
2885        if DebyePeaks['nPeaks']:
2886            print >>pFile,'\n Single peak coefficients'
2887            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
2888            line = ' names :  '
2889            for parm in parms:
2890                line += '%8s refine?'%(parm)
2891            print >>pFile,line
2892            for j,term in enumerate(DebyePeaks['peaksList']):
2893                line = ' peak'+'%2d'%(j)+':'
2894                for i in range(4):
2895                    line += '%12.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
2896                print >>pFile,line
2897       
2898    def PrintInstParms(Inst):
2899        print >>pFile,'\n Instrument Parameters:'
2900        insKeys = Inst.keys()
2901        insKeys.sort()
2902        iBeg = 0
2903        Ok = True
2904        while Ok:
2905            ptlbls = ' name  :'
2906            ptstr =  ' value :'
2907            varstr = ' refine:'
2908            iFin = min(iBeg+9,len(insKeys))
2909            for item in insKeys[iBeg:iFin]:
2910                if item not in ['Type','Source']:
2911                    ptlbls += '%12s' % (item)
2912                    ptstr += '%12.6f' % (Inst[item][1])
2913                    if item in ['Lam1','Lam2','Azimuth','fltPath','2-theta',]:
2914                        varstr += 12*' '
2915                    else:
2916                        varstr += '%12s' % (str(bool(Inst[item][2])))
2917            print >>pFile,ptlbls
2918            print >>pFile,ptstr
2919            print >>pFile,varstr
2920            iBeg = iFin
2921            if iBeg == len(insKeys):
2922                Ok = False
2923            else:
2924                print >>pFile,'\n'
2925       
2926    def PrintSampleParms(Sample):
2927        print >>pFile,'\n Sample Parameters:'
2928        print >>pFile,' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \
2929            (Sample['Omega'],Sample['Chi'],Sample['Phi'])
2930        ptlbls = ' name  :'
2931        ptstr =  ' value :'
2932        varstr = ' refine:'
2933        if 'Bragg' in Sample['Type']:
2934            for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2935                ptlbls += '%14s'%(item)
2936                ptstr += '%14.4f'%(Sample[item][0])
2937                varstr += '%14s'%(str(bool(Sample[item][1])))
2938           
2939        elif 'Debye' in Type:        #Debye-Scherrer
2940            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2941                ptlbls += '%14s'%(item)
2942                ptstr += '%14.4f'%(Sample[item][0])
2943                varstr += '%14s'%(str(bool(Sample[item][1])))
2944
2945        print >>pFile,ptlbls
2946        print >>pFile,ptstr
2947        print >>pFile,varstr
2948       
2949    histDict = {}
2950    histVary = []
2951    controlDict = {}
2952    histoList = Histograms.keys()
2953    histoList.sort()
2954    for histogram in histoList:
2955        if 'PWDR' in histogram:
2956            Histogram = Histograms[histogram]
2957            hId = Histogram['hId']
2958            pfx = ':'+str(hId)+':'
2959            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
2960            controlDict[pfx+'Limits'] = Histogram['Limits'][1]
2961            controlDict[pfx+'Exclude'] = Histogram['Limits'][2:]
2962            for excl in controlDict[pfx+'Exclude']:
2963                Histogram['Data'][0] = ma.masked_inside(Histogram['Data'][0],excl[0],excl[1])
2964            if controlDict[pfx+'Exclude']:
2965                ma.mask_rows(Histogram['Data'])
2966            Background = Histogram['Background']
2967            Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
2968            controlDict[pfx+'bakType'] = Type
2969            histDict.update(bakDict)
2970            histVary += bakVary
2971           
2972            Inst = Histogram['Instrument Parameters'][0]
2973            Type,instDict,insVary = GetInstParms(hId,Inst)
2974            controlDict[pfx+'histType'] = Type
2975            if 'XC' in Type:
2976                if pfx+'Lam1' in instDict:
2977                    controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
2978                else:
2979                    controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']           
2980            histDict.update(instDict)
2981            histVary += insVary
2982           
2983            Sample = Histogram['Sample Parameters']
2984            Type,sampDict,sampVary = GetSampleParms(hId,Sample)
2985            controlDict[pfx+'instType'] = Type
2986            histDict.update(sampDict)
2987            histVary += sampVary
2988           
2989   
2990            if Print: 
2991                print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
2992                print >>pFile,135*'-'
2993                Units = {'C':' deg','T':' msec'}
2994                units = Units[controlDict[pfx+'histType'][2]]
2995                Limits = controlDict[pfx+'Limits']
2996                print >>pFile,' Instrument type: ',Sample['Type']
2997                print >>pFile,' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)
2998                if len(controlDict[pfx+'Exclude']):
2999                    excls = controlDict[pfx+'Exclude']
3000                    for excl in excls:
3001                        print >>pFile,' Excluded region:  %8.2f%s to %8.2f%s'%(excl[0],units,excl[1],units)   
3002                PrintSampleParms(Sample)
3003                PrintInstParms(Inst)
3004                PrintBackground(Background)
3005        elif 'HKLF' in histogram:
3006            Histogram = Histograms[histogram]
3007            hId = Histogram['hId']
3008            pfx = ':'+str(hId)+':'
3009            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
3010            Inst = Histogram['Instrument Parameters'][0]
3011            controlDict[pfx+'histType'] = Inst['Type'][0]
3012            if 'X' in Inst['Type'][0]:
3013                histDict[pfx+'Lam'] = Inst['Lam'][1]
3014                controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']                   
3015    return histVary,histDict,controlDict
3016   
3017def SetHistogramData(parmDict,sigDict,Histograms,FFtables,Print=True,pFile=None):
3018    'needs a doc string'
3019   
3020    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
3021        Back = Background[0]
3022        DebyePeaks = Background[1]
3023        lenBack = len(Back[3:])
3024        backSig = [0 for i in range(lenBack+3*DebyePeaks['nDebye']+4*DebyePeaks['nPeaks'])]
3025        for i in range(lenBack):
3026            Back[3+i] = parmDict[pfx+'Back;'+str(i)]
3027            if pfx+'Back;'+str(i) in sigDict:
3028                backSig[i] = sigDict[pfx+'Back;'+str(i)]
3029        if DebyePeaks['nDebye']:
3030            for i in range(DebyePeaks['nDebye']):
3031                names = [pfx+'DebyeA;'+str(i),pfx+'DebyeR;'+str(i),pfx+'DebyeU;'+str(i)]
3032                for j,name in enumerate(names):
3033                    DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
3034                    if name in sigDict:
3035                        backSig[lenBack+3*i+j] = sigDict[name]           
3036        if DebyePeaks['nPeaks']:
3037            for i in range(DebyePeaks['nPeaks']):
3038                names = [pfx+'BkPkpos;'+str(i),pfx+'BkPkint;'+str(i),
3039                    pfx+'BkPksig;'+str(i),pfx+'BkPkgam;'+str(i)]
3040                for j,name in enumerate(names):
3041                    DebyePeaks['peaksList'][i][2*j] = parmDict[name]
3042                    if name in sigDict:
3043                        backSig[lenBack+3*DebyePeaks['nDebye']+4*i+j] = sigDict[name]
3044        return backSig
3045       
3046    def SetInstParms(pfx,Inst,parmDict,sigDict):
3047        instSig = {}
3048        insKeys = Inst.keys()
3049        insKeys.sort()
3050        for item in insKeys:
3051            insName = pfx+item
3052            Inst[item][1] = parmDict[insName]
3053            if insName in sigDict:
3054                instSig[item] = sigDict[insName]
3055            else:
3056                instSig[item] = 0
3057        return instSig
3058       
3059    def SetSampleParms(pfx,Sample,parmDict,sigDict):
3060        if 'Bragg' in Sample['Type']:             #Bragg-Brentano
3061            sampSig = [0 for i in range(5)]
3062            for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3063                Sample[item][0] = parmDict[pfx+item]
3064                if pfx+item in sigDict:
3065                    sampSig[i] = sigDict[pfx+item]
3066        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
3067            sampSig = [0 for i in range(4)]
3068            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3069                Sample[item][0] = parmDict[pfx+item]
3070                if pfx+item in sigDict:
3071                    sampSig[i] = sigDict[pfx+item]
3072        return sampSig
3073       
3074    def PrintBackgroundSig(Background,backSig):
3075        Back = Background[0]
3076        DebyePeaks = Background[1]
3077        lenBack = len(Back[3:])
3078        valstr = ' value : '
3079        sigstr = ' sig   : '
3080        refine = False
3081        for i,back in enumerate(Back[3:]):
3082            valstr += '%10.4g'%(back)
3083            if Back[1]:
3084                refine = True
3085                sigstr += '%10.4g'%(backSig[i])
3086            else:
3087                sigstr += 10*' '
3088        if refine:
3089            print >>pFile,'\n Background function: ',Back[0]
3090            print >>pFile,valstr
3091            print >>pFile,sigstr
3092        if DebyePeaks['nDebye']:
3093            ifAny = False
3094            ptfmt = "%12.3f"
3095            names =  ' names :'
3096            ptstr =  ' values:'
3097            sigstr = ' esds  :'
3098            for item in sigDict:
3099                if 'Debye' in item:
3100                    ifAny = True
3101                    names += '%12s'%(item)
3102                    ptstr += ptfmt%(parmDict[item])
3103                    sigstr += ptfmt%(sigDict[item])
3104            if ifAny:
3105                print >>pFile,'\n Debye diffuse scattering coefficients'
3106                print >>pFile,names
3107                print >>pFile,ptstr
3108                print >>pFile,sigstr
3109        if DebyePeaks['nPeaks']:
3110            print >>pFile,'\n Single peak coefficients:'
3111            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
3112            line = ' peak no. '
3113            for parm in parms:
3114                line += '%14s%12s'%(parm.center(14),'esd'.center(12))
3115            print >>pFile,line
3116            for ip in range(DebyePeaks['nPeaks']):
3117                ptstr = ' %4d '%(ip)
3118                for parm in parms:
3119                    fmt = '%14.3f'
3120                    efmt = '%12.3f'
3121                    if parm == 'BkPkpos':
3122                        fmt = '%14.4f'
3123                        efmt = '%12.4f'
3124                    name = pfx+parm+';%d'%(ip)
3125                    ptstr += fmt%(parmDict[name])
3126                    if name in sigDict:
3127                        ptstr += efmt%(sigDict[name])
3128                    else:
3129                        ptstr += 12*' '
3130                print >>pFile,ptstr
3131        sumBk = np.array(Histogram['sumBk'])
3132        print >>pFile,' Background sums: empirical %.3g, Debye %.3g, peaks %.3g, Total %.3g'    \
3133            %(sumBk[0],sumBk[1],sumBk[2],np.sum(sumBk))
3134       
3135    def PrintInstParmsSig(Inst,instSig):
3136        refine = False
3137        insKeys = instSig.keys()
3138        insKeys.sort()
3139        iBeg = 0
3140        Ok = True
3141        while Ok:
3142            ptlbls = ' names :'
3143            ptstr =  ' value :'
3144            sigstr = ' sig   :'
3145            iFin = min(iBeg+9,len(insKeys))
3146            for name in insKeys[iBeg:iFin]:
3147                if name not in  ['Type','Lam1','Lam2','Azimuth','Source','fltPath']:
3148                    ptlbls += '%12s' % (name)
3149                    ptstr += '%12.6f' % (Inst[name][1])
3150                    if instSig[name]:
3151                        refine = True
3152                        sigstr += '%12.6f' % (instSig[name])
3153                    else:
3154                        sigstr += 12*' '
3155            if refine:
3156                print >>pFile,'\n Instrument Parameters:'
3157                print >>pFile,ptlbls
3158                print >>pFile,ptstr
3159                print >>pFile,sigstr
3160            iBeg = iFin
3161            if iBeg == len(insKeys):
3162                Ok = False
3163       
3164    def PrintSampleParmsSig(Sample,sampleSig):
3165        ptlbls = ' names :'
3166        ptstr =  ' values:'
3167        sigstr = ' sig   :'
3168        refine = False
3169        if 'Bragg' in Sample['Type']:
3170            for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3171                ptlbls += '%14s'%(item)
3172                ptstr += '%14.4f'%(Sample[item][0])
3173                if sampleSig[i]:
3174                    refine = True
3175                    sigstr += '%14.4f'%(sampleSig[i])
3176                else:
3177                    sigstr += 14*' '
3178           
3179        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
3180            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3181                ptlbls += '%14s'%(item)
3182                ptstr += '%14.4f'%(Sample[item][0])
3183                if sampleSig[i]:
3184                    refine = True
3185                    sigstr += '%14.4f'%(sampleSig[i])
3186                else:
3187                    sigstr += 14*' '
3188
3189        if refine:
3190            print >>pFile,'\n Sample Parameters:'
3191            print >>pFile,ptlbls
3192            print >>pFile,ptstr
3193            print >>pFile,sigstr
3194       
3195    histoList = Histograms.keys()
3196    histoList.sort()
3197    for histogram in histoList:
3198        if 'PWDR' in histogram:
3199            Histogram = Histograms[histogram]
3200            hId = Histogram['hId']
3201            pfx = ':'+str(hId)+':'
3202            Background = Histogram['Background']
3203            backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
3204           
3205            Inst = Histogram['Instrument Parameters'][0]
3206            instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
3207       
3208            Sample = Histogram['Sample Parameters']
3209            sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
3210
3211            print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
3212            print >>pFile,135*'-'
3213            print >>pFile,' PWDR histogram weight factor = '+'%.3f'%(Histogram['wtFactor'])
3214            print >>pFile,' Final refinement wR = %.2f%% on %d observations in this histogram'%(Histogram['Residuals']['wR'],Histogram['Residuals']['Nobs'])
3215            print >>pFile,' Other residuals: R = %.2f%%, Rb = %.2f%%, wRb = %.2f%% wRmin = %.2f%%'% \
3216                (Histogram['Residuals']['R'],Histogram['Residuals']['Rb'],Histogram['Residuals']['wRb'],Histogram['Residuals']['wRmin'])
3217            if Print:
3218                print >>pFile,' Instrument type: ',Sample['Type']
3219                if FFtables != None and 'N' not in Inst['Type'][0]:
3220                    PrintFprime(FFtables,pfx,pFile)
3221                PrintSampleParmsSig(Sample,sampSig)
3222                PrintInstParmsSig(Inst,instSig)
3223                PrintBackgroundSig(Background,backSig)
3224               
Note: See TracBrowser for help on using the repository browser.