Changeset 4951 for trunk/GSASIIphsGUI.py


Ignore:
Timestamp:
Jun 11, 2021 8:42:20 AM (4 months ago)
Author:
toby
Message:

working fullrmc implementation. Much more to do (see TODO in phsGUI), but good start. Needs external fullrmc Python from Bachir until 5.0 is released.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIphsGUI.py

    r4916 r4951  
    46464646
    46474647#### RMC Data page ################################################################################
     4648# fullrmc stuff TODO:
     4649#  1) need to implement swapping in scripts
     4650#  2) need fullrmc installation instructions
     4651#  3) when GSASIIpwd.findfullrmc fails, should trigger message with link to above
     4652#  4) fullrmc tutorials
     4653#  5) better plotting when fullrmc in main Python image?
     4654
    46484655    def UpdateRMC():
    46494656        ''' Present the controls for running fullrmc or RMCProfile
     
    57125719            G2frame.dataWindow.FRMCDataEdit.Enable(G2G.wxID_RUNRMC,True)
    57135720            RMCPdict = data['RMC']['fullrmc']
     5721            if RMCPdict['Swaps']:
     5722                wx.MessageDialog(G2frame, G2G.StripIndents(
     5723                        '''GSAS-II does not yet fully support use of swapping in fullrmc.
     5724                        Edit the script by hand before using.''',True),
     5725                        'No swaps yet',wx.OK).ShowModal()
    57145726            # debug stuff
    5715             # if GSASIIpath.GetConfigValue('debug'):
    5716             #     print('reloading',G2pwd)
    5717             #     import imp
    5718             #     imp.reload(G2pwd)
     5727            #if GSASIIpath.GetConfigValue('debug'):
     5728            #    print('reloading',G2pwd)
     5729            #    import imp
     5730            #    imp.reload(G2pwd)
    57195731            # end debug stuff           
    57205732            rname = G2pwd.MakefullrmcRun(pName,data,RMCPdict)
     
    57695781        generalData = data['General']
    57705782        if G2frame.RMCchoice == 'fullrmc':
     5783            fullrmc_exec = G2pwd.findfullrmc()
     5784            if fullrmc_exec is None:
     5785                G2G.G2MessageBox(G2frame,'fullrmc Python not found. How did we get here?')
     5786                return
    57715787            pName = G2frame.GSASprojectfile.split('.')[0] + '-' + generalData['Name']
    57725788            pName = pName.replace(' ','_')
     
    57785794                OnSetupRMC(event)
    57795795            RMCPdict = data['RMC']['fullrmc']
    5780             #rmcname = pName+'.rmc'
    5781             # if os.path.isdir(rmcname) and RMCPdict['ReStart'][0]:
    5782             #     G2G.G2MessageBox(G2frame,
    5783             #         '''You have asked to restart fullrmc but have an existing
    5784             #         run as {}. You must manually delete this directory if
    5785             #         you wish to restart or change the restart checkbox to
    5786             #         continue from the previous results.
    5787             #         '''.format(rmcname),'Restart or Continue?')
    5788             #     # TODO: could do this for the user with:
    5789             #     #import shutil
    5790             #     #shutil.rmtree(rmcname)
     5796            rmcname = pName+'-fullrmc.rmc'
     5797            if os.path.isdir(rmcname) and RMCPdict['ReStart'][0]:
     5798                msg = '''You have asked to start a new fullrmc run rather than
     5799                     continue the existing {} run.
     5800                     %%Press "Yes" to continue, deleting this
     5801                     previous run or "No" to change the restart checkbox to
     5802                     continue from the previous results.'''.format(rmcname)
     5803
     5804                dlg = wx.MessageDialog(G2frame,G2G.StripIndents(msg,True),
     5805                                           'Restart or continue',
     5806                    wx.YES|wx.NO)
     5807                try:
     5808                    dlg.CenterOnParent()
     5809                    result = dlg.ShowModal()
     5810                finally:
     5811                    dlg.Destroy()
     5812                if result == wx.ID_YES:
     5813                    import shutil
     5814                    shutil.rmtree(rmcname)
     5815                else:
     5816                    return
    57915817            G2G.G2MessageBox(G2frame,
    57925818'''For use of fullrmc, please cite:
     
    58025828                logname = '%s_%d.log'%(pName,ilog)
    58035829                if os.path.isfile(logname):
     5830                    if GSASIIpath.GetConfigValue('debug'):
     5831                        print('removing',logname)
    58045832                    os.remove(logname)
    58055833                else:
     
    58095837            if sys.platform.lower().startswith('win'):
    58105838                batch = open('fullrmc.bat','w')
    5811                 batch.write('CALL '+sys.exec_prefix+'\\Scripts\\activate\n')
    5812                 batch.write(sys.exec_prefix+'\\python.exe '+rname+'\n')
     5839                #batch.write('CALL '+sys.exec_prefix+'\\Scripts\\activate\n')
     5840                batch.write(fullrmc_exec+' '+rname+'\n')
    58135841                batch.write('pause')
    58145842                batch.close()
     
    58175845                batch = open('fullrmc.sh','w')
    58185846                batch.write('#!/bin/bash\n')
    5819                 activate = os.path.split(os.environ.get('CONDA_EXE',''))[0] +'/activate'
     5847                #activate = os.path.split(os.environ.get('CONDA_EXE',''))[0] +'/activate'
    58205848                batch.write('cd ' + os.path.split(os.path.abspath(rname))[0] + '\n')
    5821                 if os.path.exists(activate):
    5822                     batch.write('source ' + activate + ' ' +
    5823                                 os.environ['CONDA_DEFAULT_ENV'] +'\n')
    5824                     batch.write('python ' + rname + '\n')
    5825                 else:
    5826                     batch.write(sys.exec_prefix+'/python ' + rname + '\n')
     5849                #if os.path.exists(activate):
     5850                #    batch.write('source ' + activate + ' ' +
     5851                #                os.environ['CONDA_DEFAULT_ENV'] +'\n')
     5852                #    batch.write('python ' + rname + '\n')
     5853                #else:
     5854                #    batch.write(sys.exec_prefix+'/python ' + rname + '\n')
     5855                batch.write(fullrmc_exec + ' ' + rname + '\n')
    58275856                batch.close()
    58285857                if sys.platform == "darwin":
     
    59075936    def OnViewRMC(event):
    59085937        if G2frame.RMCchoice == 'fullrmc':
    5909                
    59105938            RMCPdict = data['RMC']['fullrmc']
    59115939            generalData = data['General']
    59125940            pName = G2frame.GSASprojectfile.split('.')[0] + '-' + generalData['Name']
    59135941            pName = pName.replace(' ','_')
    5914             engineFilePath = pName+'.rmc'
     5942            engineFilePath = pName+'-fullrmc.rmc'
    59155943            if not os.path.exists(engineFilePath):
    5916                 #dlg = wx.DirDialog (None, "Choose fullrmc directory", "",
    5917                 #        wx.DD_DEFAULT_STYLE | wx.DD_DIR_MUST_EXIST)
    59185944                dlg = wx.FileDialog(G2frame, 'Open fullrmc directory',
    59195945                                        defaultFile='*.rmc',
     
    59285954                engineFilePath = os.path.splitext(engineFilePath)[0] + '.rmc'
    59295955                if not os.path.exists(engineFilePath): return
    5930             from fullrmc import Engine
    5931             # try to load a few times w/short delay before giving up
    5932             for i in range(10):
    5933                 try:
    5934                     ENGINE = Engine().load(engineFilePath)
    5935                     break
    5936                 except:
    5937                     time.sleep(0.1)
    5938             else:
     5956            choices = []
     5957            statFilePath = os.path.splitext(engineFilePath)[0] + '.stats'
     5958            plotFilePath = os.path.splitext(engineFilePath)[0] + '.plots'
     5959            imgDict = {}
     5960            if os.path.exists(statFilePath):
     5961                fp = open(statFilePath,'r')
     5962                vals = []
     5963                for i,line in enumerate(fp):
     5964                    v = line.strip().split(',')[:-1] # ends with comma, remove last empty element
     5965                    if i == 0:
     5966                        lbls = [i.strip() for i in v[1:]]
     5967                        continue
     5968                    try:
     5969                        vals.append([float(i) for i in v])
     5970                    except:
     5971                        print('Error reading line ',i,'in',statFilePath)
     5972                fp.close()
     5973                steps = np.array(vals)[:,0]
     5974                yvals = np.array(vals)[:,1:].T
     5975                choices = ['Constraints vs. Steps']
     5976            if os.path.exists(plotFilePath):
     5977                import pickle
     5978                fp = open(plotFilePath,'rb')
     5979                while True:
     5980                    try:
     5981                        title = pickle.load(fp)
     5982                        imgDict[title] = fp.tell()
     5983                        im = pickle.load(fp)
     5984                        choices += [title]
     5985                    except:
     5986                        break
     5987                fp.close()
     5988            if not choices:
    59395989                G2G.G2MessageBox(G2frame,
    5940                     'The fullrmc project exists but could not be loaded.',
    5941                     'Not loaded')
     5990                    'Nothing to plot. '+                 
     5991                    'No results in '+statFilePath+' or '+plotFilePath,
     5992                    'Nothing to plot')
    59425993                return
    5943 
    5944             # make a list of possible plots
    5945             plotList = []
    5946             found = False
    5947             try:
    5948                 frame = ENGINE.usedFrame
    5949                 for item in ENGINE.constraints:
    5950                     sitem = str(type(item))
    5951                     if 'PairDistribution' in sitem or 'StructureFactor' in sitem or 'PairCorrelation' in sitem:
    5952                         if 'neutron' in item.weighting:
    5953                             nameId = 'Neutron'
    5954                         else:
    5955                             nameId = 'Xray'
    5956                         if 'StructureFactor' in sitem:
    5957                             title = 'S(Q)-{} for {}'.format(nameId,pName)
    5958                         else:
    5959                             title = 'g(r)-{} for {}'.format(nameId,pName)
    5960                            
    5961                         plotList.append(title)
    5962                         plotList.append('All partials of '+title)
    5963                         plotList.append('Total partials of '+title)
    5964                         found = True
    5965                     elif 'BondConstraint' in sitem:
    5966                         try:
    5967                             bonds = item.get_constraint_value()['bondsLength']
    5968                         except Exception as msg:
    5969                             print("Error reading BondConstraint bondsLength\n  ",msg)
    5970                             continue
    5971                         atoms = ENGINE.allElements 
    5972                         bondList = item.bondsList[:2]
    5973                         bondNames = ['%s-%s'%(atoms[bondList[0][iat]],atoms[bondList[1][iat]]) for iat in range(bonds.shape[0])]
    5974                         for Bname in set(bondNames):
    5975                             title = '{} Bond lengths for {}'.format(Bname,pName)
    5976                             plotList.append(title)
    5977                         found = True
    5978                     elif 'BondsAngleConstraint' in sitem:
    5979                         angles = 180.*item.get_constraint_value()['angles']/np.pi
    5980                         angleList = item.anglesList[:3]
    5981                         atoms = ENGINE.get_original_data("allElements",frame)
    5982                         angleNames = ['%s-%s-%s'%(atoms[angleList[1][iat]],atoms[angleList[0][iat]],atoms[angleList[2][iat]]) for iat in range(angles.shape[0])]
    5983                         for Aname in set(angleNames):
    5984                             title = '{} Bond angles for {}'.format(Aname,pName)
    5985                             plotList.append(title)
    5986                         found = True
    5987                     elif 'DihedralAngleConstraint' in sitem:
    5988                         impangles = 180.*item.get_constraint_value()['angles']/np.pi
    5989                         impangleList = item.anglesList[:4]
    5990                         atoms = ENGINE.get_original_data("allElements",frame)
    5991                         impangleNames = ['%s-%s-%s-%s'%(atoms[impangleList[0][iat]],atoms[impangleList[1][iat]],
    5992                                     atoms[impangleList[2][iat]],atoms[impangleList[3][iat]]) for iat in range(impangles.shape[0])]
    5993                         for Aname in set(impangleNames):
    5994                             title = '{} Dihedral angles for {}'.format(Aname,pName)
    5995                             plotList.append(title)
    5996                         found = True
    5997                     elif hasattr(item,'plot'):
    5998                         plotList.append(item.__class__.__name__+' pyplot')
    5999                         found = True
    6000             except Exception as msg:
    6001                  print("Unexpected error reading from fullrmc engine\n  ",msg)
    6002                  return
    6003             plotList.append('fullrmc residuals for '+pName)
    6004             if not found or len(plotList) == 0:
    6005                 G2G.G2MessageBox(G2frame,'No saved information yet, wait until fullrmc does a Save',
    6006                                          'No info')
    6007                 return
    6008              
    6009             dlg = G2G.G2MultiChoiceDialog(G2frame,'Select plots to see displayed. N.B. pyplots are in a separate window.',
    6010                                               'Select plots',plotList)
     5994            dlg = G2G.G2MultiChoiceDialog(G2frame,'Select plots to see displayed.',
     5995                                              'Select plots',choices)
    60115996            try:
    60125997                result = dlg.ShowModal()
    60135998                if result == wx.ID_OK:
    6014                     selectedPlots = [plotList[i] for i in dlg.GetSelections()]
     5999                    selectedPlots = [choices[i] for i in dlg.GetSelections()]
    60156000                else:
    60166001                    return
     
    60186003                dlg.Destroy()
    60196004
    6020             eNames = ['Total',]
    6021             nObs = [0,]
    6022             frame = ENGINE.usedFrame
    6023             for item in ENGINE.constraints:
    6024                 sitem = str(type(item))
    6025                 if 'PairDistribution' in sitem or 'StructureFactor' in sitem or 'PairCorrelation' in sitem:
    6026                     if 'neutron' in item.weighting:
    6027                         nameId = 'Neutron'
    6028                     else:
    6029                         nameId = 'Xray'
    6030                     if 'StructureFactor' in sitem:
    6031                         eNames.append('S(Q)-'+nameId)
    6032                         Xlab = r'$\mathsf{Q,\AA^{-1}}$'
    6033                         Ylab = 'S(Q)'
    6034                         title = 'S(Q)-{} for {}'.format(nameId,pName)
    6035                     else:
    6036                         eNames.append('g(r)-'+nameId)
    6037                         Xlab = r'$\mathsf{r,\AA}$'
    6038                         Ylab = r'$\mathsf{g(r),\AA^{-2}}$'
    6039                         title = 'g(r)-{} for {}'.format(nameId,pName)
    6040                     dataDict= item.get_constraints_properties(frame)
    6041                     X = dataDict['frames-experimental_x'][0]
    6042                     Y = dataDict['frames-experimental_y'][0]
    6043                     nObs.append(X.shape[0])
    6044                     rdfDict = item.get_constraint_value()
    6045                     if 'total' not in rdfDict:
    6046                         print('No data yet - wait for a save')
    6047                         return
    6048                     Z = rdfDict['total']
    6049                     XY = [[X,Z],[X,Y]]
    6050                     Names = ['Calc','Obs']
    6051                     if title in selectedPlots:
    6052                         G2plt.PlotXY(G2frame,XY,labelX=Xlab,
    6053                             labelY=Ylab,newPlot=True,Title=title,lines=True,names=Names)
    6054                     PXY = []
    6055                     PXYT = []
    6056                     Names = []
    6057                     NamesT = []
    6058 
    6059                     for item in rdfDict:
    6060                         if 'va' in item:
    6061                             continue
    6062                         if 'rdf' not in item and 'g(r)' in title:
    6063                             Ylab = r'$\mathsf{G(r),\AA^{-1}}$'
    6064                             PXYT.append([X,1.+rdfDict[item]/X])
    6065                         else:
    6066                             PXYT.append([X,rdfDict[item]])
    6067                         NamesT.append(item)
    6068                         if 'total' in item:
    6069                             if 'rdf' not in item and 'g(r)' in title:
    6070                                 Ylab = r'$\mathsf{G(r),\AA^{-1}}$'
    6071                                 PXY.append([X,1.+rdfDict[item]/X])
    6072                             else:
    6073                                 PXY.append([X,rdfDict[item]])
    6074                             Names.append(item)
    6075                     if 'All partials of '+title in selectedPlots:
    6076                         G2plt.PlotXY(G2frame,PXYT,labelX=Xlab,labelY=Ylab,
    6077                                          newPlot=True,Title='All partials of '+title,
    6078                                          lines=True,names=NamesT)
    6079                     if 'Total partials of '+title in selectedPlots:
    6080                         G2plt.PlotXY(G2frame,PXY,labelX=Xlab,labelY=Ylab,
    6081                                          newPlot=True,Title='Total partials of '+title,
    6082                                          lines=True,names=Names)
    6083 
    6084                 elif 'BondConstraint' in sitem:
    6085                     try:
    6086                         bonds = item.get_constraint_value()['bondsLength']
    6087                     except Exception as msg:
    6088                         print("Error reading BondConstraint bondsLength\n  ",msg)
    6089                         continue
    6090                     atoms = ENGINE.allElements 
    6091                     bondList = item.bondsList[:2]
    6092                     bondNames = ['%s-%s'%(atoms[bondList[0][iat]],atoms[bondList[1][iat]]) for iat in range(bonds.shape[0])]
    6093                     Bonds = list(zip(bondNames,bonds))
    6094                     for Bname in set(bondNames):
    6095                         bondLens = [bond[1] for bond in Bonds if bond[0]==Bname]
    6096                         title = '{} Bond lengths for {}'.format(Bname,pName)
    6097                         if title in selectedPlots:
    6098                             G2plt.PlotBarGraph(G2frame,bondLens,Xname=r'%s $Bond,\ \AA$'%Bname,Title=title,
    6099                                                         PlotName='%s Bonds for %s'%(Bname,pName),maxBins=20)
    6100                         #print(' %d %s bonds found'%(len(bondLens),Bname))
    6101                 elif 'BondsAngleConstraint' in sitem:
    6102                     # code not tested
    6103                     angles = 180.*item.get_constraint_value()['angles']/np.pi
    6104                     angleList = item.anglesList[:3]
    6105                     atoms = ENGINE.get_original_data("allElements",frame)
    6106                     angleNames = ['%s-%s-%s'%(atoms[angleList[1][iat]],atoms[angleList[0][iat]],atoms[angleList[2][iat]]) for iat in range(angles.shape[0])]
    6107                     Angles = list(zip(angleNames,angles))
    6108                     for Aname in set(angleNames):
    6109                         bondAngs = [angle[1] for angle in Angles if angle[0]==Aname]
    6110                         title = '{} Bond angles for {}'.format(Aname,pName)
    6111                         if title in selectedPlots:
    6112                             G2plt.PlotBarGraph(G2frame,bondAngs,Xname=r'%s Angle, deg'%Aname,Title=title,
    6113                                     PlotName='%s Angles for %s'%(Aname,pName),maxBins=20)
    6114                         print(' %d %s angles found'%(len(bondAngs),Aname))
    6115                 elif 'DihedralAngleConstraint' in sitem:
    6116                     # code not tested
    6117                     impangles = 180.*item.get_constraint_value()['angles']/np.pi
    6118                     impangleList = item.anglesList[:4]
    6119                     print(' Dihedral angle chi^2 =  %2f'%item.standardError)
    6120                     atoms = ENGINE.get_original_data("allElements",frame)
    6121                     impangleNames = ['%s-%s-%s-%s'%(atoms[impangleList[0][iat]],atoms[impangleList[1][iat]],
    6122                                 atoms[impangleList[2][iat]],atoms[impangleList[3][iat]]) for iat in range(impangles.shape[0])]
    6123                     impAngles = list(zip(impangleNames,impangles))
    6124                     for Aname in set(impangleNames):
    6125                         impAngs = [angle[1] for angle in impAngles if angle[0]==Aname]
    6126                         title = '{} Dihedral angles for {}'.format(Aname,pName)
    6127                         if title in selectedPlots:
    6128                             G2plt.PlotBarGraph(G2frame,impAngs,Xname=r'%s $Dihedral Angle, deg$'%Aname,Title=title,
    6129                                     PlotName='%s Dihedral Angles for %s'%(Aname,pName),maxBins=20)
    6130                 #elif 'AtomicCoordinationNumber' in sitem or 'InterMolecular' in sitem:
    6131                 #    print(sitem+' not plotted')
    6132                 elif hasattr(item,'plot'):
    6133                     #print('skipping constraint ',sitem)
    6134                     if item.__class__.__name__+' pyplot' in selectedPlots:
    6135                         item.plot(show=True)
    6136                        
    6137             # read log files to get std err vs steps
    6138             ilog = 0
    6139             Gen = []
    6140             Err = []
    6141             while True:
    6142                 fname = '{}_{}.log'.format(pName,ilog)
    6143                 try:
    6144                     logfile = open(fname,'r')
    6145                     for line in logfile.readlines():
    6146                         if "fullrmc <STEP>" not in line: continue
    6147                         Gen.append(float(int(line.split('Gen:')[1].split()[0])))
    6148                         Err.append(float(line.split('Err:')[1].strip()))
    6149                     logfile.close()
    6150                 except FileNotFoundError:
    6151                     break
    6152                 ilog += 1
    6153             title = 'fullrmc residuals for '+pName
    6154             #GSASIIpath.IPyBreak()
    6155             if len(Gen) > 0 and title in selectedPlots:
    6156                 Gen = np.array(Gen)
    6157                 Err = np.log10(Err)
    6158                 G2plt.PlotXY(G2frame,[[Gen,Err]],labelX='generated steps',
    6159                                labelY=r'$log_{10}$ ($\mathsf{\chi^2})$',newPlot=True,Title=title,
    6160                                lines=True,names=eNames)
     6005            for plt in selectedPlots:
     6006                if plt in imgDict:
     6007                    fp = open(plotFilePath,'rb')
     6008                    fp.seek(imgDict[plt])
     6009                    im = pickle.load(fp)
     6010                    G2plt.PlotRawImage(G2frame,im,plt)
     6011                    fp.close()
     6012                else:
     6013                    plotLbls = []
     6014                    plotVals = []
     6015                    for lbl,row in zip(lbls,yvals): # deal with <=0 costs
     6016                        if sum(row**2) == 0: continue # drop if all zeros
     6017                        if min(row) <= 0:
     6018                            row = np.where(row>0,row,min(row[np.where(row>0)])/10.)
     6019                        plotLbls.append(lbl)
     6020                        plotVals.append([steps,np.log10(row)])
     6021                    title = 'fullrmc residuals for '+pName
     6022                    G2plt.PlotXY(G2frame,plotVals,
     6023                                labelX='generated steps',
     6024                                labelY=r'$log_{10}$ ($\mathsf{\chi^2})$',
     6025                                newPlot=True,Title=title,
     6026                                lines=True,names=plotLbls)
     6027            return
    61616028        else:
    61626029            generalData = data['General']
Note: See TracChangeset for help on using the changeset viewer.