source: trunk/exports/G2cif.py @ 1014

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

background peak & Debye terms added to cif

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