Changeset 4402 for trunk


Ignore:
Timestamp:
Apr 15, 2020 5:13:34 PM (3 years ago)
Author:
toby
Message:

add scattering factors to CIF; catch overwrite of 1st reflection bug(\!\!)

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIIntPDFtool.py

    r4364 r4402  
    1010########### SVN repository information ###################
    1111'''
    12 .. _GSASIIIntPDFtool:
     12.. _GSASIIautoint:
    1313
    1414*GSASIIIntPDFtool: autointegration routines*

  • trunk/exports/G2export_CIF.py

    r4361 r4402  
    5858import GSASIIstrMain as G2stMn
    5959import GSASIImapvars as G2mv
     60import GSASIIElem as G2el
    6061
    6162DEBUG = False    #True to skip printing of reflection/powder profile lists
    … …  
    376377
    377378# Refactored over here to allow access by GSASIIscriptable.py
    378 def WriteComposition(fp, phasedict, phasenam, parmDict):
     379def WriteComposition(fp, phasedict, phasenam, parmDict, quickmode=True, keV=None):
    379380    '''determine the composition for the unit cell, crudely determine Z and
    380     then compute the composition in formula units
     381    then compute the composition in formula units.
     382
     383    If quickmode is False, then scattering factors are added to the element loop.
     384
     385    If keV is specified, then resonant scattering factors are also computed and included.
    381386    '''
    382387    General = phasedict['General']
    … …  
    421426        General['cellZ'] = Z # save it
    422427
    423     # when scattering factors are included in the CIF, this needs to be
    424     # added to the loop here but only in the one-block case.
    425     # For multiblock CIFs, scattering factors go in the histogram
    426     # blocks  (for all atoms in all appropriate phases) - an example?:
    427     #loop_
    428     #    _atom_type_symbol
    429     #    _atom_type_description
    430     #    _atom_type_scat_dispersion_real
    431     #    _atom_type_scat_dispersion_imag
    432     #    _atom_type_scat_source
    433     #    'C' 'C' 0.0033 0.0016
    434     #                         'International Tables Vol C Tables 4.2.6.8 and 6.1.1.4'
    435     #    'H' 'H' 0.0000 0.0000
    436     #                         'International Tables Vol C Tables 4.2.6.8 and 6.1.1.4'
    437     #    'P' 'P' 0.1023 0.0942
    438     #                         'International Tables Vol C Tables 4.2.6.8 and 6.1.1.4'
    439     #    'Cl' 'Cl' 0.1484 0.1585
    440     #                         'International Tables Vol C Tables 4.2.6.8 and 6.1.1.4'
    441     #    'Cu' 'Cu' 0.3201 1.2651
    442     #                         'International Tables Vol C Tables 4.2.6.8 and 6.1.1.4'
    443 
    444     #if oneblock: # add scattering factors for current phase here
     428    if not quickmode:
     429        FFtable = G2el.GetFFtable(General['AtomTypes'])
     430        BLtable = G2el.GetBLtable(General)
     431
    445432    WriteCIFitem(fp, '\nloop_  _atom_type_symbol _atom_type_number_in_cell')
     433    s = '       '
     434    if not quickmode:
     435        for j in ('a1','a2','a3','a4','b1','b2','b3','b4','c',1,2):
     436            if len(s) > 80:
     437                WriteCIFitem(fp, s)
     438                s = '       '
     439            if j==1:
     440                s += ' _atom_type_scat_source'
     441            elif j==2:
     442                s += ' _atom_type_scat_length_neutron'
     443            else:
     444                s += ' _atom_type_scat_Cromer_Mann_'
     445                s += j
     446        if keV:
     447            WriteCIFitem(fp, s)
     448            s = '        _atom_type_scat_dispersion_real _atom_type_scat_dispersion_imag _atom_type_scat_dispersion_source'
     449        WriteCIFitem(fp, s)
     450
     451           
    446452    formula = ''
    447453    for elem in HillSortElements(list(compDict.keys())):
    448         WriteCIFitem(fp, '  ' + PutInCol(elem,4) +
    449                      G2mth.ValEsd(compDict[elem],-0.009,True))
     454        s = '  '
     455        s += PutInCol(elem,4)
     456        s += PutInCol(G2mth.ValEsd(compDict[elem],-0.009,True),5)
     457        if not quickmode:
     458            for i in 'fa','fb','fc':
     459                if i != 'fc':
     460                    for j in range(4):
     461                        if elem in FFtable:
     462                            val = G2mth.ValEsd(FFtable[elem][i][j],-0.0009,True)
     463                        else:
     464                            val = '?'
     465                        s += ' '
     466                        s += PutInCol(val,9)
     467                else:
     468                    if elem in FFtable:
     469                        val = G2mth.ValEsd(FFtable[elem][i],-0.0009,True)
     470                    else:
     471                        val = '?'
     472                    s += ' '
     473                    s += PutInCol(val,9)
     474            if elem in BLtable:
     475                bldata = BLtable[elem]
     476                #isotope = bldata[0]
     477                #mass = bldata[1]['Mass']
     478                if 'BW-LS' in bldata[1]:
     479                    val = 0
     480                else:
     481                    val = G2mth.ValEsd(bldata[1]['SL'][0],-0.0009,True)
     482            else:
     483                val = '?'
     484            s += ' '
     485            s += PutInCol(val,9)
     486            WriteCIFitem(fp,s.rstrip())
     487            WriteCIFitem(fp,'      https://subversion.xray.aps.anl.gov/pyGSAS/trunk/atmdata.py')
     488            if keV:
     489                Orbs = G2el.GetXsectionCoeff(elem.split('+')[0].split('-')[0])
     490                FP,FPP,Mu = G2el.FPcalc(Orbs, keV)
     491                WriteCIFitem(fp,'  {:8.3f}{:8.3f}   https://subversion.xray.aps.anl.gov/pyGSAS/trunk/atmdata.py'.format(FP,FPP))
     492        else:
     493            WriteCIFitem(fp,s.rstrip())
    450494        if formula: formula += " "
    451495        formula += elem
    … …  
    866910                labellist.append(lbl)
    867911
    868         # Factored out to above by Jack O'Donnell, so it
    869         # could be accessed by GSASIIscriptable.py
    870         #  def WriteAtomsNuclear(phasenam):
    871         #      'Write atom positions to CIF'
    872         #      phasedict = self.Phases[phasenam] # pointer to current phase info
    873         #      General = phasedict['General']
    874         #      cx,ct,cs,cia = General['AtomPtrs']
    875         #      Atoms = phasedict['Atoms']
    876         #      cfrac = cx+3
    877         #      fpfx = str(phasedict['pId'])+'::Afrac:'
    878         #      for i,at in enumerate(Atoms):
    879         #          fval = self.parmDict.get(fpfx+str(i),at[cfrac])
    880         #          if fval != 0.0:
    881         #              break
    882         #      else:
    883         #          WriteCIFitem(self.fp, '\n# PHASE HAS NO ATOMS!')
    884         #          return
    885 
    886         #      WriteCIFitem(self.fp, '\n# ATOMIC COORDINATES AND DISPLACEMENT PARAMETERS')
    887         #      WriteCIFitem(self.fp, 'loop_ '+
    888         #                   '\n   _atom_site_label'+
    889         #                   '\n   _atom_site_type_symbol'+
    890         #                   '\n   _atom_site_fract_x'+
    891         #                   '\n   _atom_site_fract_y'+
    892         #                   '\n   _atom_site_fract_z'+
    893         #                   '\n   _atom_site_occupancy'+
    894         #                   '\n   _atom_site_adp_type'+
    895         #                   '\n   _atom_site_U_iso_or_equiv'+
    896         #                   '\n   _atom_site_symmetry_multiplicity')
    897 
    898         #      varnames = {cx:'Ax',cx+1:'Ay',cx+2:'Az',cx+3:'Afrac',
    899         #                  cia+1:'AUiso',cia+2:'AU11',cia+3:'AU22',cia+4:'AU33',
    900         #                  cia+5:'AU12',cia+6:'AU13',cia+7:'AU23'}
    901         #      self.labellist = []
    902 
    903         #      pfx = str(phasedict['pId'])+'::'
    904         #      # loop over all atoms
    905         #      naniso = 0
    906         #      for i,at in enumerate(Atoms):
    907         #          if phasedict['General']['Type'] == 'macromolecular':
    908         #              label = '%s_%s_%s_%s'%(at[ct-1],at[ct-3],at[ct-4],at[ct-2])
    909         #              s = PutInCol(MakeUniqueLabel(label,self.labellist),15) # label
    910         #          else:
    911         #              s = PutInCol(MakeUniqueLabel(at[ct-1],self.labellist),6) # label
    912         #          fval = self.parmDict.get(fpfx+str(i),at[cfrac])
    913         #          if fval == 0.0: continue # ignore any atoms that have a occupancy set to 0 (exact)
    914         #          s += PutInCol(FmtAtomType(at[ct]),4) # type
    915         #          if at[cia] == 'I':
    916         #              adp = 'Uiso '
    917         #          else:
    918         #              adp = 'Uani '
    919         #              naniso += 1
    920         #              # compute Uequiv crudely
    921         #              # correct: Defined as "1/3 trace of diagonalized U matrix".
    922         #              # SEE cell2GS & Uij2Ueqv to GSASIIlattice. Former is needed to make the GS matrix used by the latter.
    923         #              t = 0.0
    924         #              for j in (2,3,4):
    925         #                  var = pfx+varnames[cia+j]+":"+str(i)
    926         #                  t += self.parmDict.get(var,at[cia+j])
    927         #          for j in (cx,cx+1,cx+2,cx+3,cia,cia+1):
    928         #              if j in (cx,cx+1,cx+2):
    929         #                  dig = 11
    930         #                  sigdig = -0.00009
    931         #              else:
    932         #                  dig = 10
    933         #                  sigdig = -0.009
    934         #              if j == cia:
    935         #                  s += adp
    936         #              else:
    937         #                  var = pfx+varnames[j]+":"+str(i)
    938         #                  dvar = pfx+"d"+varnames[j]+":"+str(i)
    939         #                  if dvar not in self.sigDict:
    940         #                      dvar = var
    941         #                  if j == cia+1 and adp == 'Uani ':
    942         #                      val = t/3.
    943         #                      sig = sigdig
    944         #                  else:
    945         #                      #print var,(var in self.parmDict),(var in self.sigDict)
    946         #                      val = self.parmDict.get(var,at[j])
    947         #                      sig = self.sigDict.get(dvar,sigdig)
    948         #                  s += PutInCol(G2mth.ValEsd(val,sig),dig)
    949         #          s += PutInCol(at[cs+1],3)
    950         #          WriteCIFitem(self.fp, s)
    951         #      if naniso == 0: return
    952         #      # now loop over aniso atoms
    953         #      WriteCIFitem(self.fp, '\nloop_' + '\n   _atom_site_aniso_label' +
    954         #                   '\n   _atom_site_aniso_U_11' + '\n   _atom_site_aniso_U_22' +
    955         #                   '\n   _atom_site_aniso_U_33' + '\n   _atom_site_aniso_U_12' +
    956         #                   '\n   _atom_site_aniso_U_13' + '\n   _atom_site_aniso_U_23')
    957         #      for i,at in enumerate(Atoms):
    958         #          fval = self.parmDict.get(fpfx+str(i),at[cfrac])
    959         #          if fval == 0.0: continue # ignore any atoms that have a occupancy set to 0 (exact)
    960         #          if at[cia] == 'I': continue
    961         #          s = PutInCol(self.labellist[i],6) # label
    962         #          for j in (2,3,4,5,6,7):
    963         #              sigdig = -0.0009
    964         #              var = pfx+varnames[cia+j]+":"+str(i)
    965         #              val = self.parmDict.get(var,at[cia+j])
    966         #              sig = self.sigDict.get(var,sigdig)
    967         #              s += PutInCol(G2mth.ValEsd(val,sig),11)
    968         #          WriteCIFitem(self.fp, s)
    969 
    970         # def HillSortElements(elmlist):
    971         #     '''Sort elements in "Hill" order: C, H, others, (where others
    972         #     are alphabetical).
    973 
    974         #     :params list elmlist: a list of element strings
    975 
    976         #     :returns: a sorted list of element strings
    977         #     '''
    978         #     newlist = []
    979         #     oldlist = elmlist[:]
    980         #     for elm in ('C','H'):
    981         #         if elm in elmlist:
    982         #             newlist.append(elm)
    983         #             oldlist.pop(oldlist.index(elm))
    984         #     return newlist+sorted(oldlist)
    985 
    986         # Factored out to above by Jackson O'Donnell
    987         # so that it can be accessed by GSASIIscriptable
    988 
    989         # def WriteComposition(phasenam):
    990         #     '''determine the composition for the unit cell, crudely determine Z and
    991         #     then compute the composition in formula units
    992         #     '''
    993         #     phasedict = self.Phases[phasenam] # pointer to current phase info
    994         #     General = phasedict['General']
    995         #     Z = General.get('cellZ',0.0)
    996         #     cx,ct,cs,cia = General['AtomPtrs']
    997         #     Atoms = phasedict['Atoms']
    998         #     fpfx = str(phasedict['pId'])+'::Afrac:'
    999         #     cfrac = cx+3
    1000         #     cmult = cs+1
    1001         #     compDict = {} # combines H,D & T
    1002         #     sitemultlist = []
    1003         #     massDict = dict(zip(General['AtomTypes'],General['AtomMass']))
    1004         #     cellmass = 0
    1005         #     for i,at in enumerate(Atoms):
    1006         #         atype = at[ct].strip()
    1007         #         if atype.find('-') != -1: atype = atype.split('-')[0]
    1008         #         if atype.find('+') != -1: atype = atype.split('+')[0]
    1009         #         atype = atype[0].upper()+atype[1:2].lower() # force case conversion
    1010         #         if atype == "D" or atype == "D": atype = "H"
    1011         #         fvar = fpfx+str(i)
    1012         #         fval = self.parmDict.get(fvar,at[cfrac])
    1013         #         mult = at[cmult]
    1014         #         if not massDict.get(at[ct]):
    1015         #             print('Error: No mass found for atom type '+at[ct])
    1016         #             print('Will not compute cell contents for phase '+phasenam)
    1017         #             return
    1018         #         cellmass += massDict[at[ct]]*mult*fval
    1019         #         compDict[atype] = compDict.get(atype,0.0) + mult*fval
    1020         #         if fval == 1: sitemultlist.append(mult)
    1021         #     if len(compDict.keys()) == 0: return # no elements!
    1022         #     if Z < 1: # Z has not been computed or set by user
    1023         #         Z = 1
    1024         #         if not sitemultlist:
    1025         #             General['cellZ'] = 1
    1026         #             return
    1027         #         for i in range(2,min(sitemultlist)+1):
    1028         #             for m in sitemultlist:
    1029         #                 if m % i != 0:
    1030         #                     break
    1031         #                 else:
    1032         #                     Z = i
    1033         #         General['cellZ'] = Z # save it
    1034 
    1035         #     # when scattering factors are included in the CIF, this needs to be
    1036         #     # added to the loop here but only in the one-block case.
    1037         #     # For multiblock CIFs, scattering factors go in the histogram
    1038         #     # blocks  (for all atoms in all appropriate phases) - an example?:
    1039 #loop_
    1040 #    _at# om_type_symbol
    1041 #    _at# om_type_description
    1042 #    _at# om_type_scat_dispersion_real
    1043 #    _at# om_type_scat_dispersion_imag
    1044 #    _at# om_type_scat_source
    1045 #    'C'#  'C' 0.0033 0.0016
    1046 #       #                   'International Tables Vol C Tables 4.2.6.8 and 6.1.1.4'
    1047 #    'H'#  'H' 0.0000 0.0000
    1048 #       #                   'International Tables Vol C Tables 4.2.6.8 and 6.1.1.4'
    1049 #    'P'#  'P' 0.1023 0.0942
    1050 #       #                   'International Tables Vol C Tables 4.2.6.8 and 6.1.1.4'
    1051 #    'Cl# ' 'Cl' 0.1484 0.1585
    1052 #       #                   'International Tables Vol C Tables 4.2.6.8 and 6.1.1.4'
    1053 #    'Cu# ' 'Cu' 0.3201 1.2651
    1054 #       #                   'International Tables Vol C Tables 4.2.6.8 and 6.1.1.4'
    1055 
    1056         #     #if oneblock: # add scattering factors for current phase here
    1057         #     WriteCIFitem(self.fp, '\nloop_  _atom_type_symbol _atom_type_number_in_cell')
    1058         #     formula = ''
    1059         #     for elem in HillSortElements(compDict.keys()):
    1060         #         WriteCIFitem(self.fp, '  ' + PutInCol(elem,4) +
    1061         #                      G2mth.ValEsd(compDict[elem],-0.009,True))
    1062         #         if formula: formula += " "
    1063         #         formula += elem
    1064         #         if compDict[elem] == Z: continue
    1065         #         formula += G2mth.ValEsd(compDict[elem]/Z,-0.009,True)
    1066         #     WriteCIFitem(self.fp,  '\n# Note that Z affects _cell_formula_sum and _weight')
    1067         #     WriteCIFitem(self.fp,  '_cell_formula_units_Z',str(Z))
    1068         #     WriteCIFitem(self.fp,  '_chemical_formula_sum',formula)
    1069         #     WriteCIFitem(self.fp,  '_chemical_formula_weight',
    1070         #                   G2mth.ValEsd(cellmass/Z,-0.09,True))
    1071 
    1072912        def WriteDistances(phasenam,SymOpList,offsetList,symOpList,G2oprList):
    1073913            '''Report bond distances and angles for the CIF
    … …  
    12141054                    WriteCIFitem(self.fp, line)
    12151055
    1216         def WritePhaseInfo(phasenam,hist=None):
     1056        def WritePhaseInfo(phasenam):
    12171057            'Write out the phase information for the selected phase'
    12181058            WriteCIFitem(self.fp, '\n# phase info for '+str(phasenam) + ' follows')
    … …  
    12751115
    12761116            # loop over histogram(s) used in this phase
    1277             if not oneblock and not self.quickmode and not hist:
     1117            if not oneblock and not self.quickmode:
    12781118                # report pointers to the histograms used in this phase
    12791119                histlist = []
    … …  
    13121152#                self.CloseFile()
    13131153#                raise Exception("no export for "+str(phasedict['General']['Type'])+" coordinates implemented")
    1314             # report cell contents
    1315             WriteComposition(self.fp, self.Phases[phasenam], phasenam, self.parmDict)
     1154            keV = None             
     1155            if oneblock: # get xray wavelength
     1156                lamlist = []
     1157                for hist in self.Histograms:
     1158                    if 'X' not in self.Histograms[hist]['Instrument Parameters'][0]['Type'][0]:
     1159                        continue
     1160                    for k in ('Lam','Lam1'):
     1161                        if k in self.Histograms[hist]['Instrument Parameters'][0]:
     1162                            lamlist.append(self.Histograms[hist]['Instrument Parameters'][0][k][0])
     1163                            break
     1164                if len(lamlist) == 1:
     1165                    keV = 12.397639/lamlist[0]
     1166
     1167            # report cell contents                       
     1168            WriteComposition(self.fp, self.Phases[phasenam], phasenam, self.parmDict, self.quickmode, keV)
    13161169            if not self.quickmode and phasedict['General']['Type'] == 'nuclear':      # report distances and angles
    13171170                WriteDistances(phasenam,SymOpList,offsetList,symOpList,G2oprList)
    … …  
    14481301                WriteCIFitem(self.fp, '_pd_meas_2theta_fixed',txt)
    14491302
    1450             # TODO: this will need help from Bob
    1451             #if not oneblock:
    1452             #WriteCIFitem(self.fp, '\n# SCATTERING FACTOR INFO')
    1453             #WriteCIFitem(self.fp, 'loop_  _atom_type_symbol')
    1454             #if histblk['Instrument Parameters'][0]['Type'][1][1] == 'X':
    1455             #    WriteCIFitem(self.fp, '      _atom_type_scat_dispersion_real')
    1456             #    WriteCIFitem(self.fp, '      _atom_type_scat_dispersion_imag')
    1457             #    for lbl in ('a1','a2','a3', 'a4', 'b1', 'b2', 'b3', 'b4', 'c'):
    1458             #        WriteCIFitem(self.fp, '      _atom_type_scat_Cromer_Mann_'+lbl)
    1459             #elif histblk['Instrument Parameters'][0]['Type'][1][1] == 'N':
    1460             #    WriteCIFitem(self.fp, '      _atom_type_scat_length_neutron')
    1461             #WriteCIFitem(self.fp, '      _atom_type_scat_source')
    1462 
    14631303            WriteCIFitem(self.fp, '_pd_proc_ls_background_function',FormatBackground(histblk['Background'],histblk['hId']))
    14641304
    … …  
    15361376                        break
    15371377                    if hklmin is None:
    1538                         hklmin = ref[0:3]
    1539                         hklmax = ref[0:3]
    1540                         dmax = dmin = ref[4]
     1378                        hklmin = copy.copy(ref[0:3])
     1379                        hklmax = copy.copy(ref[0:3])
     1380                    if dmin is None:
     1381                         dmax = dmin = ref[4]
    15411382                    if len(histblk['Reflection Lists'].keys()) > 1:
    15421383                        s = PutInCol(phaseid,2)
    … …  
    16621503                s = "  "
    16631504                if hklmin is None:
    1664                     hklmin = ref[0:3]
    1665                     hklmax = ref[0:3]
     1505                    hklmin = copy.copy(ref[0:3])
     1506                    hklmax = copy.copy(ref[0:3])
    16661507                    dmax = dmin = ref[4]
    16671508                for i,hkl in enumerate(ref[0:3]):
    … …  
    24192260            #============================================================
    24202261            # loop over histograms, exporting them
     2262            # first, get atoms across all phases
     2263            uniqueAtoms = []
     2264            for phasenam in self.Phases:
     2265                cx,ct,cs,cia = self.Phases[phasenam]['General']['AtomPtrs']
     2266                for line in self.Phases[phasenam]['Atoms']:
     2267                    atype = line[ct].strip()
     2268                    if atype.find('-') != -1: atype = atype.split('-')[0]
     2269                    if atype.find('+') != -1: atype = atype.split('+')[0]
     2270                    atype = atype[0].upper()+atype[1:2].lower() # force case conversion
     2271                    if atype == "D" or atype == "D": atype = "H"
     2272                    if atype not in uniqueAtoms:
     2273                        uniqueAtoms.append(atype)
     2274
    24212275            for i in sorted(self.powderDict.keys()):
    24222276                hist = self.powderDict[i]
    … …  
    24342288                    histprm = self.Histograms[hist]["Sample Parameters"]
    24352289                    writeCIFtemplate(histprm,'powder',histprm['InstrName']) # powder template
     2290                   
     2291                    # get xray wavelength and compute & write f' & f''
     2292                    lam = None
     2293                    if 'X' in histblk['Instrument Parameters'][0]['Type'][0]:
     2294                        for k in ('Lam','Lam1'):
     2295                            if k in histblk['Instrument Parameters'][0]:
     2296                                lam = histblk['Instrument Parameters'][0][k][0]
     2297                                break
     2298                    if lam:
     2299                        keV = 12.397639/lam
     2300                        WriteCIFitem(self.fp,'loop_')
     2301                        for item in ('_atom_type_symbol','_atom_type_scat_dispersion_real',
     2302                                         '_atom_type_scat_dispersion_imag','_atom_type_scat_dispersion_source'):
     2303                            WriteCIFitem(self.fp,'     '+item)
     2304                        for elem in HillSortElements(uniqueAtoms):
     2305                            s = '  '
     2306                            s += PutInCol(elem,4)
     2307                            Orbs = G2el.GetXsectionCoeff(elem)
     2308                            FP,FPP,Mu = G2el.FPcalc(Orbs, keV)
     2309                            s += '  {:8.3f}{:8.3f}   https://subversion.xray.aps.anl.gov/pyGSAS/trunk/atmdata.py'.format(FP,FPP)
     2310                            WriteCIFitem(self.fp,s.rstrip())
     2311                        WriteCIFitem(self.fp,'')
    24362312                    WritePowderData(hist)
    24372313            for i in sorted(self.xtalDict.keys()):
Note: See TracChangeset for help on using the changeset viewer.