source: trunk/exports/G2cif.py @ 960

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

check in CIF dev snapshot

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