source: trunk/exports/G2cif.py @ 1045

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

fix Fourier map problems - didn't like Fo<0.
use ErrorDialog? for some errors.

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