Changeset 4951


Ignore:
Timestamp:
Jun 11, 2021 8:42:20 AM (2 years 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.

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIdataGUI.py

    r4910 r4951  
    34103410    def OnPowderFPA(self,event):
    34113411        'Perform FPA simulation/peak fitting'
     3412        # if GSASIIpath.GetConfigValue('debug'):
     3413        #     print('Debug: reloading G2fpa')
     3414        #     import imp
     3415        #     imp.reload(G2fpa)
    34123416        G2fpa.GetFPAInput(self)
    34133417       

  • 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']

  • trunk/GSASIIplot.py

    r4950 r4951  
    83518351        Page.canvas.draw()
    83528352               
     8353##### PlotRawImage ################################################################################
     8354def PlotRawImage(G2frame,image,label,newPlot=False):
     8355    '''Plot an image without axes etc.
     8356    '''
     8357    new,plotNum,Page,Plot,lim = G2frame.G2plotNB.FindPlotTab(label,'mpl')
     8358    Plot.remove() # delete original axes
     8359    Plot = Page.figure.add_axes([0.0, 0.0, 1., 1.]) # fill page & don't show
     8360    Plot.axis('off')
     8361    Plot.imshow(image)
     8362    Page.canvas.draw()
     8363   
    83538364##### PlotTRImage ################################################################################
    83548365def PlotTRImage(G2frame,tax,tay,taz,newPlot=False):

  • trunk/GSASIIpwd.py

    r4919 r4951  
    576576                data['Ruland'],data['Sample Bkg.']['Mult'],res['fun'],msg))
    577577        else:
    578             G2fil.G2Print('  end:   Flat Bkg={:.1f}, BackRatio={:.3f}, Ruland={:.3f}) *** {} ***\n'.format(
     578            G2fil.G2Print('  end:   Flat Bkg={:.1f}, BackRatio={:.3f}, Ruland={:.3f} RMS:{:.4f}) *** {} ***\n'.format(
    579579                data['Flat Bkg'],data['BackRatio'],data['Ruland'],res['fun'],msg))
    580580    return res
    … …  
    28802880#     dspaces = [0.5/np.sqrt(G2lat.calc_rDsq2(H,G)) for H in np.eye(3)]
    28812881#     return min(dspaces)
    2882    
     2882
     2883def findfullrmc():
     2884    '''Find where fullrmc is installed. Tries the following:
     2885   
     2886         1. Returns the Config var 'fullrmc_exec', if defined. No check
     2887            is done that the interpreter has fullrmc
     2888         2. The current Python interpreter if fullrmc can be imported
     2889            and fullrmc is version 5+
     2890         3. The path is checked for a fullrmc image as named by Bachir
     2891
     2892    :returns: the full path to a python executable that is assumed to
     2893      have fullrmc installed or None, if it was not found.
     2894    '''
     2895    is_exe = lambda fpath: os.path.isfile(fpath) and os.access(fpath, os.X_OK)
     2896    if GSASIIpath.GetConfigValue('fullrmc_exec') is not None and is_exe(
     2897            GSASIIpath.GetConfigValue('fullrmc_exec')):
     2898        return GSASIIpath.GetConfigValue('fullrmc_exec')
     2899    try:
     2900        import fullrmc
     2901        if int(fullrmc.__version__.split('.')[0]) >= 5:
     2902            return sys.executable
     2903    except:
     2904        pass
     2905    pathlist = os.environ["PATH"].split(os.pathsep)
     2906    for p in (GSASIIpath.path2GSAS2,GSASIIpath.binaryPath,os.getcwd(),
     2907                  os.path.split(sys.executable)[0]):
     2908        if p not in pathlist: pathlist.insert(0,p)
     2909    import glob
     2910    for p in pathlist:
     2911        if sys.platform == "darwin":
     2912            lookfor = "fullrmc*macOS*i386-64bit"
     2913        elif sys.platform == "win32":
     2914            lookfor = "fullrmc*.exe"
     2915        else:
     2916            lookfor = "fullrmc*"
     2917        fl = glob.glob(lookfor)
     2918        if len(fl) > 0:
     2919            return os.path.abspath(sorted(fl)[0])
     2920       
    28832921def MakefullrmcRun(pName,Phase,RMCPdict):
     2922    '''Creates a script to run fullrmc. Returns the name of the file that was
     2923    created.
     2924    '''
    28842925    BondList = {}
    28852926    for k in RMCPdict['Pairs']:
    … …  
    28932934            angle[i] = angle[i].strip()
    28942935        AngleList.append(angle)
    2895     rmin = RMCPdict['min Contact']
     2936    # rmin = RMCPdict['min Contact']
    28962937    cell = Phase['General']['Cell'][1:7]
    28972938    SymOpList = G2spc.AllOps(Phase['General']['SGData'])[0]
    … …  
    29012942        el = ''.join([i for i in atom[ct] if i.isalpha()])
    29022943        atomsList.append([el] + atom[cx:cx+4])
    2903     rname = pName+'-fullrmc.py'
     2944    projDir,projName = os.path.split(pName)
     2945    scrname = pName+'-fullrmc.py'
    29042946    restart = '%s_restart.pdb'%pName
    29052947    Files = RMCPdict['files']
    29062948    rundata = ''
    2907     rundata += '#### fullrmc %s file; edit by hand if you so choose #####\n'%rname
     2949    rundata += '#### fullrmc %s file; edit by hand if you so choose #####\n'%scrname
    29082950    rundata += '# created in '+__file__+" v"+filversion.split()[1]
    29092951    rundata += dt.datetime.strftime(dt.datetime.now()," at %Y-%m-%dT%H:%M\n")
    … …  
    29122954import os,glob
    29132955import time
    2914 #import matplotlib.pyplot as plt
     2956import pickle
    29152957import numpy as np
    29162958from fullrmc.Core import Collection
    … …  
    29232965from fullrmc.Constraints.DihedralAngleConstraints import DihedralAngleConstraint
    29242966from fullrmc.Generators.Swaps import SwapPositionsGenerator
    2925 
    2926 ### When True, erases an existing enging to provide a fresh start
    2927 FRESH_START = {}
     2967def writeHeader(ENGINE,statFP):
     2968    statFP.write('generated-steps, total-error, ')
     2969    for c in ENGINE.constraints:
     2970        statFP.write(c.constraintName)
     2971        statFP.write(', ')
     2972    statFP.write('\\n')
     2973    statFP.flush()
     2974   
     2975def writeCurrentStatus(ENGINE,statFP,plotF):
     2976    statFP.write(str(ENGINE.generated))
     2977    statFP.write(', ')
     2978    statFP.write(str(ENGINE.totalStandardError))
     2979    statFP.write(', ')
     2980    for c in ENGINE.constraints:
     2981        statFP.write(str(c.standardError))
     2982        statFP.write(', ')
     2983    statFP.write('\\n')
     2984    statFP.flush()
     2985    mpl.use('agg')
     2986    fp = open(plotF,'wb')
     2987    for c in ENGINE.constraints:
     2988        p = c.plot(show=False)
     2989        p[0].canvas.draw()
     2990        image = p[0].canvas.buffer_rgba()
     2991        pickle.dump(c.constraintName,fp)
     2992        pickle.dump(np.array(image),fp)
     2993    fp.close()
     2994'''
     2995    rundata += '''
     2996### When True, erases an existing engine to provide a fresh start
     2997FRESH_START = {:}
     2998dirName = "{:}"
     2999prefix = "{:}"
     3000project = prefix + "-fullrmc"
    29283001time0 = time.time()
    2929 '''.format(RMCPdict['ReStart'][0])
     3002'''.format(RMCPdict['ReStart'][0],projDir,projName)
    29303003   
    29313004    rundata += '# setup structure\n'
    … …  
    29363009
    29373010    rundata += '\n# initialize engine\n'
    2938     rundata += 'engineFileName = "%s.rmc"\n'%pName
    2939 
    2940     rundata += '''\n# check Engine exists if so (and not FRESH_START) load it
     3011    rundata += '''
     3012engineFileName = os.path.join(dirName, project + '.rmc')
     3013projectStats = os.path.join(dirName, project + '.stats')
     3014projectPlots = os.path.join(dirName, project + '.plots')
     3015pdbFile = os.path.join(dirName, project + '_restart.pdb')
     3016# check Engine exists if so (and not FRESH_START) load it
    29413017# otherwise build it
    29423018ENGINE = Engine(path=None)
    … …  
    29653041            sfwt = 'atomicNumber'
    29663042        if 'G(r)' in File:
    2967             rundata += '    GR = np.loadtxt("%s").T\n'%filDat[0]
     3043            rundata += '    GR = np.loadtxt(os.path.join(dirName,"%s")).T\n'%filDat[0]
    29683044            if filDat[3] == 0:
    29693045                rundata += '''    # read and xform G(r) as defined in RMCProfile
    … …  
    29823058            rundata += '    GofR.set_limits((None, rmax))\n'
    29833059        elif '(Q)' in File:
    2984             rundata += '    SOQ = np.loadtxt("%s").T\n'%filDat[0]
     3060            rundata += '    SOQ = np.loadtxt(os.path.join(dirName,"%s")).T\n'%filDat[0]
    29853061            if filDat[3] == 0:
    29863062                rundata += '    # Read & xform F(Q) as defined in RMCProfile to S(Q)-1\n'
    … …  
    30153091               '("element","{1}","{0}","{2}",{5},{6},{5},{6},{3},{4}),\n'.format(*item))
    30163092        rundata += '             ])\n'
    3017     rundata += '    for f in glob.glob("{}_*.log"): os.remove(f)\n'.format(pName)
    30183093    rundata += '''
     3094    for f in glob.glob(os.path.join(dirName,prefix+"_*.log")): os.remove(f)
    30193095    ENGINE.save()
    30203096else:
    30213097    ENGINE = ENGINE.load(path=engineFileName)
     3098    rmax = min( [ENGINE.boundaryConditions.get_a(), ENGINE.boundaryConditions.get_b(), ENGINE.boundaryConditions.get_c()] ) /2.
     3099
     3100ENGINE.set_log_file(os.path.join(dirName,prefix))
    30223101'''
    3023     rundata += 'ENGINE.set_log_file("{}")\n'.format(pName)
    30243102    if RMCPdict['Swaps']:
    30253103        rundata += '\n#set up for site swaps\n'
    … …  
    30553133    rundata += 'for c in ENGINE.constraints:  # loop over predefined constraints\n'
    30563134    rundata += '    if type(c) is fPDF.PairDistributionConstraint:\n'
    3057     rundata += '        c.set_variance_squared(1./wtDict["Pair-"+c.weighting])\n'
     3135    # rundata += '        c.set_variance_squared(1./wtDict["Pair-"+c.weighting])\n'
    30583136    rundata += '        c.set_limits((None,rmax))\n'
    30593137    if RMCPdict['FitScale']:
    30603138        rundata += '        c.set_adjust_scale_factor((10, 0.01, 100.))\n'
    3061     rundata += '    elif type(c) is ReducedStructureFactorConstraint:\n'
    3062     rundata += '        c.set_variance_squared(1./wtDict["Struct-"+c.weighting])\n'
     3139    # rundata += '        c.set_variance_squared(1./wtDict["Struct-"+c.weighting])\n'
    30633140    if RMCPdict['FitScale']:
     3141        rundata += '    elif type(c) is ReducedStructureFactorConstraint:\n'
    30643142        rundata += '        c.set_adjust_scale_factor((10, 0.01, 100.))\n'
    30653143    # torsions difficult to implement, must be internal to cell & named with
    … …  
    30733151    #         rundata += '    %s\n'%str(tuple(torsion))
    30743152    #     rundata += '        ]})\n'           
    3075     rundata += '\n# setup runs for fullrmc\n'
    3076 
    3077     rundata += 'steps = 10000\n'
     3153    rundata += '''
     3154if FRESH_START:
     3155    # initialize engine with one step to get starting config energetics
     3156    ENGINE.run(restartPdb=pdbFile,numberOfSteps=1, saveFrequency=1)
     3157    statFP = open(projectStats,'w')
     3158    writeHeader(ENGINE,statFP)
     3159    writeCurrentStatus(ENGINE,statFP,projectPlots)
     3160else:
     3161    statFP = open(projectStats,'a')
     3162
     3163# setup runs for fullrmc
     3164'''
     3165    rundata += 'steps = {}\n'.format(RMCPdict['Steps/cycle'])
    30783166    rundata += 'for _ in range({}):\n'.format(RMCPdict['Cycles'])
    30793167    rundata += '    ENGINE.set_groups_as_atoms()\n'
    30803168    rundata += '    expected = ENGINE.generated+steps\n'
    30813169   
    3082     rundata += '    ENGINE.run(restartPdb="%s",numberOfSteps=steps, saveFrequency=1000)\n'%restart
     3170    rundata += '    ENGINE.run(restartPdb=pdbFile,numberOfSteps=steps, saveFrequency=steps)\n'
     3171    rundata += '    writeCurrentStatus(ENGINE,statFP,projectPlots)\n'
    30833172    rundata += '    if ENGINE.generated != expected: break # run was stopped\n'
     3173    rundata += 'statFP.close()\n'
    30843174    rundata += 'print("ENGINE run time %.2f s"%(time.time()-time0))\n'
    3085     rfile = open(rname,'w')
     3175    rfile = open(scrname,'w')
    30863176    rfile.writelines(rundata)
    30873177    rfile.close()
    3088     return rname
     3178    return scrname
    30893179   
    30903180def GetRMCBonds(general,RMCPdict,Atoms,bondList):
    … …  
    43244414    G2fil.G2Print ("OK")
    43254415    plotter.StartEventLoop()
     4416   
     4417#    GSASIIpath.SetBinaryPath(True,False)
     4418#    print('found',findfullrmc())
Note: See TracChangeset for help on using the changeset viewer.