source: trunk/exports/G2cif.py @ 1013

Last change on this file since 1013 was 1013, checked in by vondreele, 8 years ago

add output of profile coeff, background & size/mustrain/strain parms. to cif with esds

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