source: trunk/exports/G2cif.py @ 1057

Last change on this file since 1057 was 1057, checked in by toby, 8 years ago

reorg UserCalibrants? into ImageCalibrants?; Add trim (whitespace) routine; fix InstName? bug in GetHistogramPhaseData?; rework ValidatedTxtCtrl? for CIF; ScrolledMultiEditor? now resets after Cancel; CIF export progress

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