Changeset 4523


Ignore:
Timestamp:
Jul 17, 2020 12:51:18 PM (3 years ago)
Author:
toby
Message:

more fullrmc work, but not done

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIddataGUI.py

    r4491 r4523  
    605605               
    606606        dispSizer = wx.BoxSizer(wx.HORIZONTAL)
    607         dispRef = wx.CheckBox(DData,wx.ID_ANY,label=' Layer displacement (\xb5m): ')
     607        dispRef = wx.CheckBox(DData,wx.ID_ANY,label=u' Layer displacement (\xb5m): ')
    608608        dispRef.SetValue(UseList[G2frame.hist]['Layer Disp'][1])
    609609        dispRef.Bind(wx.EVT_CHECKBOX, OnDispRef)

  • trunk/GSASIIphsGUI.py

    r4518 r4523  
    47434743            return pairSizer
    47444744                   
    4745         def FileSizer(RMCdict):
     4745        def FileSizer(RMCdict,mainSizer):
    47464746           
    47474747            def OnFitScale(event):
    … …  
    48034803                           
    48044804            Indx = {}
    4805             titleSizer = wx.BoxSizer(wx.HORIZONTAL)
    4806             titleSizer.Add(wx.StaticText(G2frame.FRMC,label=' Select data for processing: '),0,WACV)
    4807             fitscale = wx.CheckBox(G2frame.FRMC,label=' Fit scale factors?')
    4808             fitscale.SetValue(RMCPdict['FitScale'])
    4809             fitscale.Bind(wx.EVT_CHECKBOX,OnFitScale)
    4810             titleSizer.Add(fitscale,0,WACV)
    4811             mainSizer.Add(titleSizer,0,WACV)
     4805            mainSizer.Add(wx.StaticText(G2frame.FRMC,label='Select data for processing: '),0,WACV)
    48124806            if G2frame.RMCchoice == 'fullrmc':
    4813                 mainSizer.Add(wx.StaticText(G2frame.FRMC,
    4814                     label=' NB: fullrmc data files must be 2 columns; all other lines preceeded by "#". Edit before use.'),0,WACV)
    48154807                Heads = ['Name','File','Weight','type','Plot','Delete']
    48164808                fileSizer = wx.FlexGridSizer(6,5,5)
    … …  
    48384830                            if type(RMCPdict['files'][fil][3]) is bool: RMCPdict['files'][fil][3] = 0
    48394831                            fmtTyp = G2G.G2ChoiceButton(G2frame.FRMC,choices,RMCPdict['files'][fil],3)
    4840                         elif 'F(Q)' in fil:
     4832                        elif '(Q)' in fil:
    48414833                            choices = 'F(Q)-RMCProfile','S(Q)-PDFFIT'
    48424834                            if type(RMCPdict['files'][fil][3]) is bool: RMCPdict['files'][fil][3] = 0
    … …  
    48534845                        Indx[delBtn.GetId()] = fil
    48544846                        fileSizer.Add(delBtn,0,WACV)
    4855                         if 'F(Q)' in fil:
     4847                        if '(Q)' in fil:
    48564848                            fileSizer.Add((-1,-1),0)
    48574849                            corrChk = wx.CheckBox(G2frame.FRMC,label='Apply sinc convolution? ')
    … …  
    48704862                        fileSizer.Add((-1,-1),0)
    48714863                        fileSizer.Add((-1,-1),0)
    4872                 return fileSizer
     4864                mainSizer.Add(fileSizer,0,WACV)
     4865                fitscale = wx.CheckBox(G2frame.FRMC,label=' Fit scale factors?')
     4866                fitscale.SetValue(RMCPdict['FitScale'])
     4867                fitscale.Bind(wx.EVT_CHECKBOX,OnFitScale)
     4868                mainSizer.Add(fitscale,0,WACV)
     4869                mainSizer.Add(wx.StaticText(G2frame.FRMC,
     4870                    label=' NB: fullrmc data files must be 2 columns; all other lines preceeded by "#". Edit before use.'),0,WACV)
     4871                return
    48734872            # RMCProfile
    48744873            Heads = ['Name','File','Format','Weight','Plot','Delete']
    … …  
    49074906                    fileSizer.Add((5,5),0)
    49084907                    fileSizer.Add((5,5),0)
    4909             return fileSizer
     4908            mainSizer.Add(fileSizer,0,WACV)
     4909            fitscale = wx.CheckBox(G2frame.FRMC,label=' Fit scale factors?')
     4910            fitscale.SetValue(RMCPdict['FitScale'])
     4911            fitscale.Bind(wx.EVT_CHECKBOX,OnFitScale)
     4912            mainSizer.Add(fitscale,0,WACV)
     4913            return
    49104914       
    49114915        G2frame.GetStatusBar().SetStatusText('',1)
    … …  
    49594963                Pairs = {pairs:[0.0,0.0,0.0] for pairs in Pairs}
    49604964                files = {'Neutron real space data; G(r): ':['Select',1.,'G(r)',0,True],
    4961                           'Neutron reciprocal space data; F(Q): ':['Select',1.,'F(Q)',0,True],
     4965                          'Neutron reciprocal space data; S(Q)-1: ':['Select',1.,'F(Q)',0,True],
    49624966                          'Xray real space data; G(r): ':['Select',1.,'G(r)',0,True],
    4963                           'Xray reciprocal space data; F(Q): ':['Select',1.,'F(Q)',0,True],}
     4967                          'Xray reciprocal space data; S(Q)-1: ':['Select',1.,'F(Q)',0,True],}
    49644968                data['RMC']['fullrmc'] = {'SuperCell':[1,1,1],'Box':[10.,10.,10.],'aTypes':aTypes,
    49654969                    'atSeq':atSeq,'Pairs':Pairs,'files':files,'ReStart':[False,False],'Cycles':1,
    … …  
    49704974
    49714975            def GetSuperSizer():
     4976                def ShowRmax(*args,**kwargs):
     4977                    cell = data['General']['Cell'][1:7]
     4978                    bigcell = np.array(cell)*np.array(RMCPdict['SuperCell']+[1,1,1])
     4979                    bigG = G2lat.cell2Gmat(bigcell)[0]
     4980                    rmax = min([0.5/np.sqrt(G2lat.calc_rDsq2(H,bigG)) for H in np.eye(3)])
     4981                    rmaxlbl.SetLabel('  Rmax = {:.1f}'.format(rmax))
    49724982                superSizer = wx.BoxSizer(wx.HORIZONTAL)
    49734983                axes = ['X','Y','Z']
    … …  
    49754985                    superSizer.Add(wx.StaticText(G2frame.FRMC,label=' %s-axis: '%ax),0,WACV)
    49764986                    superSizer.Add(G2G.ValidatedTxtCtrl(G2frame.FRMC,RMCPdict['SuperCell'],
    4977                         i,xmin=1,xmax=20,size=(50,25)),0,WACV)
     4987                        i,xmin=1,xmax=20,size=(50,25),OnLeave=ShowRmax),0,WACV)
     4988                rmaxlbl = wx.StaticText(G2frame.FRMC,label=' Rmax=?')
     4989                superSizer.Add(rmaxlbl,0,WACV)
     4990                ShowRmax()
    49784991                return superSizer
    49794992
    … …  
    52155228   
    52165229            G2G.HorizontalLine(mainSizer,G2frame.FRMC)
    5217             mainSizer.Add(FileSizer(RMCPdict),0,WACV)
     5230            FileSizer(RMCPdict,mainSizer)
    52185231               
    52195232        elif G2frame.RMCchoice ==  'RMCProfile':
    … …  
    56105623            mainSizer.Add(samSizer,0,WACV)
    56115624
    5612             mainSizer.Add(FileSizer(RMCPdict),0,WACV)
     5625            FileSizer(RMCPdict,mainSizer)
    56135626   
    56145627        SetPhaseWindow(G2frame.FRMC,mainSizer)
    … …  
    57975810        if G2frame.RMCchoice == 'fullrmc':
    57985811               
    5799             #dlg = wx.DirDialog (None, "Choose fullrmc directory", "",
    5800             #        wx.DD_DEFAULT_STYLE | wx.DD_DIR_MUST_EXIST)
    58015812            RMCPdict = data['RMC']['fullrmc']
    58025813            generalData = data['General']
    … …  
    58055816            engineFilePath = pName+'.rmc'
    58065817            if not os.path.exists(engineFilePath):
    5807                 dlg = wx.FileDialog(self, 'Open fullrmc directory',
     5818                #dlg = wx.DirDialog (None, "Choose fullrmc directory", "",
     5819                #        wx.DD_DEFAULT_STYLE | wx.DD_DIR_MUST_EXIST)
     5820                dlg = wx.FileDialog(G2frame, 'Open fullrmc directory',
    58085821                                        defaultFile='*.rmc',
    58095822                                        wildcard='*.rmc')
    … …  
    58175830                engineFilePath = os.path.splitext(engineFilePath)[0] + '.rmc'
    58185831                if not os.path.exists(engineFilePath): return
    5819             # import psutil
    5820             # pid = RMCPdict.get('pid',-1)
    5821             # Proc = None
    5822             # if pid and psutil.pid_exists(pid):
    5823             #     proc = psutil.Process(pid).children()
    5824             #     for child in proc:
    5825             #         if 'conhost' in child.name():       #probably very Windows specific
    5826             #             Proc = child
    58275832            from fullrmc import Engine
    5828             # load
     5833            # try to load a few times w/short delay before giving up
     5834            for i in range(10):
     5835                try:
     5836                    ENGINE = Engine().load(engineFilePath)
     5837                    break
     5838                except:
     5839                    time.sleep(0.1)
     5840            else:
     5841                G2G.G2MessageBox(G2frame,
     5842                    'The fullrmc project exists but could not be loaded.',
     5843                    'Not loaded')
     5844                return
     5845
     5846            # make a list of possible plots
     5847            plotList = []
     5848            found = False
     5849            try:
     5850                frame = ENGINE.usedFrame
     5851                for item in ENGINE.constraints:
     5852                    sitem = str(type(item))
     5853                    if 'PairDistribution' in sitem or 'StructureFactor' in sitem or 'PairCorrelation' in sitem:
     5854                        if 'neutron' in item.weighting:
     5855                            nameId = 'Neutron'
     5856                        else:
     5857                            nameId = 'Xray'
     5858                        if 'StructureFactor' in sitem:
     5859                            title = 'S(Q)-{} for {}'.format(nameId,pName)
     5860                        else:
     5861                            title = 'g(r)-{} for {}'.format(nameId,pName)
     5862                           
     5863                        plotList.append(title+pName)
     5864                        plotList.append('All partials of '+title+pName)
     5865                        plotList.append('Total partials of '+title+pName)
     5866                        found = True
     5867                    elif 'BondConstraint' in sitem:
     5868                        try:
     5869                            bonds = item.get_constraint_value()['bondsLength']
     5870                        except Exception as msg:
     5871                            print("Error reading BondConstraint bondsLength\n  ",msg)
     5872                            continue
     5873                        atoms = ENGINE.allElements 
     5874                        bondList = item.bondsList[:2]
     5875                        bondNames = ['%s-%s'%(atoms[bondList[0][iat]],atoms[bondList[1][iat]]) for iat in range(bonds.shape[0])]
     5876                        for Bname in set(bondNames):
     5877                            title = '{} Bond lengths for {}'.format(Bname,pName)
     5878                            plotList.append(title)
     5879                        found = True
     5880                    elif 'BondsAngleConstraint' in sitem:
     5881                        angles = 180.*item.get_constraint_value()['angles']/np.pi
     5882                        angleList = item.anglesList[:3]
     5883                        atoms = ENGINE.get_original_data("allElements",frame)
     5884                        angleNames = ['%s-%s-%s'%(atoms[angleList[1][iat]],atoms[angleList[0][iat]],atoms[angleList[2][iat]]) for iat in range(angles.shape[0])]
     5885                        for Aname in set(angleNames):
     5886                            title = '{} Bond angles for {}'.format(Aname,pName)
     5887                            plotList.append(title)
     5888                        found = True
     5889                    elif 'DihedralAngleConstraint' in sitem:
     5890                        impangles = 180.*item.get_constraint_value()['angles']/np.pi
     5891                        impangleList = item.anglesList[:4]
     5892                        atoms = ENGINE.get_original_data("allElements",frame)
     5893                        impangleNames = ['%s-%s-%s-%s'%(atoms[impangleList[0][iat]],atoms[impangleList[1][iat]],
     5894                                    atoms[impangleList[2][iat]],atoms[impangleList[3][iat]]) for iat in range(impangles.shape[0])]
     5895                        for Aname in set(impangleNames):
     5896                            title = '{} Dihedral angles for {}'.format(Aname,pName)
     5897                            plotList.append(title)
     5898                        found = True
     5899                    elif hasattr(item,'plot'):
     5900                        plotList.append(item.__class__.__name__+' pyplot')
     5901                        found = True
     5902            except Exception as msg:
     5903                 print("Unexpected error reading from fullrmc engine\n  ",msg)
     5904                 return
     5905            if not found or len(plotList) == 0:
     5906                G2G.G2MessageBox(G2frame,'No saved information yet, wait until fullrmc does a Save',
     5907                                         'No info')
     5908                return
     5909             
     5910            dlg = G2G.G2MultiChoiceDialog(G2frame,'Select plots to see displayed. N.B. pyplots are in a separate window.',
     5911                                              'Select plots',plotList)
     5912            try:
     5913                result = dlg.ShowModal()
     5914                if result == wx.ID_OK:
     5915                    selectedPlots = [plotList[i] for i in dlg.GetSelections()]
     5916                else:
     5917                    return
     5918            finally:
     5919                dlg.Destroy()
     5920
     5921            eNames = ['Total',]
     5922            nObs = [0,]
     5923            frame = ENGINE.usedFrame
     5924            for item in ENGINE.constraints:
     5925                sitem = str(type(item))
     5926                if 'PairDistribution' in sitem or 'StructureFactor' in sitem or 'PairCorrelation' in sitem:
     5927                    if 'neutron' in item.weighting:
     5928                        nameId = 'Neutron'
     5929                    else:
     5930                        nameId = 'Xray'
     5931                    if 'StructureFactor' in sitem:
     5932                        eNames.append('S(Q)-'+nameId)
     5933                        Xlab = r'$\mathsf{Q,\AA^{-1}}$'
     5934                        Ylab = 'S(Q)'
     5935                        title = 'S(Q)-{} for {}'.format(nameId,pName)
     5936                    else:
     5937                        eNames.append('g(r)-'+nameId)
     5938                        Xlab = r'$\mathsf{r,\AA}$'
     5939                        Ylab = r'$\mathsf{g(r),\AA^{-2}}$'
     5940                        title = 'g(r)-{} for {}'.format(nameId,pName)
     5941                    dataDict= item.get_constraints_properties(frame)
     5942                    X = dataDict['frames-experimental_x'][0]
     5943                    Y = dataDict['frames-experimental_y'][0]
     5944                    nObs.append(X.shape[0])
     5945                    rdfDict = item.get_constraint_value()
     5946                    if 'total' not in rdfDict:
     5947                        print('No data yet - wait for a save')
     5948                        return
     5949                    Z = rdfDict['total']
     5950                    XY = [[X,Z],[X,Y]]
     5951                    Names = ['Calc','Obs']
     5952                    if title in selectedPlots:
     5953                        G2plt.PlotXY(G2frame,XY,labelX=Xlab,
     5954                            labelY=Ylab,newPlot=True,Title=title,lines=True,names=Names)
     5955                    PXY = []
     5956                    PXYT = []
     5957                    Names = []
     5958                    NamesT = []
     5959
     5960                    for item in rdfDict:
     5961                        if 'va' in item:
     5962                            continue
     5963                        if 'rdf' not in item and 'g(r)' in title:
     5964                            Ylab = r'$\mathsf{G(r),\AA^{-1}}$'
     5965                            PXYT.append([X,1.+rdfDict[item]/X])
     5966                        else:
     5967                            PXYT.append([X,rdfDict[item]])
     5968                        NamesT.append(item)
     5969                        if 'total' in item:
     5970                            if 'rdf' not in item and 'g(r)' in title:
     5971                                Ylab = r'$\mathsf{G(r),\AA^{-1}}$'
     5972                                PXY.append([X,1.+rdfDict[item]/X])
     5973                            else:
     5974                                PXY.append([X,rdfDict[item]])
     5975                            Names.append(item)
     5976                    if 'All partials of '+title in selectedPlots:
     5977                        G2plt.PlotXY(G2frame,PXYT,labelX=Xlab,labelY=Ylab,
     5978                                         newPlot=True,Title='All partials of '+title,
     5979                                         lines=True,names=NamesT)
     5980                    if 'Total partials of '+title in selectedPlots:
     5981                        G2plt.PlotXY(G2frame,PXY,labelX=Xlab,labelY=Ylab,
     5982                                         newPlot=True,Title='Total partials of '+title,
     5983                                         lines=True,names=Names)
     5984
     5985                elif 'BondConstraint' in sitem:
     5986                    try:
     5987                        bonds = item.get_constraint_value()['bondsLength']
     5988                    except Exception as msg:
     5989                        print("Error reading BondConstraint bondsLength\n  ",msg)
     5990                        continue
     5991                    atoms = ENGINE.allElements 
     5992                    bondList = item.bondsList[:2]
     5993                    bondNames = ['%s-%s'%(atoms[bondList[0][iat]],atoms[bondList[1][iat]]) for iat in range(bonds.shape[0])]
     5994                    Bonds = list(zip(bondNames,bonds))
     5995                    for Bname in set(bondNames):
     5996                        bondLens = [bond[1] for bond in Bonds if bond[0]==Bname]
     5997                        title = '{} Bond lengths for {}'.format(Bname,pName)
     5998                        if title in selectedPlots:
     5999                            G2plt.PlotBarGraph(G2frame,bondLens,Xname=r'%s $Bond,\ \AA$'%Bname,Title=title,
     6000                                                        PlotName='%s Bonds for %s'%(Bname,pName),maxBins=20)
     6001                        #print(' %d %s bonds found'%(len(bondLens),Bname))
     6002                elif 'BondsAngleConstraint' in sitem:
     6003                    # code not tested
     6004                    angles = 180.*item.get_constraint_value()['angles']/np.pi
     6005                    angleList = item.anglesList[:3]
     6006                    atoms = ENGINE.get_original_data("allElements",frame)
     6007                    angleNames = ['%s-%s-%s'%(atoms[angleList[1][iat]],atoms[angleList[0][iat]],atoms[angleList[2][iat]]) for iat in range(angles.shape[0])]
     6008                    Angles = list(zip(angleNames,angles))
     6009                    for Aname in set(angleNames):
     6010                        bondAngs = [angle[1] for angle in Angles if angle[0]==Aname]
     6011                        title = '{} Bond angles for {}'.format(Aname,pName)
     6012                        if title in selectedPlots:
     6013                            G2plt.PlotBarGraph(G2frame,bondAngs,Xname=r'%s Angle, deg'%Aname,Title=title,
     6014                                    PlotName='%s Angles for %s'%(Aname,pName),maxBins=20)
     6015                        print(' %d %s angles found'%(len(bondAngs),Aname))
     6016                elif 'DihedralAngleConstraint' in sitem:
     6017                    # code not tested
     6018                    impangles = 180.*item.get_constraint_value()['angles']/np.pi
     6019                    impangleList = item.anglesList[:4]
     6020                    print(' Dihedral angle chi^2 =  %2f'%item.standardError)
     6021                    atoms = ENGINE.get_original_data("allElements",frame)
     6022                    impangleNames = ['%s-%s-%s-%s'%(atoms[impangleList[0][iat]],atoms[impangleList[1][iat]],
     6023                                atoms[impangleList[2][iat]],atoms[impangleList[3][iat]]) for iat in range(impangles.shape[0])]
     6024                    impAngles = list(zip(impangleNames,impangles))
     6025                    for Aname in set(impangleNames):
     6026                        impAngs = [angle[1] for angle in impAngles if angle[0]==Aname]
     6027                        title = '{} Dihedral angles for {}'.format(Aname,pName)
     6028                        if title in selectedPlots:
     6029                            G2plt.PlotBarGraph(G2frame,impAngs,Xname=r'%s $Dihedral Angle, deg$'%Aname,Title=title,
     6030                                    PlotName='%s Dihedral Angles for %s'%(Aname,pName),maxBins=20)
     6031                #elif 'AtomicCoordinationNumber' in sitem or 'InterMolecular' in sitem:
     6032                #    print(sitem+' not plotted')
     6033                elif hasattr(item,'plot'):
     6034                    #print('skipping constraint ',sitem)
     6035                    if item.__class__.__name__+' pyplot' in selectedPlots:
     6036                        item.plot(show=True)
    58296037            GSASIIpath.IPyBreak()
    5830             ENGINE = Engine(path=None)
    5831             ENGINE.load(engineFilePath)
    5832            
    5833             # try:
    5834             #     if engineFilePath is not None:
    5835             #         result, mes = ENGINE.is_engine(engineFilePath, mes=True)
    5836             #         if result:
    5837             #             if Proc is not None:
    5838             #                 Proc.suspend()
    5839             #             ENGINE = ENGINE.load(engineFilePath)
    5840             eNames = ['Total',]
    5841             found = False
    5842             nObs = [0,]
    5843             try:
    5844                 for frame in list(ENGINE.frames):
    5845                     ENGINE.set_used_frame(frame)
    5846                     for item in ENGINE.constraints:
    5847                         sitem = str(type(item))
    5848                         if 'PairDistribution' in sitem or 'StructureFactor' in sitem or 'PairCorrelation' in sitem:
    5849                             nameId = 'X'
    5850                             if 'neutron' in item.weighting:
    5851                                 nameId = 'N'
    5852                             found = True
    5853                             Xlab = r'$\mathsf{r,\AA}$'
    5854                             Ylab = r'$\mathsf{g(r),\AA^{-2}}$'
    5855                             title = ' g(r)%s for '%nameId
    5856                             if 'StructureFactor' in sitem:
    5857                                 eNames.append('S(Q)'+nameId)
    5858                                 Xlab = r'$\mathsf{Q,\AA^{-1}}$'
    5859                                 Ylab = 'S(Q)'
    5860                                 title = ' S(Q)%s for '%nameId
    5861                             else:
    5862                                 eNames.append('g(r)'+nameId)
    5863                             dataDict= item.get_constraints_properties(frame)
    5864                             X = dataDict['frames-experimental_x'][0]
    5865                             Y = dataDict['frames-experimental_y'][0]
    5866                             nObs.append(X.shape[0])
    5867                             rdfDict = item.get_constraint_value()
    5868                             if 'total' not in rdfDict:
    5869                                 print('No data yet - wait for a save')
    5870                                 #ENGINE.close()
    5871                                 #if Proc is not None:
    5872                                 #    Proc.resume()
    5873                                 return
    5874                             Z = rdfDict['total']
    5875                             XY = [[X,Z],[X,Y]]
    5876                             Names = ['Calc','Obs']
    5877                             G2plt.PlotXY(G2frame,XY,labelX=Xlab,
    5878                                 labelY=Ylab,newPlot=True,Title=title+pName,
    5879                                 lines=True,names=Names)
    5880                             PXY = []
    5881                             PXYT = []
    5882                             Names = []
    5883                             NamesT = []
    5884                             for item in rdfDict:
    5885                                 if 'va' in item:
    5886                                     continue
    5887                                 if 'rdf' not in item and 'g(r)' in title:
    5888                                     Ylab = r'$\mathsf{G(r),\AA^{-1}}$'
    5889                                     PXYT.append([X,1.+rdfDict[item]/X])
    5890                                 else:
    5891                                     PXYT.append([X,rdfDict[item]])
    5892                                 NamesT.append(item)
    5893                                 if 'total' in item:
    5894                                     if 'rdf' not in item and 'g(r)' in title:
    5895                                         Ylab = r'$\mathsf{G(r),\AA^{-1}}$'
    5896                                         PXY.append([X,1.+rdfDict[item]/X])
    5897                                     else:
    5898                                         PXY.append([X,rdfDict[item]])
    5899                                     Names.append(item)
    5900                             G2plt.PlotXY(G2frame,PXYT,labelX=Xlab,
    5901                                 labelY=Ylab,newPlot=True,Title=' All partials of '+title+pName,
    5902                                 lines=True,names=NamesT)
    5903                             G2plt.PlotXY(G2frame,PXY,labelX=Xlab,
    5904                                 labelY=Ylab,newPlot=True,Title=' Total partials of '+title+pName,
    5905                                 lines=True,names=Names)
    5906 
    5907                         else:
    5908                             found = True
    5909                             if 'BondConstraint' in sitem:
    5910                                 try:
    5911                                     bonds = item.get_constraint_value()['bondsLength']
    5912                                 except TypeError:
    5913                                     break
    5914                                 bondList = item.bondsList[:2]
    5915                                 atoms = ENGINE.get_original_data("allElements",frame)
    5916                                 bondNames = ['%s-%s'%(atoms[bondList[0][iat]],atoms[bondList[1][iat]]) for iat in range(bonds.shape[0])]
    5917                                 bondSet = list(set(bondNames))
    5918                                 Bonds = list(zip(bondNames,bonds))
    5919                                 for Bname in bondSet:
    5920                                     bondLens = [bond[1] for bond in Bonds if bond[0]==Bname]
    5921                                     G2plt.PlotBarGraph(G2frame,bondLens,Xname=r'%s $Bond,\ \AA$'%Bname,Title='%s Bond lengths for %s'%(Bname,pName),
    5922                                         PlotName='%s Bonds for %s'%(Bname,pName),maxBins=20)
    5923                                     print(' %d %s bonds found'%(len(bondLens),Bname))
    5924                             elif 'BondsAngleConstraint' in sitem:
    5925                                 angles = 180.*item.get_constraint_value()['angles']/np.pi
    5926                                 angleList = item.anglesList[:3]
    5927                                 atoms = ENGINE.get_original_data("allElements",frame)
    5928                                 angleNames = ['%s-%s-%s'%(atoms[angleList[1][iat]],atoms[angleList[0][iat]],atoms[angleList[2][iat]]) for iat in range(angles.shape[0])]
    5929                                 angleSet = list(set(angleNames))
    5930                                 Angles = list(zip(angleNames,angles))
    5931                                 for Aname in angleSet:
    5932                                     bondAngs = [angle[1] for angle in Angles if angle[0]==Aname]
    5933                                     G2plt.PlotBarGraph(G2frame,bondAngs,Xname=r'%s Angle, deg'%Aname,Title='%s Bond angles for %s'%(Aname,pName),
    5934                                         PlotName='%s Angles for %s'%(Aname,pName),maxBins=20)
    5935                                     print(' %d %s angles found'%(len(bondAngs),Aname))
    5936                             elif 'DihedralAngleConstraint' in sitem:
    5937                                 impangles = 180.*item.get_constraint_value()['angles']/np.pi
    5938                                 impangleList = item.anglesList[:4]
    5939                                 print(' Dihedral angle chi^2 =  %2f'%item.standardError)
    5940                                 atoms = ENGINE.get_original_data("allElements",frame)
    5941                                 impangleNames = ['%s-%s-%s-%s'%(atoms[impangleList[0][iat]],atoms[impangleList[1][iat]],
    5942                                     atoms[impangleList[2][iat]],atoms[impangleList[3][iat]]) for iat in range(impangles.shape[0])]
    5943                                 impangleSet = list(set(impangleNames))
    5944                                 impAngles = list(zip(impangleNames,impangles))
    5945                                 for Aname in impangleSet:
    5946                                     impAngs = [angle[1] for angle in impAngles if angle[0]==Aname]
    5947                                     G2plt.PlotBarGraph(G2frame,impAngs,Xname=r'%s $Dihedral Angle, deg$'%Aname,Title='%s Dihedral angles for %s'%(Aname,pName),
    5948                                         PlotName='%s Dihedral Angles for %s'%(Aname,pName),maxBins=20)
    5949                             elif 'AtomicCoordinationNumber' in sitem or 'InterMolecular' in sitem:
    5950                                 print(sitem+' not plotted')
    5951                             else:
    5952                                 print(sitem)
    5953                                 item.plot(show=True)
    5954                                 pass
    5955                     nObs[0] = np.sum(nObs[1:])
    5956                 if not found:
    5957                     print(' No saved information yet, wait until fullrmc does a Save')
    5958                     eNames = []
    5959                 #ENGINE.close()      #return lock on ENGINE repository & close it
    5960                 #if Proc is not None:
    5961                 #    Proc.resume()
    5962             except AssertionError:
    5963                  print("Can't open fullrmc engine while running")
     6038           
    59646039            # loglines = []
    59656040            # ilog = 0

  • trunk/GSASIIpwd.py

    r4520 r4523  
    27792779    rmin = RMCPdict['min Contact']
    27802780    cell = Phase['General']['Cell'][1:7]
    2781     bigcell = np.array(cell)*np.array(RMCPdict['SuperCell']+[1,1,1])
    2782     bigG = G2lat.cell2Gmat(bigcell)[0]
    2783     rmax = min([0.5/np.sqrt(G2lat.calc_rDsq2(H,bigG)) for H in np.eye(3)])
    27842781    SymOpList = G2spc.AllOps(Phase['General']['SGData'])[0]
    27852782    cx,ct,cs,cia = Phase['General']['AtomPtrs']
    … …  
    27912788    restart = '%s_restart.pdb'%pName
    27922789    Files = RMCPdict['files']
    2793     wtDict = {}
    27942790    rundata = ''
    27952791    rundata += '#### fullrmc %s file; edit by hand if you so choose #####\n'%rname
    … …  
    28552851            if filDat[3] == 0:
    28562852                rundata += '''    # read and xform G(r) as defined in RMCProfile
    2857 # see eq. 44 in Keen, J. Appl. Cryst. (2001) 34 172-177\n'''
     2853    # see eq. 44 in Keen, J. Appl. Cryst. (2001) 34 172-177\n'''
    28582854                rundata += '    GR[1] *= 4 * np.pi * GR[0] * rho0 / sumCiBi2\n'
    28592855                rundata += '    GofR = fPDF.PairDistributionConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
    … …  
    28682864            rundata += '    ENGINE.add_constraints([GofR])\n'
    28692865            rundata += '    GofR.set_limits((None, rmax))\n'
    2870             wtDict['Pair-'+sfwt] = filDat[1]
    2871         elif 'F(Q)' in File:
     2866        elif '(Q)' in File:
    28722867            rundata += '    SOQ = np.loadtxt("%s").T\n'%filDat[0]
    28732868            if filDat[3] == 0:
    2874                 rundata += '    # Read & xform F(Q) as defined in RMCProfile\n'
     2869                rundata += '    # Read & xform F(Q) as defined in RMCProfile to S(Q)-1\n'
    28752870                rundata += '    SOQ[1] *= 1 / sumCiBi2\n'
    2876                 rundata += '    SOQ[1] += 1\n'
    28772871            elif filDat[3] == 1:
    28782872                rundata += '    # This is S(Q) as defined in PDFFIT\n'
     2873                rundata += '    SOQ[1] -= 1\n'
    28792874            if filDat[4]:
    2880                 rundata += '    SOQ[1] = Collection.sinc_convolution(q=SOQ[0],sq=SOQ[1],rmax={:.3f})\n'.format(rmax)
     2875                rundata += '    SOQ[1] = Collection.sinc_convolution(q=SOQ[0],sq=SOQ[1],rmax=rmax)\n'
    28812876            rundata += '    SofQ = ReducedStructureFactorConstraint(experimentalData=SOQ.T, weighting="%s")\n'%sfwt
    28822877            rundata += '    ENGINE.add_constraints([SofQ])\n'
    28832878        else:
    28842879            print('What is this?')
    2885         wtDict['Struct-'+sfwt] = filDat[1]
    28862880    rundata += '    ENGINE.add_constraints(DistanceConstraint(defaultLowerDistance={}))\n'.format(RMCPdict['min Contact'])
    28872881    rundata += '''    B_CONSTRAINT   = BondConstraint()
    … …  
    29062900#            rundata += 'SwapB = [[idx] for idx in range(len(aN)) if aN[idx]=="%s"]\n'%swap[1]
    29072901#            rundata += 'SwapGen["%s-%s"] = [SwapPositionsGenerator(swapList=SwapA),SwapPositionsGenerator(swapList=SwapB),%.2f]\n'%(swap[0],swap[1],swap[2])
    2908     rundata += '# setup/change constraints - can be done without restart\n'   
     2902    rundata += '\n# set weights -- do this now so values can be changed without a restart\n'
     2903    rundata += 'wtDict = {}\n'
     2904    for File in Files:
     2905        filDat = RMCPdict['files'][File]
     2906        if not os.path.exists(filDat[0]): continue
     2907        if 'Xray' in File:
     2908            sfwt = 'atomicNumber'
     2909        else:
     2910            sfwt = 'neutronCohb'
     2911        if 'G(r)' in File:
     2912            typ = 'Pair'
     2913        elif '(Q)' in File:
     2914            typ = 'Struct'
     2915        rundata += 'wtDict["{}-{}"] = {}\n'.format(typ,sfwt,filDat[1])
    29092916    rundata += 'for c in ENGINE.constraints:  # loop over predefined constraints\n'
    2910     rundata += '    strcons = str(type(c))\n'
    29112917    rundata += '    if type(c) is fPDF.PairDistributionConstraint:\n'
    2912     rundata += '        c.set_variance_squared(wtDict["Pair-"+c.weighting])\n'
    2913     rundata += '        c.set_limits((None,%.3f))\n'%(rmax)
     2918    rundata += '        c.set_variance_squared(1./wtDict["Pair-"+c.weighting])\n'
     2919    rundata += '        c.set_limits((None,rmax))\n'
    29142920    if RMCPdict['FitScale']:
    29152921        rundata += '        c.set_adjust_scale_factor((10, 0.01, 100.))\n'
    2916     rundata += '    elif "StructureFactor" in strcons:\n'
    2917     rundata += '        c.set_variance_squared(wtDict["Struct-"+c.weighting])\n'
     2922    rundata += '    elif type(c) is ReducedStructureFactorConstraint:\n'
     2923    rundata += '        c.set_variance_squared(1./wtDict["Struct-"+c.weighting])\n'
    29182924    if RMCPdict['FitScale']:
    29192925        rundata += '        c.set_adjust_scale_factor((10, 0.01, 100.))\n'
    2920 #    if AngleList and not RMCPdict['Swaps']: rundata += setAngleConstraints()
    29212926    # if len(RMCPdict['Torsions']):         # Torsions currently commented out in GUI
    29222927    #     rundata += 'for c in ENGINE.constraints:  # look for Dihedral Angle Constraints\n'
    … …  
    29542959    rfile.writelines(rundata)
    29552960    rfile.close()
    2956    
    29572961    return rname
    2958 
    2959 # def MakefullrmcPDB(Name,Phase,RMCPdict):
    2960 #     generalData = Phase['General']
    2961 #     Atseq = RMCPdict['atSeq']
    2962 #     Dups,Fracs = findDup(Phase['Atoms'])
    2963 #     Sfracs = [np.cumsum(fracs) for fracs in Fracs]
    2964 #     ifSfracs = any([np.any(sfracs-1.) for sfracs in Sfracs])
    2965 #     Supercell = RMCPdict['SuperCell']
    2966 #     Cell = generalData['Cell'][1:7]
    2967 #     Trans = np.eye(3)*np.array(Supercell)
    2968 #     newPhase = copy.deepcopy(Phase)
    2969 #     newPhase['General']['SGData'] = G2spc.SpcGroup('P 1')[1]
    2970 #     newPhase['General']['Cell'][1:] = G2lat.TransformCell(Cell,Trans.T)
    2971 #     newPhase,Atcodes = G2lat.TransformPhase(Phase,newPhase,Trans,np.zeros(3),np.zeros(3),ifMag=False,Force=True)
    2972 #     Atoms = newPhase['Atoms']
    2973 
    2974 #     if ifSfracs:
    2975 #         Natm = np.core.defchararray.count(np.array(Atcodes),'+')    #no. atoms in original unit cell
    2976 #         Natm = np.count_nonzero(Natm-1)
    2977 #         Satoms = []
    2978 #         for i in range(len(Atoms)//Natm):
    2979 #             ind = i*Natm
    2980 #             Satoms.append(G2mth.sortArray(G2mth.sortArray(G2mth.sortArray(Atoms[ind:ind+Natm],5),4),3))
    2981 #         Natoms = []
    2982 #         for satoms in Satoms:
    2983 #             for idup,dup in enumerate(Dups):
    2984 #                 ldup = len(dup)
    2985 #                 natm = len(satoms)
    2986 #                 i = 0
    2987 #                 while i < natm:
    2988 #                     if satoms[i][0] in dup:
    2989 #                         atoms = satoms[i:i+ldup]
    2990 #                         try:
    2991 #                             atom = atoms[np.searchsorted(Sfracs[idup],rand.random())]
    2992 #                             Natoms.append(atom)
    2993 #                         except IndexError:      #what about vacancies?
    2994 #                             if 'Va' not in Atseq:
    2995 #                                 Atseq.append('Va')
    2996 #                                 RMCPdict['aTypes']['Va'] = 0.0
    2997 #                             atom = atoms[0]
    2998 #                             atom[1] = 'Va'
    2999 #                             Natoms.append(atom)
    3000 #                         i += ldup
    3001 #                     else:
    3002 #                        i += 1
    3003 #     else:
    3004 #         Natoms = Atoms
    3005 
    3006 #     XYZ = np.array([atom[3:6] for atom in Natoms]).T
    3007 #     XYZptp = np.array([ma.ptp(XYZ[0]),ma.ptp(XYZ[1]),ma.ptp(XYZ[2])])/2.
    3008 #     Cell = newPhase['General']['Cell'][1:7]
    3009 #     A,B = G2lat. cell2AB(Cell)
    3010 #     fname = Name+'_cbb.pdb'
    3011 #     fl = open(fname,'w')
    3012 #     fl.write('REMARK    Boundary Conditions:%6.2f  0.0  0.0  0.0%7.2f  0.0  0.0  0.0%7.2f\n'%(
    3013 #              Cell[0],Cell[1],Cell[2]))
    3014 #     fl.write('ORIGX1      1.000000  0.000000  0.000000        0.00000\n')
    3015 #     fl.write('ORIGX2      0.000000  1.000000  0.000000        0.00000\n')
    3016 #     fl.write('ORIGX3      0.000000  0.000000  1.000000        0.00000\n')
    3017 #     fl.write('CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f P 1           1\n'%(
    3018 #             Cell[0],Cell[1],Cell[2],Cell[3],Cell[4],Cell[5]))
    3019 
    3020 #     Natm = np.core.defchararray.count(np.array(Atcodes),'+')
    3021 #     Natm = np.count_nonzero(Natm-1)
    3022 #     nat = 0
    3023 #     if RMCPdict['byMolec']:
    3024 #         NPM = RMCPdict['Natoms']
    3025 #         for iat,atom in enumerate(Natoms):
    3026 #             XYZ = np.inner(A,np.array(atom[3:6])-XYZptp)    #shift origin to middle & make Cartesian;residue = 'RMC'
    3027 #             fl.write('ATOM  %5d %-4s RMC%6d%12.3f%8.3f%8.3f  1.00  0.00          %2s\n'%(       
    3028 #                     1+nat%NPM,atom[0],1+nat//NPM,XYZ[0],XYZ[1],XYZ[2],atom[1].lower()))
    3029 #             nat += 1
    3030 #     else:
    3031 #         for atm in Atseq:
    3032 #             for iat,atom in enumerate(Natoms):
    3033 #                 if atom[1] == atm:
    3034 #                     XYZ = np.inner(A,np.array(atom[3:6])-XYZptp)    #shift origin to middle & make Cartesian
    3035 #                     fl.write('ATOM  %5d %-4s RMC%6d%12.3f%8.3f%8.3f  1.00  0.00          %2s\n'%(       
    3036 #                             1+nat,atom[0],1+nat,XYZ[0],XYZ[1],XYZ[2],atom[1].lower()))
    3037 #                     nat += 1
    3038 #     fl.close()
    3039 #     return fname
    3040    
    3041 def MakePdparse(RMCPdict):
    3042     fname = 'make_pdb.py'
    3043     outName = RMCPdict['moleculePdb'].split('.')
    3044     outName[0] += '_rbb'
    3045     outName = '.'.join(outName)
    3046     RMCPdict['atomPDB'] = outName   #might be empty if pdbparser run fails
    3047     fl = open(fname,'w')
    3048     fl.write('from pdbparser.pdbparser import pdbparser\n')
    3049     fl.write('from pdbparser.Utilities.Construct import AmorphousSystem\n')
    3050     fl.write("pdb = pdbparser('%s')\n"%RMCPdict['moleculePdb'])
    3051     boxstr= 'boxsize=%s'%str(RMCPdict['Box'])
    3052     recstr = 'recursionLimit=%d'%RMCPdict['maxRecursion']
    3053     denstr = 'density=%.3f'%RMCPdict['targetDensity']
    3054     fl.write('pdb = AmorphousSystem(pdb,%s,%s,%s,\n'%(boxstr,recstr,denstr))
    3055     fl.write('    priorities={"boxSize":True, "insertionNumber":False, "density":True}).construct().get_pdb()\n')
    3056     fl.write('pdb.export_pdb("%s")\n'%outName)
    3057     fl.close
    3058     return fname
    3059 
     2962   
    30602963def GetRMCBonds(general,RMCPdict,Atoms,bondList):
    30612964    bondDist = []
Note: See TracChangeset for help on using the changeset viewer.