source: trunk/exports/G2cif.py @ 1039

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

a few more cif fix ups
an extra _pd_block_diffractogram_id removed
difficulty with sph. harm. PO being a 1-d array fixed

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