Changeset 5090


Ignore:
Timestamp:
Nov 18, 2021 4:39:00 PM (6 months ago)
Author:
vondreele
Message:

add MTZ exporter for both F2 and F options. The F style MTZ file can be read by coot.

Location:
trunk
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIdataGUI.py

    r5082 r5090  
    5353    pass
    5454import scipy as sp
    55 #import scipy.optimize as so
     55import struct as st
    5656try:
    5757    import wx
     
    8787import GSASIIseqGUI as G2seq
    8888import GSASIIddataGUI as G2ddG
     89import GSASIIspc as G2spc
    8990
    9091try:
     
    29112912                self.Bind(wx.EVT_MENU, obj.Exporter, id=item.GetId())
    29122913                self.ExportLookup[item.GetId()] = typ # lookup table for submenu item
    2913             for lbl,submenu in (('Phase',seqPhasemenu),
    2914                         ('Powder',seqHistmenu),
    2915                         ):
     2914            for lbl,submenu in (('Phase',seqPhasemenu),('Powder',seqHistmenu),):
    29162915                if lbl.lower() in obj.exporttype:
    29172916                    try:
     
    29732972        self.ExportHKL.append(item)
    29742973        self.Bind(wx.EVT_MENU, self.OnExportHKL, id=item.GetId())
     2974       
     2975        item = parent.Append(wx.ID_ANY,'Export MTZ file...','')
     2976        self.ExportMTZ.append(item)
     2977        self.Bind(wx.EVT_MENU, self.OnExportMTZ, id=item.GetId())
    29752978
    29762979        item = parent.Append(wx.ID_ANY,'Export PDF...','Select PDF item to enable')
     
    30783081        self.ExportPhase = []
    30793082        self.ExportCIF = []
     3083        self.ExportMTZ = []
    30803084        #
    30813085        self.MacroStatusList = []  # logging
     
    46134617        sys.exit()
    46144618       
     4619    def OnExportMTZ(self,event):
     4620        ''' exports MTZ file from macromoleculat Reflection Lists in multiple histograms
     4621        '''
     4622       
     4623        Histograms,Phases = self.GetUsedHistogramsAndPhasesfromTree()
     4624        MTZok = False
     4625        for phase in Phases:
     4626            Phase = Phases[phase]
     4627            General = Phase['General']
     4628            if 'macro' in General['Type']:    #really ought to be just one
     4629                SGData = General['SGData']
     4630                Cell = General['Cell'][1:7]
     4631                MTZok = True
     4632                break
     4633        if not MTZok:
     4634            wx.MessageBox('No macromolecular phase found; no MTZ file created',caption='No macromolecular phase',
     4635                style=wx.ICON_EXCLAMATION)
     4636            return
     4637        header = 'Select histograms for MTZ file:'
     4638        histList = [item for item in Histograms]
     4639        od = {'label_1':'Export Fo','value_1':False,}
     4640        dlg = G2G.G2MultiChoiceDialog(self,header,'Add to MTZ',histList,extraOpts=od)
     4641        result = None
     4642        if dlg.ShowModal() == wx.ID_OK:
     4643            result = dlg.GetSelections()
     4644        dlg.Destroy()
     4645        if not result: return
     4646        useHist = [histList[item] for item in result]
     4647        useFo = od['value_1']
     4648        nDif = len(useHist)
     4649        nCol = 4*nDif+3
     4650        refDict = {}
     4651        Hrange = [1000,-1000]
     4652        Krange = [1000,-1000]
     4653        Lrange = [1000,-1000]
     4654        dRange = [10000.,0.0]
     4655        IoRange = [[10000.0,0.0] for i in range(nDif)]
     4656        IcRange = [[10000.0,0.0] for i in range(nDif)]
     4657        WIRange = [[10000.0,0.0] for i in range(nDif)]
     4658       
     4659        for ih,hist in enumerate(useHist):
     4660            Refset = Histograms[hist]['Reflection Lists'][General['Name']]['RefList']
     4661            for ref in Refset:
     4662                Hrange[0] = min(ref[0],Hrange[0])
     4663                Hrange[1] = max(ref[0],Hrange[1])
     4664                Krange[0] = min(ref[1],Krange[0])
     4665                Krange[1] = max(ref[1],Krange[1])
     4666                Lrange[0] = min(ref[2],Lrange[0])
     4667                Lrange[1] = max(ref[2],Lrange[1])
     4668                dRange[0] = min(ref[4],dRange[0])
     4669                dRange[1] = max(ref[4],dRange[1])
     4670                if useFo:
     4671                    IoRange[ih][0] = min(np.sqrt(ref[8]),IoRange[ih][0])
     4672                    IoRange[ih][1] = max(np.sqrt(ref[8]),IoRange[ih][1])
     4673                    WIRange[ih][0] = min(1./np.sqrt(ref[8]),WIRange[ih][0])
     4674                    WIRange[ih][1] = max(1./np.sqrt(ref[8]),WIRange[ih][1])
     4675                    IcRange[ih][0] = min(np.sqrt(ref[9]),IcRange[ih][0])
     4676                    IcRange[ih][1] = max(np.sqrt(ref[9]),IcRange[ih][1])
     4677                else:
     4678                    IoRange[ih][0] = min(ref[8],IoRange[ih][0])
     4679                    IoRange[ih][1] = max(ref[8],IoRange[ih][1])
     4680                    WIRange[ih][0] = min(1./ref[8],WIRange[ih][0])
     4681                    WIRange[ih][1] = max(1./ref[8],WIRange[ih][1])
     4682                    IcRange[ih][0] = min(ref[9],IcRange[ih][0])
     4683                    IcRange[ih][1] = max(ref[9],IcRange[ih][1])
     4684                hklStr = '%5d%5d%5d'%(ref[0],ref[1],ref[2])
     4685                if useFo:
     4686                    rec = {ih:np.array([np.sqrt(ref[8]),1./np.sqrt(ref[8]),np.sqrt(ref[9]),ref[10]],dtype=np.float32)}
     4687                else:
     4688                    rec = {ih:np.array([ref[8],1./ref[8],ref[9],ref[10]],dtype=np.float32)}
     4689                if hklStr in refDict:
     4690                    refDict[hklStr].update(rec)
     4691                else:
     4692                    refDict[hklStr] = rec
     4693            nRef = len(refDict)
     4694        #set up header stuff
     4695        Header = '%s'%('VERS MTZ:V1.1'.ljust(80))
     4696        Header += ('TITLE %s'%(General['Name'])).ljust(80)
     4697        Header += ('NCOL %10d%10d%10d'%(nCol,nRef,0)).ljust(80)
     4698        Header += ('CELL %9.4f%9.4f%9.4f%9.4f%9.4f%9.4f'%(Cell[0],Cell[1],Cell[2],Cell[3],Cell[4],Cell[5])).ljust(80)
     4699        Header += 'SORT    1   2   3   0   0'.ljust(80)
     4700        Nops = len(SGData['SGOps'])
     4701        Ncent = len(SGData['SGCen'])
     4702        SGSym = SGData['SpGrp']
     4703        SGnum = G2spc.SpaceGroupNumber(SGSym)
     4704        Header += ("SYMINF %6d%6d %s %6d '%s' PG%s"%(Nops*Ncent,Nops,SGSym[0],SGnum,SGSym,SGData['SGPtGrp'])).ljust(80)
     4705        SGtext,SGTable = G2spc.SGPrint(SGData)
     4706        textOps = G2spc.TextOps(SGtext,SGTable,reverse=True)
     4707        for ops in textOps:
     4708            Header += ('SYMM  %s'%ops.upper()).ljust(80)
     4709        Header += ('RESO  %.4f  %.4f '%(1./dRange[1]*2,1./dRange[0]*2)).ljust(80)
     4710        Header += 'VALM NAN'.ljust(80)
     4711        Header += ('COLUMN %s H%18d%18d'%('H'.ljust(30),Hrange[0],Hrange[1])).ljust(79)+'0'
     4712        Header += ('COLUMN %s H%18d%18d'%('K'.ljust(30),Krange[0],Krange[1])).ljust(79)+'0'
     4713        Header += ('COLUMN %s H%18d%18d'%('L'.ljust(30),Lrange[0],Lrange[1])).ljust(79)+'0'
     4714        for ih,hist in enumerate(useHist):
     4715            hname = hist.split()[1].split('.')[0]
     4716            if useFo:
     4717                Header += ('COLUMN %s F%18.4f%18.4f'%(('F_'+hname).ljust(30),IoRange[ih][0],IoRange[ih][1])).ljust(78)+'%2d'%(ih+1)
     4718                Header += ('COLUMN %s Q%18.8f%18.8f'%(('WT_'+hname).ljust(30),WIRange[ih][0],WIRange[ih][1])).ljust(78)+'%2d'%(ih+1)
     4719                Header += ('COLUMN %s F%18.4f%18.4f'%(('FC_'+hname).ljust(30),IcRange[ih][0],IcRange[ih][1])).ljust(78)+'%2d'%(ih+1)
     4720            else:
     4721                Header += ('COLUMN %s J%18.4f%18.4f'%(('I_'+hname).ljust(30),IoRange[ih][0],IoRange[ih][1])).ljust(78)+'%2d'%(ih+1)
     4722                Header += ('COLUMN %s Q%18.8f%18.8f'%(('WT_'+hname).ljust(30),WIRange[ih][0],WIRange[ih][1])).ljust(78)+'%2d'%(ih+1)
     4723                Header += ('COLUMN %s J%18.4f%18.4f'%(('IC_'+hname).ljust(30),IcRange[ih][0],IcRange[ih][1])).ljust(78)+'%2d'%(ih+1)
     4724            Header += ('COLUMN %s P%18.4f%18.4f'%(('PHIC_'+hname).ljust(30),-180.0,180.)).ljust(78)+'%2d'%(ih+1)
     4725        Header += ('NDIF %10d'%nDif).ljust(80)
     4726        projName = os.path.split(self.GSASprojectfile)[1]
     4727        for iH,hist in enumerate(useHist):
     4728            hname = hist.split()[1].split('.')[0]
     4729            Iparm = Histograms[hist]['Instrument Parameters'][0]
     4730            wave = G2mth.getWave(Iparm)
     4731            cell = G2lat.CellDijCorr(Cell,SGData,Phase['Histograms'],hist)
     4732            Header += ('PROJECT %7d %s'%(iH+1,projName)).ljust(80)
     4733            Header += ('CRYSTAL %7d %s'%(iH+1,General['Name'])).ljust(80)
     4734            Header += ('DATASET %7d %s'%(iH+1,hname)).ljust(80)
     4735            Header += ('DCELL   %7d %10.4f%10.4f%10.4f%10.4f%10.4f%10.4f'   \
     4736                %(iH+1,cell[0],cell[1],cell[2],cell[3],cell[4],cell[5])).ljust(80) # not correct- sholud include Dijs
     4737            Header += ('DWAVE   %7d %10.5f'%(iH+1,wave)).ljust(80)
     4738        Header += ('END'.ljust(80))
     4739        Header += ('Created by GSAS-II on %s'%time.ctime()).ljust(80)
     4740        Header += ('MTZENDOFHEADER'.ljust(80))
     4741       
     4742        if useFo:
     4743            fName = General['Name'].replace(' ','_')+'F.mtz'
     4744        else:
     4745            fName = General['Name'].replace(' ','_')+'F2.mtz'
     4746        MTZ = open(fName,'wb')
     4747        MTZ.write(st.pack('4sl2s70x',b'MTZ ',21+nCol*nRef,b'DA'))
     4748       
     4749        for ref in refDict:
     4750            rec = [float(i) for i in ref.split()]
     4751            refData = refDict[ref]
     4752            for i in range(nDif):
     4753                if i in refData:
     4754                    rec += list(refData[i])
     4755                else:
     4756                    rec += [np.nan,np.nan,np.nan,np.nan]
     4757            string = b''
     4758            for item in rec:
     4759                string += st.pack('f',item)
     4760            MTZ.write(string)
     4761        MTZ.write(Header.encode())
     4762        MTZ.close()
     4763        print('MTZ file %s written'%fName)
     4764       
    46154765    def OnExportPeakList(self,event):
    46164766        pth = G2G.GetExportPath(self)
     
    68246974        tmpSizer.Add(G2G.HelpButton(G2frame.dataWindow,helpIndex='RefineType'))
    68256975        LSSizer.Add(tmpSizer,0,WACV)
    6826         Choice=['analytic Jacobian','numeric','analytic Hessian','Hessian SVD']   #TODO +'SVD refine' - what flags will it need?
     6976        Choice=['analytic Jacobian','numeric','analytic Hessian','Hessian SVD']
    68276977        derivSel = wx.ComboBox(parent=G2frame.dataWindow,value=data['deriv type'],choices=Choice,
    68286978            style=wx.CB_READONLY|wx.CB_DROPDOWN)
     
    68496999            LSSizer.Add(wx.StaticText(G2frame.dataWindow,label='SVD zero\ntolerance:'),0,WACV)
    68507000            LSSizer.Add(G2G.ValidatedTxtCtrl(G2frame.dataWindow,data,'SVDtol',nDig=(10,1,'g'),xmin=1.e-9,xmax=.01),0,WACV)
    6851         else:       #TODO what for SVD refine?
     7001        else:
    68527002            LSSizer.Add(wx.StaticText(G2frame.dataWindow,label=' Initial shift factor: '),0,WACV)
    68537003            Factr = G2G.ValidatedTxtCtrl(G2frame.dataWindow,data,'shift factor',nDig=(10,5),xmin=1.e-5,xmax=100.)
  • trunk/GSASIIlattice.py

    r5082 r5090  
    254254             parmDict[pfx+'A0']+parmDict[phfx+'D11'],0,0,0]
    255255    return A
     256
     257def CellDijCorr(Cell,SGData,Data,hist):
     258    '''Returns the cell corrected for Dij values.
     259
     260    :param list Cell: lattice parameters
     261    :param dict SGdata: a symmetry object
     262    :param dict Data: phase data structure; contains set of Dij values
     263    :param str hist: histogram name
     264
     265    :returns: cell corrected for Dij values
     266    '''
     267    A = cell2A(Cell)
     268    Dij = Data[hist]['HStrain'][0]
     269    if SGData['SGLaue'] in ['-1',]:
     270        newA = [A[0]+Dij[0],A[1]+Dij[1],A[2]+Dij[2],A[3]+Dij[3],A[4]+Dij[4],A[5]+Dij[5]]
     271    elif SGData['SGLaue'] in ['2/m',]:
     272        if SGData['SGUniq'] == 'a':
     273            newA = [A[0]+Dij[0],A[1]+Dij[1],A[2]+Dij[2],0,0,A[5]+Dij[3]]
     274        elif SGData['SGUniq'] == 'b':
     275            newA = [A[0]+Dij[0],A[1]+Dij[1],A[2]+Dij[2],0,A[4]+Dij[3]]
     276        else:
     277            newA = [A[0]+Dij[0],A[1]+Dij[1],A[2]+Dij[2],A[3]+Dij[3],0,0]
     278    elif SGData['SGLaue'] in ['mmm',]:
     279        newA = [A[0]+Dij[0],A[1]+Dij[1],A[2]+Dij[2],0,0,0]
     280    elif SGData['SGLaue'] in ['4/m','4/mmm']:
     281        newA = [A[0]+Dij[0],A[0]+Dij[0],A[2]+Dij[1],0,0,0]
     282    elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
     283        newA = [A[0]+Dij[0],A[0]+Dij[0],A[2]+Dij[1],A[0]+Dij[0],0,0]
     284    elif SGData['SGLaue'] in ['3R', '3mR']:
     285        newA = [A[0]+Dij[0],A[0]+Dij[0],A[0]+Dij[0],A[3]+Dij[1],A[3]+Dij[1],A[3]+Dij[1]]
     286    elif SGData['SGLaue'] in ['m3m','m3']:
     287        newA = [A[0]+Dij[0],A[0]+Dij[0],A[0]+Dij[0],0,0,0]
     288    return A2cell(newA)
    256289   
    257290def prodMGMT(G,Mat):
  • trunk/GSASIIphsGUI.py

    r5075 r5090  
    62716271            rname = pName+'-fullrmc.py'
    62726272            if not os.path.exists(rname):
    6273                 G2G.G2MessageBox(G2frame,
    6274                     'The fullrmc script has not been created. Running setup.',
     6273                G2G.G2MessageBox(G2frame,'The fullrmc script has not been created. Running setup.',
    62756274                    'Not setup')
    62766275                OnSetupRMC(event)
     
    64066405            rname = pName+'-PDFfit.py'
    64076406            # if not os.path.exists(rname):
    6408             #     G2G.G2MessageBox(G2frame,
    6409             #         'The PDFfit script has not been created. Running setup.',
     6407            #     G2G.G2MessageBox(G2frame,'The PDFfit script has not been created. Running setup.',
    64106408            #         'Not setup')
    64116409            #     OnSetupRMC(event)
  • trunk/GSASIIplot.py

    r5077 r5090  
    22212221            newPlot = True
    22222222        elif event.key in ['+','=']:
    2223             G2frame.plusPlot = not G2frame.plusPlot
     2223            G2frame.plusPlot = (G2frame.plusPlot+1)%3
    22242224        elif event.key == '/':
    22252225            Page.plotStyle['Normalize'] = not Page.plotStyle['Normalize']
     
    30693069        G2frame.Contour = False
    30703070        G2frame.Weight = True
    3071         G2frame.plusPlot = True
     3071        G2frame.plusPlot = 1
    30723072        G2frame.SubBack = False
    30733073        Page.plotStyle['logPlot'] = False
     
    31653165                Page.Choice += ['e: create excluded region',
    31663166                        's: toggle sqrt plot','w: toggle (Io-Ic)/sig plot',
    3167                         '+: no selection']
     3167                        '+: toggle obs line plot']
    31683168            else:
    31693169                Page.Choice += ['q: toggle q plot','s: toggle sqrt plot',
    31703170                    't: toggle d-spacing plot','w: toggle (Io-Ic)/sig plot',
    3171                     '+: no selection']
     3171                    '+: toggle obs line plot']
    31723172            if Page.plotStyle['sqrtPlot'] or Page.plotStyle['logPlot']:
    31733173                del Page.Choice[1]
     
    31853185                'b: toggle subtract background file','g: toggle grid',
    31863186                'm: toggle multidata plot','n: toggle semilog/loglog',
    3187                 'q: toggle S(q) plot','w: toggle (Io-Ic)/sig plot','+: toggle selection',]
     3187                'q: toggle S(q) plot','w: toggle (Io-Ic)/sig plot','+: toggle obs line plot',]
    31883188            if not G2frame.SinglePlot:
    31893189                Page.Choice = Page.Choice+ \
     
    34863486                    Plot.set_ylabel('Data sequence',fontsize=14)
    34873487        else:
    3488             if G2frame.plusPlot:
     3488            if not G2frame.plusPlot:
     3489                pP = ''
     3490                lW = 1.5
     3491            elif G2frame.plusPlot == 1:
    34893492                pP = '+'
    34903493                lW = 0
    34913494            else:
    3492                 pP = ''
     3495                pP = '+'
    34933496                lW = 1.5
    34943497            if plottype in ['SASD','REFD'] and Page.plotStyle['logPlot']:
  • trunk/exports/G2export_PDB.py

    r4573 r5090  
    248248            self.CloseFile()
    249249            print('Phase '+phasenam+' written to XYZ file '+self.fullpath)
    250    
     250               
Note: See TracChangeset for help on using the changeset viewer.