Changeset 5080


Ignore:
Timestamp:
Nov 13, 2021 3:45:30 PM (7 months ago)
Author:
toby
Message:

mmCIF export: histogram info

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/exports/G2export_CIF.py

    r5079 r5080  
    29222922            WriteCIFitem(self.fp, '   _pd_proc_intensity_bkg_calc')
    29232923            WriteCIFitem(self.fp, '   _pd_proc_ls_weight')
     2924            maxY = max(histblk['Data'][1].max(),histblk['Data'][3].max())
     2925            if maxY < 0: maxY *= -10 # this should never happen, but...
     2926            ndec = max(0,10-int(np.log10(maxY))-1) # 10 sig figs should be enough
     2927            maxSU = histblk['Data'][2].max()
     2928            if maxSU < 0: maxSU *= -1 # this should never happen, but...
     2929            ndecSU = max(0,8-int(np.log10(maxSU))-1) # 8 sig figs should be enough
     2930            lowlim,highlim = histblk['Limits'][1]
     2931
     2932            if DEBUG:
     2933                print('DEBUG: skipping profile list')
     2934            else:
     2935                for x,yobs,yw,ycalc,ybkg in zip(histblk['Data'][0].data,        #get the data from these masked arrays
     2936                                                histblk['Data'][1].data,
     2937                                                histblk['Data'][2].data,
     2938                                                histblk['Data'][3].data,
     2939                                                histblk['Data'][4].data):
     2940                    if lowlim <= x <= highlim:
     2941                        pass
     2942                    else:
     2943                        yw = 0.0 # show the point is not in use
     2944
     2945                    if fixedstep:
     2946                        s = ""
     2947                    elif zero:
     2948                        s = PutInCol(G2mth.ValEsd(x-zero,-0.00009),10)
     2949                    else:
     2950                        s = PutInCol(G2mth.ValEsd(x,-0.00009),10)
     2951                    s += PutInCol(Yfmt(ndec,yobs),12)
     2952                    s += PutInCol(Yfmt(ndec,ycalc),12)
     2953                    s += PutInCol(Yfmt(ndec,ybkg),11)
     2954                    s += PutInCol(Yfmt(ndecSU,yw),9)
     2955                    WriteCIFitem(self.fp, "  "+s)
     2956                   
     2957        def WritePowderDataMM(histlbl,seq=False):
     2958            'Write out the selected powder diffraction histogram info'
     2959            histblk = self.Histograms[histlbl]
     2960            inst = histblk['Instrument Parameters'][0]
     2961            hId = histblk['hId']
     2962            pfx = ':' + str(hId) + ':'
     2963
     2964            if 'Lam1' in inst:
     2965                ratio = self.parmDict.get('I(L2)/I(L1)',inst['I(L2)/I(L1)'][1])
     2966                sratio = self.sigDict.get('I(L2)/I(L1)',-0.0009)
     2967                lam1 = self.parmDict.get('Lam1',inst['Lam1'][1])
     2968                slam1 = self.sigDict.get('Lam1',-0.00009)
     2969                lam2 = self.parmDict.get('Lam2',inst['Lam2'][1])
     2970                slam2 = self.sigDict.get('Lam2',-0.00009)
     2971                # always assume Ka1 & Ka2 if two wavelengths are present
     2972                WriteCIFitem(self.fp, '_diffrn_radiation.type','K\\a~1,2~')
     2973                WriteCIFitem(self.fp, 'loop_' +
     2974                             '\n   _diffrn_radiation_wavelength.wavelength' +
     2975                             '\n   _diffrn_radiation_wavelength.wt' +
     2976                             '\n   _diffrn_radiation_wavelength.id')
     2977                WriteCIFitem(self.fp, '  ' + PutInCol(G2mth.ValEsd(lam1,slam1),15)+
     2978                             PutInCol('1.0',15) +
     2979                             PutInCol('1',5))
     2980                WriteCIFitem(self.fp, '  ' + PutInCol(G2mth.ValEsd(lam2,slam2),15)+
     2981                             PutInCol(G2mth.ValEsd(ratio,sratio),15)+
     2982                             PutInCol('2',5))
     2983            elif 'Lam' in inst:
     2984                lam1 = self.parmDict.get('Lam',inst['Lam'][1])
     2985                slam1 = self.sigDict.get('Lam',-0.00009)
     2986                WriteCIFitem(self.fp, '_diffrn_radiation_wavelength.wavelength',G2mth.ValEsd(lam1,slam1))
     2987
     2988            if not oneblock:
     2989                if seq:
     2990                    pass
     2991                elif not phasebyhistDict.get(histlbl):
     2992                    WriteCIFitem(self.fp, '\n# No phases associated with this data set')
     2993                else:
     2994                    WriteCIFitem(self.fp, '\n# PHASE TABLE')
     2995                    WriteCIFitem(self.fp, 'loop_' +
     2996                                 '\n   _pd_phase_id' +
     2997                                 '\n   _pd_phase_block_id' +
     2998                                 '\n   _pd_phase_mass_%')
     2999                    wtFrSum = 0.
     3000                    for phasenam in phasebyhistDict.get(histlbl):
     3001                        hapData = self.Phases[phasenam]['Histograms'][histlbl]
     3002                        General = self.Phases[phasenam]['General']
     3003                        wtFrSum += hapData['Scale'][0]*General['Mass']
     3004
     3005                    for phasenam in phasebyhistDict.get(histlbl):
     3006                        hapData = self.Phases[phasenam]['Histograms'][histlbl]
     3007                        General = self.Phases[phasenam]['General']
     3008                        wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum
     3009                        pfx = str(self.Phases[phasenam]['pId'])+':'+str(hId)+':'
     3010                        if pfx+'Scale' in self.sigDict:
     3011                            sig = self.sigDict[pfx+'Scale']*wtFr/hapData['Scale'][0]
     3012                        else:
     3013                            sig = -0.0001
     3014                        WriteCIFitem(self.fp,
     3015                            '  '+
     3016                            str(self.Phases[phasenam]['pId']) +
     3017                            '  '+datablockidDict[phasenam]+
     3018                            '  '+G2mth.ValEsd(wtFr,sig)
     3019                            )
     3020                    WriteCIFitem(self.fp, 'loop_' +
     3021                                 '\n   _gsas_proc_phase_R_F_factor' +
     3022                                 '\n   _gsas_proc_phase_R_Fsqd_factor' +
     3023                                 '\n   _gsas_proc_phase_id' +
     3024                                 '\n   _gsas_proc_phase_block_id')
     3025                    for phasenam in phasebyhistDict.get(histlbl):
     3026                        pfx = str(self.Phases[phasenam]['pId'])+':'+str(hId)+':'
     3027                        WriteCIFitem(self.fp,
     3028                            '  '+
     3029                            '  '+G2mth.ValEsd(histblk[pfx+'Rf']/100.,-.00009) +
     3030                            '  '+G2mth.ValEsd(histblk[pfx+'Rf^2']/100.,-.00009)+
     3031                            '  '+str(self.Phases[phasenam]['pId'])+
     3032                            '  '+datablockidDict[phasenam]
     3033                            )
     3034            elif len(self.Phases) == 1:
     3035                # single phase in this histogram
     3036                # get the phase number here
     3037                pId = self.Phases[list(self.Phases.keys())[0]]['pId']
     3038                pfx = str(pId)+':'+str(hId)+':'
     3039                WriteCIFitem(self.fp, '_refine.ls_R_factor_all    ','%.5f'%(histblk[pfx+'Rf']/100.))
     3040                WriteCIFitem(self.fp, '_refine_ls.R_Fsqd_factor   ','%.5f'%(histblk[pfx+'Rf^2']/100.))
     3041               
     3042            try:
     3043                WriteCIFitem(self.fp, '_pd_proc_ls_prof_R_factor   ','%.5f'%(histblk['R']/100.))
     3044                WriteCIFitem(self.fp, '_pd_proc_ls_prof_wR_factor  ','%.5f'%(histblk['wR']/100.))
     3045                WriteCIFitem(self.fp, '_gsas_proc_ls_prof_R_B_factor ','%.5f'%(histblk['Rb']/100.))
     3046                WriteCIFitem(self.fp, '_gsas_proc_ls_prof_wR_B_factor','%.5f'%(histblk['wRb']/100.))
     3047                WriteCIFitem(self.fp, '_pd_proc_ls_prof_wR_expected','%.5f'%(histblk['wRmin']/100.))
     3048            except KeyError:
     3049                pass
     3050
     3051            if histblk['Instrument Parameters'][0]['Type'][1][1] == 'X':
     3052                WriteCIFitem(self.fp, '_diffrn_radiation.probe','x-ray')
     3053                pola = histblk['Instrument Parameters'][0].get('Polariz.')
     3054                if pola:
     3055                    pfx = ':' + str(hId) + ':'
     3056                    sig = self.sigDict.get(pfx+'Polariz.',-0.0009)
     3057                    txt = G2mth.ValEsd(pola[1],sig)
     3058                    WriteCIFitem(self.fp, '_diffrn_radiation.polarisn_ratio',txt)
     3059            elif histblk['Instrument Parameters'][0]['Type'][1][1] == 'N':
     3060                WriteCIFitem(self.fp, '_diffrn_radiation.probe','neutron')
     3061            if 'T' in inst['Type'][0]:
     3062                txt = G2mth.ValEsd(inst['2-theta'][0],-0.009)
     3063                WriteCIFitem(self.fp, '_pd_meas_2theta_fixed',txt)
     3064
     3065            WriteCIFitem(self.fp, '_pd_proc_ls_background_function',FormatBackground(histblk['Background'],histblk['hId']))
     3066
     3067            # TODO: this will need help from Bob
     3068            #WriteCIFitem(self.fp, '_exptl_absorpt_process_details','?')
     3069            #WriteCIFitem(self.fp, '_exptl_absorpt_correction_T_min','?')
     3070            #WriteCIFitem(self.fp, '_exptl_absorpt_correction_T_max','?')
     3071            #C extinction
     3072            #WRITE(IUCIF,'(A)') '# Extinction correction'
     3073            #CALL WRVAL(IUCIF,'_gsas_exptl_extinct_corr_T_min',TEXT(1:10))
     3074            #CALL WRVAL(IUCIF,'_gsas_exptl_extinct_corr_T_max',TEXT(11:20))
     3075
     3076            # code removed because it is causing duplication in histogram block 1/26/19 BHT
     3077            #if not oneblock:                 # instrumental profile terms go here
     3078            #    WriteCIFitem(self.fp, '_pd_proc_ls_profile_function',
     3079            #        FormatInstProfile(histblk["Instrument Parameters"],histblk['hId']))
     3080
     3081            # data collection parameters for the powder dataset
     3082
     3083            temperature = histblk['Sample Parameters'].get('Temperature') # G2 uses K
     3084            if not temperature:
     3085                T = '?'
     3086            else:
     3087                T = G2mth.ValEsd(temperature,-0.009,True) # CIF uses K
     3088            WriteCIFitem(self.fp, '_diffrn_ambient.temp',T)
     3089
     3090            pressure = histblk['Sample Parameters'].get('Pressure') #G2 uses mega-Pascal
     3091            if not pressure:
     3092                P = '?'
     3093            else:
     3094                P = G2mth.ValEsd(pressure*1000,-0.09,True) # CIF uses kilopascal (G2 Mpa)
     3095            WriteCIFitem(self.fp, '_diffrn_ambient.pressure',P)
     3096
     3097            WriteCIFitem(self.fp, '\n# STRUCTURE FACTOR TABLE')
     3098            # compute maximum intensity reflection
     3099            Imax = 0
     3100            phaselist = []
     3101            for phasenam in histblk['Reflection Lists']:
     3102                try:
     3103                    scale = self.Phases[phasenam]['Histograms'][histlbl]['Scale'][0]
     3104                except KeyError: # reflection table from removed phase?
     3105                    continue
     3106                phaselist.append(phasenam)
     3107                refList = np.asarray(histblk['Reflection Lists'][phasenam]['RefList'])
     3108                I100 = scale*refList.T[8]*refList.T[11]
     3109                #Icorr = np.array([refl[13] for refl in histblk['Reflection Lists'][phasenam]])[0]
     3110                #FO2 = np.array([refl[8] for refl in histblk['Reflection Lists'][phasenam]])
     3111                #I100 = scale*FO2*Icorr
     3112#                Imax = max(Imax,max(I100))
     3113
     3114            WriteCIFitem(self.fp, 'loop_')
     3115            #refprx = '_refln.' # mm
     3116            if len(phaselist) > 1:
     3117                WriteCIFitem(self.fp, '   _pd_refln_phase_id')
     3118            WriteCIFitem(self.fp, '   _refln.index_h' +
     3119                         '\n   _refln.index_k' +
     3120                         '\n   _refln.index_l' +
     3121                         '\n   _refln.F_squared_meas' +
     3122                         '\n   _refln.F_squared_calc' +
     3123                         '\n   _refln.phase_calc' +
     3124                         '\n   _refln.d_spacing' +
     3125                         '\n   _refln.status' +
     3126                         '\n   _refln.crystal_id' +
     3127                         '\n   _refln.wavelength_id' +
     3128                         '\n   _refln.scale_group_code' +
     3129                         '\n   _refln.F_squared_sigma')
     3130           
     3131#            if Imax > 0:
     3132#                WriteCIFitem(self.fp, '   _gsas_i100_meas')
     3133
     3134            refcount = 0
     3135            hklmin = None
     3136            hklmax = None
     3137            dmax = None
     3138            dmin = None
     3139            for phasenam in phaselist:
     3140                scale = self.Phases[phasenam]['Histograms'][histlbl]['Scale'][0]
     3141                phaseid = self.Phases[phasenam]['pId']
     3142                refcount += len(histblk['Reflection Lists'][phasenam]['RefList'])
     3143                refList = np.asarray(histblk['Reflection Lists'][phasenam]['RefList'])
     3144                I100 = scale*refList.T[8]*refList.T[11]
     3145                for j,ref in enumerate(histblk['Reflection Lists'][phasenam]['RefList']):
     3146                    if DEBUG:
     3147                        print('DEBUG: skipping reflection list')
     3148                        break
     3149                    if hklmin is None:
     3150                        hklmin = copy.copy(ref[0:3])
     3151                        hklmax = copy.copy(ref[0:3])
     3152                    if dmin is None:
     3153                         dmax = dmin = ref[4]
     3154                    if len(phaselist) > 1:
     3155                        s = PutInCol(phaseid,2)
     3156                    else:
     3157                        s = ""
     3158                    for i,hkl in enumerate(ref[0:3]):
     3159                        hklmax[i] = max(hkl,hklmax[i])
     3160                        hklmin[i] = min(hkl,hklmin[i])
     3161                        s += PutInCol(int(hkl),4)
     3162                    for I in ref[8:10]:
     3163                        s += PutInCol(G2mth.ValEsd(I,-0.0009),14)
     3164                    s += PutInCol(G2mth.ValEsd(ref[10],-0.9),7)
     3165                    dmax = max(dmax,ref[4])
     3166                    dmin = min(dmin,ref[4])
     3167                    s += PutInCol(G2mth.ValEsd(ref[4],-0.00009),8)
     3168#                    if Imax > 0:
     3169#                        s += PutInCol(G2mth.ValEsd(100.*I100[j]/Imax,-0.09),6)
     3170                    s += PutInCol('o',2)
     3171                    s += PutInCol('1',2)
     3172                    s += PutInCol('1',2)
     3173                    s += PutInCol('1',2)
     3174                    s += PutInCol('.',2)
     3175                    WriteCIFitem(self.fp, "  "+s)
     3176
     3177            # Write reflection statistics
     3178            WriteCIFitem(self.fp, '_diffrn_reflns.number', str(refcount))
     3179            if hklmin is not None and len(phaselist) == 1: # hkl range has no meaning with multiple phases
     3180                WriteCIFitem(self.fp, '_diffrn_reflns.limit_h_min', str(int(hklmin[0])))
     3181                WriteCIFitem(self.fp, '_diffrn_reflns.limit_h_max', str(int(hklmax[0])))
     3182                WriteCIFitem(self.fp, '_diffrn_reflns.limit_k_min', str(int(hklmin[1])))
     3183                WriteCIFitem(self.fp, '_diffrn_reflns.limit_k_max', str(int(hklmax[1])))
     3184                WriteCIFitem(self.fp, '_diffrn_reflns.limit_l_min', str(int(hklmin[2])))
     3185                WriteCIFitem(self.fp, '_diffrn_reflns.limit_l_max', str(int(hklmax[2])))
     3186            if hklmin is not None:
     3187                WriteCIFitem(self.fp, '_reflns.d_resolution_low  ', G2mth.ValEsd(dmax,-0.009))
     3188                WriteCIFitem(self.fp, '_reflns.d_resolution_high ', G2mth.ValEsd(dmin,-0.009))
     3189               
     3190            WriteCIFitem(self.fp, '\n# POWDER DATA TABLE')
     3191
     3192            # is data fixed step? If the step varies by <0.01% treat as fixed step
     3193            fixedstep = False
     3194            zero = None
     3195            WriteCIFitem(self.fp, '_refine.pdbx_pd_meas_number_of_points', str(len(histblk['Data'][0])))
     3196            WriteCIFitem(self.fp, '\nloop_')
     3197            #            WriteCIFitem(self.fp, '   _pd_proc_d_spacing') # need easy way to get this
     3198            if 'T' in inst['Type'][0]: # and not TOF
     3199                WriteCIFitem(self.fp, '   _pd_meas_time_of_flight')
     3200            else:
     3201                WriteCIFitem(self.fp, '   _pdbx_powder_data.pd_meas_2theta_scan')
     3202            # at least for now, always report weights.
     3203            #if countsdata:
     3204            #    WriteCIFitem(self.fp, '   _pd_meas_counts_total')
     3205            #else:
     3206            WriteCIFitem(self.fp, '   _pdbx_powder_data.pd_meas_intensity_total')
     3207            WriteCIFitem(self.fp, '   _pdbx_powder_data.pd_calc_intensity_total')
     3208            WriteCIFitem(self.fp, '   _pdbx_powder_data.pd_proc_intensity_bkg_calc')
     3209            WriteCIFitem(self.fp, '   _pdbx_powder_data.pd_proc_ls_weight')
    29243210            maxY = max(histblk['Data'][1].max(),histblk['Data'][3].max())
    29253211            if maxY < 0: maxY *= -10 # this should never happen, but...
     
    34403726        elif histOnly: #====Histogram only CIF ================================
    34413727            print('Writing CIF output to file '+self.filename)
     3728            MM = False
     3729            for p in self.Phases:
     3730                if self.Phases[p]['General']['Type'] == 'macromolecular':
     3731                    MM = True
     3732                    break
    34423733            hist = histOnly
    34433734            #histname = histOnly.replace(' ','')
     
    34503741            #phasenam = self.Phases.keys()[0]
    34513742            WriteCIFitem(self.fp, 'data_'+self.CIFname)
    3452             if hist.startswith("PWDR"):
     3743            if hist.startswith("PWDR") and MM:
     3744                WritePowderDataMM(hist)
     3745            elif hist.startswith("PWDR"):
    34533746                WritePowderData(hist)
    34543747            elif hist.startswith("HKLF"):
     
    36123905                return
    36133906        self.cifdefs.Destroy()
     3907        MM = False
     3908        for p in self.Phases:
     3909            if self.Phases[p]['General']['Type'] == 'macromolecular':
     3910                MM = True
     3911                break
    36143912        #======================================================================
    36153913        #---- Start writing the CIF - single block
     
    36213919            self.CIFname = hist[5:40].replace(' ','')
    36223920            WriteCIFitem(self.fp, 'data_'+self.CIFname)
    3623             if hist.startswith("PWDR"):
     3921            if hist.startswith("PWDR") and MM:
     3922                WritePowderDataMM(hist)
     3923            elif hist.startswith("PWDR"):
    36243924                WritePowderData(hist)
    36253925            elif hist.startswith("HKLF"):
     
    36643964                histblk = self.Histograms[hist]["Sample Parameters"]
    36653965                writeCIFtemplate(histblk,'powder',histblk['InstrName']) # write powder template
    3666                 WritePowderData(hist)
     3966                if hist.startswith("PWDR") and MM:
     3967                    WritePowderDataMM(hist)
     3968                else:
     3969                    WritePowderData(hist)
    36673970            elif hist.startswith("HKLF"):
    36683971                histprm = self.Histograms[hist]["Instrument Parameters"][0]
     
    40334336                                WriteCIFitem(self.fp,s.rstrip())
    40344337                            WriteCIFitem(self.fp,'')
    4035                         WritePowderData(hist)
     4338                        if MM:
     4339                            WritePowderDataMM(hist)
     4340                        else:
     4341                            WritePowderData(hist)
    40364342                for i in sorted(self.xtalDict.keys()):
    40374343                    hist = self.xtalDict[i]
Note: See TracChangeset for help on using the changeset viewer.