Changeset 4402 for trunk/exports/G2export_CIF.py
 Timestamp:
 Apr 15, 2020 5:13:34 PM (3 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/exports/G2export_CIF.py
r4361 r4402 58 58 import GSASIIstrMain as G2stMn 59 59 import GSASIImapvars as G2mv 60 import GSASIIElem as G2el 60 61 61 62 DEBUG = False #True to skip printing of reflection/powder profile lists … … 376 377 377 378 # Refactored over here to allow access by GSASIIscriptable.py 378 def WriteComposition(fp, phasedict, phasenam, parmDict ):379 def WriteComposition(fp, phasedict, phasenam, parmDict, quickmode=True, keV=None): 379 380 '''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. 381 386 ''' 382 387 General = phasedict['General'] … … 421 426 General['cellZ'] = Z # save it 422 427 423 # when scattering factors are included in the CIF, this needs to be 424 # added to the loop here but only in the oneblock 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 445 432 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 446 452 formula = '' 447 453 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 'BWLS' 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()) 450 494 if formula: formula += " " 451 495 formula += elem … … 866 910 labellist.append(lbl) 867 911 868 # Factored out to above by Jack O'Donnell, so it869 # could be accessed by GSASIIscriptable.py870 # def WriteAtomsNuclear(phasenam):871 # 'Write atom positions to CIF'872 # phasedict = self.Phases[phasenam] # pointer to current phase info873 # General = phasedict['General']874 # cx,ct,cs,cia = General['AtomPtrs']875 # Atoms = phasedict['Atoms']876 # cfrac = cx+3877 # 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 # break882 # else:883 # WriteCIFitem(self.fp, '\n# PHASE HAS NO ATOMS!')884 # return885 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 atoms905 # naniso = 0906 # for i,at in enumerate(Atoms):907 # if phasedict['General']['Type'] == 'macromolecular':908 # label = '%s_%s_%s_%s'%(at[ct1],at[ct3],at[ct4],at[ct2])909 # s = PutInCol(MakeUniqueLabel(label,self.labellist),15) # label910 # else:911 # s = PutInCol(MakeUniqueLabel(at[ct1],self.labellist),6) # label912 # 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) # type915 # if at[cia] == 'I':916 # adp = 'Uiso '917 # else:918 # adp = 'Uani '919 # naniso += 1920 # # compute Uequiv crudely921 # # 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.0924 # 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 = 11930 # sigdig = 0.00009931 # else:932 # dig = 10933 # sigdig = 0.009934 # if j == cia:935 # s += adp936 # 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 = var941 # if j == cia+1 and adp == 'Uani ':942 # val = t/3.943 # sig = sigdig944 # 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: return952 # # now loop over aniso atoms953 # 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': continue961 # s = PutInCol(self.labellist[i],6) # label962 # for j in (2,3,4,5,6,7):963 # sigdig = 0.0009964 # 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 others972 # are alphabetical).973 974 # :params list elmlist: a list of element strings975 976 # :returns: a sorted list of element strings977 # '''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'Donnell987 # so that it can be accessed by GSASIIscriptable988 989 # def WriteComposition(phasenam):990 # '''determine the composition for the unit cell, crudely determine Z and991 # then compute the composition in formula units992 # '''993 # phasedict = self.Phases[phasenam] # pointer to current phase info994 # 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+31000 # cmult = cs+11001 # compDict = {} # combines H,D & T1002 # sitemultlist = []1003 # massDict = dict(zip(General['AtomTypes'],General['AtomMass']))1004 # cellmass = 01005 # 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 conversion1010 # 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 # return1018 # cellmass += massDict[at[ct]]*mult*fval1019 # compDict[atype] = compDict.get(atype,0.0) + mult*fval1020 # 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 user1023 # Z = 11024 # if not sitemultlist:1025 # General['cellZ'] = 11026 # return1027 # for i in range(2,min(sitemultlist)+1):1028 # for m in sitemultlist:1029 # if m % i != 0:1030 # break1031 # else:1032 # Z = i1033 # General['cellZ'] = Z # save it1034 1035 # # when scattering factors are included in the CIF, this needs to be1036 # # added to the loop here but only in the oneblock case.1037 # # For multiblock CIFs, scattering factors go in the histogram1038 # # blocks (for all atoms in all appropriate phases)  an example?:1039 #loop_1040 # _at# om_type_symbol1041 # _at# om_type_description1042 # _at# om_type_scat_dispersion_real1043 # _at# om_type_scat_dispersion_imag1044 # _at# om_type_scat_source1045 # 'C'# 'C' 0.0033 0.00161046 # # 'International Tables Vol C Tables 4.2.6.8 and 6.1.1.4'1047 # 'H'# 'H' 0.0000 0.00001048 # # 'International Tables Vol C Tables 4.2.6.8 and 6.1.1.4'1049 # 'P'# 'P' 0.1023 0.09421050 # # 'International Tables Vol C Tables 4.2.6.8 and 6.1.1.4'1051 # 'Cl# ' 'Cl' 0.1484 0.15851052 # # 'International Tables Vol C Tables 4.2.6.8 and 6.1.1.4'1053 # 'Cu# ' 'Cu' 0.3201 1.26511054 # # 'International Tables Vol C Tables 4.2.6.8 and 6.1.1.4'1055 1056 # #if oneblock: # add scattering factors for current phase here1057 # 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 += elem1064 # if compDict[elem] == Z: continue1065 # 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 1072 912 def WriteDistances(phasenam,SymOpList,offsetList,symOpList,G2oprList): 1073 913 '''Report bond distances and angles for the CIF … … 1214 1054 WriteCIFitem(self.fp, line) 1215 1055 1216 def WritePhaseInfo(phasenam ,hist=None):1056 def WritePhaseInfo(phasenam): 1217 1057 'Write out the phase information for the selected phase' 1218 1058 WriteCIFitem(self.fp, '\n# phase info for '+str(phasenam) + ' follows') … … 1275 1115 1276 1116 # 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: 1278 1118 # report pointers to the histograms used in this phase 1279 1119 histlist = [] … … 1312 1152 # self.CloseFile() 1313 1153 # 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) 1316 1169 if not self.quickmode and phasedict['General']['Type'] == 'nuclear': # report distances and angles 1317 1170 WriteDistances(phasenam,SymOpList,offsetList,symOpList,G2oprList) … … 1448 1301 WriteCIFitem(self.fp, '_pd_meas_2theta_fixed',txt) 1449 1302 1450 # TODO: this will need help from Bob1451 #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 1463 1303 WriteCIFitem(self.fp, '_pd_proc_ls_background_function',FormatBackground(histblk['Background'],histblk['hId'])) 1464 1304 … … 1536 1376 break 1537 1377 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] 1541 1382 if len(histblk['Reflection Lists'].keys()) > 1: 1542 1383 s = PutInCol(phaseid,2) … … 1662 1503 s = " " 1663 1504 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]) 1666 1507 dmax = dmin = ref[4] 1667 1508 for i,hkl in enumerate(ref[0:3]): … … 2419 2260 #============================================================ 2420 2261 # 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 2421 2275 for i in sorted(self.powderDict.keys()): 2422 2276 hist = self.powderDict[i] … … 2434 2288 histprm = self.Histograms[hist]["Sample Parameters"] 2435 2289 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,'') 2436 2312 WritePowderData(hist) 2437 2313 for i in sorted(self.xtalDict.keys()):
Note: See TracChangeset
for help on using the changeset viewer.