source: trunk/exports/G2cif.py @ 1010

Last change on this file since 1010 was 1010, checked in by toby, 10 years ago

more CIF work
  • Property svn:eol-style set to native
File size: 55.1 KB
Line 
1'''Development code to export a GSAS-II project as a CIF
2The heavy lifting is done in method export
3'''
4
5# TODO: need a mechanism for editing of instrument names, bond pub flags, templates,...
6
7import datetime as dt
8import os.path
9import numpy as np
10import GSASIIIO as G2IO
11#reload(G2IO)
12import GSASIIgrid as G2gd
13import GSASIIstrIO as G2stIO
14#reload(G2stIO)
15#import GSASIImapvars as G2mv
16import GSASIImath as G2mth
17reload(G2mth)
18import GSASIIlattice as G2lat
19import GSASIIspc as G2spg
20#reload(G2spg)
21
22def getCallerDocString(): # for development
23    "Return the calling function's doc string"
24    import inspect as ins
25    for item in ins.stack()[1][0].f_code.co_consts:
26        if type(item) is str:
27            return item
28    else:
29        return '?'
30
31class ExportCIF(G2IO.ExportBaseclass):
32    def __init__(self,G2frame):
33        super(self.__class__,self).__init__( # fancy way to say <parentclass>.__init__
34            G2frame=G2frame,
35            formatName = 'full CIF',
36            longFormatName = 'Export project as complete CIF'
37            )
38        self.author = ''
39
40    def export(self,mode='full'):
41        '''Export a CIF
42
43        :param str mode: "full" (default) to create a complete CIF of project,
44          "simple" for a simple CIF with only coordinates
45        '''
46   
47        def WriteCIFitem(name,value=''):
48            if value:
49                if "\n" in value or len(value)> 70:
50                    if name.strip(): print name
51                    print '; '+value
52                    print '; '
53                elif " " in value:
54                    if len(name)+len(value) > 65:
55                        print name,'\n   ','"' + str(value) + '"'
56                    else:
57                        print name,'  ','"' + str(value) + '"'
58                else:
59                    if len(name)+len(value) > 65:
60                        print name,'\n   ',value
61                    else:
62                        print name,'  ',value
63            else:
64                print name
65
66        def WriteAudit():
67            WriteCIFitem('_audit_creation_method',
68                         'created in GSAS-II')
69            WriteCIFitem('_audit_creation_date',self.CIFdate)
70            if self.author:
71                WriteCIFitem('_audit_author_name',self.author)
72            WriteCIFitem('_audit_update_record',
73                         self.CIFdate+'  Initial software-generated CIF')
74
75        def WriteOverall():
76            '''Write out overall refinement information.
77
78            More could be done here, but this is a good start.
79            '''
80            WriteCIFitem('_pd_proc_info_datetime', self.CIFdate)
81            WriteCIFitem('_pd_calc_method', 'Rietveld Refinement')
82            #WriteCIFitem('_refine_ls_shift/su_max',DAT1)
83            #WriteCIFitem('_refine_ls_shift/su_mean',DAT2)
84            WriteCIFitem('_computing_structure_refinement','GSAS-II')
85            try:
86                vars = str(len(self.OverallParms['Covariance']['varyList']))
87            except:
88                vars = '?'
89            WriteCIFitem('_refine_ls_number_parameters',vars)
90            try:
91                GOF = G2mth.ValEsd(self.OverallParms['Covariance']['Rvals']['GOF'],-0.009)
92            except:
93                GOF = '?'
94            WriteCIFitem('_refine_ls_goodness_of_fit_all',GOF)
95
96            # get restraint info
97            # restraintDict = self.OverallParms.get('Restraints',{})
98            # for i in  self.OverallParms['Constraints']:
99            #     print i
100            #     for j in self.OverallParms['Constraints'][i]:
101            #         print j
102            #WriteCIFitem('_refine_ls_number_restraints',TEXT)
103           
104            # other things to consider reporting
105            # _refine_ls_number_reflns
106            # _refine_ls_goodness_of_fit_obs
107            # _refine_ls_R_factor_all
108            # _refine_ls_R_factor_obs
109            # _refine_ls_wR_factor_all
110            # _refine_ls_wR_factor_obs
111            # _refine_ls_restrained_S_all
112            # _refine_ls_restrained_S_obs
113
114            # include an overall profile r-factor, if there is more than one powder histogram
115            if len(self.powderDict) > 1:
116                WriteCIFitem('\n# OVERALL POWDER R-FACTOR')
117                try:
118                    R = str(self.OverallParms['Covariance']['Rvals']['Rwp'])
119                except:
120                    R = '?'
121                WriteCIFitem('_pd_proc_ls_prof_wR_factor',R)
122                #WriteCIFitem('_pd_proc_ls_prof_R_factor',TEXT(11:20)) # who cares!
123            WriteCIFitem('_refine_ls_matrix_type','full')
124            #WriteCIFitem('_refine_ls_matrix_type','userblocks')
125
126        def WritePubTemplate():
127            '''TODO: insert the publication template ``template_publ.cif`` or some modified
128            version for this project. Store this in the GPX file?
129            '''
130            print getCallerDocString()
131
132        def WritePhaseTemplate():
133            '''TODO: insert the phase template ``template_phase.cif`` or some modified
134            version for this project
135            '''
136            print getCallerDocString()
137
138        def WritePowderTemplate():
139            '''TODO: insert the phase template ``template_instrument.cif`` or some modified
140            version for this project
141            '''
142            print getCallerDocString()
143
144        def WriteSnglXtalTemplate():
145            '''TODO: insert the single-crystal histogram template
146            for this project
147            '''
148            print getCallerDocString()
149
150        def FormatSH(phasenam):
151            'Format a full spherical harmonics texture description as a string'
152            phasedict = self.Phases[phasenam] # pointer to current phase info           
153            pfx = str(phasedict['pId'])+'::'
154            s = ""
155            textureData = phasedict['General']['SH Texture']   
156            if textureData.get('Order'):
157                s += "Spherical Harmonics correction. Order = "+str(textureData['Order'])
158                s += " Model: " + str(textureData['Model']) + "\n    Orientation angles: "
159                for name in ['omega','chi','phi']:
160                    aname = pfx+'SH '+name
161                    s += name + " = "
162                    sig = self.sigDict.get(aname,-0.09)
163                    s += G2mth.ValEsd(self.parmDict[aname],sig)
164                    s += "; "
165                s += "\n"
166                s1 = "    Coefficients:  "
167                for name in textureData['SH Coeff'][1]:
168                    aname = pfx+name
169                    if len(s1) > 60:
170                        s += s1 + "\n"
171                        s1 = "    "
172                    s1 += aname + ' = '
173                    sig = self.sigDict.get(aname,-0.0009)
174                    s1 += G2mth.ValEsd(self.parmDict[aname],sig)
175                    s1 += "; "
176                s += s1
177            return s
178
179        def FormatHAPpo(phasenam):
180            '''return the March-Dollase/SH correction for every
181            histogram in the current phase formatted into a
182            character string
183            '''
184            phasedict = self.Phases[phasenam] # pointer to current phase info           
185            s = ''
186            for histogram in sorted(phasedict['Histograms']):
187                if histogram.startswith("HKLF"): continue # powder only
188                Histogram = self.Histograms.get(histogram)
189                if not Histogram: continue
190                hapData = phasedict['Histograms'][histogram]
191                if hapData['Pref.Ori.'][0] == 'MD':
192                    aname = str(phasedict['pId'])+':'+str(Histogram['hId'])+':MD'
193                    if self.parmDict.get(aname,1.0) != 1.0: continue
194                    sig = self.sigDict.get(aname,-0.009)
195                    if s != "": s += '\n'
196                    s += 'March-Dollase correction'
197                    if len(self.powderDict) > 1:
198                        s += ', histogram '+str(Histogram['hId']+1)
199                    s += ' coef. = ' + G2mth.ValEsd(self.parmDict[aname],sig)
200                    s += ' axis = ' + str(hapData['Pref.Ori.'][3])
201                else: # must be SH
202                    if s != "": s += '\n'
203                    s += 'Simple spherical harmonic correction'
204                    if len(self.powderDict) > 1:
205                        s += ', histogram '+str(Histogram['hId']+1)
206                    s += ' Order = '+str(hapData['Pref.Ori.'][4])+'\n'
207                    s1 = "    Coefficients:  "
208                    for item in hapData['Pref.Ori.'][5]:
209                        print item
210                        aname = str(phasedict['pId'])+':'+str(Histogram['hId'])+':'+item
211                        print aname
212                        if len(s1) > 60:
213                            s += s1 + "\n"
214                            s1 = "    "
215                        s1 += aname + ' = '
216                        sig = self.sigDict.get(aname,-0.0009)
217                        s1 += G2mth.ValEsd(self.parmDict[aname],sig)
218                        s1 += "; "
219                    s += s1
220            return s
221        def FormatBackground(bkg):
222            '''Display the Background information as a descriptive text string.
223           
224            TODO: this needs to be expanded to show the diffuse peak and
225            Debye term information as well.
226
227            :returns: the text description (str)
228            '''
229            fxn, bkgdict = bkg
230            terms = fxn[2]
231            txt = 'Background function: "'+fxn[0]+'" function with '+str(terms)+' terms:\n'
232            l = "   "
233            for v in fxn[3:]:
234                if len(l) > 60:
235                    txt += l + '\n'
236                    l = '   '
237                l += G2mth.ValEsd(v,-.009)+', '
238            txt += l
239            return txt
240
241        def FormatInstProfile(instparmdict):
242            '''Format the instrumental profile parameters with a
243            string description. Will only be called on PWDR histograms
244            '''
245            #print instparmdict[0].keys()
246            return 'TODO: Instrument profile goes here'
247
248        def FormatPhaseProfile(phasenam):
249            '''Format the phase-related profile parameters (size/strain)
250            with a string description.
251            return an empty string or None there are no
252            powder histograms for this phase.
253            '''
254            s = ''
255            phasedict = self.Phases[phasenam] # pointer to current phase info           
256            for histogram in sorted(phasedict['Histograms']):
257                if histogram.startswith("HKLF"): continue # powder only
258                Histogram = self.Histograms.get(histogram)
259                if not Histogram: continue
260                hapData = phasedict['Histograms'][histogram]
261            return 'TODO: Phase profile goes here'
262       
263        def FmtAtomType(sym):
264            'Reformat a GSAS-II atom type symbol to match CIF rules'
265            sym = sym.replace('_','') # underscores are not allowed: no isotope designation?
266            # in CIF, oxidation state sign symbols come after, not before
267            if '+' in sym:
268                sym = sym.replace('+','') + '+'
269            elif '-' in sym:
270                sym = sym.replace('-','') + '-'
271            return sym
272           
273        def PutInCol(val,wid):
274            '''Pad a value to >=wid+1 columns by adding spaces at the end. Always
275            adds at least one space
276            '''
277            val = str(val).replace(' ','')
278            if not val: val = '?'
279            fmt = '{:' + str(wid) + '} '
280            return fmt.format(val)
281
282        def MakeUniqueLabel(lbl,labellist):
283            'Make sure that every atom label is unique'
284            lbl = lbl.strip()
285            if not lbl: # deal with a blank label
286                lbl = 'A_1'
287            if lbl not in labellist:
288                labellist.append(lbl)
289                return lbl
290            i = 1
291            prefix = lbl
292            if '_' in lbl:
293                prefix = lbl[:lbl.rfind('_')]
294                suffix = lbl[lbl.rfind('_')+1:]
295                try:
296                    i = int(suffix)+1
297                except:
298                    pass
299            while prefix+'_'+str(i) in labellist:
300                i += 1
301            else:
302                lbl = prefix+'_'+str(i)
303                labellist.append(lbl)
304
305        def WriteAtomsNuclear(phasenam):
306            'Write atom positions to CIF'
307            phasedict = self.Phases[phasenam] # pointer to current phase info
308            General = phasedict['General']
309            cx,ct,cs,cia = General['AtomPtrs']
310            Atoms = phasedict['Atoms']
311            cfrac = cx+3
312            fpfx = str(phasedict['pId'])+'::Afrac:'       
313            for i,at in enumerate(Atoms):
314                fval = self.parmDict.get(fpfx+str(i),at[cfrac])
315                if fval != 0.0:
316                    break
317            else:
318                WriteCIFitem('\n# PHASE HAS NO ATOMS!')
319                return
320               
321            WriteCIFitem('\n# ATOMIC COORDINATES AND DISPLACEMENT PARAMETERS')
322            WriteCIFitem('loop_ '+
323                         '\n\t_atom_site_label'+
324                         '\n\t_atom_site_type_symbol'+
325                         '\n\t_atom_site_fract_x'+
326                         '\n\t_atom_site_fract_y'+
327                         '\n\t_atom_site_fract_z'+
328                         '\n\t_atom_site_occupancy'+
329                         '\n\t_atom_site_adp_type'+
330                         '\n\t_atom_site_U_iso_or_equiv'+
331                         '\n\t_atom_site_symmetry_multiplicity')
332
333            varnames = {cx:'Ax',cx+1:'Ay',cx+2:'Az',cx+3:'Afrac',
334                        cia+1:'AUiso',cia+2:'AU11',cia+3:'AU22',cia+4:'AU33',
335                        cia+5:'AU12',cia+6:'AU13',cia+7:'AU23'}
336            self.labellist = []
337           
338            pfx = str(phasedict['pId'])+'::'
339            # loop over all atoms
340            naniso = 0
341            for i,at in enumerate(Atoms):
342                s = PutInCol(MakeUniqueLabel(at[ct-1],self.labellist),6) # label
343                fval = self.parmDict.get(fpfx+str(i),at[cfrac])
344                if fval == 0.0: continue # ignore any atoms that have a occupancy set to 0 (exact)
345                s += PutInCol(FmtAtomType(at[ct]),4) # type
346                if at[cia] == 'I':
347                    adp = 'Uiso '
348                else:
349                    adp = 'Uani '
350                    naniso += 1
351                    # compute Uequiv crudely
352                    # correct: Defined as "1/3 trace of diagonalized U matrix".
353                    # SEE cell2GS & Uij2Ueqv to GSASIIlattice. Former is needed to make the GS matrix used by the latter.
354                    t = 0.0
355                    for j in (2,3,4):
356                        var = pfx+varnames[cia+j]+":"+str(i)
357                        t += self.parmDict.get(var,at[cia+j])
358                for j in (cx,cx+1,cx+2,cx+3,cia+1):
359                    if j in (cx,cx+1,cx+2):
360                        dig = 11
361                        sigdig = -0.00009
362                    else:
363                        dig = 10
364                        sigdig = -0.009
365                    var = pfx+varnames[j]+":"+str(i)
366                    dvar = pfx+"d"+varnames[j]+":"+str(i)
367                    if dvar not in self.sigDict:
368                        dvar = var
369                    if j == cia+1 and adp == 'Uani ':
370                        val = t/3.
371                        sig = sigdig
372                    else:
373                        #print var,(var in self.parmDict),(var in self.sigDict)
374                        val = self.parmDict.get(var,at[j])
375                        sig = self.sigDict.get(dvar,sigdig)
376                    s += PutInCol(G2mth.ValEsd(val,sig),dig)
377                s += adp
378                s += PutInCol(at[cs+1],3)
379                WriteCIFitem(s)
380            if naniso == 0: return
381            # now loop over aniso atoms
382            WriteCIFitem('\nloop_' + '\n\t_atom_site_aniso_label' + 
383                         '\n\t_atom_site_aniso_U_11' + '\n\t_atom_site_aniso_U_12' +
384                         '\n\t_atom_site_aniso_U_13' + '\n\t_atom_site_aniso_U_22' +
385                         '\n\t_atom_site_aniso_U_23' + '\n\t_atom_site_aniso_U_33')
386            for i,at in enumerate(Atoms):
387                fval = self.parmDict.get(fpfx+str(i),at[cfrac])
388                if fval == 0.0: continue # ignore any atoms that have a occupancy set to 0 (exact)
389                if at[cia] == 'I': continue
390                s = PutInCol(self.labellist[i],6) # label
391                for j in (2,3,4,5,6,7):
392                    sigdig = -0.0009
393                    var = pfx+varnames[cia+j]+":"+str(i)
394                    val = self.parmDict.get(var,at[cia+j])
395                    sig = self.sigDict.get(var,sigdig)
396                    s += PutInCol(G2mth.ValEsd(val,sig),11)
397                WriteCIFitem(s)
398
399        def HillSortElements(elmlist):
400            '''Sort elements in "Hill" order: C, H, others, (where others
401            are alphabetical).
402
403            :params list elmlist: a list of element strings
404
405            :returns: a sorted list of element strings
406            '''
407            newlist = []
408            oldlist = elmlist[:]
409            for elm in ('C','H'):
410                if elm in elmlist:
411                    newlist.append(elm)
412                    oldlist.pop(oldlist.index(elm))
413            return newlist+sorted(oldlist)
414
415        def WriteComposition(phasenam):
416            '''determine the composition for the unit cell, crudely determine Z and
417            then compute the composition in formula units
418            '''
419            phasedict = self.Phases[phasenam] # pointer to current phase info
420            General = phasedict['General']
421            Z = General.get('cellZ',0.0)
422            cx,ct,cs,cia = General['AtomPtrs']
423            Atoms = phasedict['Atoms']
424            fpfx = str(phasedict['pId'])+'::Afrac:'       
425            cfrac = cx+3
426            cmult = cs+1
427            compDict = {} # combines H,D & T
428            sitemultlist = []
429            massDict = dict(zip(General['AtomTypes'],General['AtomMass']))
430            cellmass = 0
431            for i,at in enumerate(Atoms):
432                atype = at[ct].strip()
433                if atype.find('-') != -1: atype = atype.split('-')[0]
434                if atype.find('+') != -1: atype = atype.split('+')[0]
435                atype = atype[0].upper()+atype[1:2].lower() # force case conversion
436                if atype == "D" or atype == "D": atype = "H"
437                fvar = fpfx+str(i)
438                fval = self.parmDict.get(fvar,at[cfrac])
439                mult = at[cmult]
440                if not massDict.get(at[ct]):
441                    print 'No mass found for atom type '+at[ct]
442                    print 'Will not compute cell contents for phase '+phasenam
443                    return
444                cellmass += massDict[at[ct]]*mult*fval
445                compDict[atype] = compDict.get(atype,0.0) + mult*fval
446                if fval == 1: sitemultlist.append(mult)
447            if len(compDict.keys()) == 0: return # no elements!
448            if Z < 1: # Z has not been computed or set by user
449                Z = 1
450                for i in range(2,min(sitemultlist)+1):
451                    for m in sitemultlist:
452                        if m % i != 0:
453                            break
454                        else:
455                            Z = i
456                General['cellZ'] = Z # save it
457
458            # when scattering factors are included in the CIF, this needs to be
459            # added to the loop here but only in the one-block case.
460            # For multiblock CIFs, scattering factors go in the histogram
461            # blocks  (for all atoms in all appropriate phases)
462
463            #if oneblock: # add scattering factors for current phase here
464            WriteCIFitem('\nloop_  _atom_type_symbol _atom_type_number_in_cell')
465            formula = ''
466            reload(G2mth)
467            for elem in HillSortElements(compDict.keys()):
468                WriteCIFitem('  ' + PutInCol(elem,4) +
469                             G2mth.ValEsd(compDict[elem],-0.009,True))
470                if formula: formula += " "
471                formula += elem
472                if compDict[elem] == Z: continue
473                formula += G2mth.ValEsd(compDict[elem]/Z,-0.009,True)
474            WriteCIFitem( '\n# Note that Z affects _cell_formula_sum and _weight')
475            WriteCIFitem( '_cell_formula_units_Z',str(Z))
476            WriteCIFitem( '_chemical_formula_sum',formula)
477            WriteCIFitem( '_chemical_formula_weight',
478                          G2mth.ValEsd(cellmass/Z,-0.09,True))
479
480        def WriteDistances(phasenam,SymOpList,offsetList,symOpList,G2oprList):
481            '''Report bond distances and angles for the CIF
482
483            Note that _geom_*_symmetry_* fields are values of form
484            n_klm where n is the symmetry operation in SymOpList (counted
485            starting with 1) and (k-5, l-5, m-5) are translations to add
486            to (x,y,z). See
487            http://www.iucr.org/__data/iucr/cifdic_html/1/cif_core.dic/Igeom_angle_site_symmetry_.html
488
489            TODO: need a method to select publication flags for distances/angles
490            '''
491            phasedict = self.Phases[phasenam] # pointer to current phase info           
492            Atoms = phasedict['Atoms']
493            cx,ct,cs,cia = phasedict['General']['AtomPtrs']
494            fpfx = str(phasedict['pId'])+'::Afrac:'       
495            cfrac = cx+3
496            # loop over interatomic distances for this phase
497            WriteCIFitem('\n# MOLECULAR GEOMETRY')
498            WriteCIFitem('loop_' + 
499                         '\n\t_geom_bond_atom_site_label_1' +
500                         '\n\t_geom_bond_atom_site_label_2' + 
501                         '\n\t_geom_bond_distance' + 
502                         '\n\t_geom_bond_site_symmetry_1' + 
503                         '\n\t_geom_bond_site_symmetry_2' + 
504                         '\n\t_geom_bond_publ_flag')
505
506            # Note that labels should be read from self.labellist to correspond
507            # to the values reported in the atoms table and zero occupancy atoms
508            # should not be included
509            fpfx = str(phasedict['pId'])+'::Afrac:'       
510            for i,at in enumerate(Atoms):
511                if self.parmDict.get(fpfx+str(i),at[cfrac]) == 0.0: continue
512                lbl = self.labellist[i]
513
514
515            # loop over interatomic angles for this phase
516            WriteCIFitem('loop_' + 
517                         '\n\t_geom_angle_atom_site_label_1' + 
518                         '\n\t_geom_angle_atom_site_label_2' + 
519                         '\n\t_geom_angle_atom_site_label_3' + 
520                         '\n\t_geom_angle' + 
521                         '\n\t_geom_angle_site_symmetry_1' +
522                         '\n\t_geom_angle_site_symmetry_2' + 
523                         '\n\t_geom_angle_site_symmetry_3' + 
524                         '\n\t_geom_angle_publ_flag')
525
526
527        def WritePhaseInfo(phasenam):
528            WriteCIFitem('\n# phase info for '+str(phasenam) + ' follows')
529            phasedict = self.Phases[phasenam] # pointer to current phase info           
530            WriteCIFitem('_pd_phase_name', phasenam)
531            pfx = str(phasedict['pId'])+'::'
532            A,sigA = G2stIO.cellFill(pfx,phasedict['General']['SGData'],self.parmDict,self.sigDict)
533            cellSig = G2stIO.getCellEsd(pfx,
534                                       phasedict['General']['SGData'],A,
535                                       self.OverallParms['Covariance'])  # returns 7 vals, includes sigVol
536            cellList = G2lat.A2cell(A) + (G2lat.calc_V(A),)
537            defsigL = 3*[-0.00001] + 3*[-0.001] + [-0.01] # significance to use when no sigma
538            names = ['length_a','length_b','length_c',
539                     'angle_alpha','angle_beta ','angle_gamma',
540                     'volume']
541            prevsig = 0
542            for lbl,defsig,val,sig in zip(names,defsigL,cellList,cellSig):
543                if sig:
544                    txt = G2mth.ValEsd(val,sig)
545                    prevsig = -sig # use this as the significance for next value
546                else:
547                    txt = G2mth.ValEsd(val,min(defsig,prevsig),True)
548                WriteCIFitem('_cell_'+lbl,txt)
549                   
550            WriteCIFitem('_symmetry_cell_setting',
551                         phasedict['General']['SGData']['SGSys'])
552
553            spacegroup = phasedict['General']['SGData']['SpGrp'].strip()
554            # regularize capitalization and remove trailing H/R
555            spacegroup = spacegroup[0].upper() + spacegroup[1:].lower().rstrip('rh ')
556            WriteCIFitem('_symmetry_space_group_name_H-M',spacegroup)
557
558            # generate symmetry operations including centering and center of symmetry
559            SymOpList,offsetList,symOpList,G2oprList = G2spg.AllOps(
560                phasedict['General']['SGData'])
561            WriteCIFitem('loop_ _space_group_symop_id _space_group_symop_operation_xyz')
562            for i,op in enumerate(SymOpList,start=1):
563                WriteCIFitem('   {:3d}  {:}'.format(i,op.lower()))
564
565            # loop over histogram(s) used in this phase
566            if not oneblock and not self.quickmode:
567                # report pointers to the histograms used in this phase
568                histlist = []
569                for hist in self.Phases[phasenam]['Histograms']:
570                    if self.Phases[phasenam]['Histograms'][hist]['Use']:
571                        if phasebyhistDict.get(hist):
572                            phasebyhistDict[hist].append(phasenam)
573                        else:
574                            phasebyhistDict[hist] = [phasenam,]
575                        blockid = datablockidDict.get(hist)
576                        if not blockid:
577                            print "Internal error: no block for data. Phase "+str(
578                                phasenam)+" histogram "+str(hist)
579                            histlist = []
580                            break
581                        histlist.append(blockid)
582
583                if len(histlist) == 0:
584                    WriteCIFitem('# Note: phase has no associated data')
585                else:
586                    WriteCIFitem('loop_  _pd_block_diffractogram_id')
587
588            # report atom params
589            if phasedict['General']['Type'] == 'nuclear':        #this needs macromolecular variant, etc!
590                WriteAtomsNuclear(phasenam)
591            else:
592                raise Exception,"no export for mm coordinates implemented"
593            # report cell contents
594            WriteComposition(phasenam)
595            if not self.quickmode:      # report distances and angles
596                WriteDistances(phasenam,SymOpList,offsetList,symOpList,G2oprList)
597
598        def WritePowderData(histlbl):
599            text = '?'
600            histblk = self.Histograms[histlbl]
601            inst = histblk['Instrument Parameters'][0]
602            hId = histblk['hId']
603            pfx = ':' + str(hId) + ':'
604            print 'TODO: powder here data for',histblk["Sample Parameters"]['InstrName']
605            # see wrpowdhist.for & wrreflist.for
606           
607            if 'Lam1' in inst:
608                ratio = self.parmDict.get('I(L2)/I(L1)',inst['I(L2)/I(L1)'][1])
609                sratio = self.sigDict.get('I(L2)/I(L1)',-0.0009)
610                lam1 = self.parmDict.get('Lam1',inst['Lam1'][1])
611                slam1 = self.sigDict.get('Lam1',-0.00009)
612                lam2 = self.parmDict.get('Lam2',inst['Lam2'][1])
613                slam2 = self.sigDict.get('Lam2',-0.00009)
614                # always assume Ka1 & Ka2 if two wavelengths are present
615                WriteCIFitem('loop_' + 
616                             '\n\t_diffrn_radiation_wavelength' +
617                             '\n\t_diffrn_radiation_wavelength_wt' + 
618                             '\n\t_diffrn_radiation_type' + 
619                             '\n\t_diffrn_radiation_wavelength_id')
620                WriteCIFitem('  ' + PutInCol(G2mth.ValEsd(lam1,slam1),15)+
621                             PutInCol('1.0',15) + 
622                             PutInCol('K\\a~1~',10) + 
623                             PutInCol('1',5))
624                WriteCIFitem('  ' + PutInCol(G2mth.ValEsd(lam2,slam2),15)+
625                             PutInCol(G2mth.ValEsd(ratio,sratio),15)+
626                             PutInCol('K\\a~2~',10) + 
627                             PutInCol('2',5))               
628            else:
629                lam1 = self.parmDict.get('Lam',inst['Lam'][1])
630                slam1 = self.sigDict.get('Lam',-0.00009)
631                WriteCIFitem('_diffrn_radiation_wavelength',G2mth.ValEsd(lam1,slam1))
632
633
634            if not oneblock:
635                if not phasebyhistDict.get(histlbl):
636                    WriteCIFitem('\n# No phases associated with this data set')
637                else:
638                    WriteCIFitem('\n# PHASE TABLE')
639                    WriteCIFitem('loop_' +
640                                 '\n\t_pd_phase_id' + 
641                                 '\n\t_pd_phase_block_id' + 
642                                 '\n\t_pd_phase_mass_%')
643                    wtFrSum = 0.
644                    for phasenam in phasebyhistDict.get(histlbl):
645                        hapData = self.Phases[phasenam]['Histograms'][histlbl]
646                        General = self.Phases[phasenam]['General']
647                        wtFrSum += hapData['Scale'][0]*General['Mass']
648
649                    for phasenam in phasebyhistDict.get(histlbl):
650                        hapData = self.Phases[phasenam]['Histograms'][histlbl]
651                        General = self.Phases[phasenam]['General']
652                        wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum
653                        pfx = str(self.Phases[phasenam]['pId'])+':'+str(hId)+':'
654                        if pfx+'Scale' in self.sigDict:
655                            sig = self.sigDict[pfx+'Scale']*wtFr/hapData['Scale'][0]
656                        else:
657                            sig = -0.0001
658                        WriteCIFitem(
659                            '  '+
660                            str(self.Phases[phasenam]['pId']) +
661                            '  '+datablockidDict[phasenam]+
662                            '  '+G2mth.ValEsd(wtFr,sig)
663                            )
664
665            # TODO: this will need help from Bob
666            # WriteCIFitem('_pd_proc_ls_prof_R_factor','?')
667            # WriteCIFitem('_pd_proc_ls_prof_wR_factor','?')
668            # WriteCIFitem('_pd_proc_ls_prof_wR_expected','?')
669            # WriteCIFitem('_refine_ls_R_Fsqd_factor','?')
670
671            if histblk['Instrument Parameters'][0]['Type'][1][1] == 'X':
672                WriteCIFitem('_diffrn_radiation_probe','x-ray')
673                pola = histblk['Instrument Parameters'][0].get('Polariz.')
674                if pola:
675                    pfx = ':' + str(hId) + ':'
676                    sig = self.sigDict.get(pfx+'Polariz.',-0.0009)
677                    txt = G2mth.ValEsd(pola[1],sig)
678                    WriteCIFitem('_diffrn_radiation_polarisn_ratio',txt)
679            elif histblk['Instrument Parameters'][0]['Type'][1][1] == 'N':
680                WriteCIFitem('_diffrn_radiation_probe','neutron')
681            # TOF (note that this may not be defined)
682            #if histblk['Instrument Parameters'][0]['Type'][1][2] == 'T':
683            #    WriteCIFitem('_pd_meas_2theta_fixed',text)
684           
685
686            # TODO: this will need help from Bob
687            #if not oneblock:
688            #WriteCIFitem('\n# SCATTERING FACTOR INFO')
689            #WriteCIFitem('loop_  _atom_type_symbol')
690            #if histblk['Instrument Parameters'][0]['Type'][1][1] == 'X':
691            #    WriteCIFitem('      _atom_type_scat_dispersion_real')
692            #    WriteCIFitem('      _atom_type_scat_dispersion_imag')
693            #    for lbl in ('a1','a2','a3', 'a4', 'b1', 'b2', 'b3', 'b4', 'c'):
694            #        WriteCIFitem('      _atom_type_scat_Cromer_Mann_'+lbl)
695            #elif histblk['Instrument Parameters'][0]['Type'][1][1] == 'N':
696            #    WriteCIFitem('      _atom_type_scat_length_neutron')
697            #WriteCIFitem('      _atom_type_scat_source')
698
699            WriteCIFitem('_pd_proc_ls_background_function',FormatBackground(histblk['Background']))
700
701            #WriteCIFitem('_exptl_absorpt_process_details','?')
702            #WriteCIFitem('_exptl_absorpt_correction_T_min','?')
703            #WriteCIFitem('_exptl_absorpt_correction_T_max','?')
704            #C extinction
705            #WRITE(IUCIF,'(A)') '# Extinction correction'
706            #CALL WRVAL(IUCIF,'_gsas_exptl_extinct_corr_T_min',TEXT(1:10))
707            #CALL WRVAL(IUCIF,'_gsas_exptl_extinct_corr_T_max',TEXT(11:20))
708
709            if not oneblock:                 # instrumental profile terms go here
710                WriteCIFitem('_pd_proc_ls_profile_function', 
711                             FormatInstProfile(histblk["Instrument Parameters"]))
712
713            #refprx = '_refln.' # mm
714            refprx = '_refln_' # normal
715            WriteCIFitem('\n# STRUCTURE FACTOR TABLE')           
716            # compute maximum intensity reflection
717            Imax = 0
718            for phasenam in histblk['Reflection Lists']:
719                scale = self.Phases[phasenam]['Histograms'][histlbl]['Scale'][0]
720                Icorr = np.array([refl[13] for refl in histblk['Reflection Lists'][phasenam]])
721                FO2 = np.array([refl[8] for refl in histblk['Reflection Lists'][phasenam]])
722                I100 = scale*FO2*Icorr
723                Imax = max(Imax,max(I100))
724
725            WriteCIFitem('loop_')
726            if len(histblk['Reflection Lists'].keys()) > 1:
727                WriteCIFitem('\t_pd_refln_phase_id')
728            WriteCIFitem('\t' + refprx + 'index_h' + 
729                         '\n\t' + refprx + 'index_k' + 
730                         '\n\t' + refprx + 'index_l' + 
731#                         '\n\t_pd_refln_wavelength_id' +
732#                         '\n\t' + refprx + 'status' +
733                         '\n\t' + refprx + 'F_squared_meas' + 
734#                         '\n\t' + refprx + 'F_squared_sigma' +
735                         '\n\t' + refprx + 'F_squared_calc' + 
736                         '\n\t' + refprx + 'phase_calc' + 
737                         '\n\t_pd_refln_d_spacing')
738            if Imax > 0:
739                WriteCIFitem('\t_gsas_i100_meas')
740
741            refcount = 0
742            hklmin = None
743            hklmax = None
744            dmax = None
745            dmin = None
746            for phasenam in histblk['Reflection Lists']:
747                scale = self.Phases[phasenam]['Histograms'][histlbl]['Scale'][0]
748                phaseid = self.Phases[phasenam]['pId']
749                refcount += len(histblk['Reflection Lists'][phasenam])
750                for ref in histblk['Reflection Lists'][phasenam]:
751                    if hklmin is None:
752                        hklmin = ref[0:3]
753                        hklmax = ref[0:3]
754                        dmax = dmin = ref[4]
755                    if len(histblk['Reflection Lists'].keys()) > 1:
756                        s = PutInCol(phaseid,2)
757                    else:
758                        s = ""
759                    for i,hkl in enumerate(ref[0:3]):
760                        hklmax[i] = max(hkl,hklmax[i])
761                        hklmin[i] = min(hkl,hklmin[i])
762                        s += PutInCol(int(hkl),4)
763                    for I in ref[8:10]:
764                        s += PutInCol(G2mth.ValEsd(I,-0.0009),10)
765                    s += PutInCol(G2mth.ValEsd(ref[10],-0.9),7)
766                    dmax = max(dmax,ref[4])
767                    dmin = min(dmin,ref[4])
768                    s += PutInCol(G2mth.ValEsd(ref[4],-0.009),8)
769                    if Imax > 0:
770                        I100 = 100.*scale*ref[8]*ref[13]/Imax
771                        s += PutInCol(G2mth.ValEsd(I100,-0.09),6)
772                    WriteCIFitem("  "+s)
773
774            WriteCIFitem('_reflns_number_total', str(refcount))
775            if hklmin is not None and len(histblk['Reflection Lists']) == 1: # hkl range has no meaning with multiple phases
776                WriteCIFitem('_reflns_limit_h_min', str(int(hklmin[0])))
777                WriteCIFitem('_reflns_limit_h_max', str(int(hklmax[0])))
778                WriteCIFitem('_reflns_limit_k_min', str(int(hklmin[1])))
779                WriteCIFitem('_reflns_limit_k_max', str(int(hklmax[1])))
780                WriteCIFitem('_reflns_limit_l_min', str(int(hklmin[2])))
781                WriteCIFitem('_reflns_limit_l_max', str(int(hklmax[2])))
782            if hklmin is not None:
783                WriteCIFitem('_reflns_d_resolution_high', G2mth.ValEsd(dmin,-0.009))
784                WriteCIFitem('_reflns_d_resolution_low', G2mth.ValEsd(dmax,-0.0009))
785
786            WriteCIFitem('\n# POWDER DATA TABLE')
787            # is data fixed step? If the step varies by <0.01% treat as fixed step
788            steps = histblk['Data'][0][1:] - histblk['Data'][0][:-1]
789            if abs(max(steps)-min(steps)) > abs(max(steps))/10000.:
790                fixedstep = False
791            else:
792                fixedstep = True
793
794            if fixedstep: # and not TOF
795                WriteCIFitem('_pd_meas_2theta_range_min', G2mth.ValEsd(histblk['Data'][0][0],-0.00009))
796                WriteCIFitem('_pd_meas_2theta_range_max', G2mth.ValEsd(histblk['Data'][0][-1],-0.00009))
797                WriteCIFitem('_pd_meas_2theta_range_inc', G2mth.ValEsd(steps.sum()/len(steps),-0.00009))
798                # zero correct, if defined
799                zero = None
800                zerolst = histblk['Instrument Parameters'][0].get('Zero')
801                if zerolst: zero = zerolst[1]
802                zero = self.parmDict.get('Zero',zero)
803                # TODO: Bob is zero added or subtracted?
804                if zero:
805                    WriteCIFitem('_pd_proc_2theta_range_min', G2mth.ValEsd(histblk['Data'][0][0]-zero,-0.00009))
806                    WriteCIFitem('_pd_proc_2theta_range_max', G2mth.ValEsd(histblk['Data'][0][-1]-zero,-0.00009))
807                    WriteCIFitem('_pd_proc_2theta_range_inc', G2mth.ValEsd(steps.sum()/len(steps),-0.00009))
808               
809            if zero:
810                WriteCIFitem('_pd_proc_number_of_points', str(len(histblk['Data'][0])))
811            else:
812                WriteCIFitem('_pd_meas_number_of_points', str(len(histblk['Data'][0])))
813            WriteCIFitem('\nloop_')
814            #            WriteCIFitem('\t_pd_proc_d_spacing') # need easy way to get this
815            if not fixedstep:
816                if zero:
817                    WriteCIFitem('\t_pd_proc_2theta_corrected')
818                else:
819                    WriteCIFitem('\t_pd_meas_2theta_scan')
820            # at least for now, always report weights. TODO: Should they be multiplied by weights?
821            #if countsdata:
822            #    WriteCIFitem('\t_pd_meas_counts_total')
823            #else:
824            WriteCIFitem('\t_pd_meas_intensity_total')
825            WriteCIFitem('\t_pd_calc_intensity_total')
826            WriteCIFitem('\t_pd_proc_intensity_bkg_calc')
827            WriteCIFitem('\t_pd_proc_ls_weight')
828            # TODO: are data excluded?
829            for x,yobs,yw,ycalc,ybkg in zip(histblk['Data'][0],
830                                            histblk['Data'][1],
831                                            histblk['Data'][2],
832                                            histblk['Data'][3],
833                                            histblk['Data'][4]):
834                if fixedstep:
835                    s = ""
836                else:
837                    s = PutInCol(G2mth.ValEsd(x-zero,-0.00009),10)
838                s += PutInCol(G2mth.ValEsd(yobs,-0.00009),14)
839                s += PutInCol(G2mth.ValEsd(ycalc,-0.00009),14)
840                s += PutInCol(G2mth.ValEsd(ybkg,-0.00009),14)
841                s += PutInCol(G2mth.ValEsd(yw,-0.000009),14)
842                WriteCIFitem("  "+s)
843                         
844        def WriteSingleXtalData(histlbl):
845            histblk = self.Histograms[histlbl]
846            print 'TODO: single xtal here data for',histblk["Instrument Parameters"][0]['InstrName']
847
848            #refprx = '_refln.' # mm
849            refprx = '_refln_' # normal
850
851            print histblk.keys()
852           
853            WriteCIFitem('\n# STRUCTURE FACTOR TABLE')           
854            WriteCIFitem('loop_' + 
855                         '\n\t' + refprx + 'index_h' + 
856                         '\n\t' + refprx + 'index_k' + 
857                         '\n\t' + refprx + 'index_l' +
858                         '\n\t' + refprx + 'F_squared_meas' + 
859                         '\n\t' + refprx + 'F_squared_sigma' + 
860                         '\n\t' + refprx + 'F_squared_calc' + 
861                         '\n\t' + refprx + 'phase_calc' +
862                         '\n\t_pd_refln_d_spacing'
863                         )
864
865            hklmin = None
866            hklmax = None
867            dmax = None
868            dmin = None
869            refcount = len(histblk['Data'])
870            for ref in histblk['Data']:
871                s = "  "
872                if hklmin is None:
873                    hklmin = ref[0:3]
874                    hklmax = ref[0:3]
875                    dmax = dmin = ref[4]
876                for i,hkl in enumerate(ref[0:3]):
877                    hklmax[i] = max(hkl,hklmax[i])
878                    hklmin[i] = min(hkl,hklmin[i])
879                    s += PutInCol(int(hkl),4)
880                sig = ref[6] * ref[8] / ref[5]
881                s += PutInCol(G2mth.ValEsd(ref[8],-abs(sig/10)),12)
882                s += PutInCol(G2mth.ValEsd(sig,-abs(sig)/10.),10)
883                s += PutInCol(G2mth.ValEsd(ref[9],-abs(sig/10)),12)
884                s += PutInCol(G2mth.ValEsd(ref[10],-0.9),7)
885                dmax = max(dmax,ref[4])
886                dmin = min(dmin,ref[4])
887                s += PutInCol(G2mth.ValEsd(ref[4],-0.009),8)
888                WriteCIFitem(s)
889            WriteCIFitem('_reflns_number_total', str(refcount))
890            if hklmin is not None:
891                WriteCIFitem('_reflns_limit_h_min', str(int(hklmin[0])))
892                WriteCIFitem('_reflns_limit_h_max', str(int(hklmax[0])))
893                WriteCIFitem('_reflns_limit_k_min', str(int(hklmin[1])))
894                WriteCIFitem('_reflns_limit_k_max', str(int(hklmax[1])))
895                WriteCIFitem('_reflns_limit_l_min', str(int(hklmin[2])))
896                WriteCIFitem('_reflns_limit_l_max', str(int(hklmax[2])))
897                WriteCIFitem('_reflns_d_resolution_high', G2mth.ValEsd(dmin,-0.009))
898                WriteCIFitem('_reflns_d_resolution_low', G2mth.ValEsd(dmax,-0.0009))
899
900        #============================================================
901        # the export process starts here
902        # also load all of the tree into a set of dicts
903        self.loadTree()
904        #self.dumpTree()
905        # create a dict with refined values and their uncertainties
906        self.loadParmDict()
907        #
908
909        # get restraint info
910        #restraintDict = self.OverallParms.get('Restraints',{})
911        #for i in  self.OverallParms['Constraints']:
912        #    print i
913        #    for j in self.OverallParms['Constraints'][i]:
914        #        print j
915        #return
916
917        self.CIFdate = dt.datetime.strftime(dt.datetime.now(),"%Y-%m-%dT%H:%M")
918        # count phases, powder and single crystal histograms
919        self.nphase = len(self.Phases)
920        self.powderDict = {}
921        self.xtalDict = {}
922        for hist in self.Histograms:
923            i = self.Histograms[hist]['hId']
924            if hist.startswith("PWDR"): 
925                self.powderDict[i] = hist
926            elif hist.startswith("HKLF"): 
927                self.xtalDict[i] = hist
928        # is there anything to export?
929        if self.nphase + len(self.powderDict) + len(self.xtalDict) == 0: 
930            self.G2frame.ErrorDialog(
931                'Empty project',
932                'No data or phases to include in CIF')
933            return
934        # is there a file name defined?
935        self.CIFname = os.path.splitext(
936            os.path.split(self.G2frame.GSASprojectfile)[1]
937            )[0]
938        self.CIFname = self.CIFname.replace(' ','')
939        if not self.CIFname:
940            self.G2frame.ErrorDialog(
941                'No GPX name',
942                'Please save the project to provide a name')
943            return
944        # test for quick CIF mode or no data
945        self.quickmode = False
946        phasenam = phasenum = None # include all phases
947        if mode != "full" or len(self.powderDict) + len(self.xtalDict) == 0:
948            self.quickmode = True
949            oneblock = True
950            if self.nphase == 0:
951                self.G2frame.ErrorDialog(
952                    'No phase present',
953                    'Cannot create a coordinates CIF with no phases')
954                return
955            elif self.nphase > 1: # quick mode: choose one phase
956                choices = sorted(self.Phases.keys())
957                phasenum = G2gd.ItemSelector(choices,self.G2frame)
958                if phasenum is None: return
959                phasenam = choices[phasenum]
960        # will this require a multiblock CIF?
961        elif self.nphase > 1:
962            oneblock = False
963        elif len(self.powderDict) + len(self.xtalDict) > 1:
964            oneblock = False
965        else: # one phase, one dataset, Full CIF
966            oneblock = True
967
968        # make sure needed infomation is present
969        # get CIF author name -- required for full CIFs
970        try:
971            self.author = self.OverallParms['Controls'].get("Author",'').strip()
972        except KeyError:
973            pass
974        while not (self.author or self.quickmode):
975            dlg = G2gd.SingleStringDialog(self.G2frame,'Get CIF Author','Provide CIF Author name (Last, First)')
976            if not dlg.Show(): return # cancel was pressed
977            self.author = dlg.GetValue()
978            dlg.Destroy()
979        try:
980            self.OverallParms['Controls']["Author"] = self.author # save for future
981        except KeyError:
982            pass
983        self.shortauthorname = self.author.replace(',','').replace(' ','')[:20]
984
985        # check the instrument name for every histogram
986        if not self.quickmode:
987            dictlist = []
988            keylist = []
989            lbllist = []
990            invalid = 0
991            key3 = 'InstrName'
992            for hist in self.Histograms:
993                if hist.startswith("PWDR"): 
994                    key2 = "Sample Parameters"
995                    d = self.Histograms[hist][key2]
996                elif hist.startswith("HKLF"): 
997                    key2 = "Instrument Parameters"
998                    d = self.Histograms[hist][key2][0]
999                   
1000                lbllist.append(hist)
1001                dictlist.append(d)
1002                keylist.append(key3)
1003                instrname = d.get(key3)
1004                if instrname is None:
1005                    d[key3] = ''
1006                    invalid += 1
1007                elif instrname.strip() == '':
1008                    invalid += 1
1009            if invalid:
1010                msg = ""
1011                if invalid > 3: msg = (
1012                    "\n\nNote: it may be faster to set the name for\n"
1013                    "one histogram for each instrument and use the\n"
1014                    "File/Copy option to duplicate the name"
1015                    )
1016                if not G2gd.CallScrolledMultiEditor(
1017                    self.G2frame,dictlist,keylist,
1018                    prelbl=range(1,len(dictlist)+1),
1019                    postlbl=lbllist,
1020                    title='Instrument names',
1021                    header="Edit instrument names. Note that a non-blank\nname is required for all histograms"+msg,
1022                    ): return
1023
1024        if oneblock and not self.quickmode:
1025            # select a dataset to use (there should only be one set in one block,
1026            # but take whatever comes 1st
1027            for hist in self.Histograms:
1028                histblk = self.Histograms[hist]
1029                if hist.startswith("PWDR"): 
1030                    instnam = histblk["Sample Parameters"]['InstrName']
1031                    break # ignore all but 1st data histogram
1032                elif hist.startswith("HKLF"): 
1033                    instnam = histblk["Instrument Parameters"][0]['InstrName']
1034                    break # ignore all but 1st data histogram
1035        #======================================================================
1036        # Start writing the CIF - single block
1037        #======================================================================
1038        if oneblock:
1039            WriteCIFitem('data_'+self.CIFname)
1040            if phasenam is None: # if not already selected, select the first phase (should be one)
1041                phasenam = self.Phases.keys()[0]
1042            #print 'phasenam',phasenam
1043            phaseblk = self.Phases[phasenam] # pointer to current phase info
1044            if not self.quickmode:
1045                instnam = instnam.replace(' ','')
1046                WriteCIFitem('_pd_block_id',
1047                             str(self.CIFdate) + "|" + str(self.CIFname) + "|" +
1048                             str(self.shortauthorname) + "|" + instnam)
1049                WriteAudit()
1050                WritePubTemplate()
1051                WriteOverall()
1052                WritePhaseTemplate()
1053            # report the phase info
1054            WritePhaseInfo(phasenam)
1055            if hist.startswith("PWDR") and not self.quickmode:
1056                # preferred orientation
1057                SH = FormatSH(phasenam)
1058                MD = FormatHAPpo(phasenam)
1059                if SH and MD:
1060                    WriteCIFitem('_pd_proc_ls_pref_orient_corr', SH + '\n' + MD)
1061                elif SH or MD:
1062                    WriteCIFitem('_pd_proc_ls_pref_orient_corr', SH + MD)
1063                else:
1064                    WriteCIFitem('_pd_proc_ls_pref_orient_corr', 'none')
1065                    # report profile, since one-block: include both histogram and phase info
1066                WriteCIFitem('_pd_proc_ls_profile_function',
1067                             FormatInstProfile(histblk["Instrument Parameters"])
1068                             + '\n' +
1069                             FormatPhaseProfile(phasenam))
1070                WritePowderTemplate()
1071                WritePowderData(hist)
1072            elif hist.startswith("HKLF") and not self.quickmode:
1073                WriteSnglXtalTemplate()
1074                WriteSingleXtalData(hist)
1075        else:
1076        #======================================================================
1077        # Start writing the CIF - multiblock
1078        #======================================================================
1079            # publication info
1080            WriteCIFitem('\ndata_'+self.CIFname+'_publ')
1081            WriteAudit()
1082            WriteCIFitem('_pd_block_id',
1083                         str(self.CIFdate) + "|" + str(self.CIFname) + "|" +
1084                         str(self.shortauthorname) + "|Overall")
1085            WritePubTemplate()
1086            # overall info
1087            WriteCIFitem('data_'+str(self.CIFname)+'_overall')
1088            WriteOverall()
1089            #============================================================
1090            WriteCIFitem('# POINTERS TO PHASE AND HISTOGRAM BLOCKS')
1091            datablockidDict = {} # save block names here -- N.B. check for conflicts between phase & hist names (unlikely!)
1092            # loop over phase blocks
1093            if self.nphase > 1:
1094                loopprefix = ''
1095                WriteCIFitem('loop_   _pd_phase_block_id')
1096            else:
1097                loopprefix = '_pd_phase_block_id'
1098           
1099            for phasenam in sorted(self.Phases.keys()):
1100                i = self.Phases[phasenam]['pId']
1101                datablockidDict[phasenam] = (str(self.CIFdate) + "|" + str(self.CIFname) + "|" +
1102                             'phase_'+ str(i) + '|' + str(self.shortauthorname))
1103                WriteCIFitem(loopprefix,datablockidDict[phasenam])
1104            # loop over data blocks
1105            if len(self.powderDict) + len(self.xtalDict) > 1:
1106                loopprefix = ''
1107                WriteCIFitem('loop_   _pd_block_diffractogram_id')
1108            else:
1109                loopprefix = '_pd_block_diffractogram_id'
1110            for i in sorted(self.powderDict.keys()):
1111                hist = self.powderDict[i]
1112                histblk = self.Histograms[hist]
1113                instnam = histblk["Sample Parameters"]['InstrName']
1114                instnam = instnam.replace(' ','')
1115                i = histblk['hId']
1116                datablockidDict[hist] = (str(self.CIFdate) + "|" + str(self.CIFname) + "|" +
1117                                         str(self.shortauthorname) + "|" +
1118                                         instnam + "_hist_"+str(i))
1119                WriteCIFitem(loopprefix,datablockidDict[hist])
1120            for i in sorted(self.xtalDict.keys()):
1121                hist = self.xtalDict[i]
1122                histblk = self.Histograms[hist]
1123                instnam = histblk["Instrument Parameters"][0]['InstrName']
1124                instnam = instnam.replace(' ','')
1125                i = histblk['hId']
1126                datablockidDict[hist] = (str(self.CIFdate) + "|" + str(self.CIFname) + "|" +
1127                                         str(self.shortauthorname) + "|" +
1128                                         instnam + "_hist_"+str(i))
1129                WriteCIFitem(loopprefix,datablockidDict[hist])
1130            #============================================================
1131            # loop over phases, exporting them
1132            phasebyhistDict = {} # create a cross-reference to phases by histogram
1133            for j,phasenam in enumerate(sorted(self.Phases.keys())):
1134                i = self.Phases[phasenam]['pId']
1135                WriteCIFitem('\ndata_'+self.CIFname+"_phase_"+str(i))
1136                print "debug, processing ",phasenam
1137                WriteCIFitem('# Information for phase '+str(i))
1138                WriteCIFitem('_pd_block_id',datablockidDict[phasenam])
1139                # report the phase
1140                WritePhaseTemplate()
1141                WritePhaseInfo(phasenam)
1142                # preferred orientation
1143                SH = FormatSH(phasenam)
1144                MD = FormatHAPpo(phasenam)
1145                if SH and MD:
1146                    WriteCIFitem('_pd_proc_ls_pref_orient_corr', SH + '\n' + MD)
1147                elif SH or MD:
1148                    WriteCIFitem('_pd_proc_ls_pref_orient_corr', SH + MD)
1149                else:
1150                    WriteCIFitem('_pd_proc_ls_pref_orient_corr', 'none')
1151                # report sample profile terms
1152                PP = FormatPhaseProfile(phasenam)
1153                if PP:
1154                    WriteCIFitem('_pd_proc_ls_profile_function',PP)
1155                   
1156            #============================================================
1157            # loop over histograms, exporting them
1158            for i in sorted(self.powderDict.keys()):
1159                hist = self.powderDict[i]
1160                histblk = self.Histograms[hist]
1161                if hist.startswith("PWDR"): 
1162                    WriteCIFitem('\ndata_'+self.CIFname+"_pwd_"+str(i))
1163                    #instnam = histblk["Sample Parameters"]['InstrName']
1164                    # report instrumental profile terms
1165                    WriteCIFitem('_pd_proc_ls_profile_function',
1166                                 FormatInstProfile(histblk["Instrument Parameters"]))
1167                    WriteCIFitem('# Information for histogram '+str(i)+': '+
1168                                 hist)
1169                    WriteCIFitem('_pd_block_id',datablockidDict[hist])
1170                    WritePowderTemplate()
1171                    WritePowderData(hist)
1172            for i in sorted(self.xtalDict.keys()):
1173                hist = self.xtalDict[i]
1174                histblk = self.Histograms[hist]
1175                if hist.startswith("HKLF"): 
1176                    WriteCIFitem('\ndata_'+self.CIFname+"_sx_"+str(i))
1177                    #instnam = histblk["Instrument Parameters"][0]['InstrName']
1178                    WriteCIFitem('# Information for histogram '+str(i)+': '+
1179                                 hist)
1180                    WriteCIFitem('_pd_block_id',datablockidDict[hist])
1181                    WriteSnglXtalTemplate()
1182                    WriteSingleXtalData(hist)
1183
1184        WriteCIFitem('#--' + 15*'eof--' + '#')
1185
Note: See TracBrowser for help on using the repository browser.