Changeset 4518


Ignore:
Timestamp:
Jul 12, 2020 8:59:57 PM (3 years ago)
Author:
toby
Message:

start on fullrmc update for 5.0; fix for rigid bodies

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIphsGUI.py

    r4513 r4518  
    3232import platform
    3333import os
    34 import os.path
    3534import wx
    3635import wx.grid as wg
    … …  
    14241423    :param float cMin: Minimum along the *c* direction (fractional units).
    14251424       Defaults to 0.
    1426     :
    14271425    '''
    14281426
    … …  
    46284626
    46294627    def UpdateRMC():
    4630         ''' Present the controls for running fullrmc or RMCprofile
     4628        ''' Present the controls for running fullrmc or RMCProfile
    46314629        '''
    46324630        global runFile
    … …  
    47054703                wx.CallAfter(UpdateRMC)
    47064704               
    4707             def OnSwapAtSel(event):
    4708                 Obj = event.GetEventObject()
    4709                 swap,i = Indx[Obj.GetId()]
    4710                 RMCPdict['Swaps'][swap][i] = Obj.GetStringSelection()
    4711                                        
    47124705            Indx = {}
    47134706            atChoice = RMCPdict['atSeq']           
    47144707            if G2frame.RMCchoice == 'fullrmc':
    47154708                atChoice = atNames
    4716             swapSizer = wx.FlexGridSizer(4,5,5)
    4717             swapLabels = [' ','Atom-A','Atom-B',' Swap prob.']
     4709            swapSizer = wx.FlexGridSizer(6,5,5)
     4710            swapLabels = [' ','Atom-A','Atom-B',' Swap prob.',' ','delete']
    47184711            for lab in swapLabels:
    47194712                swapSizer.Add(wx.StaticText(G2frame.FRMC,label=lab),0,WACV)
    47204713            for ifx,swap in enumerate(RMCPdict['Swaps']):
    4721                 delBtn = wx.Button(G2frame.FRMC,label='Delete')
     4714                swapSizer.Add((20,-1))
     4715                for i in [0,1]:
     4716                    if swap[i] not in atChoice: swap[i] = atChoice[0]
     4717                    atmSel = G2G.EnumSelector(G2frame.FRMC,swap,i,atChoice)
     4718                    swapSizer.Add(atmSel,0,WACV)
     4719                swapSizer.Add(G2G.ValidatedTxtCtrl(G2frame.FRMC,swap,2,xmin=0.01,xmax=0.5,size=(50,25)),0,WACV)
     4720                swapSizer.Add((20,-1))
     4721                delBtn = wx.Button(G2frame.FRMC,label='Del',style=wx.BU_EXACTFIT)
    47224722                delBtn.Bind(wx.EVT_BUTTON,OnDelSwap)
    47234723                Indx[delBtn.GetId()] = ifx
    47244724                swapSizer.Add(delBtn,0,WACV)
    4725                 for i in [0,1]:
    4726                     atmSel = wx.ComboBox(G2frame.FRMC,choices=atChoice,style=wx.CB_DROPDOWN|wx.TE_READONLY)
    4727                     atmSel.SetStringSelection(swap[i])
    4728                     atmSel.Bind(wx.EVT_COMBOBOX,OnSwapAtSel)
    4729                     Indx[atmSel.GetId()] = [ifx,i]
    4730                     swapSizer.Add(atmSel,0,WACV)
    4731                 swapSizer.Add(G2G.ValidatedTxtCtrl(G2frame.FRMC,swap,2,xmin=0.01,xmax=0.5,size=(50,25)),0,WACV)
    47324725            return swapSizer
    47334726       
    … …  
    48054798                Obj = event.GetEventObject()
    48064799                fil = Indx[Obj.GetId()]
    4807                 RMCPdict['files'][fil][0] = 'Select'
     4800                RMCPdict['files'][fil][0] = 'Select file'
    48084801                RMCPdict['ReStart'][0] = True
    48094802                wx.CallAfter(UpdateRMC)
    … …  
    48174810            titleSizer.Add(fitscale,0,WACV)
    48184811            mainSizer.Add(titleSizer,0,WACV)
    4819             ncol= 6
    4820             Heads = ['Name','File','Format','Weight','Plot','Delete']
    48214812            if G2frame.RMCchoice == 'fullrmc':
    48224813                mainSizer.Add(wx.StaticText(G2frame.FRMC,
    48234814                    label=' NB: fullrmc data files must be 2 columns; all other lines preceeded by "#". Edit before use.'),0,WACV)
    4824                 Heads = ['Name','File','Weight','Plot','Corr','Delete']
    4825             fileSizer = wx.FlexGridSizer(ncol,5,5)
     4815                Heads = ['Name','File','Weight','type','Plot','Delete']
     4816                fileSizer = wx.FlexGridSizer(6,5,5)
     4817                Formats = ['RMC','GUDRUN','STOG']
     4818                for head in Heads:
     4819                    fileSizer.Add(wx.StaticText(G2frame.FRMC,label=head),0,WACV)
     4820                for fil in RMCPdict['files']:
     4821                    fileSizer.Add(wx.StaticText(G2frame.FRMC,label=fil),0,WACV)
     4822                    Rfile = RMCPdict['files'][fil][0]
     4823                    filSel = wx.Button(G2frame.FRMC,label=Rfile)
     4824                    filSel.Bind(wx.EVT_BUTTON,OnFileSel)
     4825                    Indx[filSel.GetId()] = fil
     4826                    fileSizer.Add(filSel,0,WACV)
     4827                    if Rfile and os.path.exists(Rfile): # in case .gpx file is moved away from G(R), F(Q), etc. files
     4828                        fileSizer.Add(G2G.ValidatedTxtCtrl(G2frame.FRMC,RMCPdict['files'][fil],1,
     4829                                                               size=(50,25)),0,WACV)
     4830                        #patch
     4831                        if len(RMCPdict['files'][fil]) < 4:
     4832                            RMCPdict['files'][fil].append(0)
     4833                        if len(RMCPdict['files'][fil]) < 5:
     4834                            RMCPdict['files'][fil].append(True)
     4835                        #end patch
     4836                        if 'G(r)' in fil:
     4837                            choices = 'G(r)-RMCProfile','G(r)-PDFFIT','g(r)'
     4838                            if type(RMCPdict['files'][fil][3]) is bool: RMCPdict['files'][fil][3] = 0
     4839                            fmtTyp = G2G.G2ChoiceButton(G2frame.FRMC,choices,RMCPdict['files'][fil],3)
     4840                        elif 'F(Q)' in fil:
     4841                            choices = 'F(Q)-RMCProfile','S(Q)-PDFFIT'
     4842                            if type(RMCPdict['files'][fil][3]) is bool: RMCPdict['files'][fil][3] = 0
     4843                            fmtTyp = G2G.G2ChoiceButton(G2frame.FRMC,choices,RMCPdict['files'][fil],3)
     4844                        else:
     4845                            fmtTyp = (-1,-1)
     4846                        fileSizer.Add(fmtTyp,0,WACV)
     4847                        plotBtn = wx.Button(G2frame.FRMC,label='Plot',style=wx.BU_EXACTFIT)
     4848                        plotBtn.Bind(wx.EVT_BUTTON,OnPlotBtn)
     4849                        Indx[plotBtn.GetId()] = fil
     4850                        fileSizer.Add(plotBtn,0,WACV)
     4851                        delBtn = wx.Button(G2frame.FRMC,label='Del',style=wx.BU_EXACTFIT)
     4852                        delBtn.Bind(wx.EVT_BUTTON,OnDelBtn)
     4853                        Indx[delBtn.GetId()] = fil
     4854                        fileSizer.Add(delBtn,0,WACV)
     4855                        if 'F(Q)' in fil:
     4856                            fileSizer.Add((-1,-1),0)
     4857                            corrChk = wx.CheckBox(G2frame.FRMC,label='Apply sinc convolution? ')
     4858                            corrChk.SetValue(RMCPdict['files'][fil][4])
     4859                            Indx[corrChk.GetId()] = fil
     4860                            corrChk.Bind(wx.EVT_CHECKBOX,OnCorrChk)
     4861                            fileSizer.Add(corrChk,0,WACV)
     4862                            fileSizer.Add((-1,-1),0)
     4863                            fileSizer.Add((-1,-1),0)
     4864                            fileSizer.Add((-1,-1),0)
     4865                            fileSizer.Add((-1,-1),0)
     4866                    else:
     4867                        RMCPdict['files'][fil][0] = 'Select file' # set filSel?
     4868                        fileSizer.Add((-1,-1),0)
     4869                        fileSizer.Add((-1,-1),0)
     4870                        fileSizer.Add((-1,-1),0)
     4871                        fileSizer.Add((-1,-1),0)
     4872                return fileSizer
     4873            # RMCProfile
     4874            Heads = ['Name','File','Format','Weight','Plot','Delete']
     4875            fileSizer = wx.FlexGridSizer(6,5,5)
    48264876            Formats = ['RMC','GUDRUN','STOG']
    48274877            for head in Heads:
    … …  
    48354885                fileSizer.Add(filSel,0,WACV)
    48364886                if Rfile and os.path.exists(Rfile): #incase .gpx file is moved away from G(R), F(Q), etc. files
    4837                     if G2frame.RMCchoice == 'RMCProfile':                       
    4838                         nform = 3
    4839                         if 'Xray' in fil: nform = 1
    4840                         fileFormat = wx.ComboBox(G2frame.FRMC,choices=Formats[:nform],style=wx.CB_DROPDOWN|wx.TE_READONLY)
    4841                         fileFormat.SetStringSelection(RMCPdict['files'][fil][3])
    4842                         Indx[fileFormat.GetId()] = fil
    4843                         fileFormat.Bind(wx.EVT_COMBOBOX,OnFileFormat)
    4844                         fileSizer.Add(fileFormat,0,WACV)
     4887                    nform = 3
     4888                    if 'Xray' in fil: nform = 1
     4889                    fileFormat = wx.ComboBox(G2frame.FRMC,choices=Formats[:nform],style=wx.CB_DROPDOWN|wx.TE_READONLY)
     4890                    fileFormat.SetStringSelection(RMCPdict['files'][fil][3])
     4891                    Indx[fileFormat.GetId()] = fil
     4892                    fileFormat.Bind(wx.EVT_COMBOBOX,OnFileFormat)
     4893                    fileSizer.Add(fileFormat,0,WACV)
    48454894                    fileSizer.Add(G2G.ValidatedTxtCtrl(G2frame.FRMC,RMCPdict['files'][fil],1),0,WACV)
    4846                     plotBtn = wx.Button(G2frame.FRMC,label='Plot?')
     4895                    plotBtn = wx.Button(G2frame.FRMC,label='Plot?',style=wx.BU_EXACTFIT)
    48474896                    plotBtn.Bind(wx.EVT_BUTTON,OnPlotBtn)
    48484897                    Indx[plotBtn.GetId()] = fil
    48494898                    fileSizer.Add(plotBtn,0,WACV)
    4850                     if G2frame.RMCchoice == 'fullrmc':
    4851                         lab = 'Apply sinc convolution? '
    4852                         if 'G(r)' in fil:
    4853                             lab = 'Apply r*G(r) conversion? '
    4854                         corrChk = wx.CheckBox(G2frame.FRMC,label=lab)
    4855                         #patch
    4856                         if len(RMCPdict['files'][fil]) < 4:
    4857                             RMCPdict['files'][fil].append(True)
    4858                         #end patch
    4859                         corrChk.SetValue(RMCPdict['files'][fil][3])
    4860                         Indx[corrChk.GetId()] = fil
    4861                         corrChk.Bind(wx.EVT_CHECKBOX,OnCorrChk)
    4862                         fileSizer.Add(corrChk,0,WACV)
    4863                     delBtn = wx.Button(G2frame.FRMC,label='Delete')
     4899                    delBtn = wx.Button(G2frame.FRMC,label='Del',style=wx.BU_EXACTFIT)
    48644900                    delBtn.Bind(wx.EVT_BUTTON,OnDelBtn)
    48654901                    Indx[delBtn.GetId()] = fil
    … …  
    49014937                     caption='fullrmc not installed',style=wx.ICON_INFORMATION)
    49024938                return
     4939            if int(fullrmc.__version__.split('.')[0]) < 5:
     4940                wx.MessageBox('This module requires fullrmc >= 5.0. You have {}.'.format(fullrmc.__version__) +
     4941                        ' Revert to GSAS-II version <=4517 to use older versions of fullrmc. ',
     4942                     caption='fullrmc to old',style=wx.ICON_INFORMATION)
     4943                return
    49034944            mainSizer.Add(wx.StaticText(G2frame.FRMC,label=''' "Fullrmc, a Rigid Body Reverse Monte Carlo Modeling Package Enabled with Machine Learning and Artificial Intelligence",     
    49044945 B. Aoun, Jour. Comp. Chem. 2016, 37, 1102-1111. doi: https://doi.org/10.1002/jcc.24304
    … …  
    49074948            G2frame.dataWindow.FRMCDataEdit.Enable(G2G.wxID_RUNRMC,True)
    49084949            G2frame.dataWindow.FRMCDataEdit.Enable(G2G.wxID_VIEWRMC,True)
    4909             mainSizer.Add(wx.StaticText(G2frame.FRMC,label=' fullrmc big box starting pdb file preparation:'),0,WACV)
     4950            #mainSizer.Add(wx.StaticText(G2frame.FRMC,label=' fullrmc big box starting pdb file preparation:'),0,WACV)
    49104951            if not data['RMC']['fullrmc']:
    49114952                Atypes = [atype.split('+')[0].split('-')[0] for atype in data['General']['AtomTypes']]
    … …  
    49174958                    Pairs += pair
    49184959                Pairs = {pairs:[0.0,0.0,0.0] for pairs in Pairs}
    4919                 files = {'Neutron real space data; G(r): ':['Select',0.05,'G(r)',True],
    4920                           'Neutron reciprocal space data; F(Q): ':['Select',0.05,'F(Q)',True],
    4921                           'Xray real space data; G(r): ':['Select',0.01,'G(r)',True],
    4922                           'Xray reciprocal space data; F(Q): ':['Select',0.01,'F(Q)',True],}
    4923                 data['RMC']['fullrmc'] = {'SuperCell':[1,1,1],'Box':[10.,10.,10.],'aTypes':aTypes,'byMolec':True,
    4924                     'Natoms':1,'atSeq':atSeq,'Pairs':Pairs,'files':files,'ReStart':[False,False],'Cycles':1,
     4960                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],
     4962                          '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],}
     4964                data['RMC']['fullrmc'] = {'SuperCell':[1,1,1],'Box':[10.,10.,10.],'aTypes':aTypes,
     4965                    'atSeq':atSeq,'Pairs':Pairs,'files':files,'ReStart':[False,False],'Cycles':1,
    49254966                    'Swaps':[],'useBVS':False,'FitScale':False,'AveCN':[],'FxCN':[],'Angles':[],'Angle Weight':1.e-5,
    49264967                    'moleculePdb':'Select','targetDensity':1.0,'maxRecursion':10000,'Torsions':[],'Torsion Weight':1.e-5,
    4927                     'atomPDB':'','Bond Weight':1.e-5,'min Contact':1.5}
     4968                    'Bond Weight':1.e-5,'min Contact':1.5,'periodicBound':True}
    49284969            RMCPdict = data['RMC']['fullrmc']
    49294970
    … …  
    49605001                    RMCPdict['moleculePdb'] = fName
    49615002                    pdbButton.SetLabel(fName)
    4962                    
    4963             def OnMakePDB(event):
    4964                 if not ifBox:
    4965                     generalData = data['General']
    4966                     pName = generalData['Name'].replace(' ','_')
    4967                     pdbname = G2pwd.MakefullrmcPDB(pName,data,RMCPdict)
    4968                     RMCPdict['atomPDB'] = pdbname
    4969                     RMCPdict['ReStart'][0] = True
    4970                     print(pdbname+ ' written')
    4971                     return                   
    4972                 if RMCPdict['moleculePdb'] == 'Select':
    4973                     wx.MessageDialog(G2frame,' You must select a source pdb file first','PDB generation error',
    4974                         wx.ICON_EXCLAMATION).ShowModal()
    4975                     return
    4976                 dlg = wx.MessageDialog(G2frame,'''
    4977 Warning - this step can take more than an hour; do you want to proceed?
    4978 It will be run as a separate process, and a result is required for fullrmc.
    4979 Make sure your parameters are correctly set.
    4980                     ''','Make big box pdb file',wx.YES_NO|wx.ICON_QUESTION)
    4981                 try:
    4982                     result = dlg.ShowModal()
    4983                     if result in [wx.ID_YES,]:
    4984                         pdpy = G2pwd.MakePdparse(RMCPdict)
    4985                         import subprocess as sb
    4986                         batch = open('pdbparse.bat','w')
    4987                         batch.write('CALL '+sys.exec_prefix+'\\Scripts\\activate\n')
    4988                         batch.write(sys.exec_prefix+'\\python.exe %s\n'%pdpy)
    4989                         batch.write('pause')
    4990                         batch.close()
    4991                         sb.Popen('pdbparse.bat',creationflags=sb.CREATE_NEW_CONSOLE).pid
    4992                         RMCPdict['ReStart'][0] = True
    4993                     else:
    4994                         print(' Make PDB file for fullrmc abandonded')
    4995                 finally:
    4996                     dlg.Destroy()
    4997                 wx.CallAfter(UpdateRMC)
    49985003               
    49995004            def OnAddAngle(event):
    5000                 RMCPdict['Angles'].append(['','','',0.,0.])
     5005                RMCPdict['Angles'].append(['','','',0.,0.,0.,0.])
    50015006                wx.CallAfter(UpdateRMC)
    50025007               
    … …  
    50135018                    wx.CallAfter(UpdateRMC)
    50145019                   
    5015                 def OnAngleAtSel(event):
    5016                     Obj = event.GetEventObject()
    5017                     angle,i = Indx[Obj.GetId()]
    5018                     RMCPdict['Angles'][angle][i] = Obj.GetStringSelection()
     5020                # def OnAngleAtSel(event):
     5021                #     Obj = event.GetEventObject()
     5022                #     angle,i = Indx[Obj.GetId()]
     5023                #     RMCPdict['Angles'][angle][i] = Obj.GetStringSelection()
    50195024                                           
    50205025                def SetRestart1(invalid,value,tc):
    … …  
    50235028                Indx = {}
    50245029                atChoice = [atm for atm in RMCPdict['atSeq'] if 'Va' not in atm]
    5025                 angleSizer = wx.FlexGridSizer(6,5,5)
    5026                 fxcnLabels = [' ','Atom-A','Atom-B','Atom-C',' min angle',' max angle']
    5027                 for lab in fxcnLabels:
    5028                     angleSizer.Add(wx.StaticText(G2frame.FRMC,label=lab),0,WACV)
     5030                angleSizer = wx.GridBagSizer(0,5)
     5031                fxcnLabels1 = [' ',' ','Central','',None,'angle restraint values (deg)',None,'search distance (A)']
     5032                fxcnLabels2 = [' ','Atom-A','Atom','Atom-C','min','max','from','to']
     5033                for i in range(8):
     5034                    if fxcnLabels1[i]:
     5035                        cspan=1
     5036                        coloff = 0
     5037                        if fxcnLabels1[i-1] is None:
     5038                            cspan=2
     5039                            coloff = 1
     5040                        angleSizer.Add(wx.StaticText(G2frame.FRMC,label=fxcnLabels1[i]),
     5041                                       (0,i-coloff),(1,cspan))
     5042                    if fxcnLabels2[i]:
     5043                        angleSizer.Add(wx.StaticText(G2frame.FRMC,wx.ID_ANY,
     5044                                                    label=fxcnLabels2[i],style=wx.CENTER),
     5045                                       (1,i))
     5046                row = 1
    50295047                for ifx,angle in enumerate(RMCPdict['Angles']):
    5030                     delBtn = wx.Button(G2frame.FRMC,label='Delete')
     5048                    row += 1
     5049                    angleSizer.Add((30,-1),(row,0))
     5050                    for i in range(3):
     5051                        if angle[i] not in atChoice: angle[i] = atChoice[0]
     5052                        atmSel = G2G.EnumSelector(G2frame.FRMC,angle,i,atChoice)
     5053                        angleSizer.Add(atmSel,(row,1+i))
     5054                    for i in range(4):
     5055                        if i == 0:
     5056                            xmin,xmax=0.,180.
     5057                        elif i == 2:
     5058                            xmin,xmax=0.1,6.
     5059                        angleSizer.Add(
     5060                            G2G.ValidatedTxtCtrl(G2frame.FRMC,angle,3+i,xmin=xmin,xmax=xmax,
     5061                                OnLeave=SetRestart1,size=(50,25)),(row,4+i))
     5062                    delBtn = wx.Button(G2frame.FRMC,label='Del',style=wx.BU_EXACTFIT)
    50315063                    delBtn.Bind(wx.EVT_BUTTON,OnDelAngle)
    50325064                    Indx[delBtn.GetId()] = ifx
    5033                     angleSizer.Add(delBtn,0,WACV)
    5034                     for i in [0,1,2]:
    5035                         atmSel = wx.ComboBox(G2frame.FRMC,choices=atChoice,style=wx.CB_DROPDOWN|wx.TE_READONLY)
    5036                         atmSel.SetStringSelection(angle[i])
    5037                         atmSel.Bind(wx.EVT_COMBOBOX,OnAngleAtSel)
    5038                         Indx[atmSel.GetId()] = [ifx,i]
    5039                         angleSizer.Add(atmSel,0,WACV)
    5040                     angleSizer.Add(G2G.ValidatedTxtCtrl(G2frame.FRMC,angle,3,xmin=0.,xmax=180.,OnLeave=SetRestart1,size=(50,25)),0,WACV)
    5041                     angleSizer.Add(G2G.ValidatedTxtCtrl(G2frame.FRMC,angle,4,xmin=0.,xmax=180.,OnLeave=SetRestart1,size=(50,25)),0,WACV)
     5065                    angleSizer.Add(delBtn,(row,9))
    50425066                return angleSizer
    50435067           
    … …  
    50835107            if 'moleculePdb' not in RMCPdict:
    50845108                RMCPdict.update({'moleculePdb':'Select','targetDensity':1.0,'maxRecursion':10000})
    5085             if 'byMolec' not in RMCPdict:
    5086                 RMCPdict['byMolec'] = True
    5087             if 'Natoms' not in RMCPdict:
    5088                 RMCPdict['Natoms'] = 1
    50895109            if 'FitScale' not in RMCPdict:
    50905110                RMCPdict['FitScale'] = False
    5091             if 'atomPDB' not in RMCPdict:
    5092                 RMCPdict['atomPDB'] = ''
     5111#            if 'atomPDB' not in RMCPdict:
     5112#                RMCPdict['atomPDB'] = ''
    50935113            if 'Angles' not in RMCPdict:
    50945114                RMCPdict.update({'Angles':[],'Angle Weight':1.e-5,'Bond Weight':1.e-5,'Torsions':[],'Torsion Weight':1.e-5})
    … …  
    50975117            if 'min Contact' not in RMCPdict:
    50985118                RMCPdict['min Contact'] = 1.5
     5119            if 'periodicBound' not in RMCPdict:
     5120                RMCPdict['periodicBound'] = True
    50995121#end patches
    51005122
    … …  
    51035125            atomData = data['Atoms']
    51045126            atNames = [atom[ct-1] for atom in atomData]
    5105             ifP1 = False
    5106             if generalData['SGData']['SpGrp'] == 'P 1':
    5107                 ifP1 = True               
     5127            # ifP1 = False
     5128            # if generalData['SGData']['SpGrp'] == 'P 1':
     5129            #     ifP1 = True               
    51085130            ifBox = False
    51095131            if 'macromolecular' in generalData['Type']:
    … …  
    51135135                lineSizer.Add(wx.StaticText(G2frame.FRMC,label=' Big box dimensions, %s:'%Angstr),0,WACV)               
    51145136                lineSizer.Add(GetBoxSizer(),0,WACV)
    5115             elif ifP1:
    5116                 lineSizer.Add(wx.StaticText(G2frame.FRMC,label=' Lattice multipliers:'),0,WACV)
    5117                 lineSizer.Add(GetSuperSizer(),0,WACV)
    5118                 lineSizer.Add(wx.StaticText(G2frame.FRMC,label=' Num. atoms per group '),0,WACV)
    5119                 lineSizer.Add(G2G.ValidatedTxtCtrl(G2frame.FRMC,RMCPdict,'Natoms',xmin=1,size=[40,25]),0,WACV)
    5120             else:
    5121                 lineSizer.Add(wx.StaticText(G2frame.FRMC,label=' Starting phase symmetry must be P 1; transform structure first'))
     5137#            elif ifP1:
     5138            lineSizer.Add(wx.StaticText(G2frame.FRMC,label=' Lattice multipliers:'),0,WACV)
     5139            lineSizer.Add(GetSuperSizer(),0,WACV)
     5140            lineSizer.Add((5,-1))
     5141#            lineSizer.Add(G2G.G2CheckBox(G2frame.FRMC,'Impose periodic boundaries',RMCPdict,'periodicBound'),
     5142#                              0,WACV)
     5143#            lineSizer.Add(wx.StaticText(G2frame.FRMC,label=' Num. atoms per group '),0,WACV)
     5144#            lineSizer.Add(G2G.ValidatedTxtCtrl(G2frame.FRMC,RMCPdict,'Natoms',xmin=1,size=[40,25]),0,WACV)
     5145#            else:
     5146#                lineSizer.Add(wx.StaticText(G2frame.FRMC,label=' Starting phase symmetry must be P 1; transform structure first'))
    51225147            mainSizer.Add(lineSizer,0,WACV)
    51235148            if ifBox:
    … …  
    51315156                molecSizer.Add(wx.StaticText(G2frame.FRMC,label=' max tries '),0,WACV)
    51325157                molecSizer.Add(G2G.ValidatedTxtCtrl(G2frame.FRMC,RMCPdict,'maxRecursion',xmin=1000,xmax=1000000,size=[60,25]),0,WACV)
    5133                 makePDB = wx.Button(G2frame.FRMC,label='Make big box PDB (slow!)')
    5134                 makePDB.Bind(wx.EVT_BUTTON,OnMakePDB)
    5135                 molecSizer.Add(makePDB,0,WACV)               
    51365158                mainSizer.Add(molecSizer,0,WACV)
    5137             elif ifP1:
    5138                 makePDB = wx.Button(G2frame.FRMC,label='Make big box PDB')
    5139                 makePDB.Bind(wx.EVT_BUTTON,OnMakePDB)
    5140                 mainSizer.Add(makePDB,0,WACV)
    5141             else:       #Abort because starting phase symmetry isn't P 1
    5142                 SetPhaseWindow(G2frame.FRMC,mainSizer)
    5143                 return
    5144                
    51455159            G2G.HorizontalLine(mainSizer,G2frame.FRMC)
    51465160            mainSizer.Add(wx.StaticText(G2frame.FRMC,label=' fullrmc run file preparation:'),0,WACV)
    … …  
    51595173            G2G.HorizontalLine(mainSizer,G2frame.FRMC)
    51605174            swapBox = wx.BoxSizer(wx.HORIZONTAL)
    5161             swapAdd = wx.Button(G2frame.FRMC,label='Add')
     5175            swapBox.Add(wx.StaticText(G2frame.FRMC,label='Atom swap probabiities: '),0,WACV)
     5176            swapAdd = wx.Button(G2frame.FRMC,label='Add',style=wx.BU_EXACTFIT)
    51625177            swapAdd.Bind(wx.EVT_BUTTON,OnAddSwap)
    51635178            swapBox.Add(swapAdd,0,WACV)
    5164             swapBox.Add(wx.StaticText(G2frame.FRMC,label=' Atom swap probabiities: '),0,WACV)
    51655179            mainSizer.Add(swapBox,0,WACV)       
    51665180            if len(RMCPdict['Swaps']):
    … …  
    51685182
    51695183            G2G.HorizontalLine(mainSizer,G2frame.FRMC)
    5170             mainSizer.Add(wx.StaticText(G2frame.FRMC,label=' Enter constraints && restraints:'),0,WACV)
     5184            mainSizer.Add(wx.StaticText(G2frame.FRMC,label='Geometry constraints && restraints'),0,WACV)
    51715185            distBox = wx.BoxSizer(wx.HORIZONTAL)
    51725186            distBox.Add(wx.StaticText(G2frame.FRMC,label=' Distance constraints, weight: :'),0,WACV)       
    … …  
    51755189            distBox.Add(G2G.ValidatedTxtCtrl(G2frame.FRMC,RMCPdict,'min Contact',xmin=0.,xmax=4.,size=(50,25)),0,WACV)           
    51765190            mainSizer.Add(distBox,0,WACV)
    5177             if RMCPdict['byMolec']:
    5178                 mainSizer.Add(GetPairSizer(RMCPdict),0,WACV)
    5179                
    5180                 angBox = wx.BoxSizer(wx.HORIZONTAL)
    5181                 angAdd = wx.Button(G2frame.FRMC,label='Add')
    5182                 angAdd.Bind(wx.EVT_BUTTON,OnAddAngle)
    5183                 angBox.Add(angAdd,0,WACV)
    5184                 angBox.Add(wx.StaticText(G2frame.FRMC,label=' A-B-C angle restraints, weight: '),0,WACV)
    5185                 angBox.Add(G2G.ValidatedTxtCtrl(G2frame.FRMC,RMCPdict,'Angle Weight',xmin=0.,xmax=100.,size=(50,25)),0,WACV)
    5186                 mainSizer.Add(angBox,0,WACV)
    5187                 if len(RMCPdict['Angles']):
    5188                     mainSizer.Add(GetAngleSizer(),0,WACV)
    5189                    
    5190                 torBox = wx.BoxSizer(wx.HORIZONTAL)
    5191                 torAdd = wx.Button(G2frame.FRMC,label='Add')
    5192                 torAdd.Bind(wx.EVT_BUTTON,OnAddTorsion)
    5193                 torBox.Add(torAdd,0,WACV)
    5194                 torBox.Add(wx.StaticText(G2frame.FRMC,label=' A-B-C-D torsion angle restraints, weight: '),0,WACV)
    5195                 torBox.Add(G2G.ValidatedTxtCtrl(G2frame.FRMC,RMCPdict,'Torsion Weight',xmin=0.,xmax=100.,size=(50,25)),0,WACV)
    5196                 mainSizer.Add(torBox,0,WACV)
    5197                 if len(RMCPdict['Torsions']):
    5198                     mainSizer.Add(GetTorsionSizer(),0,WACV)
     5191            mainSizer.Add(GetPairSizer(RMCPdict),0,WACV)
     5192
     5193            mainSizer.Add((-1,10))
     5194            angBox = wx.BoxSizer(wx.HORIZONTAL)
     5195            angBox.Add(wx.StaticText(G2frame.FRMC,label=' A-B-C angle restraints, weight: '),0,WACV)
     5196            angBox.Add(G2G.ValidatedTxtCtrl(G2frame.FRMC,RMCPdict,'Angle Weight',xmin=0.,xmax=100.,size=(50,25)),0,WACV)
     5197            angBox.Add((20,-1))
     5198            angAdd = wx.Button(G2frame.FRMC,label='Add',style=wx.BU_EXACTFIT)
     5199            angAdd.Bind(wx.EVT_BUTTON,OnAddAngle)
     5200            angBox.Add(angAdd,0,WACV)
     5201            mainSizer.Add(angBox,0,WACV)
     5202            if len(RMCPdict['Angles']):
     5203                mainSizer.Add(GetAngleSizer(),0,WACV)
     5204
     5205            # Torsions are probably not implemented correctly, hide them for now
     5206            # torBox = wx.BoxSizer(wx.HORIZONTAL)
     5207            # torAdd = wx.Button(G2frame.FRMC,label='Add')
     5208            # torAdd.Bind(wx.EVT_BUTTON,OnAddTorsion)
     5209            # torBox.Add(torAdd,0,WACV)
     5210            # torBox.Add(wx.StaticText(G2frame.FRMC,label=' A-B-C-D torsion angle restraints (intracell only), weight: '),0,WACV)
     5211            # torBox.Add(G2G.ValidatedTxtCtrl(G2frame.FRMC,RMCPdict,'Torsion Weight',xmin=0.,xmax=100.,size=(50,25)),0,WACV)
     5212            # mainSizer.Add(torBox,0,WACV)
     5213            # if len(RMCPdict['Torsions']):
     5214            #     mainSizer.Add(GetTorsionSizer(),0,WACV)
    51995215   
    52005216            G2G.HorizontalLine(mainSizer,G2frame.FRMC)
    … …  
    56005616    def OnSetupRMC(event):
    56015617        generalData = data['General']
    5602         pName = generalData['Name'].replace(' ','_')
    56035618        if not G2frame.GSASprojectfile:     #force a project save
    56045619            G2frame.OnFileSaveas(event)
    56055620        if G2frame.RMCchoice == 'fullrmc':
    5606             DisAglCtls = {}
    5607             if 'DisAglCtls' in generalData:
    5608                 DisAglCtls = generalData['DisAglCtls']
    5609             dlg = G2G.DisAglDialog(G2frame,DisAglCtls,generalData,Reset=False)
    5610             if dlg.ShowModal() == wx.ID_OK:
    5611                 DisAglCtls = dlg.GetData()
    5612                 if 'H' not in DisAglCtls['AtomTypes']:
    5613                     DisAglCtls['AtomTypes'].append('H')
    5614                     DisAglCtls['AngleRadii'].append(0.5)
    5615                     DisAglCtls['BondRadii'].append(0.5)
    5616             dlg.Destroy()
    5617             generalData['DisAglCtls'] = DisAglCtls
     5621            pName = G2frame.GSASprojectfile.split('.')[0] + '-' + generalData['Name']
     5622            pName = pName.replace(' ','_')
    56185623            G2frame.dataWindow.FRMCDataEdit.Enable(G2G.wxID_RUNRMC,True)
    56195624            RMCPdict = data['RMC']['fullrmc']
    5620             fname = G2pwd.MakefullrmcRun(pName,data,RMCPdict)
    5621             if fname is None:
    5622                 wx.MessageDialog(G2frame,' Big box pdb file is missing; fullrmc will not run','Missing pdb file',wx.OK).ShowModal()
    5623             else:   
    5624                 print('fullrmc file %s build completed'%fname)
     5625            # debug stuff
     5626            import imp
     5627            imp.reload(G2pwd)
     5628            # end debug stuff           
     5629            rname = G2pwd.MakefullrmcRun(pName,data,RMCPdict)
     5630            print('build of fullrmc file {} completed'.format(rname))
    56255631        else:
     5632            pName = generalData['Name'].replace(' ','_')
    56265633            G2frame.dataWindow.FRMCDataEdit.Enable(G2G.wxID_RUNRMC,True)
    56275634            RMCPdict = data['RMC']['RMCProfile']
    … …  
    56675674           
    56685675    def OnRunRMC(event):
    5669        
     5676        '''Run a previously created RMCProfile/fullrmc script
     5677        '''
    56705678        generalData = data['General']
    5671         pName = generalData['Name'].replace(' ','_')
    56725679        if G2frame.RMCchoice == 'fullrmc':
     5680            pName = G2frame.GSASprojectfile.split('.')[0] + '-' + generalData['Name']
     5681            pName = pName.replace(' ','_')
     5682            rname = pName+'-fullrmc.py'
     5683            if not os.path.exists(rname):
     5684                G2G.G2MessageBox(G2frame,
     5685                    'The fullrmc script has not been created. Running setup.',
     5686                    'Not setup')
     5687                OnSetupRMC(event)
    56735688            RMCPdict = data['RMC']['fullrmc']
    5674             wx.MessageBox(''' For use of fullrmc, please cite:
    5675       Fullrmc, a Rigid Body Reverse Monte Carlo Modeling Package Enabled with
    5676       Machine Learning and Artificial Intelligence,
     5689            rmcname = pName+'.rmc'
     5690            if os.path.isdir(rmcname) and RMCPdict['ReStart'][0]:
     5691                G2G.G2MessageBox(G2frame,
     5692                    '''You have asked to restart fullrmc but have an existing
     5693                    run as {}. You must manually delete this directory if
     5694                    you wish to restart or change the restart checkbox to
     5695                    continue from the previous results.
     5696                    '''.format(rmcname),'Restart or Continue?')
     5697                # TODO: could do this for the user with:
     5698                #import shutil
     5699                #shutil.rmtree(rmcname)
     5700            G2G.G2MessageBox(G2frame,'''For use of fullrmc, please cite:
     5701
     5702      "Fullrmc, a Rigid Body Reverse Monte Carlo
     5703      Modeling Package Enabled with Machine Learning
     5704      and Artificial Intelligence",
    56775705      B. Aoun, Jour. Comp. Chem. 2016, 37, 1102-1111.
    5678       doi: https://doi.org/10.1002/jcc.24304''',caption='fullrmc',style=wx.ICON_INFORMATION)
    5679             rmcname = pName+'.rmc'           
    5680             while os.path.isdir(rmcname) and RMCPdict['ReStart'][0]:
    5681                 msg = wx.MessageBox(''' fullrmc will fail to restart if %s exists. You must delete %s by hand now.'''%(rmcname,rmcname),
    5682                     caption='fullrmc file error',style=wx.ICON_EXCLAMATION|wx.OK|wx.CANCEL)
    5683                 if msg != wx.OK:
    5684                     return
    5685             if os.path.isfile('pdbparser_0.log'):
    5686                 os.remove('pdbparser_0.log')           
     5706      DOI: https://doi.org/10.1002/jcc.24304''','Please cite fullrmc')
    56875707            ilog = 0
    56885708            while True:
    … …  
    56935713                    break
    56945714                ilog += 1
    5695 # TBD - remove filedialog & use defined run.py file name here
    5696             rname = pName+'-run.py'
    5697             # dlg = wx.FileDialog(G2frame, 'Choose fullrmc python file to execute', G2G.GetImportPath(G2frame),
    5698             #     wildcard='fullrmc python file (*.py)|*.py',style=wx.FD_CHANGE_DIR)
    5699             # try:
    5700             #     if dlg.ShowModal() == wx.ID_OK:
    5701             #         rname = dlg.GetPath()
    5702             #     else:
    5703             #         return
    5704             # finally:
    5705             #     dlg.Destroy()
    57065715            import subprocess as sb
    5707             batch = open('fullrmc.bat','w')
    5708             batch.write('CALL '+sys.exec_prefix+'\\Scripts\\activate\n')
    5709             batch.write(sys.exec_prefix+'\\python.exe '+rname+'\n')
    5710             batch.write('pause')
    5711             batch.close()
    5712             RMCPdict['pid'] = sb.Popen('fullrmc.bat',creationflags=sb.CREATE_NEW_CONSOLE).pid
     5716            if sys.platform.lower().startswith('win'):
     5717                batch = open('fullrmc.bat','w')
     5718                batch.write('CALL '+sys.exec_prefix+'\\Scripts\\activate\n')
     5719                batch.write(sys.exec_prefix+'\\python.exe '+rname+'\n')
     5720                batch.write('pause')
     5721                batch.close()
     5722                sb.Popen('fullrmc.bat',creationflags=sb.CREATE_NEW_CONSOLE)
     5723            else:
     5724                batch = open('fullrmc.sh','w')
     5725                batch.write('#!/bin/bash\n')
     5726                activate = os.path.split(os.environ.get('CONDA_EXE',''))[0] +'/activate'
     5727                batch.write('cd ' + os.path.split(os.path.abspath(rname))[0] + '\n')
     5728                if os.path.exists(activate):
     5729                    batch.write('source ' + activate + ' ' +
     5730                                os.environ['CONDA_DEFAULT_ENV'] +'\n')
     5731                    batch.write('python ' + rname + '\n')
     5732                else:
     5733                    batch.write(sys.exec_prefix+'/python ' + rname + '\n')
     5734                batch.close()
     5735                if sys.platform == "darwin":
     5736                    GSASIIpath.MacRunScript(os.path.abspath('fullrmc.sh'))
     5737                else:
     5738                    # TODO: better to create this in a new terminal on Linux
     5739                    sb.Popen(['/bin/bash','fullrmc.sh'])
    57135740        else:
     5741            pName = generalData['Name'].replace(' ','_')
    57145742            RMCPdict = data['RMC']['RMCProfile']
    57155743            rmcfile = G2fl.find('rmcprofile.exe',GSASIIpath.path2GSAS2)
    … …  
    57685796    def OnViewRMC(event):
    57695797        if G2frame.RMCchoice == 'fullrmc':
     5798               
     5799            #dlg = wx.DirDialog (None, "Choose fullrmc directory", "",
     5800            #        wx.DD_DEFAULT_STYLE | wx.DD_DIR_MUST_EXIST)
    57705801            RMCPdict = data['RMC']['fullrmc']
    57715802            generalData = data['General']
    5772             pName = generalData['Name'].replace(' ','_')
    5773             import psutil
    5774             pid = RMCPdict.get('pid',-1)
    5775             Proc = None
    5776             if pid and psutil.pid_exists(pid):
    5777                 proc = psutil.Process(pid).children()
    5778                 for child in proc:
    5779                     if 'conhost' in child.name():       #probably very Windows specific
    5780                         Proc = child
     5803            pName = G2frame.GSASprojectfile.split('.')[0] + '-' + generalData['Name']
     5804            pName = pName.replace(' ','_')
     5805            engineFilePath = pName+'.rmc'
     5806            if not os.path.exists(engineFilePath):
     5807                dlg = wx.FileDialog(self, 'Open fullrmc directory',
     5808                                        defaultFile='*.rmc',
     5809                                        wildcard='*.rmc')
     5810                try:
     5811                    if dlg.ShowModal() == wx.ID_OK:
     5812                        engineFilePath = dlg.GetPath()
     5813                    else:
     5814                        return
     5815                finally:
     5816                    dlg.Destroy()
     5817                engineFilePath = os.path.splitext(engineFilePath)[0] + '.rmc'
     5818                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
    57815827            from fullrmc import Engine
    57825828            # load
    5783             engineFilePath = pName+'.rmc'
     5829            GSASIIpath.IPyBreak()
    57845830            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)
    57855840            eNames = ['Total',]
     5841            found = False
    57865842            nObs = [0,]
    57875843            try:
    5788                 if engineFilePath is not None:
    5789                     result, mes = ENGINE.is_engine(engineFilePath, mes=True)
    5790                     if result:
    5791                         if Proc is not None:
    5792                             Proc.suspend()
    5793                         ENGINE = ENGINE.load(engineFilePath)
    5794                         found = False
    5795                         for frame in list(ENGINE.frames):
    5796                             ENGINE.set_used_frame(frame)
    5797                             FRitems = ENGINE.constraints
    5798                             for item in FRitems:
    5799                                 sitem = str(type(item))
    5800                                 if 'PairDistribution' in sitem or 'StructureFactor' in sitem or 'PairCorrelation' in sitem:
    5801                                     nameId = 'X'
    5802                                     if 'neutron' in item.weighting:
    5803                                         nameId = 'N'
    5804                                     found = True
    5805                                     Xlab = r'$\mathsf{r,\AA}$'
    5806                                     Ylab = r'$\mathsf{g(r),\AA^{-2}}$'
    5807                                     title = ' g(r)%s for '%nameId
    5808                                     if 'StructureFactor' in sitem:
    5809                                         eNames.append('S(Q)'+nameId)
    5810                                         Xlab = r'$\mathsf{Q,\AA^{-1}}$'
    5811                                         Ylab = 'S(Q)'
    5812                                         title = ' S(Q)%s for '%nameId
     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])
    58135897                                    else:
    5814                                         eNames.append('g(r)'+nameId)
    5815                                     dataDict= item.get_constraints_properties(frame)
    5816                                     X = dataDict['frames-experimental_x'][0]
    5817                                     Y = dataDict['frames-experimental_y'][0]
    5818                                     nObs.append(X.shape[0])
    5819                                     rdfDict = item.get_constraint_value()
    5820                                     if 'total' not in rdfDict:
    5821                                         print('No data yet - wait for a save')
    5822                                         ENGINE.close()
    5823                                         if Proc is not None:
    5824                                             Proc.resume()
    5825                                         return
    5826                                     Z = rdfDict['total']
    5827                                     XY = [[X,Z],[X,Y]]
    5828                                     Names = ['Calc','Obs']
    5829                                     G2plt.PlotXY(G2frame,XY,labelX=Xlab,
    5830                                         labelY=Ylab,newPlot=True,Title=title+pName,
    5831                                         lines=True,names=Names)
    5832                                     PXY = []
    5833                                     PXYT = []
    5834                                     Names = []
    5835                                     NamesT = []
    5836                                     for item in rdfDict:
    5837                                         if 'va' in item:
    5838                                             continue
    5839                                         if 'rdf' not in item and 'g(r)' in title:
    5840                                             Ylab = r'$\mathsf{G(r),\AA^{-1}}$'
    5841                                             PXYT.append([X,1.+rdfDict[item]/X])
    5842                                         else:
    5843                                             PXYT.append([X,rdfDict[item]])
    5844                                         NamesT.append(item)
    5845                                         if 'total' in item:
    5846                                             if 'rdf' not in item and 'g(r)' in title:
    5847                                                 Ylab = r'$\mathsf{G(r),\AA^{-1}}$'
    5848                                                 PXY.append([X,1.+rdfDict[item]/X])
    5849                                             else:
    5850                                                 PXY.append([X,rdfDict[item]])
    5851                                             Names.append(item)
    5852                                     G2plt.PlotXY(G2frame,PXYT,labelX=Xlab,
    5853                                         labelY=Ylab,newPlot=True,Title=' All partials of '+title+pName,
    5854                                         lines=True,names=NamesT)
    5855                                     G2plt.PlotXY(G2frame,PXY,labelX=Xlab,
    5856                                         labelY=Ylab,newPlot=True,Title=' Total partials of '+title+pName,
    5857                                         lines=True,names=Names)
    5858                                    
    5859                                 else:
    5860                                     found = True
    5861                                     if 'BondConstraint' in sitem:
    5862                                         try:
    5863                                             bonds = item.get_constraint_value()['bondsLength']
    5864                                         except TypeError:
    5865                                             break
    5866                                         bondList = item.bondsList[:2]
    5867                                         atoms = ENGINE.get_original_data("allElements",frame)
    5868                                         bondNames = ['%s-%s'%(atoms[bondList[0][iat]],atoms[bondList[1][iat]]) for iat in range(bonds.shape[0])]
    5869                                         bondSet = list(set(bondNames))
    5870                                         Bonds = list(zip(bondNames,bonds))
    5871                                         for Bname in bondSet:
    5872                                             bondLens = [bond[1] for bond in Bonds if bond[0]==Bname]
    5873                                             G2plt.PlotBarGraph(G2frame,bondLens,Xname=r'%s $Bond,\ \AA$'%Bname,Title='%s Bond lengths for %s'%(Bname,pName),
    5874                                                 PlotName='%s Bonds for %s'%(Bname,pName),maxBins=20)
    5875                                             print(' %d %s bonds found'%(len(bondLens),Bname))
    5876                                     elif 'BondsAngleConstraint' in sitem:
    5877                                         angles = 180.*item.get_constraint_value()['angles']/np.pi
    5878                                         angleList = item.anglesList[:3]
    5879                                         atoms = ENGINE.get_original_data("allElements",frame)
    5880                                         angleNames = ['%s-%s-%s'%(atoms[angleList[1][iat]],atoms[angleList[0][iat]],atoms[angleList[2][iat]]) for iat in range(angles.shape[0])]
    5881                                         angleSet = list(set(angleNames))
    5882                                         Angles = list(zip(angleNames,angles))
    5883                                         for Aname in angleSet:
    5884                                             bondAngs = [angle[1] for angle in Angles if angle[0]==Aname]
    5885                                             G2plt.PlotBarGraph(G2frame,bondAngs,Xname=r'%s Angle, deg'%Aname,Title='%s Bond angles for %s'%(Aname,pName),
    5886                                                 PlotName='%s Angles for %s'%(Aname,pName),maxBins=20)
    5887                                             print(' %d %s angles found'%(len(bondAngs),Aname))
    5888                                     elif 'DihedralAngleConstraint' in sitem:
    5889                                         impangles = 180.*item.get_constraint_value()['angles']/np.pi
    5890                                         impangleList = item.anglesList[:4]
    5891                                         print(' Dihedral angle chi^2 =  %2f'%item.standardError)
    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                                         impangleSet = list(set(impangleNames))
    5896                                         impAngles = list(zip(impangleNames,impangles))
    5897                                         for Aname in impangleSet:
    5898                                             impAngs = [angle[1] for angle in impAngles if angle[0]==Aname]
    5899                                             G2plt.PlotBarGraph(G2frame,impAngs,Xname=r'%s $Dihedral Angle, deg$'%Aname,Title='%s Dihedral angles for %s'%(Aname,pName),
    5900                                                 PlotName='%s Dihedral Angles for %s'%(Aname,pName),maxBins=20)
    5901                                     elif 'AtomicCoordinationNumber' in sitem or 'InterMolecular' in sitem:
    5902                                         print(sitem+' not plotted')
    5903                                     else:
    5904                                         print(sitem)
    5905                                         item.plot(show=True)
    5906                                         pass
    5907                             nObs[0] = np.sum(nObs[1:])
    5908                         if not found:
    5909                             print(' No saved information yet, wait until fullrmc does a Save')
    5910                             eNames = []
    5911                         ENGINE.close()      #return lock on ENGINE repository & close it
    5912                         if Proc is not None:
    5913                             Proc.resume()
     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()
    59145962            except AssertionError:
    59155963                 print("Can't open fullrmc engine while running")
    5916             loglines = []
    5917             ilog = 0
    5918             while True:
    5919                 fname = '%s_%d.log'%(pName,ilog)
    5920                 try:
    5921                     logfile = open(fname,'r')
    5922                     loglines += logfile.readlines()
    5923                     logfile.close()
    5924                 except FileNotFoundError:
    5925                     break
    5926                 ilog += 1
    5927             if not len(loglines):
    5928                 print('no log file found')
    5929                 return
    5930             start = 0
    5931             while True:
    5932                 if start == len(loglines):
    5933                     print('No log info to plot')
    5934                     return
    5935                 line = loglines[start]               
    5936                 if 'Err:' in line:
    5937                     break
    5938                 else:
    5939                     start += 1
    5940             Gen = []
    5941             Err = []
    5942             start -= 1
    5943             while True:
    5944                 start += 1
    5945                 try:
    5946                     line = loglines[start]
    5947                 except:
    5948                     break
    5949                 if 'Err' not in line:
    5950                     continue
    5951                 items = line.split(' - ')
    5952                 try:    # could be a trashed line at end
    5953                     errStr = items[5][:-1].split('Err:')[1]
    5954                     Err.append([float(val) for val in errStr.split(',')])
    5955                 except ValueError:
    5956                     break
    5957                 Gen.append(int(items[1].split('Gen:')[1]))
    5958            
    5959             Gen = np.array(Gen)
    5960             Err = np.array(Err)
    5961             nObs = np.array(nObs)
    5962             if np.any(nObs):
    5963                 Err /= nObs[:Err.shape[1]]
    5964                 ptstr1 = ''
    5965                 ptstr2 = ''
    5966                 for it,item in enumerate(eNames):
    5967                     ptstr1 += ' %s obs: %d'%(item,nObs[it])
    5968                     ptstr2 += ' %s reduced chi^2: %.5f'%(item,Err[-1][it])
    5969                 print(ptstr1)
    5970                 print(ptstr2)
    5971                 Err = np.log10(Err)
    5972                 XY = [[Gen,Erri] for Erri in Err.T]
    5973                 G2plt.PlotXY(G2frame,XY,labelX='no. generated',
    5974                     labelY=r'$log_{10}$ (reduced $\mathsf{\chi^2})$',newPlot=True,Title='fullrmc residuals for '+pName,
    5975                     lines=True,names=eNames)
     5964            # loglines = []
     5965            # ilog = 0
     5966            # while True:
     5967            #     fname = '%s_%d.log'%(pName,ilog)
     5968            #     try:
     5969            #         logfile = open(fname,'r')
     5970            #         loglines += logfile.readlines()
     5971            #         logfile.close()
     5972            #     except FileNotFoundError:
     5973            #         break
     5974            #     ilog += 1
     5975            # if not len(loglines):
     5976            #     print('no log file found')
     5977            #     return
     5978            # start = 0
     5979            # while True:
     5980            #     if start == len(loglines):
     5981            #         print('No log info to plot')
     5982            #         return
     5983            #     line = loglines[start]               
     5984            #     if 'Err:' in line:
     5985            #         break
     5986            #     else:
     5987            #         start += 1
     5988            # Gen = []
     5989            # Err = []
     5990            # start -= 1
     5991            # while True:
     5992            #     start += 1
     5993            #     try:
     5994            #         line = loglines[start]
     5995            #     except:
     5996            #         break
     5997            #     if 'Err' not in line:
     5998            #         continue
     5999            #     items = line.split(' - ')
     6000            #     try:    # could be a trashed line at end
     6001            #         errStr = items[5][:-1].split('Err:')[1]
     6002            #         Err.append([float(val) for val in errStr.split(',')])
     6003            #     except ValueError:
     6004            #         break
     6005            #     Gen.append(int(items[1].split('Gen:')[1]))
     6006           
     6007            # Gen = np.array(Gen)
     6008            # Err = np.array(Err)
     6009            # nObs = np.array(nObs)
     6010            # if np.any(nObs):
     6011            #     Err /= nObs[:Err.shape[1]]
     6012            #     ptstr1 = ''
     6013            #     ptstr2 = ''
     6014            #     for it,item in enumerate(eNames):
     6015            #         ptstr1 += ' %s obs: %d'%(item,nObs[it])
     6016            #         ptstr2 += ' %s reduced chi^2: %.5f'%(item,Err[-1][it])
     6017            #     print(ptstr1)
     6018            #     print(ptstr2)
     6019            #     Err = np.log10(Err)
     6020            #     XY = [[Gen,Erri] for Erri in Err.T]
     6021            #     G2plt.PlotXY(G2frame,XY,labelX='no. generated',
     6022            #         labelY=r'$log_{10}$ (reduced $\mathsf{\chi^2})$',newPlot=True,Title='fullrmc residuals for '+pName,
     6023            #         lines=True,names=eNames)
    59766024                     
    59776025        else:
    … …  
    1023310281                dups = []
    1023410282                assigned = []
     10283                atomsLabels = [a[0] for a in data['Atoms']]
    1023510284                for r in range(tbl.GetRowsCount()):
    1023610285                    sel = tbl.GetValue(r,4).strip()
    1023710286                    if sel == 'Create new': continue # ignore positions of new atoms
    10238                     if sel not in labelsChoices: continue
    10239                     atmNum = labelsChoices.index(sel)-1
     10287                    if sel not in atomsLabels: continue
     10288                    atmNum = atomsLabels.index(sel)
    1024010289                    if atmNum < 0: continue
    1024110290                    if atmNum in assigned:
    … …  
    1041810467                    Indx[ang.GetId()] = [t,torSlide]
    1041910468                    TorSizer.Add(ang,0,WACV)                           
    10420                 mainSizer.Add(TorSizer,1,wx.EXPAND|wx.RIGHT)
     10469                mainSizer.Add(TorSizer,0,wx.EXPAND|wx.RIGHT)
    1042110470            else:
    1042210471                mainSizer.Add(wx.StaticText(RigidBodies,wx.ID_ANY,'No side chain torsions'),0,WACV)

  • trunk/GSASIIpwd.py

    r4483 r4518  
    2020import os.path
    2121import subprocess as subp
     22import datetime as dt
    2223import copy
    2324
    … …  
    3334
    3435import GSASIIpath
     36filversion = "$Revision$"
    3537GSASIIpath.SetVersionNumber("$Revision$")
    3638import GSASIIlattice as G2lat
    … …  
    25112513    return fname
    25122514
    2513 def FindBonds(Phase,RMCPdict):
    2514     generalData = Phase['General']
    2515     cx,ct,cs,cia = generalData['AtomPtrs']
    2516     atomData = Phase['Atoms']
    2517     Res = 'RMC'
    2518     if 'macro' in generalData['Type']:
    2519         Res = atomData[0][ct-3]
    2520     AtDict = {atom[ct-1]:atom[ct] for atom in atomData}
    2521     Pairs = RMCPdict['Pairs']   #dict!
    2522     BondList = []
    2523     notNames = []
    2524     for FrstName in AtDict:
    2525         nbrs = G2mth.FindAllNeighbors(Phase,FrstName,list(AtDict.keys()),notName=notNames,Short=True)[0]
    2526         Atyp1 = AtDict[FrstName]
    2527         if 'Va' in Atyp1:
     2515# def FindBonds(Phase,RMCPdict):
     2516#     generalData = Phase['General']
     2517#     cx,ct,cs,cia = generalData['AtomPtrs']
     2518#     atomData = Phase['Atoms']
     2519#     Res = 'RMC'
     2520#     if 'macro' in generalData['Type']:
     2521#         Res = atomData[0][ct-3]
     2522#     AtDict = {atom[ct-1]:atom[ct] for atom in atomData}
     2523#     Pairs = RMCPdict['Pairs']   #dict!
     2524#     BondList = []
     2525#     notNames = []
     2526#     for FrstName in AtDict:
     2527#         nbrs = G2mth.FindAllNeighbors(Phase,FrstName,list(AtDict.keys()),notName=notNames,Short=True)[0]
     2528#         Atyp1 = AtDict[FrstName]
     2529#         if 'Va' in Atyp1:
     2530#             continue
     2531#         for nbr in nbrs:
     2532#             Atyp2 = AtDict[nbr[0]]
     2533#             if 'Va' in Atyp2:
     2534#                 continue
     2535#             try:
     2536#                 bndData = Pairs[' %s-%s'%(Atyp1,Atyp2)][1:]
     2537#             except KeyError:
     2538#                 bndData = Pairs[' %s-%s'%(Atyp2,Atyp1)][1:]
     2539#             if any(bndData):
     2540#                 if bndData[0] <= nbr[1] <= bndData[1]:
     2541#                     bondStr = str((FrstName,nbr[0])+tuple(bndData))+',\n'
     2542#                     revbondStr = str((nbr[0],FrstName)+tuple(bndData))+',\n'
     2543#                     if bondStr not in BondList and revbondStr not in BondList:
     2544#                         BondList.append(bondStr)
     2545#         notNames.append(FrstName)
     2546#     return Res,BondList
     2547
     2548# def FindAngles(Phase,RMCPdict):
     2549#     generalData = Phase['General']
     2550#     Cell = generalData['Cell'][1:7]
     2551#     Amat = G2lat.cell2AB(Cell)[0]
     2552#     cx,ct,cs,cia = generalData['AtomPtrs']
     2553#     atomData = Phase['Atoms']
     2554#     AtLookup = G2mth.FillAtomLookUp(atomData,cia+8)
     2555#     AtDict = {atom[ct-1]:atom[ct] for atom in atomData}
     2556#     Angles = RMCPdict['Angles']
     2557#     AngDict = {'%s-%s-%s'%(angle[0],angle[1],angle[2]):angle[3:] for angle in Angles}
     2558#     AngleList = []
     2559#     for MidName in AtDict:
     2560#         nbrs,nbrIds = G2mth.FindAllNeighbors(Phase,MidName,list(AtDict.keys()),Short=True)
     2561#         if len(nbrs) < 2: #need 2 neighbors to make an angle
     2562#             continue
     2563#         Atyp2 = AtDict[MidName]
     2564#         for i,nbr1 in enumerate(nbrs):
     2565#             Atyp1 = AtDict[nbr1[0]]
     2566#             for j,nbr3 in enumerate(nbrs[i+1:]):
     2567#                 Atyp3 = AtDict[nbr3[0]]
     2568#                 IdList = [nbrIds[1][i],nbrIds[0],nbrIds[1][i+j+1]]
     2569#                 try:
     2570#                     angData = AngDict['%s-%s-%s'%(Atyp1,Atyp2,Atyp3)]
     2571#                 except KeyError:
     2572#                     try:
     2573#                         angData = AngDict['%s-%s-%s'%(Atyp3,Atyp2,Atyp1)]
     2574#                     except KeyError:
     2575#                         continue
     2576#                 XYZ = np.array(G2mth.GetAtomItemsById(atomData,AtLookup,IdList,cx,numItems=3))
     2577#                 calAngle = G2mth.getRestAngle(XYZ,Amat)
     2578#                 if angData[0] <= calAngle <= angData[1]:
     2579#                     angStr = str((MidName,nbr1[0],nbr3[0])+tuple(angData))+',\n'
     2580#                     revangStr = str((MidName,nbr3[0],nbr1[0])+tuple(angData))+',\n'
     2581#                     if angStr not in AngleList and revangStr not in AngleList:
     2582#                         AngleList.append(angStr)
     2583#     return AngleList
     2584
     2585# def GetSqConvolution(XY,d):
     2586
     2587#     n = XY.shape[1]
     2588#     snew = np.zeros(n)
     2589#     dq = np.zeros(n)
     2590#     sold = XY[1]
     2591#     q = XY[0]
     2592#     dq[1:] = np.diff(q)
     2593#     dq[0] = dq[1]
     2594   
     2595#     for j in range(n):
     2596#         for i in range(n):
     2597#             b = abs(q[i]-q[j])
     2598#             t = q[i]+q[j]
     2599#             if j == i:
     2600#                 snew[j] += q[i]*sold[i]*(d-np.sin(t*d)/t)*dq[i]
     2601#             else:
     2602#                 snew[j] += q[i]*sold[i]*(np.sin(b*d)/b-np.sin(t*d)/t)*dq[i]
     2603#         snew[j] /= np.pi*q[j]
     2604   
     2605#     snew[0] = snew[1]
     2606#     return snew
     2607
     2608# def GetMaxSphere(pdbName):
     2609#     try:
     2610#         pFil = open(pdbName,'r')
     2611#     except FileNotFoundError:
     2612#         return None
     2613#     while True:
     2614#         line = pFil.readline()
     2615#         if 'Boundary' in line:
     2616#             line = line.split()[3:]
     2617#             G = np.array([float(item) for item in line])
     2618#             G = np.reshape(G,(3,3))**2
     2619#             G = nl.inv(G)
     2620#             pFil.close()
     2621#             break
     2622#     dspaces = [0.5/np.sqrt(G2lat.calc_rDsq2(H,G)) for H in np.eye(3)]
     2623#     return min(dspaces)
     2624   
     2625def MakefullrmcRun(pName,Phase,RMCPdict):
     2626    BondList = {}
     2627    for k in RMCPdict['Pairs']:
     2628        if RMCPdict['Pairs'][k][1]+RMCPdict['Pairs'][k][2]>0:
     2629            BondList[k] = (RMCPdict['Pairs'][k][1],RMCPdict['Pairs'][k][2])
     2630    AngleList = []
     2631    for angle in RMCPdict['Angles']:
     2632        if angle[3] == angle[4] or angle[5] >= angle[6] or angle[6] <= 0:
    25282633            continue
    2529         for nbr in nbrs:
    2530             Atyp2 = AtDict[nbr[0]]
    2531             if 'Va' in Atyp2:
    2532                 continue
    2533             try:
    2534                 bndData = Pairs[' %s-%s'%(Atyp1,Atyp2)][1:]
    2535             except KeyError:
    2536                 bndData = Pairs[' %s-%s'%(Atyp2,Atyp1)][1:]
    2537             if any(bndData):
    2538                 if bndData[0] <= nbr[1] <= bndData[1]:
    2539                     bondStr = str((FrstName,nbr[0])+tuple(bndData))+',\n'
    2540                     revbondStr = str((nbr[0],FrstName)+tuple(bndData))+',\n'
    2541                     if bondStr not in BondList and revbondStr not in BondList:
    2542                         BondList.append(bondStr)
    2543         notNames.append(FrstName)
    2544     return Res,BondList
    2545 
    2546 def FindAngles(Phase,RMCPdict):
    2547     generalData = Phase['General']
    2548     Cell = generalData['Cell'][1:7]
    2549     Amat = G2lat.cell2AB(Cell)[0]
    2550     cx,ct,cs,cia = generalData['AtomPtrs']
    2551     atomData = Phase['Atoms']
    2552     AtLookup = G2mth.FillAtomLookUp(atomData,cia+8)
    2553     AtDict = {atom[ct-1]:atom[ct] for atom in atomData}
    2554     Angles = RMCPdict['Angles']
    2555     AngDict = {'%s-%s-%s'%(angle[0],angle[1],angle[2]):angle[3:] for angle in Angles}
    2556     AngleList = []
    2557     for MidName in AtDict:
    2558         nbrs,nbrIds = G2mth.FindAllNeighbors(Phase,MidName,list(AtDict.keys()),Short=True)
    2559         if len(nbrs) < 2: #need 2 neighbors to make an angle
    2560             continue
    2561         Atyp2 = AtDict[MidName]
    2562         for i,nbr1 in enumerate(nbrs):
    2563             Atyp1 = AtDict[nbr1[0]]
    2564             for j,nbr3 in enumerate(nbrs[i+1:]):
    2565                 Atyp3 = AtDict[nbr3[0]]
    2566                 IdList = [nbrIds[1][i],nbrIds[0],nbrIds[1][i+j+1]]
    2567                 try:
    2568                     angData = AngDict['%s-%s-%s'%(Atyp1,Atyp2,Atyp3)]
    2569                 except KeyError:
    2570                     try:
    2571                         angData = AngDict['%s-%s-%s'%(Atyp3,Atyp2,Atyp1)]
    2572                     except KeyError:
    2573                         continue
    2574                 XYZ = np.array(G2mth.GetAtomItemsById(atomData,AtLookup,IdList,cx,numItems=3))
    2575                 calAngle = G2mth.getRestAngle(XYZ,Amat)
    2576                 if angData[0] <= calAngle <= angData[1]:
    2577                     angStr = str((MidName,nbr1[0],nbr3[0])+tuple(angData))+',\n'
    2578                     revangStr = str((MidName,nbr3[0],nbr1[0])+tuple(angData))+',\n'
    2579                     if angStr not in AngleList and revangStr not in AngleList:
    2580                         AngleList.append(angStr)
    2581     return AngleList
    2582 
    2583 def GetSqConvolution(XY,d):
    2584 
    2585     n = XY.shape[1]
    2586     snew = np.zeros(n)
    2587     dq = np.zeros(n)
    2588     sold = XY[1]
    2589     q = XY[0]
    2590     dq[1:] = np.diff(q)
    2591     dq[0] = dq[1]
    2592    
    2593     for j in range(n):
    2594         for i in range(n):
    2595             b = abs(q[i]-q[j])
    2596             t = q[i]+q[j]
    2597             if j == i:
    2598                 snew[j] += q[i]*sold[i]*(d-np.sin(t*d)/t)*dq[i]
    2599             else:
    2600                 snew[j] += q[i]*sold[i]*(np.sin(b*d)/b-np.sin(t*d)/t)*dq[i]
    2601         snew[j] /= np.pi*q[j]
    2602    
    2603     snew[0] = snew[1]
    2604     return snew
    2605 
    2606 def GetMaxSphere(pdbName):
    2607     try:
    2608         pFil = open(pdbName,'r')
    2609     except FileNotFoundError:
    2610         return None
    2611     while True:
    2612         line = pFil.readline()
    2613         if 'Boundary' in line:
    2614             line = line.split()[3:]
    2615             G = np.array([float(item) for item in line])
    2616             G = np.reshape(G,(3,3))**2
    2617             G = nl.inv(G)
    2618             pFil.close()
    2619             break
    2620     dspaces = [0.5/np.sqrt(G2lat.calc_rDsq2(H,G)) for H in np.eye(3)]
    2621     return min(dspaces)
    2622    
    2623 def MakefullrmcRun(pName,Phase,RMCPdict):
    2624     Res,BondList = FindBonds(Phase,RMCPdict)
    2625     AngleList = FindAngles(Phase,RMCPdict)
     2634        AngleList.append(angle)
    26262635    rmin = RMCPdict['min Contact']
    2627     rmax = GetMaxSphere(RMCPdict['atomPDB'])
    2628     if rmax is None:
    2629         return None
    2630     rname = pName+'-run.py'
     2636    cell = Phase['General']['Cell'][1:7]
     2637    bigcell = np.array(cell)*np.array(RMCPdict['SuperCell']+[1,1,1])
     2638    bigG = G2lat.cell2Gmat(bigcell)[0]
     2639    rmax = min([0.5/np.sqrt(G2lat.calc_rDsq2(H,bigG)) for H in np.eye(3)])
     2640    SymOpList = G2spc.AllOps(Phase['General']['SGData'])[0]
     2641    cx,ct,cs,cia = Phase['General']['AtomPtrs']
     2642    atomsList = []
     2643    for atom in Phase['Atoms']:
     2644        el = ''.join([i for i in atom[ct] if i.isalpha()])
     2645        atomsList.append([el] + atom[cx:cx+4])
     2646    rname = pName+'-fullrmc.py'
    26312647    restart = '%s_restart.pdb'%pName
    26322648    Files = RMCPdict['files']
    … …  
    26342650    rundata = ''
    26352651    rundata += '#### fullrmc %s file; edit by hand if you so choose #####\n'%rname
     2652    rundata += '# created in '+__file__+" v"+filversion.split()[1]
     2653    rundata += dt.datetime.strftime(dt.datetime.now()," at %Y-%m-%dT%H:%M\n")
    26362654    rundata += '''
    26372655# fullrmc imports (all that are potentially useful)
     2656#import matplotlib.pyplot as plt
    26382657import numpy as np
    26392658import time
    2640 from fullrmc.sincConvolution import sincConvolution
    2641 from fullrmc.Globals import LOGGER
     2659from fullrmc.Core import Collection
    26422660from fullrmc.Engine import Engine
    2643 from fullrmc.Constraints.PairDistributionConstraints import PairDistributionConstraint
    2644 from fullrmc.Constraints.StructureFactorConstraints import ReducedStructureFactorConstraint
    2645 from fullrmc.Constraints.DistanceConstraints import InterMolecularDistanceConstraint
     2661import fullrmc.Constraints.PairDistributionConstraints as fPDF
     2662from fullrmc.Constraints.StructureFactorConstraints import ReducedStructureFactorConstraint, StructureFactorConstraint
     2663from fullrmc.Constraints.DistanceConstraints import DistanceConstraint
    26462664from fullrmc.Constraints.BondConstraints import BondConstraint
    26472665from fullrmc.Constraints.AngleConstraints import BondsAngleConstraint
    26482666from fullrmc.Constraints.DihedralAngleConstraints import DihedralAngleConstraint
    26492667from fullrmc.Generators.Swaps import SwapPositionsGenerator
    2650 from fullrmc.debugStuff import *
    2651 InvokeDebugOpts()
     2668
     2669### When True, erases an existing enging to provide a fresh start
     2670FRESH_START = {}
    26522671time0 = time.time()
    2653 SwapGen = {}
    2654 # engine setup\n'''
    2655 #Unused imports
    2656 # from fullrmc.Constraints.PairCorrelationConstraints import PairCorrelationConstraint
    2657 # from fullrmc.Core.MoveGenerator import MoveGeneratorCollector
    2658 # from fullrmc.Core.GroupSelector import RecursiveGroupSelector
    2659 # from fullrmc.Selectors.RandomSelectors import RandomSelector
    2660 # from fullrmc.Selectors.OrderedSelectors import DefinedOrderSelector
    2661 # from fullrmc.Generators.Translations import TranslationGenerator, TranslationAlongSymmetryAxisGenerator
    2662 # from fullrmc.Generators.Agitations import DistanceAgitationGenerator, AngleAgitationGenerator
    2663 # from fullrmc.Generators.Rotations import RotationGenerator, RotationAboutAxisGenerator
    2664 # from fullrmc.Core.Collection import get_principal_axis
    2665 #End unused imports
    2666     rundata += 'LOGGER.set_log_file_basename("%s")\n'%pName
     2672'''.format(RMCPdict['ReStart'][0])
     2673   
     2674    rundata += '# setup structure\n'
     2675    rundata += 'cell = ' + str(cell) + '\n'
     2676    rundata += "SymOpList = "+str([i.lower() for i in SymOpList]) + '\n'
     2677    rundata += 'atomList = ' + str(atomsList).replace('],','],\n  ') + '\n'
     2678    rundata += 'supercell = ' + str(RMCPdict['SuperCell']) + '\n'
     2679
     2680    rundata += '\n# initialize engine\n'
    26672681    rundata += 'engineFileName = "%s.rmc"\n'%pName
    2668     rundata += 'ENGINE = Engine(path=None)\n'
    2669     rundata += 'if not ENGINE.is_engine(engineFileName):\n'
    2670     rundata += '# create engine & set atomic (pdb) model\n'
    2671     rundata += '    ENGINE = Engine(path=engineFileName)\n'
    2672     rundata += '    ENGINE.set_pdb("%s")\n'%RMCPdict['atomPDB']
    2673     rundata += '# create experimental constraints must restart to change these\n'
     2682
     2683    rundata += '''\n# check Engine exists if so (and not FRESH_START) load it
     2684# otherwise build it
     2685ENGINE = Engine(path=None)
     2686if not ENGINE.is_engine(engineFileName) or FRESH_START:
     2687    ## create structure
     2688    ENGINE = Engine(path=engineFileName, freshStart=True)
     2689    ENGINE.build_crystal_set_pdb(symOps     = SymOpList,
     2690                                 atoms      = atomList,
     2691                                 unitcellBC = cell,
     2692                                 supercell  = supercell)
     2693    rmax = min( [ENGINE.boundaryConditions.get_a(), ENGINE.boundaryConditions.get_b(), ENGINE.boundaryConditions.get_c()] ) /2.
     2694'''   
     2695    import atmdata
     2696    rundata += '# conversion factors (may be needed)\n'
     2697    rundata += '    sumCiBi2 = 0.\n'
     2698    for elem in Phase['General']['AtomTypes']:
     2699        rundata += '    Ci = ENGINE.numberOfAtomsPerElement["{}"]/len(ENGINE.allElements)\n'.format(elem)
     2700        rundata += '    sumCiBi2 += (Ci*{})**2\n'.format(atmdata.AtmBlens[elem+'_']['SL'][0])
     2701    rundata += '    rho0 = len(ENGINE.allNames)/ENGINE.volume\n'
     2702    # settings that require a new Engine
    26742703    for File in Files:
    26752704        filDat = RMCPdict['files'][File]
    2676         if filDat[0] != 'Select':
    2677             sfwt = 'neutronCohb'
    2678             if 'Xray' in File:
    2679                 sfwt = 'atomicNumber'
    2680             if 'G(r)' in File:
    2681                 rundata += '    RGR = np.loadtxt("%s").T\n'%filDat[0]
    2682                 if filDat[3]:
    2683                     rundata += '    RGR[1] *= RGR[0]\n'
    2684                 rundata += '    GofR = PairDistributionConstraint(experimentalData=RGR.T, weighting="%s")\n'%sfwt
    2685                 rundata += '    ENGINE.add_constraints([GofR])\n'
    2686                 wtDict['Pair-'+sfwt] = filDat[1]
     2705        if not os.path.exists(filDat[0]): continue
     2706        sfwt = 'neutronCohb'
     2707        if 'Xray' in File:
     2708            sfwt = 'atomicNumber'
     2709        if 'G(r)' in File:
     2710            rundata += '    GR = np.loadtxt("%s").T\n'%filDat[0]
     2711            if filDat[3] == 0:
     2712                rundata += '''    # read and xform G(r) as defined in RMCProfile
     2713# see eq. 44 in Keen, J. Appl. Cryst. (2001) 34 172-177\n'''
     2714                rundata += '    GR[1] *= 4 * np.pi * GR[0] * rho0 / sumCiBi2\n'
     2715                rundata += '    GofR = fPDF.PairDistributionConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
     2716            elif filDat[3] == 1:
     2717                rundata += '    # This is G(r) as defined in PDFFIT\n'
     2718                rundata += '    GofR = fPDF.PairDistributionConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
     2719            elif filDat[3] == 2:
     2720                rundata += '    # This is g(r)\n'
     2721                rundata += '    GofR = fPDF.PairCorrelationConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
    26872722            else:
    2688                 rundata += '    SOQ = np.loadtxt("%s").T\n'%filDat[0]
    2689                 if filDat[3]:
    2690                     rundata += '    SOQ[1] = sincConvolution(SOQ,%.3f)\n'%rmax
    2691                 rundata += '    FofQ = ReducedStructureFactorConstraint(experimentalData=SOQ.T, weighting="%s")\n'%sfwt
    2692                 rundata += '    ENGINE.add_constraints([FofQ])\n'
    2693                 wtDict['Struct-'+sfwt] = filDat[1]
    2694     rundata += '    ENGINE.add_constraints(InterMolecularDistanceConstraint())\n'
    2695     if RMCPdict['byMolec']:
    2696         if len(BondList):
    2697             rundata += '    B_CONSTRAINT   = BondConstraint()\n'
    2698             rundata += '    ENGINE.add_constraints(B_CONSTRAINT)\n'
    2699         if len(AngleList):
    2700             rundata += '    A_CONSTRAINT   = BondsAngleConstraint()\n'
    2701             rundata += '    ENGINE.add_constraints(A_CONSTRAINT)\n'
    2702         if len(RMCPdict['Torsions']):
    2703             rundata += '    T_CONSTRAINT   = DihedralAngleConstraint()\n'
    2704             rundata += '    ENGINE.add_constraints(T_CONSTRAINT)\n'
    2705     rundata += '    ENGINE.save()\n'
    2706     rundata += 'else:\n'
    2707     rundata += '    ENGINE = ENGINE.load(path=engineFileName)\n'
    2708     rundata += '#fill & change constraints - can be done without restart\n'
    2709     rundata += 'wtDict = %s\n'%str(wtDict)
    2710     rundata += 'Constraints = ENGINE.constraints\n'
    2711     rundata += 'for constraint in Constraints:\n'
    2712     rundata += '    strcons = str(type(constraint))\n'
    2713     rundata += '    if "InterMolecular" in strcons:\n'
    2714     rundata += '        constraint.set_default_distance(%f)\n'%RMCPdict['min Contact']
    2715     rundata += '    elif "PairDistribution" in strcons:\n'
    2716     rundata += '        constraint.set_variance_squared(wtDict["Pair-"+constraint.weighting])\n'
    2717     rundata += '        constraint.set_limits((None,%.3f))\n'%(rmax)
    2718 #    rundata += '        constraint.set_limits((%.3f,%.3f))\n'%(rmin,rmax)
     2723                raise ValueError('Invalid G(r) type: '+str(filDat[3]))
     2724            rundata += '    ENGINE.add_constraints([GofR])\n'
     2725            rundata += '    GofR.set_limits((None, rmax))\n'
     2726            wtDict['Pair-'+sfwt] = filDat[1]
     2727        elif 'F(Q)' in File:
     2728            rundata += '    SOQ = np.loadtxt("%s").T\n'%filDat[0]
     2729            if filDat[3] == 0:
     2730                rundata += '    # Read & xform F(Q) as defined in RMCProfile\n'
     2731                rundata += '    SOQ[1] *= 1 / sumCiBi2\n'
     2732                rundata += '    SOQ[1] += 1\n'
     2733            elif filDat[3] == 1:
     2734                rundata += '    # This is S(Q) as defined in PDFFIT\n'
     2735            if filDat[4]:
     2736                rundata += '    SOQ[1] = Collection.sinc_convolution(q=SOQ[0],sq=SOQ[1],rmax={:.3f})\n'.format(rmax)
     2737            rundata += '    SofQ = ReducedStructureFactorConstraint(experimentalData=SOQ.T, weighting="%s")\n'%sfwt
     2738            rundata += '    ENGINE.add_constraints([SofQ])\n'
     2739        else:
     2740            print('What is this?')
     2741        wtDict['Struct-'+sfwt] = filDat[1]
     2742    rundata += '    ENGINE.add_constraints(DistanceConstraint(defaultLowerDistance={}))\n'.format(RMCPdict['min Contact'])
     2743    rundata += '''    B_CONSTRAINT   = BondConstraint()
     2744    ENGINE.add_constraints(B_CONSTRAINT)
     2745    B_CONSTRAINT.create_supercell_bonds(bondsDefinition=[
     2746'''
     2747    for pair in BondList:
     2748        e1,e2 = pair.split('-')
     2749        rundata += '            ("element","{}","{}",{},{}),\n'.format(
     2750                                    e1.strip(),e2.strip(),*BondList[pair])
     2751    rundata += '''             ])
     2752    ENGINE.save()
     2753else:
     2754    ENGINE = ENGINE.load(path=engineFileName)
     2755'''               
     2756#    if RMCPdict['Swaps']:
     2757#        rundata += '#set up for site swaps\n'
     2758#        rundata += 'aN = ENGINE.allNames\n'
     2759#        rundata += 'SwapGen = {}\n'
     2760#        for swap in RMCPdict['Swaps']:
     2761#            rundata += 'SwapA = [[idx] for idx in range(len(aN)) if aN[idx]=="%s"]\n'%swap[0]
     2762#            rundata += 'SwapB = [[idx] for idx in range(len(aN)) if aN[idx]=="%s"]\n'%swap[1]
     2763#            rundata += 'SwapGen["%s-%s"] = [SwapPositionsGenerator(swapList=SwapA),SwapPositionsGenerator(swapList=SwapB),%.2f]\n'%(swap[0],swap[1],swap[2])
     2764    rundata += '# setup/change constraints - can be done without restart\n'   
     2765    rundata += 'for c in ENGINE.constraints:  # loop over predefined constraints\n'
     2766    rundata += '    strcons = str(type(c))\n'
     2767    rundata += '    if type(c) is fPDF.PairDistributionConstraint:\n'
     2768    rundata += '        c.set_variance_squared(wtDict["Pair-"+c.weighting])\n'
     2769    rundata += '        c.set_limits((None,%.3f))\n'%(rmax)
    27192770    if RMCPdict['FitScale']:
    2720         rundata += '        constraint.set_adjust_scale_factor((10, 0.01, 100.))\n'
     2771        rundata += '        c.set_adjust_scale_factor((10, 0.01, 100.))\n'
    27212772    rundata += '    elif "StructureFactor" in strcons:\n'
    2722     rundata += '        constraint.set_variance_squared(wtDict["Struct-"+constraint.weighting])\n'
     2773    rundata += '        c.set_variance_squared(wtDict["Struct-"+c.weighting])\n'
    27232774    if RMCPdict['FitScale']:
    2724         rundata += '        constraint.set_adjust_scale_factor((10, 0.01, 100.))\n'
    2725     if RMCPdict['byMolec']:
    2726         if len(BondList):
    2727             rundata += '    elif "BondConstraint" in strcons:\n'
    2728             rundata += '        constraint.set_variance_squared(%f)\n'%RMCPdict['Bond Weight']
    2729             rundata += '        constraint.create_bonds_by_definition(bondsDefinition={"%s":[\n'%Res
    2730             for bond in BondList:
    2731                 rundata += '        %s'%bond
    2732             rundata += '        ]})\n'
    2733         if len(AngleList):
    2734             rundata += '    elif "BondsAngleConstraint" in strcons:\n'
    2735             rundata += '        constraint.set_variance_squared(%f)\n'%RMCPdict['Angle Weight']
    2736             rundata += '        constraint.create_angles_by_definition(anglesDefinition={"%s":[\n'%Res
    2737             for angle in AngleList:
    2738                 rundata += '        %s'%angle
    2739             rundata += '        ]})\n'
    2740         if len(RMCPdict['Torsions']):
    2741             rundata += '    elif "DihedralAngleConstraint" in strcons:\n'
    2742             rundata += '        constraint.set_variance_squared(%f)\n'%RMCPdict['Torsion Weight']
    2743             rundata += '        constraint.create_angles_by_definition(anglesDefinition={"%s":[\n'%Res
    2744             for torsion in RMCPdict['Torsions']:
    2745                 rundata += '    %s\n'%str(tuple(torsion))
    2746             rundata += '        ]})\n'
    2747     if len(RMCPdict['Swaps']):
    2748         rundata += '    allNames = ENGINE.allNames\n'
    2749         for swap in RMCPdict['Swaps']:
    2750             rundata += '    SwapA = [[idx] for idx in range(len(allNames)) if allNames[idx]=="%s"]\n'%swap[0]
    2751             rundata += '    SwapB = [[idx] for idx in range(len(allNames)) if allNames[idx]=="%s"]\n'%swap[1]
    2752             rundata += '    SwapGen["%s-%s"] = [SwapPositionsGenerator(swapList=SwapA),SwapPositionsGenerator(swapList=SwapB),%.2f]\n'%(swap[0],swap[1],swap[2])
    2753     rundata += 'ENGINE.save()\n'
    2754     rundata += '#setup runs for fullrmc\n'
     2775        rundata += '        c.set_adjust_scale_factor((10, 0.01, 100.))\n'
     2776#    if AngleList and not RMCPdict['Swaps']: rundata += setAngleConstraints()
     2777    # if len(RMCPdict['Torsions']):         # Torsions currently commented out in GUI
     2778    #     rundata += 'for c in ENGINE.constraints:  # look for Dihedral Angle Constraints\n'
     2779    #     rundata += '    if type(c) is DihedralAngleConstraint:\n'
     2780    #     rundata += '        c.set_variance_squared(%f)\n'%RMCPdict['Torsion Weight']
     2781    #     rundata += '        c.create_angles_by_definition(anglesDefinition={"%s":[\n'%Res
     2782    #     for torsion in RMCPdict['Torsions']:
     2783    #         rundata += '    %s\n'%str(tuple(torsion))
     2784    #     rundata += '        ]})\n'           
     2785    rundata += '\n# setup runs for fullrmc\n'
    27552786
    27562787    rundata += 'for _ in range(%d):\n'%RMCPdict['Cycles']
     2788    if BondList and RMCPdict['Swaps']: rundata += setBondConstraints('    ')
     2789    if AngleList and RMCPdict['Swaps']: rundata += setAngleConstraints('    ')
    27572790    rundata += '    ENGINE.set_groups_as_atoms()\n'
    27582791    rundata += '    ENGINE.run(restartPdb="%s",numberOfSteps=10000, saveFrequency=1000)\n'%restart
    2759     if len(RMCPdict['Swaps']):
     2792    if RMCPdict['Swaps']:
     2793        if BondList: rundata += setBondConstraints('    ',clear=True)
     2794        if AngleList: rundata += setAngleConstraints('    ',clear=True)
    27602795        rundata += '    for swaps in SwapGen:\n'
    27612796        rundata += '        AB = swaps.split("-")\n'
    27622797        rundata += '        ENGINE.set_groups_as_atoms()\n'
    27632798        rundata += '        for g in ENGINE.groups:\n'
    2764         rundata += '            if allNames[g.indexes[0]]==AB[0]:\n'
     2799        rundata += '            if aN[g.indexes[0]]==AB[0]:\n'
    27652800        rundata += '                g.set_move_generator(SwapGen[swaps][0])\n'
    2766         rundata += '            elif allNames[g.indexes[0]]==AB[1]:\n'
     2801        rundata += '            elif aN[g.indexes[0]]==AB[1]:\n'
    27672802        rundata += '                g.set_move_generator(SwapGen[swaps][1])\n'
    27682803        rundata += '            sProb = SwapGen[swaps][2]\n'
    … …  
    27702805        rundata += '        ENGINE.set_groups_as_atoms()\n'
    27712806        rundata += '        ENGINE.run(restartPdb="%s",numberOfSteps=10000*(1.-sProb), saveFrequency=1000)\n'%restart
    2772     rundata += 'ENGINE.close()\n'
     2807    #rundata += 'ENGINE.close()\n'
    27732808    rundata += 'print("ENGINE run time %.2f s"%(time.time()-time0))\n'
    27742809    rfile = open(rname,'w')
    … …  
    27782813    return rname
    27792814
    2780 def MakefullrmcPDB(Name,Phase,RMCPdict):
    2781     generalData = Phase['General']
    2782     Atseq = RMCPdict['atSeq']
    2783     Dups,Fracs = findDup(Phase['Atoms'])
    2784     Sfracs = [np.cumsum(fracs) for fracs in Fracs]
    2785     ifSfracs = any([np.any(sfracs-1.) for sfracs in Sfracs])
    2786     Supercell = RMCPdict['SuperCell']
    2787     Cell = generalData['Cell'][1:7]
    2788     Trans = np.eye(3)*np.array(Supercell)
    2789     newPhase = copy.deepcopy(Phase)
    2790     newPhase['General']['SGData'] = G2spc.SpcGroup('P 1')[1]
    2791     newPhase['General']['Cell'][1:] = G2lat.TransformCell(Cell,Trans.T)
    2792     newPhase,Atcodes = G2lat.TransformPhase(Phase,newPhase,Trans,np.zeros(3),np.zeros(3),ifMag=False,Force=True)
    2793     Atoms = newPhase['Atoms']
    2794 
    2795     if ifSfracs:
    2796         Natm = np.core.defchararray.count(np.array(Atcodes),'+')    #no. atoms in original unit cell
    2797         Natm = np.count_nonzero(Natm-1)
    2798         Satoms = []
    2799         for i in range(len(Atoms)//Natm):
    2800             ind = i*Natm
    2801             Satoms.append(G2mth.sortArray(G2mth.sortArray(G2mth.sortArray(Atoms[ind:ind+Natm],5),4),3))
    2802         Natoms = []
    2803         for satoms in Satoms:
    2804             for idup,dup in enumerate(Dups):
    2805                 ldup = len(dup)
    2806                 natm = len(satoms)
    2807                 i = 0
    2808                 while i < natm:
    2809                     if satoms[i][0] in dup:
    2810                         atoms = satoms[i:i+ldup]
    2811                         try:
    2812                             atom = atoms[np.searchsorted(Sfracs[idup],rand.random())]
    2813                             Natoms.append(atom)
    2814                         except IndexError:      #what about vacancies?
    2815                             if 'Va' not in Atseq:
    2816                                 Atseq.append('Va')
    2817                                 RMCPdict['aTypes']['Va'] = 0.0
    2818                             atom = atoms[0]
    2819                             atom[1] = 'Va'
    2820                             Natoms.append(atom)
    2821                         i += ldup
    2822                     else:
    2823                        i += 1
    2824     else:
    2825         Natoms = Atoms
    2826 
    2827     XYZ = np.array([atom[3:6] for atom in Natoms]).T
    2828     XYZptp = np.array([ma.ptp(XYZ[0]),ma.ptp(XYZ[1]),ma.ptp(XYZ[2])])/2.
    2829     Cell = newPhase['General']['Cell'][1:7]
    2830     A,B = G2lat. cell2AB(Cell)
    2831     fname = Name+'_cbb.pdb'
    2832     fl = open(fname,'w')
    2833     fl.write('REMARK    Boundary Conditions:%6.2f  0.0  0.0  0.0%7.2f  0.0  0.0  0.0%7.2f\n'%(
    2834              Cell[0],Cell[1],Cell[2]))
    2835     fl.write('ORIGX1      1.000000  0.000000  0.000000        0.00000\n')
    2836     fl.write('ORIGX2      0.000000  1.000000  0.000000        0.00000\n')
    2837     fl.write('ORIGX3      0.000000  0.000000  1.000000        0.00000\n')
    2838     fl.write('CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f P 1           1\n'%(
    2839             Cell[0],Cell[1],Cell[2],Cell[3],Cell[4],Cell[5]))
    2840 
    2841     Natm = np.core.defchararray.count(np.array(Atcodes),'+')
    2842     Natm = np.count_nonzero(Natm-1)
    2843     nat = 0
    2844     if RMCPdict['byMolec']:
    2845         NPM = RMCPdict['Natoms']
    2846         for iat,atom in enumerate(Natoms):
    2847             XYZ = np.inner(A,np.array(atom[3:6])-XYZptp)    #shift origin to middle & make Cartesian;residue = 'RMC'
    2848             fl.write('ATOM  %5d %-4s RMC%6d%12.3f%8.3f%8.3f  1.00  0.00          %2s\n'%(       
    2849                     1+nat%NPM,atom[0],1+nat//NPM,XYZ[0],XYZ[1],XYZ[2],atom[1].lower()))
    2850             nat += 1
    2851     else:
    2852         for atm in Atseq:
    2853             for iat,atom in enumerate(Natoms):
    2854                 if atom[1] == atm:
    2855                     XYZ = np.inner(A,np.array(atom[3:6])-XYZptp)    #shift origin to middle & make Cartesian
    2856                     fl.write('ATOM  %5d %-4s RMC%6d%12.3f%8.3f%8.3f  1.00  0.00          %2s\n'%(       
    2857                             1+nat,atom[0],1+nat,XYZ[0],XYZ[1],XYZ[2],atom[1].lower()))
    2858                     nat += 1
    2859     fl.close()
    2860     return fname
     2815# def MakefullrmcPDB(Name,Phase,RMCPdict):
     2816#     generalData = Phase['General']
     2817#     Atseq = RMCPdict['atSeq']
     2818#     Dups,Fracs = findDup(Phase['Atoms'])
     2819#     Sfracs = [np.cumsum(fracs) for fracs in Fracs]
     2820#     ifSfracs = any([np.any(sfracs-1.) for sfracs in Sfracs])
     2821#     Supercell = RMCPdict['SuperCell']
     2822#     Cell = generalData['Cell'][1:7]
     2823#     Trans = np.eye(3)*np.array(Supercell)
     2824#     newPhase = copy.deepcopy(Phase)
     2825#     newPhase['General']['SGData'] = G2spc.SpcGroup('P 1')[1]
     2826#     newPhase['General']['Cell'][1:] = G2lat.TransformCell(Cell,Trans.T)
     2827#     newPhase,Atcodes = G2lat.TransformPhase(Phase,newPhase,Trans,np.zeros(3),np.zeros(3),ifMag=False,Force=True)
     2828#     Atoms = newPhase['Atoms']
     2829
     2830#     if ifSfracs:
     2831#         Natm = np.core.defchararray.count(np.array(Atcodes),'+')    #no. atoms in original unit cell
     2832#         Natm = np.count_nonzero(Natm-1)
     2833#         Satoms = []
     2834#         for i in range(len(Atoms)//Natm):
     2835#             ind = i*Natm
     2836#             Satoms.append(G2mth.sortArray(G2mth.sortArray(G2mth.sortArray(Atoms[ind:ind+Natm],5),4),3))
     2837#         Natoms = []
     2838#         for satoms in Satoms:
     2839#             for idup,dup in enumerate(Dups):
     2840#                 ldup = len(dup)
     2841#                 natm = len(satoms)
     2842#                 i = 0
     2843#                 while i < natm:
     2844#                     if satoms[i][0] in dup:
     2845#                         atoms = satoms[i:i+ldup]
     2846#                         try:
     2847#                             atom = atoms[np.searchsorted(Sfracs[idup],rand.random())]
     2848#                             Natoms.append(atom)
     2849#                         except IndexError:      #what about vacancies?
     2850#                             if 'Va' not in Atseq:
     2851#                                 Atseq.append('Va')
     2852#                                 RMCPdict['aTypes']['Va'] = 0.0
     2853#                             atom = atoms[0]
     2854#                             atom[1] = 'Va'
     2855#                             Natoms.append(atom)
     2856#                         i += ldup
     2857#                     else:
     2858#                        i += 1
     2859#     else:
     2860#         Natoms = Atoms
     2861
     2862#     XYZ = np.array([atom[3:6] for atom in Natoms]).T
     2863#     XYZptp = np.array([ma.ptp(XYZ[0]),ma.ptp(XYZ[1]),ma.ptp(XYZ[2])])/2.
     2864#     Cell = newPhase['General']['Cell'][1:7]
     2865#     A,B = G2lat. cell2AB(Cell)
     2866#     fname = Name+'_cbb.pdb'
     2867#     fl = open(fname,'w')
     2868#     fl.write('REMARK    Boundary Conditions:%6.2f  0.0  0.0  0.0%7.2f  0.0  0.0  0.0%7.2f\n'%(
     2869#              Cell[0],Cell[1],Cell[2]))
     2870#     fl.write('ORIGX1      1.000000  0.000000  0.000000        0.00000\n')
     2871#     fl.write('ORIGX2      0.000000  1.000000  0.000000        0.00000\n')
     2872#     fl.write('ORIGX3      0.000000  0.000000  1.000000        0.00000\n')
     2873#     fl.write('CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f P 1           1\n'%(
     2874#             Cell[0],Cell[1],Cell[2],Cell[3],Cell[4],Cell[5]))
     2875
     2876#     Natm = np.core.defchararray.count(np.array(Atcodes),'+')
     2877#     Natm = np.count_nonzero(Natm-1)
     2878#     nat = 0
     2879#     if RMCPdict['byMolec']:
     2880#         NPM = RMCPdict['Natoms']
     2881#         for iat,atom in enumerate(Natoms):
     2882#             XYZ = np.inner(A,np.array(atom[3:6])-XYZptp)    #shift origin to middle & make Cartesian;residue = 'RMC'
     2883#             fl.write('ATOM  %5d %-4s RMC%6d%12.3f%8.3f%8.3f  1.00  0.00          %2s\n'%(       
     2884#                     1+nat%NPM,atom[0],1+nat//NPM,XYZ[0],XYZ[1],XYZ[2],atom[1].lower()))
     2885#             nat += 1
     2886#     else:
     2887#         for atm in Atseq:
     2888#             for iat,atom in enumerate(Natoms):
     2889#                 if atom[1] == atm:
     2890#                     XYZ = np.inner(A,np.array(atom[3:6])-XYZptp)    #shift origin to middle & make Cartesian
     2891#                     fl.write('ATOM  %5d %-4s RMC%6d%12.3f%8.3f%8.3f  1.00  0.00          %2s\n'%(       
     2892#                             1+nat,atom[0],1+nat,XYZ[0],XYZ[1],XYZ[2],atom[1].lower()))
     2893#                     nat += 1
     2894#     fl.close()
     2895#     return fname
    28612896   
    28622897def MakePdparse(RMCPdict):
    … …  
    32523287    Reflectometry as a function of kz for a set of slabs.
    32533288
    3254     :param kz: float[n] (1/Ang). Scattering vector, :math:`2\pi\sin(\\theta)/\lambda`.
     3289    :param kz: float[n] (1/Ang). Scattering vector, :math:`2\\pi\\sin(\\theta)/\\lambda`.
    32553290        This is :math:`\\tfrac12 Q_z`.       
    32563291    :param depth:  float[m] (Ang).

  • trunk/GSASIIscriptable.py

    r4517 r4518  
    676676subset have direct mechanism implemented for change from the GSASIIscriptable
    677677API. In practice all parameters in a .gpx file can be edited via scripting,
    678 but sometimes determining what should be changed to implement a change can be
    679 complex.
    680 Several routines, :meth:`G2Phase.getHAPentryValue`,
    681 :meth:`G2Phase.getPhaseEntryValue` and :meth:`G2PwdrData.getHistEntryList`
     678but sometimes determining what should be set to implement a parameter
     679change can be complex.
     680Several routines, :meth:`G2Phase.getHAPentryList`,
     681:meth:`G2Phase.getPhaseEntryList` and :meth:`G2PwdrData.getHistEntryList`
     682(and their related get...Value and set.Value entries),
    682683provide a mechanism to discover what the GUI is changing inside a .gpx file.
    683684
    … …  
    43794380        See :meth:`G2Phase.getHAPentryList` for a related example.
    43804381
     4382        .. seealso::
     4383            :meth:`getHistEntryValue`
     4384            :meth:`setHistEntryValue`
     4385
    43814386        """
    43824387        return [i for i in dictDive(self.data,keyname) if i[0] != ['Histograms']]
    … …  
    51805185        See :meth:`getHAPentryList` for a related example.
    51815186
     5187        .. seealso::
     5188            :meth:`getPhaseEntryValue`
     5189            :meth:`setPhaseEntryValue`
     5190
    51825191        """
    51835192        return [i for i in dictDive(self.data,keyname) if i[0] != ['Histograms']]
    … …  
    51875196        Where the value returned is a list, it may be used as the target of
    51885197        an assignment (as in
    5189         ``getHistEntryValue(...)[...] = val``)
     5198        ``getPhaseEntryValue(...)[...] = val``)
    51905199        to set a value inside a list.       
    51915200
    … …  
    52455254        [(['PWDR test Bank 1', 'Scale'], list, [1.0, False])]
    52465255
     5256        .. seealso::
     5257            :meth:`getHAPentryValue`
     5258            :meth:`setHAPentryValue`
     5259
    52475260        """
    52485261        if histname:
    … …  
    52565269        value returned is a list, it may be used as the target of
    52575270        an assignment (as in
    5258         ``getHAPentryValue(``...``)[``...``] = ``val)
     5271        ``getHAPentryValue(...)[...] = val``)
    52595272        to set a value inside a list.
    52605273
Note: See TracChangeset for help on using the changeset viewer.