source: trunk/GSASIIstrIO.py @ 2481

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

force Transform to delete nonmagnetic atoms if phase made magnetic & add 'mag' to new phase name
fix TOF cosine background function
magnetic structure refinement works (in numeric mode only)
magnetic structure factor calculation correct

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 145.1 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2016-09-30 15:33:12 +0000 (Fri, 30 Sep 2016) $
4# $Author: vondreele $
5# $Revision: 2481 $
6# $URL: trunk/GSASIIstrIO.py $
7# $Id: GSASIIstrIO.py 2481 2016-09-30 15:33:12Z 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: 2481 $")
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 = G2spc.SytSym(at[cx:cx+3],SGData)
1212                    if SpnFlp[Nop] > 0.:    #black use p
1213                        mId,mCoef = G2spc.GetCSpqinel(SytSym)[0]
1214                    else:                   #red use q
1215                        mId,mCoef = G2spc.GetCSpqinel(SytSym)[1]
1216                    names = [pfx+'AMx:'+str(i),pfx+'AMy:'+str(i),pfx+'AMz:'+str(i)]
1217                    equivs = [[],[],[]]
1218                    for j in range(3):
1219                        if mId[j] > 0:
1220                            phaseVary.append(names[j])
1221                            equivs[mId[j]-1].append([names[j],mCoef[j]])
1222                    for equiv in equivs:
1223                        if len(equiv) > 1:
1224                            name = equiv[0][0]
1225                            coef = equiv[0][1]
1226                            for eqv in equiv[1:]:
1227                                eqv[1] /= coef
1228                                G2mv.StoreEquivalence(name,(eqv,))
1229                if General.get('Modulated',False):
1230                    AtomSS = at[-1]['SS1']
1231                    waveType = AtomSS['waveType']
1232                    phaseDict[pfx+'waveType:'+str(i)] = waveType
1233                    for Stype in ['Sfrac','Spos','Sadp','Smag']:
1234                        Waves = AtomSS[Stype]
1235                        nx = 0
1236                        for iw,wave in enumerate(Waves):
1237                            if not iw:
1238                                if waveType in ['ZigZag','Block']:
1239                                    nx = 1
1240                                CSI = G2spc.GetSSfxuinel(waveType,1,at[cx:cx+3],SGData,SSGData)
1241                            else:
1242                                CSI = G2spc.GetSSfxuinel('Fourier',iw+1-nx,at[cx:cx+3],SGData,SSGData)
1243                            uId,uCoef = CSI[Stype]
1244                            stiw = str(i)+':'+str(iw)
1245                            if Stype == 'Spos':
1246                                if waveType in ['ZigZag','Block',] and not iw:
1247                                    names = [pfx+'Tmin:'+stiw,pfx+'Tmax:'+stiw,pfx+'Xmax:'+stiw,pfx+'Ymax:'+stiw,pfx+'Zmax:'+stiw]
1248                                    equivs = [[],[], [],[],[]]
1249                                else:
1250                                    names = [pfx+'Xsin:'+stiw,pfx+'Ysin:'+stiw,pfx+'Zsin:'+stiw,
1251                                        pfx+'Xcos:'+stiw,pfx+'Ycos:'+stiw,pfx+'Zcos:'+stiw]
1252                                    equivs = [[],[],[], [],[],[]]
1253                            elif Stype == 'Sadp':
1254                                names = [pfx+'U11sin:'+stiw,pfx+'U22sin:'+stiw,pfx+'U33sin:'+stiw,
1255                                    pfx+'U12sin:'+stiw,pfx+'U13sin:'+stiw,pfx+'U23sin:'+stiw,
1256                                    pfx+'U11cos:'+stiw,pfx+'U22cos:'+stiw,pfx+'U33cos:'+stiw,
1257                                    pfx+'U12cos:'+stiw,pfx+'U13cos:'+stiw,pfx+'U23cos:'+stiw]
1258                                equivs = [[],[],[],[],[],[], [],[],[],[],[],[]]
1259                            elif Stype == 'Sfrac':
1260                                equivs = [[],[]]
1261                                if 'Crenel' in waveType and not iw:
1262                                    names = [pfx+'Fzero:'+stiw,pfx+'Fwid:'+stiw]
1263                                else:
1264                                    names = [pfx+'Fsin:'+stiw,pfx+'Fcos:'+stiw]
1265                            elif Stype == 'Smag':
1266                                equivs = [[],[],[], [],[],[]]
1267                                names = [pfx+'MXsin:'+stiw,pfx+'MYsin:'+stiw,pfx+'MZsin:'+stiw,
1268                                    pfx+'MXcos:'+stiw,pfx+'MYcos:'+stiw,pfx+'MZcos:'+stiw]
1269                            phaseDict.update(dict(zip(names,wave[0])))
1270                            if wave[1]: #what do we do here for multiple terms in modulation constraints?
1271                                for j in range(len(equivs)):
1272                                    if uId[j][0] > 0:                               
1273                                        phaseVary.append(names[j])
1274                                        equivs[uId[j][0]-1].append([names[j],uCoef[j][0]])
1275                                for equiv in equivs:
1276                                    if len(equiv) > 1:
1277                                        name = equiv[0][0]
1278                                        coef = equiv[0][1]
1279                                        for eqv in equiv[1:]:
1280                                            eqv[1] /= coef
1281                                        G2mv.StoreEquivalence(name,equiv[1:])
1282                            maxSSwave[pfx][Stype] = max(maxSSwave[pfx][Stype],iw+1)
1283            textureData = General['SH Texture']
1284            if textureData['Order'] and not seqRef:
1285                phaseDict[pfx+'SHorder'] = textureData['Order']
1286                phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']]
1287                for item in ['omega','chi','phi']:
1288                    phaseDict[pfx+'SH '+item] = textureData['Sample '+item][1]
1289                    if textureData['Sample '+item][0]:
1290                        phaseVary.append(pfx+'SH '+item)
1291                for item in textureData['SH Coeff'][1]:
1292                    phaseDict[pfx+item] = textureData['SH Coeff'][1][item]
1293                    if textureData['SH Coeff'][0]:
1294                        phaseVary.append(pfx+item)
1295               
1296            if Print:
1297                print >>pFile,'\n Phase name: ',General['Name']
1298                print >>pFile,135*'-'
1299                PrintFFtable(FFtable)
1300                PrintBLtable(BLtable)
1301                if General['Type'] == 'magnetic':
1302                    PrintMFtable(MFtable)
1303                print >>pFile,''
1304                #how do we print magnetic symmetry table? TBD
1305                if len(SSGtext):    #if superstructure
1306                    for line in SSGtext: print >>pFile,line
1307                    if len(SSGtable):
1308                        for item in SSGtable:
1309                            line = ' %s '%(item)
1310                            print >>pFile,line   
1311                    else:
1312                        print >>pFile,' ( 1)    %s'%(SSGtable[0])
1313                else:
1314                    for line in SGtext: print >>pFile,line
1315                    if len(SGtable):
1316                        for item in SGtable:
1317                            line = ' %s '%(item)
1318                            print >>pFile,line   
1319                    else:
1320                        print >>pFile,' ( 1)    %s'%(SGtable[0])
1321                PrintRBObjects(resRBData,vecRBData)
1322                PrintAtoms(General,Atoms)
1323                if General['Type'] == 'magnetic':
1324                    PrintMoments(General,Atoms)
1325                if General.get('Modulated',False):
1326                    PrintWaves(General,Atoms)
1327                print >>pFile,'\n Unit cell: a = %.5f'%(cell[1]),' b = %.5f'%(cell[2]),' c = %.5f'%(cell[3]), \
1328                    ' alpha = %.3f'%(cell[4]),' beta = %.3f'%(cell[5]),' gamma = %.3f'%(cell[6]), \
1329                    ' volume = %.3f'%(cell[7]),' Refine?',cell[0]
1330                if len(SSGtext):    #if superstructure
1331                    print >>pFile,'\n Modulation vector: mV0 = %.4f'%(Vec[0]),' mV1 = %.4f'%(Vec[1]),   \
1332                        ' mV2 = %.4f'%(Vec[2]),' max mod. index = %d'%(maxH),' Refine?',vRef
1333                if not seqRef:
1334                    PrintTexture(textureData)
1335                if name in RestraintDict:
1336                    PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
1337                        textureData,RestraintDict[name],pFile)
1338                   
1339        elif PawleyRef:
1340            if Print:
1341                print >>pFile,'\n Phase name: ',General['Name']
1342                print >>pFile,135*'-'
1343                print >>pFile,''
1344                if len(SSGtext):    #if superstructure
1345                    for line in SSGtext: print >>pFile,line
1346                    if len(SSGtable):
1347                        for item in SSGtable:
1348                            line = ' %s '%(item)
1349                            print >>pFile,line   
1350                    else:
1351                        print >>pFile,' ( 1)    %s'%(SSGtable[0])
1352                else:
1353                    for line in SGtext: print >>pFile,line
1354                    if len(SGtable):
1355                        for item in SGtable:
1356                            line = ' %s '%(item)
1357                            print >>pFile,line   
1358                    else:
1359                        print >>pFile,' ( 1)    %s'%(SGtable[0])
1360                print >>pFile,'\n Unit cell: a = %.5f'%(cell[1]),' b = %.5f'%(cell[2]),' c = %.5f'%(cell[3]), \
1361                    ' alpha = %.3f'%(cell[4]),' beta = %.3f'%(cell[5]),' gamma = %.3f'%(cell[6]), \
1362                    ' volume = %.3f'%(cell[7]),' Refine?',cell[0]
1363                if len(SSGtext):    #if superstructure
1364                    print >>pFile,'\n Modulation vector: mV0 = %.4f'%(Vec[0]),' mV1 = %.4f'%(Vec[1]),   \
1365                        ' mV2 = %.4f'%(Vec[2]),' max mod. index = %d'%(maxH),' Refine?',vRef
1366            pawleyVary = []
1367            for i,refl in enumerate(PawleyRef):
1368                phaseDict[pfx+'PWLref:'+str(i)] = refl[6+im]
1369                if im:
1370                    pawleyLookup[pfx+'%d,%d,%d,%d'%(refl[0],refl[1],refl[2],refl[3])] = i
1371                else:
1372                    pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
1373                if refl[5+im]:
1374                    pawleyVary.append(pfx+'PWLref:'+str(i))
1375            GetPawleyConstr(SGData['SGLaue'],PawleyRef,im,pawleyVary)      #does G2mv.StoreEquivalence
1376            phaseVary += pawleyVary
1377               
1378    return Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,maxSSwave
1379   
1380def cellFill(pfx,SGData,parmDict,sigDict): 
1381    '''Returns the filled-out reciprocal cell (A) terms and their uncertainties
1382    from the parameter and sig dictionaries.
1383
1384    :param str pfx: parameter prefix ("n::", where n is a phase number)
1385    :param dict SGdata: a symmetry object
1386    :param dict parmDict: a dictionary of parameters
1387    :param dict sigDict:  a dictionary of uncertainties on parameters
1388
1389    :returns: A,sigA where each is a list of six terms with the A terms
1390    '''
1391    if SGData['SGLaue'] in ['-1',]:
1392        A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1393            parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
1394    elif SGData['SGLaue'] in ['2/m',]:
1395        if SGData['SGUniq'] == 'a':
1396            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1397                0,0,parmDict[pfx+'A5']]
1398        elif SGData['SGUniq'] == 'b':
1399            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1400                0,parmDict[pfx+'A4'],0]
1401        else:
1402            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1403                parmDict[pfx+'A3'],0,0]
1404    elif SGData['SGLaue'] in ['mmm',]:
1405        A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
1406    elif SGData['SGLaue'] in ['4/m','4/mmm']:
1407        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
1408    elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1409        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
1410            parmDict[pfx+'A0'],0,0]
1411    elif SGData['SGLaue'] in ['3R', '3mR']:
1412        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
1413            parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
1414    elif SGData['SGLaue'] in ['m3m','m3']:
1415        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
1416
1417    try:
1418        if SGData['SGLaue'] in ['-1',]:
1419            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1420                sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
1421        elif SGData['SGLaue'] in ['2/m',]:
1422            if SGData['SGUniq'] == 'a':
1423                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1424                    0,0,sigDict[pfx+'A5']]
1425            elif SGData['SGUniq'] == 'b':
1426                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1427                    0,sigDict[pfx+'A4'],0]
1428            else:
1429                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1430                    sigDict[pfx+'A3'],0,0]
1431        elif SGData['SGLaue'] in ['mmm',]:
1432            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
1433        elif SGData['SGLaue'] in ['4/m','4/mmm']:
1434            sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1435        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1436            sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1437        elif SGData['SGLaue'] in ['3R', '3mR']:
1438            sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
1439        elif SGData['SGLaue'] in ['m3m','m3']:
1440            sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
1441    except KeyError:
1442        sigA = [0,0,0,0,0,0]
1443    return A,sigA
1444       
1445def PrintRestraints(cell,SGData,AtPtrs,Atoms,AtLookup,textureData,phaseRest,pFile):
1446    'needs a doc string'
1447    if phaseRest:
1448        Amat = G2lat.cell2AB(cell)[0]
1449        cx,ct,cs = AtPtrs[:3]
1450        names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
1451            ['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
1452            ['ChemComp','Sites'],['Texture','HKLs']]
1453        for name,rest in names:
1454            itemRest = phaseRest[name]
1455            if itemRest[rest] and itemRest['Use']:
1456                print >>pFile,'\n  %s %10.3f Use: %s'%(name+' restraint weight factor',itemRest['wtFactor'],str(itemRest['Use']))
1457                if name in ['Bond','Angle','Plane','Chiral']:
1458                    print >>pFile,'     calc       obs      sig   delt/sig  atoms(symOp): '
1459                    for indx,ops,obs,esd in itemRest[rest]:
1460                        try:
1461                            AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1462                            AtName = ''
1463                            for i,Aname in enumerate(AtNames):
1464                                AtName += Aname
1465                                if ops[i] == '1':
1466                                    AtName += '-'
1467                                else:
1468                                    AtName += '+('+ops[i]+')-'
1469                            XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1470                            XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1471                            if name == 'Bond':
1472                                calc = G2mth.getRestDist(XYZ,Amat)
1473                            elif name == 'Angle':
1474                                calc = G2mth.getRestAngle(XYZ,Amat)
1475                            elif name == 'Plane':
1476                                calc = G2mth.getRestPlane(XYZ,Amat)
1477                            elif name == 'Chiral':
1478                                calc = G2mth.getRestChiral(XYZ,Amat)
1479                            print >>pFile,' %9.3f %9.3f %8.3f %8.3f   %s'%(calc,obs,esd,(obs-calc)/esd,AtName[:-1])
1480                        except KeyError:
1481                            del itemRest[rest]
1482                elif name in ['Torsion','Rama']:
1483                    print >>pFile,'  atoms(symOp)  calc  obs  sig  delt/sig  torsions: '
1484                    coeffDict = itemRest['Coeff']
1485                    for indx,ops,cofName,esd in itemRest[rest]:
1486                        AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1487                        AtName = ''
1488                        for i,Aname in enumerate(AtNames):
1489                            AtName += Aname+'+('+ops[i]+')-'
1490                        XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1491                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1492                        if name == 'Torsion':
1493                            tor = G2mth.getRestTorsion(XYZ,Amat)
1494                            restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
1495                            print >>pFile,' %8.3f %8.3f %.3f %8.3f %8.3f %s'%(calc,obs,esd,(obs-calc)/esd,tor,AtName[:-1])
1496                        else:
1497                            phi,psi = G2mth.getRestRama(XYZ,Amat)
1498                            restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])                               
1499                            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])
1500                elif name == 'ChemComp':
1501                    print >>pFile,'     atoms   mul*frac  factor     prod'
1502                    for indx,factors,obs,esd in itemRest[rest]:
1503                        try:
1504                            atoms = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1505                            mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs+1))
1506                            frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs-1))
1507                            mulfrac = mul*frac
1508                            calcs = mul*frac*factors
1509                            for iatm,[atom,mf,fr,clc] in enumerate(zip(atoms,mulfrac,factors,calcs)):
1510                                print >>pFile,' %10s %8.3f %8.3f %8.3f'%(atom,mf,fr,clc)
1511                            print >>pFile,' Sum:                   calc: %8.3f obs: %8.3f esd: %8.3f'%(np.sum(calcs),obs,esd)
1512                        except KeyError:
1513                            del itemRest[rest]
1514                elif name == 'Texture' and textureData['Order']:
1515                    Start = False
1516                    SHCoef = textureData['SH Coeff'][1]
1517                    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1518                    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
1519                    print '    HKL  grid  neg esd  sum neg  num neg use unit?  unit esd  '
1520                    for hkl,grid,esd1,ifesd2,esd2 in itemRest[rest]:
1521                        PH = np.array(hkl)
1522                        phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
1523                        ODFln = G2lat.Flnh(Start,SHCoef,phi,beta,SGData)
1524                        R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
1525                        Z = ma.masked_greater(Z,0.0)
1526                        num = ma.count(Z)
1527                        sum = 0
1528                        if num:
1529                            sum = np.sum(Z)
1530                        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)
1531       
1532def getCellEsd(pfx,SGData,A,covData):
1533    'needs a doc string'
1534    dpr = 180./np.pi
1535    rVsq = G2lat.calc_rVsq(A)
1536    G,g = G2lat.A2Gmat(A)       #get recip. & real metric tensors
1537    cell = np.array(G2lat.Gmat2cell(g))   #real cell
1538    cellst = np.array(G2lat.Gmat2cell(G)) #recip. cell
1539    scos = cosd(cellst[3:6])
1540    ssin = sind(cellst[3:6])
1541    scot = scos/ssin
1542    rcos = cosd(cell[3:6])
1543    rsin = sind(cell[3:6])
1544    rcot = rcos/rsin
1545    RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
1546    varyList = covData['varyList']
1547    covMatrix = covData['covMatrix']
1548    vcov = G2mth.getVCov(RMnames,varyList,covMatrix)
1549    Ax = np.array(A)
1550    Ax[3:] /= 2.
1551    drVdA = np.array([
1552        Ax[1]*Ax[2]-Ax[5]**2,
1553        Ax[0]*Ax[2]-Ax[4]**2,
1554        Ax[0]*Ax[1]-Ax[3]**2,
1555        Ax[4]*Ax[5]-Ax[2]*Ax[3],
1556        Ax[3]*Ax[5]-Ax[1]*Ax[4],
1557        Ax[3]*Ax[4]-Ax[0]*Ax[5]])
1558    srcvlsq = np.inner(drVdA,np.inner(drVdA,vcov))
1559    Vol = 1/np.sqrt(rVsq)
1560    sigVol = Vol**3*np.sqrt(srcvlsq)/2.         #ok - checks with GSAS
1561    R123 = Ax[0]*Ax[1]*Ax[2]
1562    dsasdg = np.zeros((3,6))
1563    dadg = np.zeros((6,6))
1564    for i0 in range(3):         #0  1   2
1565        i1 = (i0+1)%3           #1  2   0
1566        i2 = (i1+1)%3           #2  0   1
1567        i3 = 5-i2               #3  5   4
1568        i4 = 5-i1               #4  3   5
1569        i5 = 5-i0               #5  4   3
1570        dsasdg[i0][i1] = 0.5*scot[i0]*scos[i0]/Ax[i1]
1571        dsasdg[i0][i2] = 0.5*scot[i0]*scos[i0]/Ax[i2]
1572        dsasdg[i0][i5] = -scot[i0]/np.sqrt(Ax[i1]*Ax[i2])
1573        denmsq = Ax[i0]*(R123-Ax[i1]*Ax[i4]**2-Ax[i2]*Ax[i3]**2+(Ax[i4]*Ax[i3])**2)
1574        denom = np.sqrt(denmsq)
1575        dadg[i5][i0] = -Ax[i5]/denom-rcos[i0]/denmsq*(R123-0.5*Ax[i1]*Ax[i4]**2-0.5*Ax[i2]*Ax[i3]**2)
1576        dadg[i5][i1] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i2]-Ax[i0]*Ax[i4]**2)
1577        dadg[i5][i2] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i1]-Ax[i0]*Ax[i3]**2)
1578        dadg[i5][i3] = Ax[i4]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i2]*Ax[i3]-Ax[i3]*Ax[i4]**2)
1579        dadg[i5][i4] = Ax[i3]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i1]*Ax[i4]-Ax[i3]**2*Ax[i4])
1580        dadg[i5][i5] = -Ax[i0]/denom
1581    for i0 in range(3):
1582        i1 = (i0+1)%3
1583        i2 = (i1+1)%3
1584        i3 = 5-i2
1585        for ij in range(6):
1586            dadg[i0][ij] = cell[i0]*(rcot[i2]*dadg[i3][ij]/rsin[i2]-dsasdg[i1][ij]/ssin[i1])
1587            if ij == i0:
1588                dadg[i0][ij] = dadg[i0][ij]-0.5*cell[i0]/Ax[i0]
1589            dadg[i3][ij] = -dadg[i3][ij]*rsin[2-i0]*dpr
1590    sigMat = np.inner(dadg,np.inner(dadg,vcov))
1591    var = np.diag(sigMat)
1592    CS = np.where(var>0.,np.sqrt(var),0.)
1593    return [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol]
1594   
1595def SetPhaseData(parmDict,sigDict,Phases,RBIds,covData,RestraintDict=None,pFile=None):
1596    'needs a doc string'
1597   
1598    def PrintAtomsAndSig(General,Atoms,atomsSig):
1599        print >>pFile,'\n Atoms:'
1600        line = '   name      x         y         z      frac   Uiso     U11     U22     U33     U12     U13     U23'
1601        if General['Type'] == 'macromolecular':
1602            line = ' res no residue chain '+line
1603        cx,ct,cs,cia = General['AtomPtrs']
1604        print >>pFile,line
1605        print >>pFile,135*'-'
1606        fmt = {0:'%7s',ct:'%7s',cx:'%10.5f',cx+1:'%10.5f',cx+2:'%10.5f',cx+3:'%8.3f',cia+1:'%8.5f',
1607            cia+2:'%8.5f',cia+3:'%8.5f',cia+4:'%8.5f',cia+5:'%8.5f',cia+6:'%8.5f',cia+7:'%8.5f'}
1608        noFXsig = {cx:[10*' ','%10s'],cx+1:[10*' ','%10s'],cx+2:[10*' ','%10s'],cx+3:[8*' ','%8s']}
1609        for atyp in General['AtomTypes']:       #zero composition
1610            General['NoAtoms'][atyp] = 0.
1611        for i,at in enumerate(Atoms):
1612            General['NoAtoms'][at[ct]] += at[cx+3]*float(at[cx+5])     #new composition
1613            if General['Type'] == 'macromolecular':
1614                name = ' %s %s %s %s:'%(at[0],at[1],at[2],at[3])
1615                valstr = ' values:          '
1616                sigstr = ' sig   :          '
1617            else:
1618                name = fmt[0]%(at[ct-1])+fmt[1]%(at[ct])+':'
1619                valstr = ' values:'
1620                sigstr = ' sig   :'
1621            for ind in range(cx,cx+4):
1622                sigind = str(i)+':'+str(ind)
1623                valstr += fmt[ind]%(at[ind])                   
1624                if sigind in atomsSig:
1625                    sigstr += fmt[ind]%(atomsSig[sigind])
1626                else:
1627                    sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
1628            if at[cia] == 'I':
1629                valstr += fmt[cia+1]%(at[cia+1])
1630                if '%d:%d'%(i,cia+1) in atomsSig:
1631                    sigstr += fmt[cia+1]%(atomsSig['%d:%d'%(i,cia+1)])
1632                else:
1633                    sigstr += 8*' '
1634            else:
1635                valstr += 8*' '
1636                sigstr += 8*' '
1637                for ind in range(cia+2,cia+8):
1638                    sigind = str(i)+':'+str(ind)
1639                    valstr += fmt[ind]%(at[ind])
1640                    if sigind in atomsSig:                       
1641                        sigstr += fmt[ind]%(atomsSig[sigind])
1642                    else:
1643                        sigstr += 8*' '
1644            print >>pFile,name
1645            print >>pFile,valstr
1646            print >>pFile,sigstr
1647           
1648    def PrintMomentsAndSig(General,Atoms,atomsSig):
1649        print >>pFile,'\n Magnetic Moments:'    #add magnitude & angle, etc.? TBD
1650        line = '   name      Mx        My        Mz'
1651        cx,ct,cs,cia = General['AtomPtrs']
1652        cmx = cx+4
1653        AtInfo = dict(zip(General['AtomTypes'],General['Lande g']))
1654        print >>pFile,line
1655        print >>pFile,135*'-'
1656        fmt = {0:'%7s',ct:'%7s',cmx:'%10.5f',cmx+1:'%10.5f',cmx+2:'%10.5f'}
1657        noFXsig = {cmx:[10*' ','%10s'],cmx+1:[10*' ','%10s'],cmx+2:[10*' ','%10s']}
1658        for i,at in enumerate(Atoms):
1659            if AtInfo[at[ct]]:
1660                name = fmt[0]%(at[ct-1])+fmt[1]%(at[ct])+':'
1661                valstr = ' values:'
1662                sigstr = ' sig   :'
1663                for ind in range(cmx,cmx+3):
1664                    sigind = str(i)+':'+str(ind)
1665                    valstr += fmt[ind]%(at[ind])                   
1666                    if sigind in atomsSig:
1667                        sigstr += fmt[ind]%(atomsSig[sigind])
1668                    else:
1669                        sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
1670                print >>pFile,name
1671                print >>pFile,valstr
1672                print >>pFile,sigstr
1673           
1674    def PrintWavesAndSig(General,Atoms,wavesSig):
1675        cx,ct,cs,cia = General['AtomPtrs']
1676        print >>pFile,'\n Modulation waves'
1677        names = {'Sfrac':['Fsin','Fcos','Fzero','Fwid'],'Spos':['Xsin','Ysin','Zsin','Xcos','Ycos','Zcos','Tmin','Tmax','Xmax','Ymax','Zmax'],
1678            'Sadp':['U11sin','U22sin','U33sin','U12sin','U13sin','U23sin','U11cos','U22cos',
1679            'U33cos','U12cos','U13cos','U23cos'],'Smag':['MXsin','MYsin','MZsin','MXcos','MYcos','MZcos']}
1680        print >>pFile,135*'-'
1681        for i,at in enumerate(Atoms):
1682            AtomSS = at[-1]['SS1']
1683            waveType = AtomSS['waveType']
1684            for Stype in ['Sfrac','Spos','Sadp','Smag']:
1685                Waves = AtomSS[Stype]
1686                if len(Waves):
1687                    print >>pFile,' atom: %s, site sym: %s, type: %s wave type: %s:'    \
1688                        %(at[ct-1],at[cs],Stype,waveType)
1689                    line = ''
1690                    for iw,wave in enumerate(Waves):
1691                        stiw = ':'+str(i)+':'+str(iw)
1692                        namstr = '  names :'
1693                        valstr = '  values:'
1694                        sigstr = '  esds  :'
1695                        if Stype == 'Spos':
1696                            nt = 6
1697                            ot = 0
1698                            if waveType in ['ZigZag','Block',] and not iw:
1699                                nt = 5
1700                                ot = 6
1701                            for j in range(nt):
1702                                name = names['Spos'][j+ot]
1703                                namstr += '%12s'%(name)
1704                                valstr += '%12.4f'%(wave[0][j])
1705                                if name+stiw in wavesSig:
1706                                    sigstr += '%12.4f'%(wavesSig[name+stiw])
1707                                else:
1708                                    sigstr += 12*' '
1709                        elif Stype == 'Sfrac':
1710                            ot = 0
1711                            if 'Crenel' in waveType and not iw:
1712                                ot = 2
1713                            for j in range(2):
1714                                name = names['Sfrac'][j+ot]
1715                                namstr += '%12s'%(names['Sfrac'][j+ot])
1716                                valstr += '%12.4f'%(wave[0][j])
1717                                if name+stiw in wavesSig:
1718                                    sigstr += '%12.4f'%(wavesSig[name+stiw])
1719                                else:
1720                                    sigstr += 12*' '
1721                        elif Stype == 'Sadp':
1722                            for j in range(12):
1723                                name = names['Sadp'][j]
1724                                namstr += '%10s'%(names['Sadp'][j])
1725                                valstr += '%10.6f'%(wave[0][j])
1726                                if name+stiw in wavesSig:
1727                                    sigstr += '%10.6f'%(wavesSig[name+stiw])
1728                                else:
1729                                    sigstr += 10*' '
1730                        elif Stype == 'Smag':
1731                            for j in range(6):
1732                                name = names['Smag'][j]
1733                                namstr += '%12s'%(names['Smag'][j])
1734                                valstr += '%12.4f'%(wave[0][j])
1735                                if name+stiw in wavesSig:
1736                                    sigstr += '%12.4f'%(wavesSig[name+stiw])
1737                                else:
1738                                    sigstr += 12*' '
1739                               
1740                    print >>pFile,namstr
1741                    print >>pFile,valstr
1742                    print >>pFile,sigstr
1743       
1744               
1745    def PrintRBObjPOAndSig(rbfx,rbsx):
1746        namstr = '  names :'
1747        valstr = '  values:'
1748        sigstr = '  esds  :'
1749        for i,px in enumerate(['Px:','Py:','Pz:']):
1750            name = pfx+rbfx+px+rbsx
1751            namstr += '%12s'%('Pos '+px[1])
1752            valstr += '%12.5f'%(parmDict[name])
1753            if name in sigDict:
1754                sigstr += '%12.5f'%(sigDict[name])
1755            else:
1756                sigstr += 12*' '
1757        for i,po in enumerate(['Oa:','Oi:','Oj:','Ok:']):
1758            name = pfx+rbfx+po+rbsx
1759            namstr += '%12s'%('Ori '+po[1])
1760            valstr += '%12.5f'%(parmDict[name])
1761            if name in sigDict:
1762                sigstr += '%12.5f'%(sigDict[name])
1763            else:
1764                sigstr += 12*' '
1765        print >>pFile,namstr
1766        print >>pFile,valstr
1767        print >>pFile,sigstr
1768       
1769    def PrintRBObjTLSAndSig(rbfx,rbsx,TLS):
1770        namstr = '  names :'
1771        valstr = '  values:'
1772        sigstr = '  esds  :'
1773        if 'N' not in TLS:
1774            print >>pFile,' Thermal motion:'
1775        if 'T' in TLS:
1776            for i,pt in enumerate(['T11:','T22:','T33:','T12:','T13:','T23:']):
1777                name = pfx+rbfx+pt+rbsx
1778                namstr += '%12s'%(pt[:3])
1779                valstr += '%12.5f'%(parmDict[name])
1780                if name in sigDict:
1781                    sigstr += '%12.5f'%(sigDict[name])
1782                else:
1783                    sigstr += 12*' '
1784            print >>pFile,namstr
1785            print >>pFile,valstr
1786            print >>pFile,sigstr
1787        if 'L' in TLS:
1788            namstr = '  names :'
1789            valstr = '  values:'
1790            sigstr = '  esds  :'
1791            for i,pt in enumerate(['L11:','L22:','L33:','L12:','L13:','L23:']):
1792                name = pfx+rbfx+pt+rbsx
1793                namstr += '%12s'%(pt[:3])
1794                valstr += '%12.3f'%(parmDict[name])
1795                if name in sigDict:
1796                    sigstr += '%12.3f'%(sigDict[name])
1797                else:
1798                    sigstr += 12*' '
1799            print >>pFile,namstr
1800            print >>pFile,valstr
1801            print >>pFile,sigstr
1802        if 'S' in TLS:
1803            namstr = '  names :'
1804            valstr = '  values:'
1805            sigstr = '  esds  :'
1806            for i,pt in enumerate(['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']):
1807                name = pfx+rbfx+pt+rbsx
1808                namstr += '%12s'%(pt[:3])
1809                valstr += '%12.4f'%(parmDict[name])
1810                if name in sigDict:
1811                    sigstr += '%12.4f'%(sigDict[name])
1812                else:
1813                    sigstr += 12*' '
1814            print >>pFile,namstr
1815            print >>pFile,valstr
1816            print >>pFile,sigstr
1817        if 'U' in TLS:
1818            name = pfx+rbfx+'U:'+rbsx
1819            namstr = '  names :'+'%12s'%('Uiso')
1820            valstr = '  values:'+'%12.5f'%(parmDict[name])
1821            if name in sigDict:
1822                sigstr = '  esds  :'+'%12.5f'%(sigDict[name])
1823            else:
1824                sigstr = '  esds  :'+12*' '
1825            print >>pFile,namstr
1826            print >>pFile,valstr
1827            print >>pFile,sigstr
1828       
1829    def PrintRBObjTorAndSig(rbsx):
1830        namstr = '  names :'
1831        valstr = '  values:'
1832        sigstr = '  esds  :'
1833        nTors = len(RBObj['Torsions'])
1834        if nTors:
1835            print >>pFile,' Torsions:'
1836            for it in range(nTors):
1837                name = pfx+'RBRTr;'+str(it)+':'+rbsx
1838                namstr += '%12s'%('Tor'+str(it))
1839                valstr += '%12.4f'%(parmDict[name])
1840                if name in sigDict:
1841                    sigstr += '%12.4f'%(sigDict[name])
1842            print >>pFile,namstr
1843            print >>pFile,valstr
1844            print >>pFile,sigstr
1845               
1846    def PrintSHtextureAndSig(textureData,SHtextureSig):
1847        print >>pFile,'\n Spherical harmonics texture: Order:' + str(textureData['Order'])
1848        names = ['omega','chi','phi']
1849        namstr = '  names :'
1850        ptstr =  '  values:'
1851        sigstr = '  esds  :'
1852        for name in names:
1853            namstr += '%12s'%(name)
1854            ptstr += '%12.3f'%(textureData['Sample '+name][1])
1855            if 'Sample '+name in SHtextureSig:
1856                sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
1857            else:
1858                sigstr += 12*' '
1859        print >>pFile,namstr
1860        print >>pFile,ptstr
1861        print >>pFile,sigstr
1862        print >>pFile,'\n Texture coefficients:'
1863        SHcoeff = textureData['SH Coeff'][1]
1864        SHkeys = SHcoeff.keys()
1865        nCoeff = len(SHcoeff)
1866        nBlock = nCoeff/10+1
1867        iBeg = 0
1868        iFin = min(iBeg+10,nCoeff)
1869        for block in range(nBlock):
1870            namstr = '  names :'
1871            ptstr =  '  values:'
1872            sigstr = '  esds  :'
1873            for name in SHkeys[iBeg:iFin]:
1874                namstr += '%12s'%(name)
1875                ptstr += '%12.3f'%(SHcoeff[name])
1876                if name in SHtextureSig:
1877                    sigstr += '%12.3f'%(SHtextureSig[name])
1878                else:
1879                    sigstr += 12*' '
1880            print >>pFile,namstr
1881            print >>pFile,ptstr
1882            print >>pFile,sigstr
1883            iBeg += 10
1884            iFin = min(iBeg+10,nCoeff)
1885           
1886    print >>pFile,'\n Phases:'
1887    for phase in Phases:
1888        print >>pFile,' Result for phase: ',phase
1889        Phase = Phases[phase]
1890        General = Phase['General']
1891        SGData = General['SGData']
1892        Atoms = Phase['Atoms']
1893        if Atoms and not General.get('doPawley'):
1894            cx,ct,cs,cia = General['AtomPtrs']
1895            AtLookup = G2mth.FillAtomLookUp(Atoms,cia+8)
1896        cell = General['Cell']
1897        pId = Phase['pId']
1898        pfx = str(pId)+'::'
1899        if cell[0]:
1900            A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
1901            cellSig = getCellEsd(pfx,SGData,A,covData)  #includes sigVol
1902            print >>pFile,' Reciprocal metric tensor: '
1903            ptfmt = "%15.9f"
1904            names = ['A11','A22','A33','A12','A13','A23']
1905            namstr = '  names :'
1906            ptstr =  '  values:'
1907            sigstr = '  esds  :'
1908            for name,a,siga in zip(names,A,sigA):
1909                namstr += '%15s'%(name)
1910                ptstr += ptfmt%(a)
1911                if siga:
1912                    sigstr += ptfmt%(siga)
1913                else:
1914                    sigstr += 15*' '
1915            print >>pFile,namstr
1916            print >>pFile,ptstr
1917            print >>pFile,sigstr
1918            cell[1:7] = G2lat.A2cell(A)
1919            cell[7] = G2lat.calc_V(A)
1920            print >>pFile,' New unit cell:'
1921            ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
1922            names = ['a','b','c','alpha','beta','gamma','Volume']
1923            namstr = '  names :'
1924            ptstr =  '  values:'
1925            sigstr = '  esds  :'
1926            for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
1927                namstr += '%12s'%(name)
1928                ptstr += fmt%(a)
1929                if siga:
1930                    sigstr += fmt%(siga)
1931                else:
1932                    sigstr += 12*' '
1933            print >>pFile,namstr
1934            print >>pFile,ptstr
1935            print >>pFile,sigstr
1936        ik = 6  #for Pawley stuff below
1937        if General.get('Modulated',False):
1938            ik = 7
1939            Vec,vRef,maxH = General['SuperVec']
1940            if vRef:
1941                print >>pFile,' New modulation vector:'
1942                namstr = '  names :'
1943                ptstr =  '  values:'
1944                sigstr = '  esds  :'
1945                for var in ['mV0','mV1','mV2']:
1946                    namstr += '%12s'%(pfx+var)
1947                    ptstr += '%12.6f'%(parmDict[pfx+var])
1948                    if pfx+var in sigDict:
1949                        sigstr += '%12.6f'%(sigDict[pfx+var])
1950                    else:
1951                        sigstr += 12*' '
1952                print >>pFile,namstr
1953                print >>pFile,ptstr
1954                print >>pFile,sigstr
1955           
1956        General['Mass'] = 0.
1957        if Phase['General'].get('doPawley'):
1958            pawleyRef = Phase['Pawley ref']
1959            for i,refl in enumerate(pawleyRef):
1960                key = pfx+'PWLref:'+str(i)
1961                refl[ik] = parmDict[key]
1962                if key in sigDict:
1963                    refl[ik+1] = sigDict[key]
1964                else:
1965                    refl[ik+1] = 0
1966        else:
1967            VRBIds = RBIds['Vector']
1968            RRBIds = RBIds['Residue']
1969            RBModels = Phase['RBModels']
1970            for irb,RBObj in enumerate(RBModels.get('Vector',[])):
1971                jrb = VRBIds.index(RBObj['RBId'])
1972                rbsx = str(irb)+':'+str(jrb)
1973                print >>pFile,' Vector rigid body parameters:'
1974                PrintRBObjPOAndSig('RBV',rbsx)
1975                PrintRBObjTLSAndSig('RBV',rbsx,RBObj['ThermalMotion'][0])
1976            for irb,RBObj in enumerate(RBModels.get('Residue',[])):
1977                jrb = RRBIds.index(RBObj['RBId'])
1978                rbsx = str(irb)+':'+str(jrb)
1979                print >>pFile,' Residue rigid body parameters:'
1980                PrintRBObjPOAndSig('RBR',rbsx)
1981                PrintRBObjTLSAndSig('RBR',rbsx,RBObj['ThermalMotion'][0])
1982                PrintRBObjTorAndSig(rbsx)
1983            atomsSig = {}
1984            wavesSig = {}
1985            cx,ct,cs,cia = General['AtomPtrs']
1986            for i,at in enumerate(Atoms):
1987                names = {cx:pfx+'Ax:'+str(i),cx+1:pfx+'Ay:'+str(i),cx+2:pfx+'Az:'+str(i),cx+3:pfx+'Afrac:'+str(i),
1988                    cia+1:pfx+'AUiso:'+str(i),cia+2:pfx+'AU11:'+str(i),cia+3:pfx+'AU22:'+str(i),cia+4:pfx+'AU33:'+str(i),
1989                    cia+5:pfx+'AU12:'+str(i),cia+6:pfx+'AU13:'+str(i),cia+7:pfx+'AU23:'+str(i),
1990                    cx+4:pfx+'AMx:'+str(i),cx+5:pfx+'AMy:'+str(i),cx+6:pfx+'AMz:'+str(i)}
1991                for ind in range(cx,cx+4):
1992                    at[ind] = parmDict[names[ind]]
1993                    if ind in range(cx,cx+3):
1994                        name = names[ind].replace('A','dA')
1995                    else:
1996                        name = names[ind]
1997                    if name in sigDict:
1998                        atomsSig[str(i)+':'+str(ind)] = sigDict[name]
1999                if at[cia] == 'I':
2000                    at[cia+1] = parmDict[names[cia+1]]
2001                    if names[cia+1] in sigDict:
2002                        atomsSig['%d:%d'%(i,cia+1)] = sigDict[names[cia+1]]
2003                else:
2004                    for ind in range(cia+2,cia+8):
2005                        at[ind] = parmDict[names[ind]]
2006                        if names[ind] in sigDict:
2007                            atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
2008                if General['Type'] == 'magnetic':
2009                    for ind in range(cx+4,cx+7):
2010                        at[ind] = parmDict[names[ind]]
2011                        if names[ind] in sigDict:
2012                            atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
2013                ind = General['AtomTypes'].index(at[ct])
2014                General['Mass'] += General['AtomMass'][ind]*at[cx+3]*at[cx+5]
2015                if General.get('Modulated',False):
2016                    AtomSS = at[-1]['SS1']
2017                    waveType = AtomSS['waveType']
2018                    for Stype in ['Sfrac','Spos','Sadp','Smag']:
2019                        Waves = AtomSS[Stype]
2020                        for iw,wave in enumerate(Waves):
2021                            stiw = str(i)+':'+str(iw)
2022                            if Stype == 'Spos':
2023                                if waveType in ['ZigZag','Block',] and not iw:
2024                                    names = ['Tmin:'+stiw,'Tmax:'+stiw,'Xmax:'+stiw,'Ymax:'+stiw,'Zmax:'+stiw]
2025                                else:
2026                                    names = ['Xsin:'+stiw,'Ysin:'+stiw,'Zsin:'+stiw,
2027                                        'Xcos:'+stiw,'Ycos:'+stiw,'Zcos:'+stiw]
2028                            elif Stype == 'Sadp':
2029                                names = ['U11sin:'+stiw,'U22sin:'+stiw,'U33sin:'+stiw,
2030                                    'U12sin:'+stiw,'U13sin:'+stiw,'U23sin:'+stiw,
2031                                    'U11cos:'+stiw,'U22cos:'+stiw,'U33cos:'+stiw,
2032                                    'U12cos:'+stiw,'U13cos:'+stiw,'U23cos:'+stiw]
2033                            elif Stype == 'Sfrac':
2034                                if 'Crenel' in waveType and not iw:
2035                                    names = ['Fzero:'+stiw,'Fwid:'+stiw]
2036                                else:
2037                                    names = ['Fsin:'+stiw,'Fcos:'+stiw]
2038                            elif Stype == 'Smag':
2039                                names = ['MXsin:'+stiw,'MYsin:'+stiw,'MZsin:'+stiw,
2040                                    'MXcos:'+stiw,'MYcos:'+stiw,'MZcos:'+stiw]
2041                            for iname,name in enumerate(names):
2042                                AtomSS[Stype][iw][0][iname] = parmDict[pfx+name]
2043                                if pfx+name in sigDict:
2044                                    wavesSig[name] = sigDict[pfx+name]
2045                   
2046            PrintAtomsAndSig(General,Atoms,atomsSig)
2047            if General['Type'] == 'magnetic':
2048                PrintMomentsAndSig(General,Atoms,atomsSig)
2049            if General.get('Modulated',False):
2050                PrintWavesAndSig(General,Atoms,wavesSig)
2051           
2052       
2053        textureData = General['SH Texture']   
2054        if textureData['Order']:
2055            SHtextureSig = {}
2056            for name in ['omega','chi','phi']:
2057                aname = pfx+'SH '+name
2058                textureData['Sample '+name][1] = parmDict[aname]
2059                if aname in sigDict:
2060                    SHtextureSig['Sample '+name] = sigDict[aname]
2061            for name in textureData['SH Coeff'][1]:
2062                aname = pfx+name
2063                textureData['SH Coeff'][1][name] = parmDict[aname]
2064                if aname in sigDict:
2065                    SHtextureSig[name] = sigDict[aname]
2066            PrintSHtextureAndSig(textureData,SHtextureSig)
2067        if phase in RestraintDict and not Phase['General'].get('doPawley'):
2068            PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
2069                textureData,RestraintDict[phase],pFile)
2070                   
2071################################################################################
2072##### Histogram & Phase data
2073################################################################################       
2074                   
2075def GetHistogramPhaseData(Phases,Histograms,Print=True,pFile=None,resetRefList=True):
2076    '''Loads the HAP histogram/phase information into dicts
2077
2078    :param dict Phases: phase information
2079    :param dict Histograms: Histogram information
2080    :param bool Print: prints information as it is read
2081    :param file pFile: file object to print to (the default, None causes printing to the console)
2082    :param bool resetRefList: Should the contents of the Reflection List be initialized
2083      on loading. The default, True, initializes the Reflection List as it is loaded.
2084
2085    :returns: (hapVary,hapDict,controlDict)
2086      * hapVary: list of refined variables
2087      * hapDict: dict with refined variables and their values
2088      * controlDict: dict with fixed parameters
2089    '''
2090   
2091    def PrintSize(hapData):
2092        if hapData[0] in ['isotropic','uniaxial']:
2093            line = '\n Size model    : %9s'%(hapData[0])
2094            line += ' equatorial:'+'%12.5f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
2095            if hapData[0] == 'uniaxial':
2096                line += ' axial:'+'%12.5f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
2097            line += '\n\t LG mixing coeff.: %12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2098            print >>pFile,line
2099        else:
2100            print >>pFile,'\n Size model    : %s'%(hapData[0])+ \
2101                '\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2102            Snames = ['S11','S22','S33','S12','S13','S23']
2103            ptlbls = ' names :'
2104            ptstr =  ' values:'
2105            varstr = ' refine:'
2106            for i,name in enumerate(Snames):
2107                ptlbls += '%12s' % (name)
2108                ptstr += '%12.6f' % (hapData[4][i])
2109                varstr += '%12s' % (str(hapData[5][i]))
2110            print >>pFile,ptlbls
2111            print >>pFile,ptstr
2112            print >>pFile,varstr
2113       
2114    def PrintMuStrain(hapData,SGData):
2115        if hapData[0] in ['isotropic','uniaxial']:
2116            line = '\n Mustrain model: %9s'%(hapData[0])
2117            line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
2118            if hapData[0] == 'uniaxial':
2119                line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
2120            line +='\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2121            print >>pFile,line
2122        else:
2123            print >>pFile,'\n Mustrain model: %s'%(hapData[0])+ \
2124                '\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2125            Snames = G2spc.MustrainNames(SGData)
2126            ptlbls = ' names :'
2127            ptstr =  ' values:'
2128            varstr = ' refine:'
2129            for i,name in enumerate(Snames):
2130                ptlbls += '%12s' % (name)
2131                ptstr += '%12.6f' % (hapData[4][i])
2132                varstr += '%12s' % (str(hapData[5][i]))
2133            print >>pFile,ptlbls
2134            print >>pFile,ptstr
2135            print >>pFile,varstr
2136
2137    def PrintHStrain(hapData,SGData):
2138        print >>pFile,'\n Hydrostatic/elastic strain: '
2139        Hsnames = G2spc.HStrainNames(SGData)
2140        ptlbls = ' names :'
2141        ptstr =  ' values:'
2142        varstr = ' refine:'
2143        for i,name in enumerate(Hsnames):
2144            ptlbls += '%12s' % (name)
2145            ptstr += '%12.4g' % (hapData[0][i])
2146            varstr += '%12s' % (str(hapData[1][i]))
2147        print >>pFile,ptlbls
2148        print >>pFile,ptstr
2149        print >>pFile,varstr
2150
2151    def PrintSHPO(hapData):
2152        print >>pFile,'\n Spherical harmonics preferred orientation: Order:' + \
2153            str(hapData[4])+' Refine? '+str(hapData[2])
2154        ptlbls = ' names :'
2155        ptstr =  ' values:'
2156        for item in hapData[5]:
2157            ptlbls += '%12s'%(item)
2158            ptstr += '%12.3f'%(hapData[5][item]) 
2159        print >>pFile,ptlbls
2160        print >>pFile,ptstr
2161   
2162    def PrintBabinet(hapData):
2163        print >>pFile,'\n Babinet form factor modification: '
2164        ptlbls = ' names :'
2165        ptstr =  ' values:'
2166        varstr = ' refine:'
2167        for name in ['BabA','BabU']:
2168            ptlbls += '%12s' % (name)
2169            ptstr += '%12.6f' % (hapData[name][0])
2170            varstr += '%12s' % (str(hapData[name][1]))
2171        print >>pFile,ptlbls
2172        print >>pFile,ptstr
2173        print >>pFile,varstr
2174       
2175    hapDict = {}
2176    hapVary = []
2177    controlDict = {}
2178    poType = {}
2179    poAxes = {}
2180    spAxes = {}
2181    spType = {}
2182   
2183    for phase in Phases:
2184        HistoPhase = Phases[phase]['Histograms']
2185        SGData = Phases[phase]['General']['SGData']
2186        cell = Phases[phase]['General']['Cell'][1:7]
2187        A = G2lat.cell2A(cell)
2188        if Phases[phase]['General'].get('Modulated',False):
2189            SSGData = Phases[phase]['General']['SSGData']
2190            Vec,x,maxH = Phases[phase]['General']['SuperVec']
2191        pId = Phases[phase]['pId']
2192        histoList = HistoPhase.keys()
2193        histoList.sort()
2194        for histogram in histoList:
2195            try:
2196                Histogram = Histograms[histogram]
2197            except KeyError:                       
2198                #skip if histogram not included e.g. in a sequential refinement
2199                continue
2200            hapData = HistoPhase[histogram]
2201            hId = Histogram['hId']
2202            if 'PWDR' in histogram:
2203                limits = Histogram['Limits'][1]
2204                inst = Histogram['Instrument Parameters'][0]
2205                Zero = inst['Zero'][0]
2206                if 'C' in inst['Type'][1]:
2207                    try:
2208                        wave = inst['Lam'][1]
2209                    except KeyError:
2210                        wave = inst['Lam1'][1]
2211                    dmin = wave/(2.0*sind(limits[1]/2.0))
2212                elif 'T' in inst['Type'][0]:
2213                    dmin = limits[0]/inst['difC'][1]
2214                if Phases[phase]['General']['Type'] == 'magnetic':
2215                    dmin = max(dmin,Phases[phase]['General']['MagDmin'])
2216                pfx = str(pId)+':'+str(hId)+':'
2217                for item in ['Scale','Extinction']:
2218                    hapDict[pfx+item] = hapData[item][0]
2219                    if hapData[item][1]:
2220                        hapVary.append(pfx+item)
2221                names = G2spc.HStrainNames(SGData)
2222                HSvals = []
2223                for i,name in enumerate(names):
2224                    hapDict[pfx+name] = hapData['HStrain'][0][i]
2225                    HSvals.append(hapDict[pfx+name])
2226                    if hapData['HStrain'][1][i]:
2227                        hapVary.append(pfx+name)
2228                DIJS = G2spc.HStrainVals(HSvals,SGData)
2229                controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
2230                if hapData['Pref.Ori.'][0] == 'MD':
2231                    hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
2232                    controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
2233                    if hapData['Pref.Ori.'][2]:
2234                        hapVary.append(pfx+'MD')
2235                else:                           #'SH' spherical harmonics
2236                    controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
2237                    controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
2238                    controlDict[pfx+'SHnames'] = G2lat.GenSHCoeff(SGData['SGLaue'],'0',controlDict[pfx+'SHord'],False)
2239                    controlDict[pfx+'SHhkl'] = []
2240                    try: #patch for old Pref.Ori. items
2241                        controlDict[pfx+'SHtoler'] = 0.1
2242                        if hapData['Pref.Ori.'][6][0] != '':
2243                            controlDict[pfx+'SHhkl'] = [eval(a.replace(' ',',')) for a in hapData['Pref.Ori.'][6]]
2244                        controlDict[pfx+'SHtoler'] = hapData['Pref.Ori.'][7]
2245                    except IndexError:
2246                        pass
2247                    for item in hapData['Pref.Ori.'][5]:
2248                        hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
2249                        if hapData['Pref.Ori.'][2]:
2250                            hapVary.append(pfx+item)
2251                for item in ['Mustrain','Size']:
2252                    controlDict[pfx+item+'Type'] = hapData[item][0]
2253                    hapDict[pfx+item+';mx'] = hapData[item][1][2]
2254                    if hapData[item][2][2]:
2255                        hapVary.append(pfx+item+';mx')
2256                    if hapData[item][0] in ['isotropic','uniaxial']:
2257                        hapDict[pfx+item+';i'] = hapData[item][1][0]
2258                        if hapData[item][2][0]:
2259                            hapVary.append(pfx+item+';i')
2260                        if hapData[item][0] == 'uniaxial':
2261                            controlDict[pfx+item+'Axis'] = hapData[item][3]
2262                            hapDict[pfx+item+';a'] = hapData[item][1][1]
2263                            if hapData[item][2][1]:
2264                                hapVary.append(pfx+item+';a')
2265                    else:       #generalized for mustrain or ellipsoidal for size
2266                        Nterms = len(hapData[item][4])
2267                        if item == 'Mustrain':
2268                            names = G2spc.MustrainNames(SGData)
2269                            pwrs = []
2270                            for name in names:
2271                                h,k,l = name[1:]
2272                                pwrs.append([int(h),int(k),int(l)])
2273                            controlDict[pfx+'MuPwrs'] = pwrs
2274                        for i in range(Nterms):
2275                            sfx = ';'+str(i)
2276                            hapDict[pfx+item+sfx] = hapData[item][4][i]
2277                            if hapData[item][5][i]:
2278                                hapVary.append(pfx+item+sfx)
2279                for bab in ['BabA','BabU']:
2280                    hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2281                    if hapData['Babinet'][bab][1]:
2282                        hapVary.append(pfx+bab)
2283                               
2284                if Print: 
2285                    print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
2286                    print >>pFile,135*'-'
2287                    print >>pFile,' Phase fraction  : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
2288                    print >>pFile,' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
2289                    if hapData['Pref.Ori.'][0] == 'MD':
2290                        Ax = hapData['Pref.Ori.'][3]
2291                        print >>pFile,' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
2292                            ' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
2293                    else: #'SH' for spherical harmonics
2294                        PrintSHPO(hapData['Pref.Ori.'])
2295                        print >>pFile,' Penalty hkl list: '+str(controlDict[pfx+'SHhkl'])+' tolerance: %.2f'%(controlDict[pfx+'SHtoler'])
2296                    PrintSize(hapData['Size'])
2297                    PrintMuStrain(hapData['Mustrain'],SGData)
2298                    PrintHStrain(hapData['HStrain'],SGData)
2299                    if hapData['Babinet']['BabA'][0]:
2300                        PrintBabinet(hapData['Babinet'])                       
2301                if resetRefList:
2302                    refList = []
2303                    Uniq = []
2304                    Phi = []
2305                    useExt = 'magnetic' in Phases[phase]['General']['Type'] and 'N' in inst['Type'][0]
2306                    if Phases[phase]['General'].get('Modulated',False):
2307                        ifSuper = True
2308                        HKLd = np.array(G2lat.GenSSHLaue(dmin,SGData,SSGData,Vec,maxH,A))
2309                        HKLd = G2mth.sortArray(HKLd,4,reverse=True)
2310                        for h,k,l,m,d in HKLd:
2311                            ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2312                            mul *= 2      # for powder overlap of Friedel pairs
2313                            if m or not ext or useExt:
2314                                if 'C' in inst['Type'][0]:
2315                                    pos = G2lat.Dsp2pos(inst,d)
2316                                    if limits[0] < pos < limits[1]:
2317                                        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])
2318                                        #... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2319                                        Uniq.append(uniq)
2320                                        Phi.append(phi)
2321                                elif 'T' in inst['Type'][0]:
2322                                    pos = G2lat.Dsp2pos(inst,d)
2323                                    if limits[0] < pos < limits[1]:
2324                                        wave = inst['difC'][1]*d/(252.816*inst['fltPath'][0])
2325                                        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])
2326                                        # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2327                                        Uniq.append(uniq)
2328                                        Phi.append(phi)
2329                    else:
2330                        ifSuper = False
2331                        HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
2332                        HKLd = G2mth.sortArray(HKLd,3,reverse=True)
2333                        for h,k,l,d in HKLd:
2334                            ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2335                            mul *= 2      # for powder overlap of Friedel pairs
2336                            if ext and not useExt:
2337                                continue
2338                            if 'C' in inst['Type'][0]:
2339                                pos = G2lat.Dsp2pos(inst,d)
2340                                if limits[0] < pos < limits[1]:
2341                                    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])
2342                                    #... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2343                                    Uniq.append(uniq)
2344                                    Phi.append(phi)
2345                            elif 'T' in inst['Type'][0]:
2346                                pos = G2lat.Dsp2pos(inst,d)
2347                                if limits[0] < pos < limits[1]:
2348                                    wave = inst['difC'][1]*d/(252.816*inst['fltPath'][0])
2349                                    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])
2350                                    # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2351                                    Uniq.append(uniq)
2352                                    Phi.append(phi)
2353                    Histogram['Reflection Lists'][phase] = {'RefList':np.array(refList),'FF':{},'Type':inst['Type'][0],'Super':ifSuper}
2354            elif 'HKLF' in histogram:
2355                inst = Histogram['Instrument Parameters'][0]
2356                hId = Histogram['hId']
2357                hfx = ':%d:'%(hId)
2358                for item in inst:
2359                    if type(inst) is not list and item != 'Type': continue # skip over non-refined items (such as InstName)
2360                    hapDict[hfx+item] = inst[item][1]
2361                pfx = str(pId)+':'+str(hId)+':'
2362                hapDict[pfx+'Scale'] = hapData['Scale'][0]
2363                if hapData['Scale'][1]:
2364                    hapVary.append(pfx+'Scale')
2365                               
2366                extApprox,extType,extParms = hapData['Extinction']
2367                controlDict[pfx+'EType'] = extType
2368                controlDict[pfx+'EApprox'] = extApprox
2369                if 'C' in inst['Type'][0]:
2370                    controlDict[pfx+'Tbar'] = extParms['Tbar']
2371                    controlDict[pfx+'Cos2TM'] = extParms['Cos2TM']
2372                if 'Primary' in extType:
2373                    Ekey = ['Ep',]
2374                elif 'I & II' in extType:
2375                    Ekey = ['Eg','Es']
2376                elif 'Secondary Type II' == extType:
2377                    Ekey = ['Es',]
2378                elif 'Secondary Type I' == extType:
2379                    Ekey = ['Eg',]
2380                else:   #'None'
2381                    Ekey = []
2382                for eKey in Ekey:
2383                    hapDict[pfx+eKey] = extParms[eKey][0]
2384                    if extParms[eKey][1]:
2385                        hapVary.append(pfx+eKey)
2386                for bab in ['BabA','BabU']:
2387                    hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2388                    if hapData['Babinet'][bab][1]:
2389                        hapVary.append(pfx+bab)
2390                Twins = hapData.get('Twins',[[np.array([[1,0,0],[0,1,0],[0,0,1]]),[1.0,False,0]],])
2391                if len(Twins) == 1:
2392                    hapDict[pfx+'Flack'] = hapData.get('Flack',[0.,False])[0]
2393                    if hapData.get('Flack',[0,False])[1]:
2394                        hapVary.append(pfx+'Flack')
2395                sumTwFr = 0.
2396                controlDict[pfx+'TwinLaw'] = []
2397                controlDict[pfx+'TwinInv'] = []
2398                NTL = 0           
2399                for it,twin in enumerate(Twins):
2400                    if 'bool' in str(type(twin[0])):
2401                        controlDict[pfx+'TwinInv'].append(twin[0])
2402                        controlDict[pfx+'TwinLaw'].append(np.zeros((3,3)))
2403                    else:
2404                        NTL += 1
2405                        controlDict[pfx+'TwinInv'].append(False)
2406                        controlDict[pfx+'TwinLaw'].append(twin[0])
2407                    if it:
2408                        hapDict[pfx+'TwinFr:'+str(it)] = twin[1]
2409                        sumTwFr += twin[1]
2410                    else:
2411                        hapDict[pfx+'TwinFr:0'] = twin[1][0]
2412                        controlDict[pfx+'TwinNMN'] = twin[1][1]
2413                    if Twins[0][1][1]:
2414                        hapVary.append(pfx+'TwinFr:'+str(it))
2415                controlDict[pfx+'NTL'] = NTL
2416                #need to make constraint on TwinFr
2417                controlDict[pfx+'TwinLaw'] = np.array(controlDict[pfx+'TwinLaw'])
2418                if len(Twins) > 1:    #force sum to unity
2419                    hapDict[pfx+'TwinFr:0'] = 1.-sumTwFr
2420                if Print: 
2421                    print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
2422                    print >>pFile,135*'-'
2423                    print >>pFile,' Scale factor     : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
2424                    if extType != 'None':
2425                        print >>pFile,' Extinction  Type: %15s'%(extType),' approx: %10s'%(extApprox)
2426                        text = ' Parameters       :'
2427                        for eKey in Ekey:
2428                            text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
2429                        print >>pFile,text
2430                    if hapData['Babinet']['BabA'][0]:
2431                        PrintBabinet(hapData['Babinet'])
2432                    if not SGData['SGInv'] and len(Twins) == 1:
2433                        print >>pFile,' Flack parameter: %10.3f'%(hapData['Flack'][0]),' Refine?',hapData['Flack'][1]
2434                    if len(Twins) > 1:
2435                        for it,twin in enumerate(Twins):
2436                            if 'bool' in str(type(twin[0])):
2437                                print >>pFile,' Nonmerohedral twin fr.: %5.3f Inv? %s Refine? '%(hapDict[pfx+'TwinFr:'+str(it)],str(controlDict[pfx+'TwinInv'][it])),Twins[0][1][1] 
2438                            else:
2439                                print >>pFile,' Twin law: %s'%(str(twin[0]).replace('\n',',')),' Twin fr.: %5.3f Refine? '%(hapDict[pfx+'TwinFr:'+str(it)]),Twins[0][1][1] 
2440                       
2441                Histogram['Reflection Lists'] = phase       
2442               
2443    return hapVary,hapDict,controlDict
2444   
2445def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,FFtables,Print=True,pFile=None):
2446    'needs a doc string'
2447   
2448    def PrintSizeAndSig(hapData,sizeSig):
2449        line = '\n Size model:     %9s'%(hapData[0])
2450        refine = False
2451        if hapData[0] in ['isotropic','uniaxial']:
2452            line += ' equatorial:%12.5f'%(hapData[1][0])
2453            if sizeSig[0][0]:
2454                line += ', sig:%8.4f'%(sizeSig[0][0])
2455                refine = True
2456            if hapData[0] == 'uniaxial':
2457                line += ' axial:%12.4f'%(hapData[1][1])
2458                if sizeSig[0][1]:
2459                    refine = True
2460                    line += ', sig:%8.4f'%(sizeSig[0][1])
2461            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2462            if sizeSig[0][2]:
2463                refine = True
2464                line += ', sig:%8.4f'%(sizeSig[0][2])
2465            if refine:
2466                print >>pFile,line
2467        else:
2468            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2469            if sizeSig[0][2]:
2470                refine = True
2471                line += ', sig:%8.4f'%(sizeSig[0][2])
2472            Snames = ['S11','S22','S33','S12','S13','S23']
2473            ptlbls = ' name  :'
2474            ptstr =  ' value :'
2475            sigstr = ' sig   :'
2476            for i,name in enumerate(Snames):
2477                ptlbls += '%12s' % (name)
2478                ptstr += '%12.6f' % (hapData[4][i])
2479                if sizeSig[1][i]:
2480                    refine = True
2481                    sigstr += '%12.6f' % (sizeSig[1][i])
2482                else:
2483                    sigstr += 12*' '
2484            if refine:
2485                print >>pFile,line
2486                print >>pFile,ptlbls
2487                print >>pFile,ptstr
2488                print >>pFile,sigstr
2489       
2490    def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
2491        line = '\n Mustrain model: %9s'%(hapData[0])
2492        refine = False
2493        if hapData[0] in ['isotropic','uniaxial']:
2494            line += ' equatorial:%12.1f'%(hapData[1][0])
2495            if mustrainSig[0][0]:
2496                line += ', sig:%8.1f'%(mustrainSig[0][0])
2497                refine = True
2498            if hapData[0] == 'uniaxial':
2499                line += ' axial:%12.1f'%(hapData[1][1])
2500                if mustrainSig[0][1]:
2501                     line += ', sig:%8.1f'%(mustrainSig[0][1])
2502            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2503            if mustrainSig[0][2]:
2504                refine = True
2505                line += ', sig:%8.3f'%(mustrainSig[0][2])
2506            if refine:
2507                print >>pFile,line
2508        else:
2509            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2510            if mustrainSig[0][2]:
2511                refine = True
2512                line += ', sig:%8.3f'%(mustrainSig[0][2])
2513            Snames = G2spc.MustrainNames(SGData)
2514            ptlbls = ' name  :'
2515            ptstr =  ' value :'
2516            sigstr = ' sig   :'
2517            for i,name in enumerate(Snames):
2518                ptlbls += '%12s' % (name)
2519                ptstr += '%12.1f' % (hapData[4][i])
2520                if mustrainSig[1][i]:
2521                    refine = True
2522                    sigstr += '%12.1f' % (mustrainSig[1][i])
2523                else:
2524                    sigstr += 12*' '
2525            if refine:
2526                print >>pFile,line
2527                print >>pFile,ptlbls
2528                print >>pFile,ptstr
2529                print >>pFile,sigstr
2530           
2531    def PrintHStrainAndSig(hapData,strainSig,SGData):
2532        Hsnames = G2spc.HStrainNames(SGData)
2533        ptlbls = ' name  :'
2534        ptstr =  ' value :'
2535        sigstr = ' sig   :'
2536        refine = False
2537        for i,name in enumerate(Hsnames):
2538            ptlbls += '%12s' % (name)
2539            ptstr += '%12.4g' % (hapData[0][i])
2540            if name in strainSig:
2541                refine = True
2542                sigstr += '%12.4g' % (strainSig[name])
2543            else:
2544                sigstr += 12*' '
2545        if refine:
2546            print >>pFile,'\n Hydrostatic/elastic strain: '
2547            print >>pFile,ptlbls
2548            print >>pFile,ptstr
2549            print >>pFile,sigstr
2550       
2551    def PrintSHPOAndSig(pfx,hapData,POsig):
2552        print >>pFile,'\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
2553        ptlbls = ' names :'
2554        ptstr =  ' values:'
2555        sigstr = ' sig   :'
2556        for item in hapData[5]:
2557            ptlbls += '%12s'%(item)
2558            ptstr += '%12.3f'%(hapData[5][item])
2559            if pfx+item in POsig:
2560                sigstr += '%12.3f'%(POsig[pfx+item])
2561            else:
2562                sigstr += 12*' ' 
2563        print >>pFile,ptlbls
2564        print >>pFile,ptstr
2565        print >>pFile,sigstr
2566       
2567    def PrintExtAndSig(pfx,hapData,ScalExtSig):
2568        print >>pFile,'\n Single crystal extinction: Type: ',hapData[0],' Approx: ',hapData[1]
2569        text = ''
2570        for item in hapData[2]:
2571            if pfx+item in ScalExtSig:
2572                text += '       %s: '%(item)
2573                text += '%12.2e'%(hapData[2][item][0])
2574                if pfx+item in ScalExtSig:
2575                    text += ' sig: %12.2e'%(ScalExtSig[pfx+item])
2576        print >>pFile,text       
2577       
2578    def PrintBabinetAndSig(pfx,hapData,BabSig):
2579        print >>pFile,'\n Babinet form factor modification: '
2580        ptlbls = ' names :'
2581        ptstr =  ' values:'
2582        sigstr = ' sig   :'
2583        for item in hapData:
2584            ptlbls += '%12s'%(item)
2585            ptstr += '%12.3f'%(hapData[item][0])
2586            if pfx+item in BabSig:
2587                sigstr += '%12.3f'%(BabSig[pfx+item])
2588            else:
2589                sigstr += 12*' ' 
2590        print >>pFile,ptlbls
2591        print >>pFile,ptstr
2592        print >>pFile,sigstr
2593       
2594    def PrintTwinsAndSig(pfx,twinData,TwinSig):
2595        print >>pFile,'\n Twin Law fractions : '
2596        ptlbls = ' names :'
2597        ptstr =  ' values:'
2598        sigstr = ' sig   :'
2599        for it,item in enumerate(twinData):
2600            ptlbls += '     twin #%d'%(it)
2601            if it:
2602                ptstr += '%12.3f'%(item[1])
2603            else:
2604                ptstr += '%12.3f'%(item[1][0])
2605            if pfx+'TwinFr:'+str(it) in TwinSig:
2606                sigstr += '%12.3f'%(TwinSig[pfx+'TwinFr:'+str(it)])
2607            else:
2608                sigstr += 12*' ' 
2609        print >>pFile,ptlbls
2610        print >>pFile,ptstr
2611        print >>pFile,sigstr
2612       
2613   
2614    PhFrExtPOSig = {}
2615    SizeMuStrSig = {}
2616    ScalExtSig = {}
2617    BabSig = {}
2618    TwinFrSig = {}
2619    wtFrSum = {}
2620    for phase in Phases:
2621        HistoPhase = Phases[phase]['Histograms']
2622        General = Phases[phase]['General']
2623        SGData = General['SGData']
2624        pId = Phases[phase]['pId']
2625        histoList = HistoPhase.keys()
2626        histoList.sort()
2627        for histogram in histoList:
2628            try:
2629                Histogram = Histograms[histogram]
2630            except KeyError:                       
2631                #skip if histogram not included e.g. in a sequential refinement
2632                continue
2633            hapData = HistoPhase[histogram]
2634            hId = Histogram['hId']
2635            pfx = str(pId)+':'+str(hId)+':'
2636            if hId not in wtFrSum:
2637                wtFrSum[hId] = 0.
2638            if 'PWDR' in histogram:
2639                for item in ['Scale','Extinction']:
2640                    hapData[item][0] = parmDict[pfx+item]
2641                    if pfx+item in sigDict:
2642                        PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2643                wtFrSum[hId] += hapData['Scale'][0]*General['Mass']
2644                if hapData['Pref.Ori.'][0] == 'MD':
2645                    hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
2646                    if pfx+'MD' in sigDict:
2647                        PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],})
2648                else:                           #'SH' spherical harmonics
2649                    for item in hapData['Pref.Ori.'][5]:
2650                        hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
2651                        if pfx+item in sigDict:
2652                            PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2653                SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
2654                    pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
2655                    pfx+'HStrain':{}})                 
2656                for item in ['Mustrain','Size']:
2657                    hapData[item][1][2] = parmDict[pfx+item+';mx']
2658#                    hapData[item][1][2] = min(1.,max(0.,hapData[item][1][2]))
2659                    if pfx+item+';mx' in sigDict:
2660                        SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
2661                    if hapData[item][0] in ['isotropic','uniaxial']:                   
2662                        hapData[item][1][0] = parmDict[pfx+item+';i']
2663                        if item == 'Size':
2664                            hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
2665                        if pfx+item+';i' in sigDict: 
2666                            SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
2667                        if hapData[item][0] == 'uniaxial':
2668                            hapData[item][1][1] = parmDict[pfx+item+';a']
2669                            if item == 'Size':
2670                                hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))                       
2671                            if pfx+item+';a' in sigDict:
2672                                SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
2673                    else:       #generalized for mustrain or ellipsoidal for size
2674                        Nterms = len(hapData[item][4])
2675                        for i in range(Nterms):
2676                            sfx = ';'+str(i)
2677                            hapData[item][4][i] = parmDict[pfx+item+sfx]
2678                            if pfx+item+sfx in sigDict:
2679                                SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx]
2680                names = G2spc.HStrainNames(SGData)
2681                for i,name in enumerate(names):
2682                    hapData['HStrain'][0][i] = parmDict[pfx+name]
2683                    if pfx+name in sigDict:
2684                        SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
2685                for name in ['BabA','BabU']:
2686                    hapData['Babinet'][name][0] = parmDict[pfx+name]
2687                    if pfx+name in sigDict:
2688                        BabSig[pfx+name] = sigDict[pfx+name]               
2689               
2690            elif 'HKLF' in histogram:
2691                for item in ['Scale','Flack']:
2692                    if parmDict.get(pfx+item):
2693                        hapData[item][0] = parmDict[pfx+item]
2694                        if pfx+item in sigDict:
2695                            ScalExtSig[pfx+item] = sigDict[pfx+item]
2696                for item in ['Ep','Eg','Es']:
2697                    if parmDict.get(pfx+item):
2698                        hapData['Extinction'][2][item][0] = parmDict[pfx+item]
2699                        if pfx+item in sigDict:
2700                            ScalExtSig[pfx+item] = sigDict[pfx+item]
2701                for item in ['BabA','BabU']:
2702                    hapData['Babinet'][item][0] = parmDict[pfx+item]
2703                    if pfx+item in sigDict:
2704                        BabSig[pfx+item] = sigDict[pfx+item]
2705                if 'Twins' in hapData:
2706                    it = 1
2707                    sumTwFr = 0.
2708                    while True:
2709                        try:
2710                            hapData['Twins'][it][1] = parmDict[pfx+'TwinFr:'+str(it)]
2711                            if pfx+'TwinFr:'+str(it) in sigDict:
2712                                TwinFrSig[pfx+'TwinFr:'+str(it)] = sigDict[pfx+'TwinFr:'+str(it)]
2713                            if it:
2714                                sumTwFr += hapData['Twins'][it][1]
2715                            it += 1
2716                        except KeyError:
2717                            break
2718                    hapData['Twins'][0][1][0] = 1.-sumTwFr
2719
2720    if Print:
2721        for phase in Phases:
2722            HistoPhase = Phases[phase]['Histograms']
2723            General = Phases[phase]['General']
2724            SGData = General['SGData']
2725            pId = Phases[phase]['pId']
2726            histoList = HistoPhase.keys()
2727            histoList.sort()
2728            for histogram in histoList:
2729                try:
2730                    Histogram = Histograms[histogram]
2731                except KeyError:                       
2732                    #skip if histogram not included e.g. in a sequential refinement
2733                    continue
2734                print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
2735                print >>pFile,130*'-'
2736                hapData = HistoPhase[histogram]
2737                hId = Histogram['hId']
2738                Histogram['Residuals'][str(pId)+'::Name'] = phase
2739                pfx = str(pId)+':'+str(hId)+':'
2740                hfx = ':%s:'%(hId)
2741                if 'PWDR' in histogram:
2742                    print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections'   \
2743                        %(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'])
2744                    print >>pFile,' Bragg intensity sum = %.3g'%(Histogram['Residuals'][pfx+'sumInt'])
2745               
2746                    if pfx+'Scale' in PhFrExtPOSig:
2747                        wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum[hId]
2748                        sigwtFr = PhFrExtPOSig[pfx+'Scale']*wtFr/hapData['Scale'][0]
2749                        print >>pFile,' Phase fraction  : %10.5f, sig %10.5f Weight fraction  : %8.5f, sig %10.5f' \
2750                            %(hapData['Scale'][0],PhFrExtPOSig[pfx+'Scale'],wtFr,sigwtFr)
2751                    if pfx+'Extinction' in PhFrExtPOSig:
2752                        print >>pFile,' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction'])
2753                    if hapData['Pref.Ori.'][0] == 'MD':
2754                        if pfx+'MD' in PhFrExtPOSig:
2755                            print >>pFile,' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD'])
2756                    else:
2757                        PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig)
2758                    PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size'])
2759                    PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
2760                    PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
2761                    if len(BabSig):
2762                        PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2763                   
2764                elif 'HKLF' in histogram:
2765                    Inst = Histogram['Instrument Parameters'][0]
2766                    print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections (%d user rejected, %d sp.gp.extinct)'   \
2767                        %(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'],
2768                        Histogram['Residuals'][pfx+'Nrej'],Histogram['Residuals'][pfx+'Next'])
2769                    if FFtables != None and 'N' not in Inst['Type'][0]:
2770                        PrintFprime(FFtables,hfx,pFile)
2771                    print >>pFile,' HKLF histogram weight factor = ','%.3f'%(Histogram['wtFactor'])
2772                    if pfx+'Scale' in ScalExtSig:
2773                        print >>pFile,' Scale factor : %10.4f, sig %10.4f'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale'])
2774                    if hapData['Extinction'][0] != 'None':
2775                        PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig)
2776                    if len(BabSig):
2777                        PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2778                    if pfx+'Flack' in ScalExtSig:
2779                        print >>pFile,' Flack parameter : %10.3f, sig %10.3f'%(hapData['Flack'][0],ScalExtSig[pfx+'Flack'])
2780                    if len(TwinFrSig):
2781                        PrintTwinsAndSig(pfx,hapData['Twins'],TwinFrSig)
2782
2783################################################################################
2784##### Histogram data
2785################################################################################       
2786                   
2787def GetHistogramData(Histograms,Print=True,pFile=None):
2788    'needs a doc string'
2789   
2790    def GetBackgroundParms(hId,Background):
2791        Back = Background[0]
2792        DebyePeaks = Background[1]
2793        bakType,bakFlag = Back[:2]
2794        backVals = Back[3:]
2795        backNames = [':'+str(hId)+':Back;'+str(i) for i in range(len(backVals))]
2796        backDict = dict(zip(backNames,backVals))
2797        backVary = []
2798        if bakFlag:
2799            backVary = backNames
2800        backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
2801        backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
2802        debyeDict = {}
2803        debyeList = []
2804        for i in range(DebyePeaks['nDebye']):
2805            debyeNames = [':'+str(hId)+':DebyeA;'+str(i),':'+str(hId)+':DebyeR;'+str(i),':'+str(hId)+':DebyeU;'+str(i)]
2806            debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
2807            debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
2808        debyeVary = []
2809        for item in debyeList:
2810            if item[1]:
2811                debyeVary.append(item[0])
2812        backDict.update(debyeDict)
2813        backVary += debyeVary
2814        peakDict = {}
2815        peakList = []
2816        for i in range(DebyePeaks['nPeaks']):
2817            peakNames = [':'+str(hId)+':BkPkpos;'+str(i),':'+str(hId)+ \
2818                ':BkPkint;'+str(i),':'+str(hId)+':BkPksig;'+str(i),':'+str(hId)+':BkPkgam;'+str(i)]
2819            peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
2820            peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
2821        peakVary = []
2822        for item in peakList:
2823            if item[1]:
2824                peakVary.append(item[0])
2825        backDict.update(peakDict)
2826        backVary += peakVary
2827        return bakType,backDict,backVary           
2828       
2829    def GetInstParms(hId,Inst):     
2830        dataType = Inst['Type'][0]
2831        instDict = {}
2832        insVary = []
2833        pfx = ':'+str(hId)+':'
2834        insKeys = Inst.keys()
2835        insKeys.sort()
2836        for item in insKeys:
2837            insName = pfx+item
2838            instDict[insName] = Inst[item][1]
2839            if len(Inst[item]) > 2 and Inst[item][2]:
2840                insVary.append(insName)
2841        if 'C' in dataType:
2842            instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
2843        return dataType,instDict,insVary
2844       
2845    def GetSampleParms(hId,Sample):
2846        sampVary = []
2847        hfx = ':'+str(hId)+':'       
2848        sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
2849            hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
2850        for key in ('Temperature','Pressure','FreePrm1','FreePrm2','FreePrm3'):
2851            if key in Sample:
2852                sampDict[hfx+key] = Sample[key]
2853        Type = Sample['Type']
2854        if 'Bragg' in Type:             #Bragg-Brentano
2855            for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2856                sampDict[hfx+item] = Sample[item][0]
2857                if Sample[item][1]:
2858                    sampVary.append(hfx+item)
2859        elif 'Debye' in Type:        #Debye-Scherrer
2860            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2861                sampDict[hfx+item] = Sample[item][0]
2862                if Sample[item][1]:
2863                    sampVary.append(hfx+item)
2864        return Type,sampDict,sampVary
2865       
2866    def PrintBackground(Background):
2867        Back = Background[0]
2868        DebyePeaks = Background[1]
2869        print >>pFile,'\n Background function: ',Back[0],' Refine?',bool(Back[1])
2870        line = ' Coefficients: '
2871        for i,back in enumerate(Back[3:]):
2872            line += '%10.3f'%(back)
2873            if i and not i%10:
2874                line += '\n'+15*' '
2875        print >>pFile,line
2876        if DebyePeaks['nDebye']:
2877            print >>pFile,'\n Debye diffuse scattering coefficients'
2878            parms = ['DebyeA','DebyeR','DebyeU']
2879            line = ' names :  '
2880            for parm in parms:
2881                line += '%8s refine?'%(parm)
2882            print >>pFile,line
2883            for j,term in enumerate(DebyePeaks['debyeTerms']):
2884                line = ' term'+'%2d'%(j)+':'
2885                for i in range(3):
2886                    line += '%10.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
2887                print >>pFile,line
2888        if DebyePeaks['nPeaks']:
2889            print >>pFile,'\n Single peak coefficients'
2890            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
2891            line = ' names :  '
2892            for parm in parms:
2893                line += '%8s refine?'%(parm)
2894            print >>pFile,line
2895            for j,term in enumerate(DebyePeaks['peaksList']):
2896                line = ' peak'+'%2d'%(j)+':'
2897                for i in range(4):
2898                    line += '%12.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
2899                print >>pFile,line
2900       
2901    def PrintInstParms(Inst):
2902        print >>pFile,'\n Instrument Parameters:'
2903        insKeys = Inst.keys()
2904        insKeys.sort()
2905        iBeg = 0
2906        Ok = True
2907        while Ok:
2908            ptlbls = ' name  :'
2909            ptstr =  ' value :'
2910            varstr = ' refine:'
2911            iFin = min(iBeg+9,len(insKeys))
2912            for item in insKeys[iBeg:iFin]:
2913                if item not in ['Type','Source']:
2914                    ptlbls += '%12s' % (item)
2915                    ptstr += '%12.6f' % (Inst[item][1])
2916                    if item in ['Lam1','Lam2','Azimuth','fltPath','2-theta',]:
2917                        varstr += 12*' '
2918                    else:
2919                        varstr += '%12s' % (str(bool(Inst[item][2])))
2920            print >>pFile,ptlbls
2921            print >>pFile,ptstr
2922            print >>pFile,varstr
2923            iBeg = iFin
2924            if iBeg == len(insKeys):
2925                Ok = False
2926            else:
2927                print >>pFile,'\n'
2928       
2929    def PrintSampleParms(Sample):
2930        print >>pFile,'\n Sample Parameters:'
2931        print >>pFile,' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \
2932            (Sample['Omega'],Sample['Chi'],Sample['Phi'])
2933        ptlbls = ' name  :'
2934        ptstr =  ' value :'
2935        varstr = ' refine:'
2936        if 'Bragg' in Sample['Type']:
2937            for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2938                ptlbls += '%14s'%(item)
2939                ptstr += '%14.4f'%(Sample[item][0])
2940                varstr += '%14s'%(str(bool(Sample[item][1])))
2941           
2942        elif 'Debye' in Type:        #Debye-Scherrer
2943            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2944                ptlbls += '%14s'%(item)
2945                ptstr += '%14.4f'%(Sample[item][0])
2946                varstr += '%14s'%(str(bool(Sample[item][1])))
2947
2948        print >>pFile,ptlbls
2949        print >>pFile,ptstr
2950        print >>pFile,varstr
2951       
2952    histDict = {}
2953    histVary = []
2954    controlDict = {}
2955    histoList = Histograms.keys()
2956    histoList.sort()
2957    for histogram in histoList:
2958        if 'PWDR' in histogram:
2959            Histogram = Histograms[histogram]
2960            hId = Histogram['hId']
2961            pfx = ':'+str(hId)+':'
2962            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
2963            controlDict[pfx+'Limits'] = Histogram['Limits'][1]
2964            controlDict[pfx+'Exclude'] = Histogram['Limits'][2:]
2965            for excl in controlDict[pfx+'Exclude']:
2966                Histogram['Data'][0] = ma.masked_inside(Histogram['Data'][0],excl[0],excl[1])
2967            if controlDict[pfx+'Exclude']:
2968                ma.mask_rows(Histogram['Data'])
2969            Background = Histogram['Background']
2970            Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
2971            controlDict[pfx+'bakType'] = Type
2972            histDict.update(bakDict)
2973            histVary += bakVary
2974           
2975            Inst = Histogram['Instrument Parameters'][0]
2976            Type,instDict,insVary = GetInstParms(hId,Inst)
2977            controlDict[pfx+'histType'] = Type
2978            if 'XC' in Type:
2979                if pfx+'Lam1' in instDict:
2980                    controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
2981                else:
2982                    controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']           
2983            histDict.update(instDict)
2984            histVary += insVary
2985           
2986            Sample = Histogram['Sample Parameters']
2987            Type,sampDict,sampVary = GetSampleParms(hId,Sample)
2988            controlDict[pfx+'instType'] = Type
2989            histDict.update(sampDict)
2990            histVary += sampVary
2991           
2992   
2993            if Print: 
2994                print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
2995                print >>pFile,135*'-'
2996                Units = {'C':' deg','T':' msec'}
2997                units = Units[controlDict[pfx+'histType'][2]]
2998                Limits = controlDict[pfx+'Limits']
2999                print >>pFile,' Instrument type: ',Sample['Type']
3000                print >>pFile,' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)
3001                if len(controlDict[pfx+'Exclude']):
3002                    excls = controlDict[pfx+'Exclude']
3003                    for excl in excls:
3004                        print >>pFile,' Excluded region:  %8.2f%s to %8.2f%s'%(excl[0],units,excl[1],units)   
3005                PrintSampleParms(Sample)
3006                PrintInstParms(Inst)
3007                PrintBackground(Background)
3008        elif 'HKLF' in histogram:
3009            Histogram = Histograms[histogram]
3010            hId = Histogram['hId']
3011            pfx = ':'+str(hId)+':'
3012            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
3013            Inst = Histogram['Instrument Parameters'][0]
3014            controlDict[pfx+'histType'] = Inst['Type'][0]
3015            if 'X' in Inst['Type'][0]:
3016                histDict[pfx+'Lam'] = Inst['Lam'][1]
3017                controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']                   
3018    return histVary,histDict,controlDict
3019   
3020def SetHistogramData(parmDict,sigDict,Histograms,FFtables,Print=True,pFile=None):
3021    'needs a doc string'
3022   
3023    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
3024        Back = Background[0]
3025        DebyePeaks = Background[1]
3026        lenBack = len(Back[3:])
3027        backSig = [0 for i in range(lenBack+3*DebyePeaks['nDebye']+4*DebyePeaks['nPeaks'])]
3028        for i in range(lenBack):
3029            Back[3+i] = parmDict[pfx+'Back;'+str(i)]
3030            if pfx+'Back;'+str(i) in sigDict:
3031                backSig[i] = sigDict[pfx+'Back;'+str(i)]
3032        if DebyePeaks['nDebye']:
3033            for i in range(DebyePeaks['nDebye']):
3034                names = [pfx+'DebyeA;'+str(i),pfx+'DebyeR;'+str(i),pfx+'DebyeU;'+str(i)]
3035                for j,name in enumerate(names):
3036                    DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
3037                    if name in sigDict:
3038                        backSig[lenBack+3*i+j] = sigDict[name]           
3039        if DebyePeaks['nPeaks']:
3040            for i in range(DebyePeaks['nPeaks']):
3041                names = [pfx+'BkPkpos;'+str(i),pfx+'BkPkint;'+str(i),
3042                    pfx+'BkPksig;'+str(i),pfx+'BkPkgam;'+str(i)]
3043                for j,name in enumerate(names):
3044                    DebyePeaks['peaksList'][i][2*j] = parmDict[name]
3045                    if name in sigDict:
3046                        backSig[lenBack+3*DebyePeaks['nDebye']+4*i+j] = sigDict[name]
3047        return backSig
3048       
3049    def SetInstParms(pfx,Inst,parmDict,sigDict):
3050        instSig = {}
3051        insKeys = Inst.keys()
3052        insKeys.sort()
3053        for item in insKeys:
3054            insName = pfx+item
3055            Inst[item][1] = parmDict[insName]
3056            if insName in sigDict:
3057                instSig[item] = sigDict[insName]
3058            else:
3059                instSig[item] = 0
3060        return instSig
3061       
3062    def SetSampleParms(pfx,Sample,parmDict,sigDict):
3063        if 'Bragg' in Sample['Type']:             #Bragg-Brentano
3064            sampSig = [0 for i in range(5)]
3065            for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3066                Sample[item][0] = parmDict[pfx+item]
3067                if pfx+item in sigDict:
3068                    sampSig[i] = sigDict[pfx+item]
3069        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
3070            sampSig = [0 for i in range(4)]
3071            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3072                Sample[item][0] = parmDict[pfx+item]
3073                if pfx+item in sigDict:
3074                    sampSig[i] = sigDict[pfx+item]
3075        return sampSig
3076       
3077    def PrintBackgroundSig(Background,backSig):
3078        Back = Background[0]
3079        DebyePeaks = Background[1]
3080        lenBack = len(Back[3:])
3081        valstr = ' value : '
3082        sigstr = ' sig   : '
3083        refine = False
3084        for i,back in enumerate(Back[3:]):
3085            valstr += '%10.4g'%(back)
3086            if Back[1]:
3087                refine = True
3088                sigstr += '%10.4g'%(backSig[i])
3089            else:
3090                sigstr += 10*' '
3091        if refine:
3092            print >>pFile,'\n Background function: ',Back[0]
3093            print >>pFile,valstr
3094            print >>pFile,sigstr
3095        if DebyePeaks['nDebye']:
3096            ifAny = False
3097            ptfmt = "%12.3f"
3098            names =  ' names :'
3099            ptstr =  ' values:'
3100            sigstr = ' esds  :'
3101            for item in sigDict:
3102                if 'Debye' in item:
3103                    ifAny = True
3104                    names += '%12s'%(item)
3105                    ptstr += ptfmt%(parmDict[item])
3106                    sigstr += ptfmt%(sigDict[item])
3107            if ifAny:
3108                print >>pFile,'\n Debye diffuse scattering coefficients'
3109                print >>pFile,names
3110                print >>pFile,ptstr
3111                print >>pFile,sigstr
3112        if DebyePeaks['nPeaks']:
3113            print >>pFile,'\n Single peak coefficients:'
3114            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
3115            line = ' peak no. '
3116            for parm in parms:
3117                line += '%14s%12s'%(parm.center(14),'esd'.center(12))
3118            print >>pFile,line
3119            for ip in range(DebyePeaks['nPeaks']):
3120                ptstr = ' %4d '%(ip)
3121                for parm in parms:
3122                    fmt = '%14.3f'
3123                    efmt = '%12.3f'
3124                    if parm == 'BkPkpos':
3125                        fmt = '%14.4f'
3126                        efmt = '%12.4f'
3127                    name = pfx+parm+';%d'%(ip)
3128                    ptstr += fmt%(parmDict[name])
3129                    if name in sigDict:
3130                        ptstr += efmt%(sigDict[name])
3131                    else:
3132                        ptstr += 12*' '
3133                print >>pFile,ptstr
3134        sumBk = np.array(Histogram['sumBk'])
3135        print >>pFile,' Background sums: empirical %.3g, Debye %.3g, peaks %.3g, Total %.3g'    \
3136            %(sumBk[0],sumBk[1],sumBk[2],np.sum(sumBk))
3137       
3138    def PrintInstParmsSig(Inst,instSig):
3139        refine = False
3140        insKeys = instSig.keys()
3141        insKeys.sort()
3142        iBeg = 0
3143        Ok = True
3144        while Ok:
3145            ptlbls = ' names :'
3146            ptstr =  ' value :'
3147            sigstr = ' sig   :'
3148            iFin = min(iBeg+9,len(insKeys))
3149            for name in insKeys[iBeg:iFin]:
3150                if name not in  ['Type','Lam1','Lam2','Azimuth','Source','fltPath']:
3151                    ptlbls += '%12s' % (name)
3152                    ptstr += '%12.6f' % (Inst[name][1])
3153                    if instSig[name]:
3154                        refine = True
3155                        sigstr += '%12.6f' % (instSig[name])
3156                    else:
3157                        sigstr += 12*' '
3158            if refine:
3159                print >>pFile,'\n Instrument Parameters:'
3160                print >>pFile,ptlbls
3161                print >>pFile,ptstr
3162                print >>pFile,sigstr
3163            iBeg = iFin
3164            if iBeg == len(insKeys):
3165                Ok = False
3166       
3167    def PrintSampleParmsSig(Sample,sampleSig):
3168        ptlbls = ' names :'
3169        ptstr =  ' values:'
3170        sigstr = ' sig   :'
3171        refine = False
3172        if 'Bragg' in Sample['Type']:
3173            for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3174                ptlbls += '%14s'%(item)
3175                ptstr += '%14.4f'%(Sample[item][0])
3176                if sampleSig[i]:
3177                    refine = True
3178                    sigstr += '%14.4f'%(sampleSig[i])
3179                else:
3180                    sigstr += 14*' '
3181           
3182        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
3183            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3184                ptlbls += '%14s'%(item)
3185                ptstr += '%14.4f'%(Sample[item][0])
3186                if sampleSig[i]:
3187                    refine = True
3188                    sigstr += '%14.4f'%(sampleSig[i])
3189                else:
3190                    sigstr += 14*' '
3191
3192        if refine:
3193            print >>pFile,'\n Sample Parameters:'
3194            print >>pFile,ptlbls
3195            print >>pFile,ptstr
3196            print >>pFile,sigstr
3197       
3198    histoList = Histograms.keys()
3199    histoList.sort()
3200    for histogram in histoList:
3201        if 'PWDR' in histogram:
3202            Histogram = Histograms[histogram]
3203            hId = Histogram['hId']
3204            pfx = ':'+str(hId)+':'
3205            Background = Histogram['Background']
3206            backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
3207           
3208            Inst = Histogram['Instrument Parameters'][0]
3209            instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
3210       
3211            Sample = Histogram['Sample Parameters']
3212            sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
3213
3214            print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
3215            print >>pFile,135*'-'
3216            print >>pFile,' PWDR histogram weight factor = '+'%.3f'%(Histogram['wtFactor'])
3217            print >>pFile,' Final refinement wR = %.2f%% on %d observations in this histogram'%(Histogram['Residuals']['wR'],Histogram['Residuals']['Nobs'])
3218            print >>pFile,' Other residuals: R = %.2f%%, Rb = %.2f%%, wRb = %.2f%% wRmin = %.2f%%'% \
3219                (Histogram['Residuals']['R'],Histogram['Residuals']['Rb'],Histogram['Residuals']['wRb'],Histogram['Residuals']['wRmin'])
3220            if Print:
3221                print >>pFile,' Instrument type: ',Sample['Type']
3222                if FFtables != None and 'N' not in Inst['Type'][0]:
3223                    PrintFprime(FFtables,pfx,pFile)
3224                PrintSampleParmsSig(Sample,sampSig)
3225                PrintInstParmsSig(Inst,instSig)
3226                PrintBackgroundSig(Background,backSig)
3227               
Note: See TracBrowser for help on using the repository browser.