source: trunk/exports/G2cif.py @ 981

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

introduce regress option; fix esd printing; more docs; new Mac app with drag & drop for open; control reset of ref list on load

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