source: trunk/exports/G2cif.py @ 1026

Last change on this file since 1026 was 1026, checked in by toby, 9 years ago

now with distances & angles

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