Changeset 960 for trunk/exports


Ignore:
Timestamp:
Jun 19, 2013 5:51:24 PM (8 years ago)
Author:
toby
Message:

check in CIF dev snapshot

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/exports/G2cif.py

    r956 r960  
    8080            WriteCIFitem('_refine_ls_number_parameters',vars)
    8181            try:
    82                 GOF = str(self.OverallParms['Covariance']['Rvals']['GOF'])
     82                GOF = G2mth.ValEsd(self.OverallParms['Covariance']['Rvals']['GOF'],-0.009)
    8383            except:
    8484                GOF = '?'
    8585            WriteCIFitem('_refine_ls_goodness_of_fit_all',GOF)
    86             #WriteCIFitem('_refine_ls_number_restraints',TEXT(1:7))
     86
     87            #WriteCIFitem('_refine_ls_number_restraints',TEXT)
     88           
    8789            # other things to consider reporting
    8890            # _refine_ls_number_reflns
     
    131133            print getCallerDocString()
    132134
    133         def FormatSH(phasedict):
     135        def FormatSH(phasenam):
    134136            'Format a full spherical harmonics texture description as a string'
    135             # SH Texture
     137            phasedict = self.Phases[phasenam] # pointer to current phase info           
    136138            pfx = str(phasedict['pId'])+'::'
    137139            s = ""
     
    160162            return s
    161163
    162         def FormatHAPpo(phasedict):
    163             'return the March-Dollase/SH correction for every histogram in the current phase'
     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           
    164170            s = ''
    165171            for histogram in sorted(phasedict['Histograms']):
     172                if histogram.startswith("HKLF"): continue # powder only
    166173                Histogram = self.Histograms.get(histogram)
    167174                if not Histogram: continue
     
    197204                    s += s1
    198205            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'
    199228       
    200229        def FmtAtomType(sym):
    201230            'Reformat a GSAS-II atom type symbol to match CIF rules'
    202231            sym = sym.replace('_','') # underscores are not allowed: no isotope designation?
    203             # in CIF oxidation state symbols come after, not before
     232            # in CIF, oxidation state sign symbols come after, not before
    204233            if '+' in sym:
    205234                sym = sym.replace('+','') + '+'
     
    241270
    242271        def WriteAtomsNuclear(phasenam):
    243             phasedict = self.Phases[phasenam] # pointer to current phase info           
     272            'Write atom positions to CIF'
     273            phasedict = self.Phases[phasenam] # pointer to current phase info
    244274            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               
    245287            WriteCIFitem('\n# ATOMIC COORDINATES AND DISPLACEMENT PARAMETERS')
    246             Atoms = phasedict['Atoms']
    247288            WriteCIFitem('loop_ '+
    248289                         '\n\t_atom_site_label'+
     
    256297                         '\n\t_atom_site_symmetry_multiplicity')
    257298
    258             varnames = {3:'Ax',4:'Ay',5:'Az',6:'Afrac',
    259                         10:'AUiso',11:'AU11',12:'AU22',13:'AU33',
    260                         14:'AU12',15:'AU13',16:'AU23'}
     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'}
    261302            labellist = []
    262303           
    263304            pfx = str(phasedict['pId'])+'::'
     305            # loop over all atoms
     306            naniso = 0
    264307            for i,at in enumerate(Atoms):
    265                 print at
    266             for i,at in enumerate(Atoms):
    267                 s = " "
    268                 s += PutInCol(MakeUniqueLabel(at[0],labellist),6) # label
    269                 #s += PutInCol(MakeUniqueLabel('A',labellist),6) # label
    270                 s += PutInCol(FmtAtomType(at[1]),6) # type
    271                 if at[9] == 'I':
     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':
    272313                    adp = 'Uiso '
    273314                else:
    274315                    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.
    275320                    t = 0.0
    276                     for j in (11,12,13):
    277                         var = pfx+varnames[j]+":"+str(i)
    278                         t += self.parmDict.get(var,at[j])
    279                 for j in (3,4,5,6,10):
    280                     if j in (3,4,5):
     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):
    281326                        dig = 11
    282327                        sigdig = -0.00009
    283328                    else:
    284                         dig = 9
     329                        dig = 10
    285330                        sigdig = -0.009
    286331                    var = pfx+varnames[j]+":"+str(i)
     
    288333                    if dvar not in self.sigDict:
    289334                        dvar = var
    290                     if j == 10 and adp == 'Uani ':
    291                         # compute Uequiv crudely
    292                         t = 0.0
    293                         for k in (11,12,13):
    294                             var = pfx+varnames[j]+":"+str(i)
    295                             t += self.parmDict.get(var,at[k])
     335                    if j == cia+1 and adp == 'Uani ':
    296336                        val = t/3.
    297337                        sig = sigdig
    298338                    else:
     339                        #print var,(var in self.parmDict),(var in self.sigDict)
    299340                        val = self.parmDict.get(var,at[j])
    300341                        sig = self.sigDict.get(dvar,sigdig)
    301342                    s += PutInCol(G2mth.ValEsd(val,sig),dig)
    302343                s += adp
    303                 print s
    304        
     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
    305446        def WritePhaseInfo(phasenam):
    306447            # see writepha.for
     
    309450            phasedict = self.Phases[phasenam] # pointer to current phase info           
    310451            WriteCIFitem('_pd_phase_name', phasenam)
    311             # replace some of these later
    312             #General = phasedict['General']
    313             #SGData = phasedict['General']['SGData']
    314             #cell = General['Cell']
    315452            pfx = str(phasedict['pId'])+'::'
    316             #covData = self.OverallParms['Covariance']
    317453            A,sigA = G2stIO.cellFill(pfx,phasedict['General']['SGData'],self.parmDict,self.sigDict)
    318454            cellSig = G2stIO.getCellEsd(pfx,
     
    346482                WriteCIFitem('   {:3d}  {:}'.format(i,op.lower()))
    347483
    348             # preferred orientation (always reported by phase)
    349             SH = FormatSH(phasedict)
    350             MD = FormatHAPpo(phasedict)
    351             if SH and MD:
    352                 WriteCIFitem('_pd_proc_ls_pref_orient_corr', SH + '\n' + MD)
    353             elif SH or MD:
    354                 WriteCIFitem('_pd_proc_ls_pref_orient_corr', SH + MD)
    355             else:
    356                 WriteCIFitem('_pd_proc_ls_pref_orient_corr', 'none')
    357 
    358484            # loop over histogram(s) used in this phase
    359             if oneblock:
    360                 # report all profile information here (include histogram & phase)
    361                 # _pd_proc_ls_profile_function
    362                 pass
    363             else:
    364                 # pointers to histograms used in this phase
     485            if not oneblock and not self.quickmode:
     486                # report pointers to the histograms used in this phase
    365487                histlist = []
    366488                for hist in self.Phases[phasenam]['Histograms']:
     
    383505                    WriteCIFitem('loop_  _pd_block_diffractogram_id')
    384506
    385                 # report sample profile information here (phase only)
    386                 # _pd_proc_ls_profile_function
    387            
     507            # report atom params
    388508            if phasedict['General']['Type'] == 'nuclear':        #this needs macromolecular variant, etc!
    389509                WriteAtomsNuclear(phasenam)
    390510            else:
    391511                raise Exception,"no export for mm coordinates implemented"
    392      
    393 
    394             raise Exception,'Testing'
    395 
    396             WriteCIFitem('loop_' + '\n\t_atom_site_aniso_label' +
    397                          '\n\t_atom_site_aniso_U_11' + '\n\t_atom_site_aniso_U_12' +
    398                          '\n\t_atom_site_aniso_U_13' + '\n\t_atom_site_aniso_U_22' +
    399                          '\n\t_atom_site_aniso_U_23' + '\n\t_atom_site_aniso_U_33')
    400 
    401             # process the chemical formula: pick a Z value & generate molecular weight
    402             # find the maximum possible Z value
    403 
    404             # order the elements in "Hill" order: C, H, D, T & alphabetical or alphabetical
    405             if not oneblock: # in single block, this is combined with the scattering factors
    406                 WriteCIFitem('loop_  _atom_type_symbol _atom_type_number_in_cell')
    407 
    408             WriteCIFitem('# If you change Z, be sure to change all 3 of the following')
    409             WriteCIFitem( '_chemical_formula_sum',text)
    410             #WRITE(TEXT,'(F15.2)') ATMASS
    411             WriteCIFitem( '_chemical_formula_weight',text)
    412             #WRITE(TEXT,'(I4)') Z
    413             WriteCIFitem( '_cell_formula_units_Z',text)
     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'
    414519
    415520            #C now loop over interatomic distances for this phase
     
    714819                    ): return
    715820
     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
    716832        #======================================================================
    717833        # Start writing the CIF - single block
     
    724840            phaseblk = self.Phases[phasenam] # pointer to current phase info
    725841            if not self.quickmode:
    726                 # select data, should only be one set in oneblock, but take whatever comes 1st
    727                 for hist in self.Histograms:
    728                     histblk = self.Histograms[hist]
    729                     if hist.startswith("PWDR"):
    730                         instnam = histblk["Sample Parameters"]['InstrName']
    731                         break # ignore all but 1st data histogram
    732                     elif hist.startswith("HKLF"):
    733                         instnam = histblk["Instrument Parameters"][0]['InstrName']
    734                         break # ignore all but 1st data histogram
    735842                instnam = instnam.replace(' ','')
    736843                WriteCIFitem('_pd_block_id',
     
    743850            WritePhaseTemplate()
    744851            WritePhaseInfo(phasenam)
    745             if not self.quickmode:
    746                 if hist.startswith("PWDR"):
    747                     WritePowderTemplate()
    748                     WritePowderData(hist)
    749                 elif hist.startswith("HKLF"):
    750                     WriteSnglXtalTemplate()
    751                     WriteSingleXtalData(hist)
     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)
    752872        else:
    753873        #======================================================================
     
    809929                WritePhaseTemplate()
    810930                WritePhaseInfo(phasenam)
    811 
     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                   
    812945            #============================================================
    813946            # loop over histograms, exporting them
     
    818951                    WriteCIFitem('\ndata_'+self.CIFname+"_pwd_"+str(i))
    819952                    #instnam = histblk["Sample Parameters"]['InstrName']
     953                    # report instrumental profile terms
     954                    WriteCIFitem('_pd_proc_ls_profile_function',
     955                                 FormatInstProfile(histblk["Instrument Parameters"]))
    820956                    WriteCIFitem('# Information for histogram '+str(i)+': '+
    821957                                 hist)
Note: See TracChangeset for help on using the changeset viewer.