source: trunk/exports/G2cif.py @ 963

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

more CIF work, more docs

  • Property svn:eol-style set to native
File size: 45.5 KB
Line 
1'''Development code to export a GSAS-II project as a CIF
2The heavy lifting is done in method export
3'''
4
5import datetime as dt
6import os.path
7import GSASIIIO as G2IO
8#reload(G2IO)
9import GSASIIgrid as G2gd
10import GSASIIstrIO as G2stIO
11#import GSASIImapvars as G2mv
12import GSASIImath as G2mth
13import GSASIIlattice as G2lat
14import GSASIIspc as G2spg
15reload(G2spg)
16
17def getCallerDocString(): # for development
18    "Return the calling function's doc string"
19    import inspect as ins
20    for item in ins.stack()[1][0].f_code.co_consts:
21        if type(item) is str:
22            return item
23    else:
24        return '?'
25
26class ExportCIF(G2IO.ExportBaseclass):
27    def __init__(self,G2frame):
28        super(self.__class__,self).__init__( # fancy way to say <parentclass>.__init__
29            G2frame=G2frame,
30            formatName = 'full CIF',
31            longFormatName = 'Export project as complete CIF'
32            )
33        self.author = ''
34
35    def export(self,mode='full'):
36        '''Export a CIF
37
38        :param str mode: "full" (default) to create a complete CIF of project,
39          "simple" for a simple CIF with only coordinates
40        '''
41        def WriteCIFitem(name,value=''):
42            if value:
43                if "\n" in value or len(value)> 70:
44                    if name.strip(): print name
45                    print '; '+value
46                    print '; '
47                elif " " in value:
48                    if len(name)+len(value) > 65:
49                        print name,'\n   ','"' + str(value) + '"'
50                    else:
51                        print name,'  ','"' + str(value) + '"'
52                else:
53                    if len(name)+len(value) > 65:
54                        print name,'\n   ',value
55                    else:
56                        print name,'  ',value
57            else:
58                print name
59
60        def WriteAudit():
61            WriteCIFitem('_audit_creation_method',
62                         'created in GSAS-II')
63            WriteCIFitem('_audit_creation_date',self.CIFdate)
64            if self.author:
65                WriteCIFitem('_audit_author_name',self.author)
66            WriteCIFitem('_audit_update_record',
67                         self.CIFdate+'  Initial software-generated CIF')
68
69        def WriteOverall():
70            '''Write out overall refinement information.
71
72            More could be done here, but this is a good start.
73            '''
74            WriteCIFitem('_pd_proc_info_datetime', self.CIFdate)
75            WriteCIFitem('_pd_calc_method', 'Rietveld Refinement')
76            #WriteCIFitem('_refine_ls_shift/su_max',DAT1)
77            #WriteCIFitem('_refine_ls_shift/su_mean',DAT2)
78            WriteCIFitem('_computing_structure_refinement','GSAS-II')
79            try:
80                vars = str(len(self.OverallParms['Covariance']['varyList']))
81            except:
82                vars = '?'
83            WriteCIFitem('_refine_ls_number_parameters',vars)
84            try:
85                GOF = G2mth.ValEsd(self.OverallParms['Covariance']['Rvals']['GOF'],-0.009)
86            except:
87                GOF = '?'
88            WriteCIFitem('_refine_ls_goodness_of_fit_all',GOF)
89
90            # get restraint info
91            # restraintDict = self.OverallParms.get('Restraints',{})
92            # for i in  self.OverallParms['Constraints']:
93            #     print i
94            #     for j in self.OverallParms['Constraints'][i]:
95            #         print j
96            #WriteCIFitem('_refine_ls_number_restraints',TEXT)
97           
98            # other things to consider reporting
99            # _refine_ls_number_reflns
100            # _refine_ls_goodness_of_fit_obs
101            # _refine_ls_R_factor_all
102            # _refine_ls_R_factor_obs
103            # _refine_ls_wR_factor_all
104            # _refine_ls_wR_factor_obs
105            # _refine_ls_restrained_S_all
106            # _refine_ls_restrained_S_obs
107
108            # include an overall profile r-factor, if there is more than one powder histogram
109            if self.npowder > 1:
110                WriteCIFitem('\n# OVERALL POWDER R-FACTOR')
111                try:
112                    R = str(self.OverallParms['Covariance']['Rvals']['Rwp'])
113                except:
114                    R = '?'
115                WriteCIFitem('_pd_proc_ls_prof_wR_factor',R)
116                #WriteCIFitem('_pd_proc_ls_prof_R_factor',TEXT(11:20)) # who cares!
117            WriteCIFitem('_refine_ls_matrix_type','full')
118            #WriteCIFitem('_refine_ls_matrix_type','userblocks')
119
120        def WritePubTemplate():
121            '''TODO: insert the publication template ``template_publ.cif`` or some modified
122            version for this project. Store this in the GPX file?
123            '''
124            print getCallerDocString()
125
126        def WritePhaseTemplate():
127            '''TODO: insert the phase template ``template_phase.cif`` or some modified
128            version for this project
129            '''
130            print getCallerDocString()
131
132        def WritePowderTemplate():
133            '''TODO: insert the phase template ``template_instrument.cif`` or some modified
134            version for this project
135            '''
136            print getCallerDocString()
137
138        def WriteSnglXtalTemplate():
139            '''TODO: insert the single-crystal histogram template
140            for this project
141            '''
142            print getCallerDocString()
143
144        def FormatSH(phasenam):
145            'Format a full spherical harmonics texture description as a string'
146            phasedict = self.Phases[phasenam] # pointer to current phase info           
147            pfx = str(phasedict['pId'])+'::'
148            s = ""
149            textureData = phasedict['General']['SH Texture']   
150            if textureData.get('Order'):
151                s += "Spherical Harmonics correction. Order = "+str(textureData['Order'])
152                s += " Model: " + str(textureData['Model']) + "\n    Orientation angles: "
153                for name in ['omega','chi','phi']:
154                    aname = pfx+'SH '+name
155                    s += name + " = "
156                    sig = self.sigDict.get(aname,-0.09)
157                    s += G2mth.ValEsd(self.parmDict[aname],sig)
158                    s += "; "
159                s += "\n"
160                s1 = "    Coefficients:  "
161                for name in textureData['SH Coeff'][1]:
162                    aname = pfx+name
163                    if len(s1) > 60:
164                        s += s1 + "\n"
165                        s1 = "    "
166                    s1 += aname + ' = '
167                    sig = self.sigDict.get(aname,-0.0009)
168                    s1 += G2mth.ValEsd(self.parmDict[aname],sig)
169                    s1 += "; "
170                s += s1
171            return s
172
173        def FormatHAPpo(phasenam):
174            '''return the March-Dollase/SH correction for every
175            histogram in the current phase formatted into a
176            character string
177            '''
178            phasedict = self.Phases[phasenam] # pointer to current phase info           
179            s = ''
180            for histogram in sorted(phasedict['Histograms']):
181                if histogram.startswith("HKLF"): continue # powder only
182                Histogram = self.Histograms.get(histogram)
183                if not Histogram: continue
184                hapData = phasedict['Histograms'][histogram]
185                if hapData['Pref.Ori.'][0] == 'MD':
186                    aname = str(phasedict['pId'])+':'+str(Histogram['hId'])+':MD'
187                    if self.parmDict.get(aname,1.0) != 1.0: continue
188                    sig = self.sigDict.get(aname,-0.009)
189                    if s != "": s += '\n'
190                    s += 'March-Dollase correction'
191                    if self.npowder > 1:
192                        s += ', histogram '+str(Histogram['hId']+1)
193                    s += ' coef. = ' + G2mth.ValEsd(self.parmDict[aname],sig)
194                    s += ' axis = ' + str(hapData['Pref.Ori.'][3])
195                else: # must be SH
196                    if s != "": s += '\n'
197                    s += 'Simple spherical harmonic correction'
198                    if self.npowder > 1:
199                        s += ', histogram '+str(Histogram['hId']+1)
200                    s += ' Order = '+str(hapData['Pref.Ori.'][4])+'\n'
201                    s1 = "    Coefficients:  "
202                    for item in hapData['Pref.Ori.'][5]:
203                        print item
204                        aname = str(phasedict['pId'])+':'+str(Histogram['hId'])+':'+item
205                        print aname
206                        if len(s1) > 60:
207                            s += s1 + "\n"
208                            s1 = "    "
209                        s1 += aname + ' = '
210                        sig = self.sigDict.get(aname,-0.0009)
211                        s1 += G2mth.ValEsd(self.parmDict[aname],sig)
212                        s1 += "; "
213                    s += s1
214            return s
215
216        def FormatInstProfile(instparmdict):
217            '''Format the instrumental profile parameters with a
218            string description. Will only be called on PWDR histograms
219            '''
220            #print instparmdict[0].keys()
221            return 'TODO: Instrument profile goes here'
222
223        def FormatPhaseProfile(phasenam):
224            '''Format the phase-related profile parameters (size/strain)
225            with a string description.
226            return an empty string or None there are no
227            powder histograms for this phase.
228            '''
229            s = ''
230            phasedict = self.Phases[phasenam] # pointer to current phase info           
231            for histogram in sorted(phasedict['Histograms']):
232                if histogram.startswith("HKLF"): continue # powder only
233                Histogram = self.Histograms.get(histogram)
234                if not Histogram: continue
235                hapData = phasedict['Histograms'][histogram]
236            return 'TODO: Phase profile goes here'
237       
238        def FmtAtomType(sym):
239            'Reformat a GSAS-II atom type symbol to match CIF rules'
240            sym = sym.replace('_','') # underscores are not allowed: no isotope designation?
241            # in CIF, oxidation state sign symbols come after, not before
242            if '+' in sym:
243                sym = sym.replace('+','') + '+'
244            elif '-' in sym:
245                sym = sym.replace('-','') + '-'
246            return sym
247           
248        def PutInCol(val,wid):
249            '''Pad a value to >=wid+1 columns by adding spaces at the end. Always
250            adds at least one space
251            '''
252            val = str(val).replace(' ','')
253            if not val: val = '?'
254            fmt = '{:' + str(wid) + '} '
255            return fmt.format(val)
256
257        def MakeUniqueLabel(lbl,labellist):
258            'Make sure that every atom label is unique'
259            lbl = lbl.strip()
260            if not lbl: # deal with a blank label
261                lbl = 'A_1'
262            if lbl not in labellist:
263                labellist.append(lbl)
264                return lbl
265            i = 1
266            prefix = lbl
267            if '_' in lbl:
268                prefix = lbl[:lbl.rfind('_')]
269                suffix = lbl[lbl.rfind('_')+1:]
270                try:
271                    i = int(suffix)+1
272                except:
273                    pass
274            while prefix+'_'+str(i) in labellist:
275                i += 1
276            else:
277                lbl = prefix+'_'+str(i)
278                labellist.append(lbl)
279
280        def WriteAtomsNuclear(phasenam):
281            'Write atom positions to CIF'
282            phasedict = self.Phases[phasenam] # pointer to current phase info
283            General = phasedict['General']
284            cx,ct,cs,cia = General['AtomPtrs']
285            Atoms = phasedict['Atoms']
286            cfrac = cx+3
287            fpfx = str(phasedict['pId'])+'::Afrac:'       
288            for i,at in enumerate(Atoms):
289                fval = self.parmDict.get(fpfx+str(i),at[cfrac])
290                if fval != 0.0:
291                    break
292            else:
293                WriteCIFitem('\n# PHASE HAS NO ATOMS!')
294                return
295               
296            WriteCIFitem('\n# ATOMIC COORDINATES AND DISPLACEMENT PARAMETERS')
297            WriteCIFitem('loop_ '+
298                         '\n\t_atom_site_label'+
299                         '\n\t_atom_site_type_symbol'+
300                         '\n\t_atom_site_fract_x'+
301                         '\n\t_atom_site_fract_y'+
302                         '\n\t_atom_site_fract_z'+
303                         '\n\t_atom_site_occupancy'+
304                         '\n\t_atom_site_adp_type'+
305                         '\n\t_atom_site_U_iso_or_equiv'+
306                         '\n\t_atom_site_symmetry_multiplicity')
307
308            varnames = {cx:'Ax',cx+1:'Ay',cx+2:'Az',cx+3:'Afrac',
309                        cia+1:'AUiso',cia+2:'AU11',cia+3:'AU22',cia+4:'AU33',
310                        cia+5:'AU12',cia+6:'AU13',cia+7:'AU23'}
311            self.labellist = []
312           
313            pfx = str(phasedict['pId'])+'::'
314            # loop over all atoms
315            naniso = 0
316            for i,at in enumerate(Atoms):
317                s = PutInCol(MakeUniqueLabel(at[ct-1],self.labellist),6) # label
318                fval = self.parmDict.get(fpfx+str(i),at[cfrac])
319                if fval == 0.0: continue # ignore any atoms that have a occupancy set to 0 (exact)
320                s += PutInCol(FmtAtomType(at[ct]),4) # type
321                if at[cia] == 'I':
322                    adp = 'Uiso '
323                else:
324                    adp = 'Uani '
325                    naniso += 1
326                    # compute Uequiv crudely
327                    # correct: Defined as "1/3 trace of diagonalized U matrix".
328                    # SEE cell2GS & Uij2Ueqv to GSASIIlattice. Former is needed to make the GS matrix used by the latter.
329                    t = 0.0
330                    for j in (2,3,4):
331                        var = pfx+varnames[cia+j]+":"+str(i)
332                        t += self.parmDict.get(var,at[cia+j])
333                for j in (cx,cx+1,cx+2,cx+3,cia+1):
334                    if j in (cx,cx+1,cx+2):
335                        dig = 11
336                        sigdig = -0.00009
337                    else:
338                        dig = 10
339                        sigdig = -0.009
340                    var = pfx+varnames[j]+":"+str(i)
341                    dvar = pfx+"d"+varnames[j]+":"+str(i)
342                    if dvar not in self.sigDict:
343                        dvar = var
344                    if j == cia+1 and adp == 'Uani ':
345                        val = t/3.
346                        sig = sigdig
347                    else:
348                        #print var,(var in self.parmDict),(var in self.sigDict)
349                        val = self.parmDict.get(var,at[j])
350                        sig = self.sigDict.get(dvar,sigdig)
351                    s += PutInCol(G2mth.ValEsd(val,sig),dig)
352                s += adp
353                s += PutInCol(at[cs+1],3)
354                WriteCIFitem(s)
355            if naniso == 0: return
356            # now loop over aniso atoms
357            WriteCIFitem('\nloop_' + '\n\t_atom_site_aniso_label' + 
358                         '\n\t_atom_site_aniso_U_11' + '\n\t_atom_site_aniso_U_12' +
359                         '\n\t_atom_site_aniso_U_13' + '\n\t_atom_site_aniso_U_22' +
360                         '\n\t_atom_site_aniso_U_23' + '\n\t_atom_site_aniso_U_33')
361            for i,at in enumerate(Atoms):
362                fval = self.parmDict.get(fpfx+str(i),at[cfrac])
363                if fval == 0.0: continue # ignore any atoms that have a occupancy set to 0 (exact)
364                if at[cia] == 'I': continue
365                s = PutInCol(self.labellist[i],6) # label
366                for j in (2,3,4,5,6,7):
367                    sigdig = -0.0009
368                    var = pfx+varnames[cia+j]+":"+str(i)
369                    val = self.parmDict.get(var,at[cia+j])
370                    sig = self.sigDict.get(var,sigdig)
371                    s += PutInCol(G2mth.ValEsd(val,sig),11)
372                WriteCIFitem(s)
373
374        def HillSortElements(elmlist):
375            '''Sort elements in "Hill" order: C, H, others, (where others
376            are alphabetical).
377
378            :params list elmlist: a list of element strings
379
380            :returns: a sorted list of element strings
381            '''
382            newlist = []
383            oldlist = elmlist[:]
384            for elm in ('C','H'):
385                if elm in elmlist:
386                    newlist.append(elm)
387                    oldlist.pop(oldlist.index(elm))
388            return newlist+sorted(oldlist)
389
390        def WriteComposition(phasenam):
391            '''determine the composition for the unit cell, crudely determine Z and
392            then compute the composition in formula units
393            '''
394            phasedict = self.Phases[phasenam] # pointer to current phase info
395            General = phasedict['General']
396            Z = General.get('cellZ',0.0)
397            cx,ct,cs,cia = General['AtomPtrs']
398            Atoms = phasedict['Atoms']
399            fpfx = str(phasedict['pId'])+'::Afrac:'       
400            cfrac = cx+3
401            cmult = cs+1
402            compDict = {} # combines H,D & T
403            sitemultlist = []
404            massDict = dict(zip(General['AtomTypes'],General['AtomMass']))
405            cellmass = 0
406            for i,at in enumerate(Atoms):
407                atype = at[ct].strip()
408                if atype.find('-') != -1: atype = atype.split('-')[0]
409                if atype.find('+') != -1: atype = atype.split('+')[0]
410                atype = atype[0].upper()+atype[1:2].lower() # force case conversion
411                if atype == "D" or atype == "D": atype = "H"
412                fvar = fpfx+str(i)
413                fval = self.parmDict.get(fvar,at[cfrac])
414                mult = at[cmult]
415                if not massDict.get(at[ct]):
416                    print 'No mass found for atom type '+at[ct]
417                    print 'Will not compute cell contents for phase '+phasenam
418                    return
419                cellmass += massDict[at[ct]]*mult*fval
420                compDict[atype] = compDict.get(atype,0.0) + mult*fval
421                if fval == 1: sitemultlist.append(mult)
422            if len(compDict.keys()) == 0: return # no elements!
423            if Z < 1: # Z has not been computed or set by user
424                Z = 1
425                for i in range(2,min(sitemultlist)+1):
426                    for m in sitemultlist:
427                        if m % i != 0:
428                            break
429                        else:
430                            Z = i
431                General['cellZ'] = Z # save it
432
433            # when scattering factors are included in the CIF, this needs to be
434            # added to the loop here but only in the one-block case.
435            # For multiblock CIFs, scattering factors go in the histogram
436            # blocks  (for all atoms in all appropriate phases)
437
438            #if oneblock: # add scattering factors for current phase here
439            WriteCIFitem('\nloop_  _atom_type_symbol _atom_type_number_in_cell')
440            formula = ''
441            reload(G2mth)
442            for elem in HillSortElements(compDict.keys()):
443                WriteCIFitem('  ' + PutInCol(elem,4) +
444                             G2mth.ValEsd(compDict[elem],-0.009,True))
445                if formula: formula += " "
446                formula += elem
447                if compDict[elem] == Z: continue
448                formula += G2mth.ValEsd(compDict[elem]/Z,-0.009,True)
449            WriteCIFitem( '\n# Note that Z affects _cell_formula_sum and _weight')
450            WriteCIFitem( '_cell_formula_units_Z',str(Z))
451            WriteCIFitem( '_chemical_formula_sum',formula)
452            WriteCIFitem( '_chemical_formula_weight',
453                          G2mth.ValEsd(cellmass/Z,-0.09,True))
454
455        def WriteDistances(phasenam,SymOpList,offsetList,symOpList,G2oprList):
456            '''Report bond distances and angles for the CIF
457
458            Note that _geom_*_symmetry_* fields are values of form
459            n_klm where n is the symmetry operation in SymOpList (counted
460            starting with 1) and (k-5, l-5, m-5) are translations to add
461            to (x,y,z). See
462            http://www.iucr.org/__data/iucr/cifdic_html/1/cif_core.dic/Igeom_angle_site_symmetry_.html
463
464            TODO: need a method to select publication flags for distances/angles
465            '''
466            phasedict = self.Phases[phasenam] # pointer to current phase info           
467            Atoms = phasedict['Atoms']
468            cx,ct,cs,cia = phasedict['General']['AtomPtrs']
469            fpfx = str(phasedict['pId'])+'::Afrac:'       
470            cfrac = cx+3
471            # loop over interatomic distances for this phase
472            WriteCIFitem('\n# MOLECULAR GEOMETRY')
473            WriteCIFitem('loop_' + 
474                         '\n\t_geom_bond_atom_site_label_1' +
475                         '\n\t_geom_bond_atom_site_label_2' + 
476                         '\n\t_geom_bond_distance' + 
477                         '\n\t_geom_bond_site_symmetry_1' + 
478                         '\n\t_geom_bond_site_symmetry_2' + 
479                         '\n\t_geom_bond_publ_flag')
480
481            # Note that labels should be read from self.labellist to correspond
482            # to the values reported in the atoms table and zero occupancy atoms
483            # should not be included
484            fpfx = str(phasedict['pId'])+'::Afrac:'       
485            for i,at in enumerate(Atoms):
486                if self.parmDict.get(fpfx+str(i),at[cfrac]) == 0.0: continue
487                lbl = self.labellist[i]
488
489
490            # loop over interatomic angles for this phase
491            WriteCIFitem('loop_' + 
492                         '\n\t_geom_angle_atom_site_label_1' + 
493                         '\n\t_geom_angle_atom_site_label_2' + 
494                         '\n\t_geom_angle_atom_site_label_3' + 
495                         '\n\t_geom_angle' + 
496                         '\n\t_geom_angle_site_symmetry_1' +
497                         '\n\t_geom_angle_site_symmetry_2' + 
498                         '\n\t_geom_angle_site_symmetry_3' + 
499                         '\n\t_geom_angle_publ_flag')
500
501
502        def WritePhaseInfo(phasenam):
503            WriteCIFitem('\n# phase info for '+str(phasenam) + ' follows')
504            phasedict = self.Phases[phasenam] # pointer to current phase info           
505            WriteCIFitem('_pd_phase_name', phasenam)
506            pfx = str(phasedict['pId'])+'::'
507            A,sigA = G2stIO.cellFill(pfx,phasedict['General']['SGData'],self.parmDict,self.sigDict)
508            cellSig = G2stIO.getCellEsd(pfx,
509                                       phasedict['General']['SGData'],A,
510                                       self.OverallParms['Covariance'])  # returns 7 vals, includes sigVol
511            cellList = G2lat.A2cell(A) + (G2lat.calc_V(A),)
512            defsigL = 3*[-0.00001] + 3*[-0.001] + [-0.01] # significance to use when no sigma
513            names = ['length_a','length_b','length_c',
514                     'angle_alpha','angle_beta ','angle_gamma',
515                     'volume']
516            prevsig = 0
517            for lbl,defsig,val,sig in zip(names,defsigL,cellList,cellSig):
518                if sig:
519                    txt = G2mth.ValEsd(val,sig)
520                    prevsig = -sig # use this as the significance for next value
521                else:
522                    txt = G2mth.ValEsd(val,min(defsig,prevsig),True)
523                WriteCIFitem('_cell_'+lbl,txt)
524                   
525            WriteCIFitem('_symmetry_cell_setting',
526                         phasedict['General']['SGData']['SGSys'])
527
528            spacegroup = phasedict['General']['SGData']['SpGrp'].strip()
529            # regularize capitalization and remove trailing H/R
530            spacegroup = spacegroup[0].upper() + spacegroup[1:].lower().rstrip('rh ')
531            WriteCIFitem('_symmetry_space_group_name_H-M',spacegroup)
532
533            # generate symmetry operations including centering and center of symmetry
534            SymOpList,offsetList,symOpList,G2oprList = G2spg.AllOps(
535                phasedict['General']['SGData'])
536            WriteCIFitem('loop_ _space_group_symop_id _space_group_symop_operation_xyz')
537            for i,op in enumerate(SymOpList,start=1):
538                WriteCIFitem('   {:3d}  {:}'.format(i,op.lower()))
539
540            # loop over histogram(s) used in this phase
541            if not oneblock and not self.quickmode:
542                # report pointers to the histograms used in this phase
543                histlist = []
544                for hist in self.Phases[phasenam]['Histograms']:
545                    if self.Phases[phasenam]['Histograms'][hist]['Use']:
546                        if phasebyhistDict.get(hist):
547                            phasebyhistDict[hist].append(phasenam)
548                        else:
549                            phasebyhistDict[hist] = [phasenam,]
550                        blockid = datablockidDict.get(hist)
551                        if not blockid:
552                            print "Internal error: no block for data. Phase "+str(
553                                phasenam)+" histogram "+str(hist)
554                            histlist = []
555                            break
556                        histlist.append(blockid)
557
558                if len(histlist) == 0:
559                    WriteCIFitem('# Note: phase has no associated data')
560                else:
561                    WriteCIFitem('loop_  _pd_block_diffractogram_id')
562
563            # report atom params
564            if phasedict['General']['Type'] == 'nuclear':        #this needs macromolecular variant, etc!
565                WriteAtomsNuclear(phasenam)
566            else:
567                raise Exception,"no export for mm coordinates implemented"
568            # report cell contents
569            WriteComposition(phasenam)
570            if not self.quickmode:      # report distances and angles
571                WriteDistances(phasenam,SymOpList,offsetList,symOpList,G2oprList)
572
573            #raise Exception,'Testing'
574
575        def WritePowderData(histlbl):
576            text = '?'
577            histblk = self.Histograms[histlbl]
578            print 'TODO: powder here data for',histblk["Sample Parameters"]['InstrName']
579            # see wrpowdhist.for & wrreflist.for
580           
581            refprx = '_refln.' # mm
582            refprx = '_refln_' # normal
583
584            if not oneblock:
585                if not phasebyhistDict.get(histlbl):
586                    WriteCIFitem('\n# No phases associated with this data set')
587                else:
588                    WriteCIFitem('\n# PHASE TABLE')
589                    WriteCIFitem('loop_' +
590                                 '\n\t_pd_phase_id' + 
591                                 '\n\t_pd_phase_block_id' + 
592                                 '\n\t_pd_phase_mass_%')
593                    for phasenam in phasebyhistDict.get(histlbl):
594                        pass
595
596            WriteCIFitem('\n# SCATTERING FACTOR INFO')
597            WriteCIFitem('_diffrn_radiation_wavelength' ,text)
598            #WriteCIFitem('_diffrn_radiation_type',text)
599            #C always assume Ka1 & Ka2 if two wavelengths are present
600            #WriteCIFitem('loop_' +
601            #             '\n\t_diffrn_radiation_wavelength' +
602            #             '\n\t_diffrn_radiation_wavelength_wt' +
603            #             '\n\t_diffrn_radiation_type' +
604            #             '\n\t_diffrn_radiation_wavelength_id')
605            #WRITE LAM1,1.0,'K\\a~1~',1, LAM2,ratio,'K\\a~2~',2
606
607            WriteCIFitem('_pd_proc_ls_prof_R_factor','?')
608            WriteCIFitem('_pd_proc_ls_prof_wR_factor','?')
609            WriteCIFitem('_pd_proc_ls_prof_wR_expected','?')
610            WriteCIFitem('_refine_ls_R_Fsqd_factor','?')
611           
612            #WriteCIFitem('_pd_meas_2theta_fixed',text)
613            WriteCIFitem('_diffrn_radiation_probe','x-ray')
614            WriteCIFitem('_diffrn_radiation_probe','neutron')
615            WriteCIFitem('_diffrn_radiation_polarisn_ratio','?')
616           
617            WriteCIFitem('loop_  _atom_type_symbol')
618            if oneblock:
619                WriteCIFitem('       _atom_type_number_in_cell')
620            #IF (HTYP(2:2) .eq. 'X' .AND. HTYP(3:3) .ne. 'E') THEN
621            WriteCIFitem('      _atom_type_scat_dispersion_real')
622            WriteCIFitem('      _atom_type_scat_dispersion_imag')
623            for lbl in ('a1','a2','a3', 'a4', 'b1', 'b2', 'b3', 'b4', 'c'):
624                WriteCIFitem('      _atom_type_scat_Cromer_Mann_'+lbl)
625            #ELSEIF (HTYP(2:2) .eq. 'N') THEN
626            WriteCIFitem('      _atom_type_scat_length_neutron')
627            #ENDIF
628            WriteCIFitem('      _atom_type_scat_source')
629
630            #C document the background function used
631            WriteCIFitem('_pd_proc_ls_background_function','?')
632
633            WriteCIFitem('_exptl_absorpt_process_details','?')
634            WriteCIFitem('_exptl_absorpt_correction_T_min','?')
635            WriteCIFitem('_exptl_absorpt_correction_T_max','?')
636            #C extinction
637            #WRITE(IUCIF,'(A)') '# Extinction correction'
638            #CALL WRVAL(IUCIF,'_gsas_exptl_extinct_corr_T_min',TEXT(1:10))
639            #CALL WRVAL(IUCIF,'_gsas_exptl_extinct_corr_T_max',TEXT(11:20))
640
641            if not oneblock:
642                # instrumental profile terms go here
643                WriteCIFitem('_pd_proc_ls_profile_function','?')
644
645            WriteCIFitem('\n# STRUCTURE FACTOR TABLE')           
646            WriteCIFitem('loop_' + 
647                         '\n\t' + refprx + 'index_h' + 
648                         '\n\t' + refprx + 'index_k' + 
649                         '\n\t' + refprx + 'index_l' + 
650                         '\n\t_pd_refln_wavelength_id' +
651                         '\n\t_pd_refln_phase_id' + 
652                         '\n\t' + refprx + 'status' + 
653                         '\n\t' + refprx + 'F_squared_meas' + 
654                         '\n\t' + refprx + 'F_squared_sigma' + 
655                         '\n\t' + refprx + 'F_squared_calc' + 
656                         '\n\t' + refprx + 'phase_calc' + 
657                         '\n\t_pd_refln_d_spacing' + 
658                         '\n\t_gsas_i100_meas')
659
660            WriteCIFitem('_reflns_number_total', text)
661            WriteCIFitem('_reflns_limit_h_min', text)
662            WriteCIFitem('_reflns_limit_h_max', text)
663            WriteCIFitem('_reflns_limit_k_min', text)
664            WriteCIFitem('_reflns_limit_k_max', text)
665            WriteCIFitem('_reflns_limit_l_min', text)
666            WriteCIFitem('_reflns_limit_l_max', text)
667            WriteCIFitem('_reflns_d_resolution_high', text)
668            WriteCIFitem('_reflns_d_resolution_low', text)
669           
670            WriteCIFitem('\n# POWDER DATA TABLE')
671            # is data fixed step?
672            fixedstep = False
673            # are ESDs sqrt(I)
674            countsdata = False
675            zero = 0.01
676            if fixedstep:
677                WriteCIFitem('_pd_meas_2theta_range_min', text)
678                WriteCIFitem('_pd_meas_2theta_range_max', text)
679                WriteCIFitem('_pd_meas_2theta_range_inc', text)
680                # zero correct
681                if zero != 0.0:
682                    WriteCIFitem('_pd_proc_2theta_range_min', text)
683                    WriteCIFitem('_pd_proc_2theta_range_max', text)
684                    WriteCIFitem('_pd_proc_2theta_range_inc', text)
685            WriteCIFitem('loop_' +
686                         '\n\t_pd_proc_d_spacing')
687                         #'_pd_meas_time_of_flight'
688            if not fixedstep:
689                if zero != 0.0:
690                    WriteCIFitem('\t_pd_proc_2theta_corrected')
691                else:
692                    WriteCIFitem('\t_pd_meas_2theta_scan')
693            if countsdata:
694                WriteCIFitem('\t_pd_meas_counts_total')
695            else:
696                WriteCIFitem('\t_pd_meas_intensity_total')
697            WriteCIFitem('\t_pd_proc_ls_weight')
698            WriteCIFitem('\t_pd_proc_intensity_bkg_calc')
699            WriteCIFitem('\t_pd_calc_intensity_total')
700            if zero != 0.0:
701                WriteCIFitem('_pd_proc_number_of_points', text)
702            else:
703                WriteCIFitem('_pd_meas_number_of_points', text)
704                         
705        def WriteSingleXtalData(histlbl):
706            histblk = self.Histograms[histlbl]
707            print 'TODO: single xtal here data for',histblk["Instrument Parameters"][0]['InstrName']
708            # see wrreflist.for
709            refprx = '_refln.' # mm
710            refprx = '_refln_' # normal
711           
712            WriteCIFitem('\n# STRUCTURE FACTOR TABLE')           
713            WriteCIFitem('loop_' + 
714                         '\n\t' + refprx + 'index_h' + 
715                         '\n\t' + refprx + 'index_k' + 
716                         '\n\t' + refprx + 'index_l' + 
717                         '\n\t' + refprx + 'status' + 
718                         '\n\t' + refprx + 'F_squared_meas' + 
719                         '\n\t' + refprx + 'F_squared_sigma' + 
720                         '\n\t' + refprx + 'F_squared_calc' + 
721                         '\n\t' + refprx + 'phase_calc')
722            WriteCIFitem('_reflns_number_total', text)
723            WriteCIFitem('_reflns_number_observed', text)
724            WriteCIFitem('_reflns_limit_h_min', text)
725            WriteCIFitem('_reflns_limit_h_max', text)
726            WriteCIFitem('_reflns_limit_k_min', text)
727            WriteCIFitem('_reflns_limit_k_max', text)
728            WriteCIFitem('_reflns_limit_l_min', text)
729            WriteCIFitem('_reflns_limit_l_max', text)
730            WriteCIFitem('_reflns_d_resolution_high', text)
731            WriteCIFitem('_reflns_d_resolution_low', text)
732
733        #============================================================
734        # the export process starts here
735        # create a dict with refined values and their uncertainties
736        self.loadParmDict()
737        # also load all of the tree into a set of dicts
738        self.loadTree()
739        #self.dumpTree()
740
741        # get restraint info
742        #restraintDict = self.OverallParms.get('Restraints',{})
743        #for i in  self.OverallParms['Constraints']:
744        #    print i
745        #    for j in self.OverallParms['Constraints'][i]:
746        #        print j
747        #return
748
749        self.CIFdate = dt.datetime.strftime(dt.datetime.now(),"%Y-%m-%dT%H:%M")
750        # count phases, powder and single crystal histograms
751        self.nphase = len(self.Phases)
752        self.npowder = 0
753        self.nsingle = 0
754        for hist in self.Histograms:
755            if hist.startswith("PWDR"): 
756                self.npowder += 1
757            elif hist.startswith("HKLF"): 
758                self.nsingle += 1
759        # is there anything to export?
760        if self.nphase + self.npowder + self.nsingle == 0: 
761            self.G2frame.ErrorDialog(
762                'Empty project',
763                'No data or phases to include in CIF')
764            return
765        # is there a file name defined?
766        self.CIFname = os.path.splitext(
767            os.path.split(self.G2frame.GSASprojectfile)[1]
768            )[0]
769        self.CIFname = self.CIFname.replace(' ','')
770        if not self.CIFname:
771            self.G2frame.ErrorDialog(
772                'No GPX name',
773                'Please save the project to provide a name')
774            return
775        # test for quick CIF mode or no data
776        self.quickmode = False
777        phasenam = phasenum = None # include all phases
778        if mode != "full" or self.npowder + self.nsingle == 0:
779            self.quickmode = True
780            oneblock = True
781            if self.nphase == 0:
782                self.G2frame.ErrorDialog(
783                    'No phase present',
784                    'Cannot create a coordinates CIF with no phases')
785                return
786            elif self.nphase > 1: # quick mode: choose one phase
787                choices = sorted(self.Phases.keys())
788                phasenum = G2gd.ItemSelector(choices,self.G2frame)
789                if phasenum is None: return
790                phasenam = choices[phasenum]
791        # will this require a multiblock CIF?
792        elif self.nphase > 1:
793            oneblock = False
794        elif self.npowder + self.nsingle > 1:
795            oneblock = False
796        else: # one phase, one dataset, Full CIF
797            oneblock = True
798
799        # make sure needed infomation is present
800        # get CIF author name -- required for full CIFs
801        try:
802            self.author = self.OverallParms['Controls'].get("Author",'').strip()
803        except KeyError:
804            pass
805        while not (self.author or self.quickmode):
806            dlg = G2gd.SingleStringDialog(self.G2frame,'Get CIF Author','Provide CIF Author name (Last, First)')
807            if not dlg.Show(): return # cancel was pressed
808            self.author = dlg.GetValue()
809            dlg.Destroy()
810        try:
811            self.OverallParms['Controls']["Author"] = self.author # save for future
812        except KeyError:
813            pass
814        self.shortauthorname = self.author.replace(',','').replace(' ','')[:20]
815
816        # check the instrument name for every histogram
817        if not self.quickmode:
818            dictlist = []
819            keylist = []
820            lbllist = []
821            invalid = 0
822            key3 = 'InstrName'
823            for hist in self.Histograms:
824                if hist.startswith("PWDR"): 
825                    key2 = "Sample Parameters"
826                    d = self.Histograms[hist][key2]
827                elif hist.startswith("HKLF"): 
828                    key2 = "Instrument Parameters"
829                    d = self.Histograms[hist][key2][0]
830                   
831                lbllist.append(hist)
832                dictlist.append(d)
833                keylist.append(key3)
834                instrname = d.get(key3)
835                if instrname is None:
836                    d[key3] = ''
837                    invalid += 1
838                elif instrname.strip() == '':
839                    invalid += 1
840            if invalid:
841                msg = ""
842                if invalid > 3: msg = (
843                    "\n\nNote: it may be faster to set the name for\n"
844                    "one histogram for each instrument and use the\n"
845                    "File/Copy option to duplicate the name"
846                    )
847                if not G2gd.CallScrolledMultiEditor(
848                    self.G2frame,dictlist,keylist,
849                    prelbl=range(1,len(dictlist)+1),
850                    postlbl=lbllist,
851                    title='Instrument names',
852                    header="Edit instrument names. Note that a non-blank\nname is required for all histograms"+msg,
853                    ): return
854
855        if oneblock and not self.quickmode:
856            # select a dataset to use (there should only be one set in one block,
857            # but take whatever comes 1st
858            for hist in self.Histograms:
859                histblk = self.Histograms[hist]
860                if hist.startswith("PWDR"): 
861                    instnam = histblk["Sample Parameters"]['InstrName']
862                    break # ignore all but 1st data histogram
863                elif hist.startswith("HKLF"): 
864                    instnam = histblk["Instrument Parameters"][0]['InstrName']
865                    break # ignore all but 1st data histogram
866        #======================================================================
867        # Start writing the CIF - single block
868        #======================================================================
869        if oneblock:
870            WriteCIFitem('data_'+self.CIFname)
871            if phasenam is None: # if not already selected, select the first phase (should be one)
872                phasenam = self.Phases.keys()[0]
873            #print 'phasenam',phasenam
874            phaseblk = self.Phases[phasenam] # pointer to current phase info
875            if not self.quickmode:
876                instnam = instnam.replace(' ','')
877                WriteCIFitem('_pd_block_id',
878                             str(self.CIFdate) + "|" + str(self.CIFname) + "|" +
879                             str(self.shortauthorname) + "|" + instnam)
880                WriteAudit()
881                WritePubTemplate()
882                WriteOverall()
883            # report the phase info
884            WritePhaseTemplate()
885            WritePhaseInfo(phasenam)
886            if hist.startswith("PWDR") and not self.quickmode:
887                # preferred orientation
888                SH = FormatSH(phasenam)
889                MD = FormatHAPpo(phasenam)
890                if SH and MD:
891                    WriteCIFitem('_pd_proc_ls_pref_orient_corr', SH + '\n' + MD)
892                elif SH or MD:
893                    WriteCIFitem('_pd_proc_ls_pref_orient_corr', SH + MD)
894                else:
895                    WriteCIFitem('_pd_proc_ls_pref_orient_corr', 'none')
896                    # report profile, since one-block: include both histogram and phase info
897                WriteCIFitem('_pd_proc_ls_profile_function',
898                             FormatInstProfile(histblk["Instrument Parameters"])
899                             + '\n' +
900                             FormatPhaseProfile(phasenam))
901                WritePowderTemplate()
902                WritePowderData(hist)
903            elif hist.startswith("HKLF") and not self.quickmode:
904                WriteSnglXtalTemplate()
905                WriteSingleXtalData(hist)
906        else:
907        #======================================================================
908        # Start writing the CIF - multiblock
909        #======================================================================
910            # publication info
911            WriteCIFitem('\ndata_'+self.CIFname+'_publ')
912            WriteAudit()
913            WriteCIFitem('_pd_block_id',
914                         str(self.CIFdate) + "|" + str(self.CIFname) + "|" +
915                         str(self.shortauthorname) + "|Overall")
916            WritePubTemplate()
917            # overall info
918            WriteCIFitem('data_'+str(self.CIFname)+'_overall')
919            WriteOverall()
920            #============================================================
921            WriteCIFitem('# POINTERS TO PHASE AND HISTOGRAM BLOCKS')
922            datablockidDict = {} # save block names here -- N.B. check for conflicts between phase & hist names (unlikely!)
923            # loop over phase blocks
924            if self.nphase > 1:
925                loopprefix = ''
926                WriteCIFitem('loop_   _pd_phase_block_id')
927            else:
928                loopprefix = '_pd_phase_block_id'
929           
930            for phasenam in sorted(self.Phases.keys()):
931                i = self.Phases[phasenam]['pId']
932                datablockidDict[phasenam] = (str(self.CIFdate) + "|" + str(self.CIFname) + "|" +
933                             'phase_'+ str(i) + '|' + str(self.shortauthorname))
934                WriteCIFitem(loopprefix,datablockidDict[phasenam])
935            # loop over data blocks
936            if self.npowder + self.nsingle > 1:
937                loopprefix = ''
938                WriteCIFitem('loop_   _pd_block_diffractogram_id')
939            else:
940                loopprefix = '_pd_block_diffractogram_id'
941            for hist in self.Histograms:
942                histblk = self.Histograms[hist]
943                if hist.startswith("PWDR"): 
944                    instnam = histblk["Sample Parameters"]['InstrName']
945                elif hist.startswith("HKLF"): 
946                    instnam = histblk["Instrument Parameters"][0]['InstrName']
947                instnam = instnam.replace(' ','')
948                i = histblk['hId']
949                datablockidDict[hist] = (str(self.CIFdate) + "|" + str(self.CIFname) + "|" +
950                                         str(self.shortauthorname) + "|" +
951                                         instnam + "_hist_"+str(i))
952                WriteCIFitem(loopprefix,datablockidDict[hist])
953            #============================================================
954            # loop over phases, exporting them
955            phasebyhistDict = {} # create a cross-reference to phases by histogram
956            for j,phasenam in enumerate(sorted(self.Phases.keys())):
957                i = self.Phases[phasenam]['pId']
958                WriteCIFitem('\ndata_'+self.CIFname+"_phase_"+str(i))
959                print "debug, processing ",phasenam
960                WriteCIFitem('# Information for phase '+str(i))
961                WriteCIFitem('_pd_block_id',datablockidDict[phasenam])
962                # report the phase
963                WritePhaseTemplate()
964                WritePhaseInfo(phasenam)
965                # preferred orientation
966                SH = FormatSH(phasenam)
967                MD = FormatHAPpo(phasenam)
968                if SH and MD:
969                    WriteCIFitem('_pd_proc_ls_pref_orient_corr', SH + '\n' + MD)
970                elif SH or MD:
971                    WriteCIFitem('_pd_proc_ls_pref_orient_corr', SH + MD)
972                else:
973                    WriteCIFitem('_pd_proc_ls_pref_orient_corr', 'none')
974                # report sample profile terms
975                PP = FormatPhaseProfile(phasenam)
976                if PP:
977                    WriteCIFitem('_pd_proc_ls_profile_function',PP)
978                   
979            #============================================================
980            # loop over histograms, exporting them
981            for hist in self.Histograms:
982                histblk = self.Histograms[hist]
983                i = histblk['hId']
984                if hist.startswith("PWDR"): 
985                    WriteCIFitem('\ndata_'+self.CIFname+"_pwd_"+str(i))
986                    #instnam = histblk["Sample Parameters"]['InstrName']
987                    # report instrumental profile terms
988                    WriteCIFitem('_pd_proc_ls_profile_function',
989                                 FormatInstProfile(histblk["Instrument Parameters"]))
990                    WriteCIFitem('# Information for histogram '+str(i)+': '+
991                                 hist)
992                    WriteCIFitem('_pd_block_id',datablockidDict[hist])
993                    WritePowderTemplate()
994                    WritePowderData(key1)
995                elif hist.startswith("HKLF"): 
996                    WriteCIFitem('\ndata_'+self.CIFname+"_sx_"+str(i))
997                    #instnam = histblk["Instrument Parameters"][0]['InstrName']
998                    WriteCIFitem('# Information for histogram '+str(i)+': '+
999                                 hist)
1000                    WriteCIFitem('_pd_block_id',datablockidDict[hist])
1001                    WriteSnglXtalTemplate()
1002                    WriteSingleXtalData(key1)
1003
1004        WriteCIFitem('#--' + 15*'eof--' + '#')
1005
Note: See TracBrowser for help on using the repository browser.