Changeset 1010 for trunk/exports/G2cif.py
 Timestamp:
 Jul 24, 2013 9:04:01 PM (9 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/exports/G2cif.py
r990 r1010 3 3 ''' 4 4 5 # TODO: need a mechanism for editing of instrument names, bond pub flags, templates,... 6 5 7 import datetime as dt 6 8 import os.path 9 import numpy as np 7 10 import GSASIIIO as G2IO 8 11 #reload(G2IO) … … 110 113 111 114 # include an overall profile rfactor, if there is more than one powder histogram 112 if self.npowder> 1:115 if len(self.powderDict) > 1: 113 116 WriteCIFitem('\n# OVERALL POWDER RFACTOR') 114 117 try: … … 192 195 if s != "": s += '\n' 193 196 s += 'MarchDollase correction' 194 if self.npowder> 1:197 if len(self.powderDict) > 1: 195 198 s += ', histogram '+str(Histogram['hId']+1) 196 199 s += ' coef. = ' + G2mth.ValEsd(self.parmDict[aname],sig) … … 199 202 if s != "": s += '\n' 200 203 s += 'Simple spherical harmonic correction' 201 if self.npowder> 1:204 if len(self.powderDict) > 1: 202 205 s += ', histogram '+str(Histogram['hId']+1) 203 206 s += ' Order = '+str(hapData['Pref.Ori.'][4])+'\n' … … 602 605 # see wrpowdhist.for & wrreflist.for 603 606 604 #refprx = '_refln.' # mm605 refprx = '_refln_' # normal606 607 607 if 'Lam1' in inst: 608 608 ratio = self.parmDict.get('I(L2)/I(L1)',inst['I(L2)/I(L1)'][1]) … … 627 627 PutInCol('2',5)) 628 628 else: 629 lam1 = self.parmDict.get('Lam',inst['Lam'] )629 lam1 = self.parmDict.get('Lam',inst['Lam'][1]) 630 630 slam1 = self.sigDict.get('Lam',0.00009) 631 631 WriteCIFitem('_diffrn_radiation_wavelength',G2mth.ValEsd(lam1,slam1)) … … 669 669 # WriteCIFitem('_refine_ls_R_Fsqd_factor','?') 670 670 671 print histblk['Instrument Parameters'][0]['Type']672 673 671 if histblk['Instrument Parameters'][0]['Type'][1][1] == 'X': 674 672 WriteCIFitem('_diffrn_radiation_probe','xray') … … 709 707 #CALL WRVAL(IUCIF,'_gsas_exptl_extinct_corr_T_max',TEXT(11:20)) 710 708 711 if not oneblock: 712 # instrumental profile terms go here 713 WriteCIFitem('_pd_proc_ls_profile_function','?') 714 715 raise Exception, "testing" 716 phasenam = self.Phases.keys()[0] 717 for key in self.Phases[phasenam]['Histograms']: 718 print key 719 print '' 720 print self.Phases[phasenam]['Histograms'][key] 721 print histblk.keys() 722 for key in histblk: 723 print key,histblk[key] 724 print inst 725 #print self.parmDict.keys() 726 #print self.sigDict.keys() 727 728 #print 'Data' 729 #for item in histblk['Data']: 730 # print item 731 #try: 732 # print key,histblk[key].keys() 733 #except: 734 # print key 735 # print histblk[key] 736 #print 'Background' 737 print histblk['Reflection Lists']['Garnet'][1] 738 for i in range(0,80): 739 for j in [0,1,2,13]: 740 print histblk['Reflection Lists']['Garnet'][i][j], 741 print 742 #print histblk['Reflection Lists']['Garnet'][i][12].shape 743 #print histblk['Reflection Lists']['Garnet'][i][14] 744 #print histblk['Background'][0] 745 #print histblk['Background'][1] 746 import numpy as np 747 refList = np.array([refl[:11] for refl in histblk['Reflection Lists']['Garnet']]) 748 #refList = histblk['Reflection Lists']['Garnet'] 749 Icorr = np.array([refl[13] for refl in histblk['Reflection Lists']['Garnet']]) 750 FO2 = np.array([refl[8] for refl in histblk['Reflection Lists']['Garnet']]) 751 print Icorr 752 I100 = refList.T[8]*Icorr 753 print I100 754 print I100/max(I100) 755 Icorr = np.array([refl[13] for refl in histblk['Reflection Lists']['Garnet']]) * np.array([refl[8] for refl in histblk['Reflection Lists']['Garnet']]) 756 print I100/max(I100) 757 709 if not oneblock: # instrumental profile terms go here 710 WriteCIFitem('_pd_proc_ls_profile_function', 711 FormatInstProfile(histblk["Instrument Parameters"])) 712 713 #refprx = '_refln.' # mm 714 refprx = '_refln_' # normal 758 715 WriteCIFitem('\n# STRUCTURE FACTOR TABLE') 759 WriteCIFitem('loop_' + 760 '\n\t' + refprx + 'index_h' + 716 # compute maximum intensity reflection 717 Imax = 0 718 for phasenam in histblk['Reflection Lists']: 719 scale = self.Phases[phasenam]['Histograms'][histlbl]['Scale'][0] 720 Icorr = np.array([refl[13] for refl in histblk['Reflection Lists'][phasenam]]) 721 FO2 = np.array([refl[8] for refl in histblk['Reflection Lists'][phasenam]]) 722 I100 = scale*FO2*Icorr 723 Imax = max(Imax,max(I100)) 724 725 WriteCIFitem('loop_') 726 if len(histblk['Reflection Lists'].keys()) > 1: 727 WriteCIFitem('\t_pd_refln_phase_id') 728 WriteCIFitem('\t' + refprx + 'index_h' + 761 729 '\n\t' + refprx + 'index_k' + 762 730 '\n\t' + refprx + 'index_l' + 763 '\n\t_pd_refln_wavelength_id' + 764 '\n\t_pd_refln_phase_id' + 765 '\n\t' + refprx + 'status' + 731 # '\n\t_pd_refln_wavelength_id' + 732 # '\n\t' + refprx + 'status' + 766 733 '\n\t' + refprx + 'F_squared_meas' + 767 '\n\t' + refprx + 'F_squared_sigma' +734 # '\n\t' + refprx + 'F_squared_sigma' + 768 735 '\n\t' + refprx + 'F_squared_calc' + 769 736 '\n\t' + refprx + 'phase_calc' + 770 '\n\t_pd_refln_d_spacing' + 771 '\n\t_gsas_i100_meas') 772 773 WriteCIFitem('_reflns_number_total', text) 774 WriteCIFitem('_reflns_limit_h_min', text) 775 WriteCIFitem('_reflns_limit_h_max', text) 776 WriteCIFitem('_reflns_limit_k_min', text) 777 WriteCIFitem('_reflns_limit_k_max', text) 778 WriteCIFitem('_reflns_limit_l_min', text) 779 WriteCIFitem('_reflns_limit_l_max', text) 780 WriteCIFitem('_reflns_d_resolution_high', text) 781 WriteCIFitem('_reflns_d_resolution_low', text) 782 737 '\n\t_pd_refln_d_spacing') 738 if Imax > 0: 739 WriteCIFitem('\t_gsas_i100_meas') 740 741 refcount = 0 742 hklmin = None 743 hklmax = None 744 dmax = None 745 dmin = None 746 for phasenam in histblk['Reflection Lists']: 747 scale = self.Phases[phasenam]['Histograms'][histlbl]['Scale'][0] 748 phaseid = self.Phases[phasenam]['pId'] 749 refcount += len(histblk['Reflection Lists'][phasenam]) 750 for ref in histblk['Reflection Lists'][phasenam]: 751 if hklmin is None: 752 hklmin = ref[0:3] 753 hklmax = ref[0:3] 754 dmax = dmin = ref[4] 755 if len(histblk['Reflection Lists'].keys()) > 1: 756 s = PutInCol(phaseid,2) 757 else: 758 s = "" 759 for i,hkl in enumerate(ref[0:3]): 760 hklmax[i] = max(hkl,hklmax[i]) 761 hklmin[i] = min(hkl,hklmin[i]) 762 s += PutInCol(int(hkl),4) 763 for I in ref[8:10]: 764 s += PutInCol(G2mth.ValEsd(I,0.0009),10) 765 s += PutInCol(G2mth.ValEsd(ref[10],0.9),7) 766 dmax = max(dmax,ref[4]) 767 dmin = min(dmin,ref[4]) 768 s += PutInCol(G2mth.ValEsd(ref[4],0.009),8) 769 if Imax > 0: 770 I100 = 100.*scale*ref[8]*ref[13]/Imax 771 s += PutInCol(G2mth.ValEsd(I100,0.09),6) 772 WriteCIFitem(" "+s) 773 774 WriteCIFitem('_reflns_number_total', str(refcount)) 775 if hklmin is not None and len(histblk['Reflection Lists']) == 1: # hkl range has no meaning with multiple phases 776 WriteCIFitem('_reflns_limit_h_min', str(int(hklmin[0]))) 777 WriteCIFitem('_reflns_limit_h_max', str(int(hklmax[0]))) 778 WriteCIFitem('_reflns_limit_k_min', str(int(hklmin[1]))) 779 WriteCIFitem('_reflns_limit_k_max', str(int(hklmax[1]))) 780 WriteCIFitem('_reflns_limit_l_min', str(int(hklmin[2]))) 781 WriteCIFitem('_reflns_limit_l_max', str(int(hklmax[2]))) 782 if hklmin is not None: 783 WriteCIFitem('_reflns_d_resolution_high', G2mth.ValEsd(dmin,0.009)) 784 WriteCIFitem('_reflns_d_resolution_low', G2mth.ValEsd(dmax,0.0009)) 785 783 786 WriteCIFitem('\n# POWDER DATA TABLE') 784 # is data fixed step? 785 fixedstep = False 786 # are ESDs sqrt(I) 787 countsdata = False 788 zero = 0.01 789 if fixedstep: 790 WriteCIFitem('_pd_meas_2theta_range_min', text) 791 WriteCIFitem('_pd_meas_2theta_range_max', text) 792 WriteCIFitem('_pd_meas_2theta_range_inc', text) 793 # zero correct 794 if zero != 0.0: 795 WriteCIFitem('_pd_proc_2theta_range_min', text) 796 WriteCIFitem('_pd_proc_2theta_range_max', text) 797 WriteCIFitem('_pd_proc_2theta_range_inc', text) 798 WriteCIFitem('loop_' + 799 '\n\t_pd_proc_d_spacing') 800 #'_pd_meas_time_of_flight' 787 # is data fixed step? If the step varies by <0.01% treat as fixed step 788 steps = histblk['Data'][0][1:]  histblk['Data'][0][:1] 789 if abs(max(steps)min(steps)) > abs(max(steps))/10000.: 790 fixedstep = False 791 else: 792 fixedstep = True 793 794 if fixedstep: # and not TOF 795 WriteCIFitem('_pd_meas_2theta_range_min', G2mth.ValEsd(histblk['Data'][0][0],0.00009)) 796 WriteCIFitem('_pd_meas_2theta_range_max', G2mth.ValEsd(histblk['Data'][0][1],0.00009)) 797 WriteCIFitem('_pd_meas_2theta_range_inc', G2mth.ValEsd(steps.sum()/len(steps),0.00009)) 798 # zero correct, if defined 799 zero = None 800 zerolst = histblk['Instrument Parameters'][0].get('Zero') 801 if zerolst: zero = zerolst[1] 802 zero = self.parmDict.get('Zero',zero) 803 # TODO: Bob is zero added or subtracted? 804 if zero: 805 WriteCIFitem('_pd_proc_2theta_range_min', G2mth.ValEsd(histblk['Data'][0][0]zero,0.00009)) 806 WriteCIFitem('_pd_proc_2theta_range_max', G2mth.ValEsd(histblk['Data'][0][1]zero,0.00009)) 807 WriteCIFitem('_pd_proc_2theta_range_inc', G2mth.ValEsd(steps.sum()/len(steps),0.00009)) 808 809 if zero: 810 WriteCIFitem('_pd_proc_number_of_points', str(len(histblk['Data'][0]))) 811 else: 812 WriteCIFitem('_pd_meas_number_of_points', str(len(histblk['Data'][0]))) 813 WriteCIFitem('\nloop_') 814 # WriteCIFitem('\t_pd_proc_d_spacing') # need easy way to get this 801 815 if not fixedstep: 802 if zero != 0.0:816 if zero: 803 817 WriteCIFitem('\t_pd_proc_2theta_corrected') 804 818 else: 805 819 WriteCIFitem('\t_pd_meas_2theta_scan') 806 if countsdata: 807 WriteCIFitem('\t_pd_meas_counts_total') 808 else: 809 WriteCIFitem('\t_pd_meas_intensity_total') 820 # at least for now, always report weights. TODO: Should they be multiplied by weights? 821 #if countsdata: 822 # WriteCIFitem('\t_pd_meas_counts_total') 823 #else: 824 WriteCIFitem('\t_pd_meas_intensity_total') 825 WriteCIFitem('\t_pd_calc_intensity_total') 826 WriteCIFitem('\t_pd_proc_intensity_bkg_calc') 810 827 WriteCIFitem('\t_pd_proc_ls_weight') 811 WriteCIFitem('\t_pd_proc_intensity_bkg_calc') 812 WriteCIFitem('\t_pd_calc_intensity_total') 813 if zero != 0.0: 814 WriteCIFitem('_pd_proc_number_of_points', text) 815 else: 816 WriteCIFitem('_pd_meas_number_of_points', text) 828 # TODO: are data excluded? 829 for x,yobs,yw,ycalc,ybkg in zip(histblk['Data'][0], 830 histblk['Data'][1], 831 histblk['Data'][2], 832 histblk['Data'][3], 833 histblk['Data'][4]): 834 if fixedstep: 835 s = "" 836 else: 837 s = PutInCol(G2mth.ValEsd(xzero,0.00009),10) 838 s += PutInCol(G2mth.ValEsd(yobs,0.00009),14) 839 s += PutInCol(G2mth.ValEsd(ycalc,0.00009),14) 840 s += PutInCol(G2mth.ValEsd(ybkg,0.00009),14) 841 s += PutInCol(G2mth.ValEsd(yw,0.000009),14) 842 WriteCIFitem(" "+s) 817 843 818 844 def WriteSingleXtalData(histlbl): 819 845 histblk = self.Histograms[histlbl] 820 846 print 'TODO: single xtal here data for',histblk["Instrument Parameters"][0]['InstrName'] 821 # see wrreflist.for 822 refprx = '_refln.' # mm847 848 #refprx = '_refln.' # mm 823 849 refprx = '_refln_' # normal 850 851 print histblk.keys() 824 852 825 853 WriteCIFitem('\n# STRUCTURE FACTOR TABLE') … … 827 855 '\n\t' + refprx + 'index_h' + 828 856 '\n\t' + refprx + 'index_k' + 829 '\n\t' + refprx + 'index_l' + 830 '\n\t' + refprx + 'status' + 857 '\n\t' + refprx + 'index_l' + 831 858 '\n\t' + refprx + 'F_squared_meas' + 832 859 '\n\t' + refprx + 'F_squared_sigma' + 833 860 '\n\t' + refprx + 'F_squared_calc' + 834 '\n\t' + refprx + 'phase_calc') 835 WriteCIFitem('_reflns_number_total', text) 836 WriteCIFitem('_reflns_number_observed', text) 837 WriteCIFitem('_reflns_limit_h_min', text) 838 WriteCIFitem('_reflns_limit_h_max', text) 839 WriteCIFitem('_reflns_limit_k_min', text) 840 WriteCIFitem('_reflns_limit_k_max', text) 841 WriteCIFitem('_reflns_limit_l_min', text) 842 WriteCIFitem('_reflns_limit_l_max', text) 843 WriteCIFitem('_reflns_d_resolution_high', text) 844 WriteCIFitem('_reflns_d_resolution_low', text) 861 '\n\t' + refprx + 'phase_calc' + 862 '\n\t_pd_refln_d_spacing' 863 ) 864 865 hklmin = None 866 hklmax = None 867 dmax = None 868 dmin = None 869 refcount = len(histblk['Data']) 870 for ref in histblk['Data']: 871 s = " " 872 if hklmin is None: 873 hklmin = ref[0:3] 874 hklmax = ref[0:3] 875 dmax = dmin = ref[4] 876 for i,hkl in enumerate(ref[0:3]): 877 hklmax[i] = max(hkl,hklmax[i]) 878 hklmin[i] = min(hkl,hklmin[i]) 879 s += PutInCol(int(hkl),4) 880 sig = ref[6] * ref[8] / ref[5] 881 s += PutInCol(G2mth.ValEsd(ref[8],abs(sig/10)),12) 882 s += PutInCol(G2mth.ValEsd(sig,abs(sig)/10.),10) 883 s += PutInCol(G2mth.ValEsd(ref[9],abs(sig/10)),12) 884 s += PutInCol(G2mth.ValEsd(ref[10],0.9),7) 885 dmax = max(dmax,ref[4]) 886 dmin = min(dmin,ref[4]) 887 s += PutInCol(G2mth.ValEsd(ref[4],0.009),8) 888 WriteCIFitem(s) 889 WriteCIFitem('_reflns_number_total', str(refcount)) 890 if hklmin is not None: 891 WriteCIFitem('_reflns_limit_h_min', str(int(hklmin[0]))) 892 WriteCIFitem('_reflns_limit_h_max', str(int(hklmax[0]))) 893 WriteCIFitem('_reflns_limit_k_min', str(int(hklmin[1]))) 894 WriteCIFitem('_reflns_limit_k_max', str(int(hklmax[1]))) 895 WriteCIFitem('_reflns_limit_l_min', str(int(hklmin[2]))) 896 WriteCIFitem('_reflns_limit_l_max', str(int(hklmax[2]))) 897 WriteCIFitem('_reflns_d_resolution_high', G2mth.ValEsd(dmin,0.009)) 898 WriteCIFitem('_reflns_d_resolution_low', G2mth.ValEsd(dmax,0.0009)) 845 899 846 900 #============================================================ … … 864 918 # count phases, powder and single crystal histograms 865 919 self.nphase = len(self.Phases) 866 self. npowder = 0867 self. nsingle = 0920 self.powderDict = {} 921 self.xtalDict = {} 868 922 for hist in self.Histograms: 923 i = self.Histograms[hist]['hId'] 869 924 if hist.startswith("PWDR"): 870 self. npowder += 1925 self.powderDict[i] = hist 871 926 elif hist.startswith("HKLF"): 872 self. nsingle += 1927 self.xtalDict[i] = hist 873 928 # is there anything to export? 874 if self.nphase + self.npowder + self.nsingle== 0:929 if self.nphase + len(self.powderDict) + len(self.xtalDict) == 0: 875 930 self.G2frame.ErrorDialog( 876 931 'Empty project', … … 890 945 self.quickmode = False 891 946 phasenam = phasenum = None # include all phases 892 if mode != "full" or self.npowder + self.nsingle== 0:947 if mode != "full" or len(self.powderDict) + len(self.xtalDict) == 0: 893 948 self.quickmode = True 894 949 oneblock = True … … 906 961 elif self.nphase > 1: 907 962 oneblock = False 908 elif self.npowder + self.nsingle> 1:963 elif len(self.powderDict) + len(self.xtalDict) > 1: 909 964 oneblock = False 910 965 else: # one phase, one dataset, Full CIF … … 1048 1103 WriteCIFitem(loopprefix,datablockidDict[phasenam]) 1049 1104 # loop over data blocks 1050 if self.npowder + self.nsingle> 1:1105 if len(self.powderDict) + len(self.xtalDict) > 1: 1051 1106 loopprefix = '' 1052 1107 WriteCIFitem('loop_ _pd_block_diffractogram_id') 1053 1108 else: 1054 1109 loopprefix = '_pd_block_diffractogram_id' 1055 for hist in self.Histograms: 1110 for i in sorted(self.powderDict.keys()): 1111 hist = self.powderDict[i] 1056 1112 histblk = self.Histograms[hist] 1057 if hist.startswith("PWDR"): 1058 instnam = histblk["Sample Parameters"]['InstrName'] 1059 elif hist.startswith("HKLF"): 1060 instnam = histblk["Instrument Parameters"][0]['InstrName'] 1113 instnam = histblk["Sample Parameters"]['InstrName'] 1114 instnam = instnam.replace(' ','') 1115 i = histblk['hId'] 1116 datablockidDict[hist] = (str(self.CIFdate) + "" + str(self.CIFname) + "" + 1117 str(self.shortauthorname) + "" + 1118 instnam + "_hist_"+str(i)) 1119 WriteCIFitem(loopprefix,datablockidDict[hist]) 1120 for i in sorted(self.xtalDict.keys()): 1121 hist = self.xtalDict[i] 1122 histblk = self.Histograms[hist] 1123 instnam = histblk["Instrument Parameters"][0]['InstrName'] 1061 1124 instnam = instnam.replace(' ','') 1062 1125 i = histblk['hId'] … … 1093 1156 #============================================================ 1094 1157 # loop over histograms, exporting them 1095 for hist in self.Histograms: 1158 for i in sorted(self.powderDict.keys()): 1159 hist = self.powderDict[i] 1096 1160 histblk = self.Histograms[hist] 1097 i = histblk['hId']1098 1161 if hist.startswith("PWDR"): 1099 1162 WriteCIFitem('\ndata_'+self.CIFname+"_pwd_"+str(i)) … … 1107 1170 WritePowderTemplate() 1108 1171 WritePowderData(hist) 1109 elif hist.startswith("HKLF"): 1172 for i in sorted(self.xtalDict.keys()): 1173 hist = self.xtalDict[i] 1174 histblk = self.Histograms[hist] 1175 if hist.startswith("HKLF"): 1110 1176 WriteCIFitem('\ndata_'+self.CIFname+"_sx_"+str(i)) 1111 1177 #instnam = histblk["Instrument Parameters"][0]['InstrName']
Note: See TracChangeset
for help on using the changeset viewer.