source: trunk/GSASIIconstrGUI.py @ 4980

Last change on this file since 4980 was 4980, checked in by toby, 17 months ago

Fix coeff. display for EQUIV

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 169.7 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASIIconstrGUI - constraint GUI routines
3########### SVN repository information ###################
4# $Date: 2021-06-28 16:38:26 +0000 (Mon, 28 Jun 2021) $
5# $Author: toby $
6# $Revision: 4980 $
7# $URL: trunk/GSASIIconstrGUI.py $
8# $Id: GSASIIconstrGUI.py 4980 2021-06-28 16:38:26Z toby $
9########### SVN repository information ###################
10'''
11*GSASIIconstrGUI: Constraint GUI routines*
12------------------------------------------
13
14Used to define constraints and rigid bodies.
15
16'''
17from __future__ import division, print_function
18import platform
19import sys
20import copy
21import os.path
22import wx
23import wx.grid as wg
24import wx.lib.scrolledpanel as wxscroll
25import wx.lib.gridmovers as gridmovers
26import random as ran
27import numpy as np
28import numpy.ma as ma
29import numpy.linalg as nl
30import GSASIIpath
31GSASIIpath.SetVersionNumber("$Revision: 4980 $")
32import GSASIIElem as G2elem
33import GSASIIElemGUI as G2elemGUI
34import GSASIIstrIO as G2stIO
35import GSASIImapvars as G2mv
36import GSASIImath as G2mth
37import GSASIIlattice as G2lat
38import GSASIIdataGUI as G2gd
39import GSASIIctrlGUI as G2G
40import GSASIIfiles as G2fl
41import GSASIIplot as G2plt
42import GSASIIobj as G2obj
43import GSASIIspc as G2spc
44import GSASIIpy3 as G2py3
45import GSASIIphsGUI as G2phG
46import GSASIIIO as G2IO
47import GSASIIscriptable as G2sc
48VERY_LIGHT_GREY = wx.Colour(235,235,235)
49WACV = wx.ALIGN_CENTER_VERTICAL
50
51class G2BoolEditor(wg.GridCellBoolEditor):
52    '''Substitute for wx.grid.GridCellBoolEditor except toggles
53    grid items immediately when opened, updates grid & table contents after every
54    item change
55    '''
56    def __init__(self):
57        self.saveVals = None
58        wx.grid.GridCellBoolEditor.__init__(self)
59
60
61    def Create(self, parent, id, evtHandler):
62        '''Create the editing control (wx.CheckBox) when cell is opened
63        for edit
64        '''
65        self._tc = wx.CheckBox(parent, -1, "")
66        self._tc.Bind(wx.EVT_CHECKBOX, self.onCheckSet)
67        self.SetControl(self._tc)
68        if evtHandler:
69            self._tc.PushEventHandler(evtHandler)
70
71    def onCheckSet(self, event):
72        '''Callback used when checkbox is toggled.
73        Makes change to table immediately (creating event)
74        '''
75        if self.saveVals:
76            self.ApplyEdit(*self.saveVals)
77
78       
79    def SetSize(self, rect):
80        '''Set position/size the edit control within the cell's rectangle.
81        '''
82#        self._tc.SetDimensions(rect.x, rect.y, rect.width+2, rect.height+2, # older
83        self._tc.SetSize(rect.x, rect.y, rect.width+2, rect.height+2,
84                               wx.SIZE_ALLOW_MINUS_ONE)
85
86    def BeginEdit(self, row, col, grid):
87        '''Prepares the edit control by loading the initial
88        value from the table (toggles it since you would not
89        click on it if you were not planning to change it),
90        buts saves the original, pre-change value.
91        Makes change to table immediately.
92        Saves the info needed to make updates in self.saveVals.
93        Sets the focus.
94        '''
95        if grid.GetTable().GetValue(row,col) not in [True,False]:
96            return
97        self.startValue = int(grid.GetTable().GetValue(row, col))
98        self.saveVals = row, col, grid
99        # invert state and set in editor
100        if self.startValue:
101            grid.GetTable().SetValue(row, col, 0)
102            self._tc.SetValue(0)
103        else:
104            grid.GetTable().SetValue(row, col, 1)
105            self._tc.SetValue(1)
106        self._tc.SetFocus()
107        self.ApplyEdit(*self.saveVals)
108
109    def EndEdit(self, row, col, grid, oldVal=None):
110        '''End editing the cell.  This is supposed to
111        return None if the value has not changed, but I am not
112        sure that actually works.
113        '''
114        val = int(self._tc.GetValue())
115        if val != oldVal:   #self.startValue:?
116            return val
117        else:
118            return None
119
120    def ApplyEdit(self, row, col, grid):
121        '''Save the value into the table, and create event.
122        Called after EndEdit(), BeginEdit and onCheckSet.
123        '''
124        val = int(self._tc.GetValue())
125        grid.GetTable().SetValue(row, col, val) # update the table
126
127    def Reset(self):
128        '''Reset the value in the control back to its starting value.
129        '''
130        self._tc.SetValue(self.startValue)
131
132    def StartingClick(self):
133        '''This seems to be needed for BeginEdit to work properly'''
134        pass
135
136    def Destroy(self):
137        "final cleanup"
138        super(G2BoolEditor, self).Destroy()
139
140    def Clone(self):
141        'required'
142        return G2BoolEditor()
143   
144class DragableRBGrid(wg.Grid):
145    '''Simple grid implentation for display of rigid body positions.
146
147    :param parent: frame or panel where grid will be placed
148    :param dict rb: dict with atom labels, types and positions
149    :param function onChange: a callback used every time a value in
150      rb is changed.
151    '''
152    def __init__(self, parent, rb, onChange=None):
153        #wg.Grid.__init__(self, parent, wx.ID_ANY,size=(-1,200))
154        wg.Grid.__init__(self, parent, wx.ID_ANY)
155        self.SetTable(RBDataTable(rb,onChange), True)
156        # Enable Row moving
157        gridmovers.GridRowMover(self)
158        self.Bind(gridmovers.EVT_GRID_ROW_MOVE, self.OnRowMove, self)
159        self.SetColSize(0, 60)
160        self.SetColSize(1, 40)
161        self.SetColSize(2, 35)
162        for r in range(len(rb['RBlbls'])):
163            self.SetReadOnly(r,0,isReadOnly=True)
164            self.SetCellEditor(r, 1, G2BoolEditor())           
165            self.SetCellRenderer(r, 1, wg.GridCellBoolRenderer())
166            self.SetReadOnly(r,2,isReadOnly=True)
167            self.SetCellEditor(r,3, wg.GridCellFloatEditor())
168            self.SetCellEditor(r,4, wg.GridCellFloatEditor())
169            self.SetCellEditor(r,6, wg.GridCellFloatEditor())
170
171    def OnRowMove(self,evt):
172        'called when a row move needs to take place'
173        frm = evt.GetMoveRow()          # Row being moved
174        to = evt.GetBeforeRow()         # Before which row to insert
175        self.GetTable().MoveRow(frm,to)
176       
177    def completeEdits(self):
178        'complete any outstanding edits'
179        if self.IsCellEditControlEnabled(): # complete any grid edits in progress
180            #if GSASIIpath.GetConfigValue('debug'): print ('Completing grid edit')
181            self.SaveEditControlValue()
182            self.HideCellEditControl()
183            self.DisableCellEditControl()
184           
185class RBDataTable(wg.GridTableBase):
186    '''A Table to support :class:`DragableRBGrid`
187    '''
188    def __init__(self,rb,onChange):
189        wg.GridTableBase.__init__(self)
190        self.colLabels = ['Label','Select','Type','x','y','z']
191        self.coords = rb['RBcoords']
192        self.labels = rb['RBlbls']
193        self.types = rb['RBtypes']
194        self.index = rb['RBindex']
195        self.select = rb['RBselection']
196        self.onChange = onChange
197
198    # required methods
199    def GetNumberRows(self):
200        return len(self.labels)
201    def GetNumberCols(self):
202        return len(self.colLabels)
203    def IsEmptyCell(self, row, col):
204        return False
205    def GetValue(self, row, col):
206        row = self.index[row]
207        if col == 0:
208            return self.labels[row]
209        elif col == 1:
210            if self.select[row]:
211                return '1'
212            else:
213                return ''
214        elif col == 2:
215            return self.types[row]
216        else:
217            return '{:.5f}'.format(self.coords[row][col-3])
218    def SetValue(self, row, col, value):
219        row = self.index[row]
220        try:
221            if col == 0:
222                self.labels[row] = value
223            elif col == 1:
224                self.select[row] = bool(value)
225            elif col == 2:
226                self.types[row] = value
227            else:
228                self.coords[row][col-3] = float(value)
229        except:
230            pass
231        if self.onChange:
232            self.onChange()
233    # Display column & row labels
234    def GetColLabelValue(self, col):
235        return self.colLabels[col]
236    def GetRowLabelValue(self,row):
237        return str(row)
238
239    # Implement "row movement" by updating the pointer array
240    def MoveRow(self,frm,to):
241        grid = self.GetView()
242        if grid:
243            move = self.index[frm]
244            del self.index[frm]
245            if frm > to:
246                self.index.insert(to,move)
247            else:
248                self.index.insert(to-1,move)
249           
250            # Notify the grid
251            grid.BeginBatch()
252            msg = wg.GridTableMessage(
253                    self, wg.GRIDTABLE_NOTIFY_ROWS_DELETED, frm, 1
254                    )
255            grid.ProcessTableMessage(msg)
256            msg = wg.GridTableMessage(
257                    self, wg.GRIDTABLE_NOTIFY_ROWS_INSERTED, to, 1
258                    )
259            grid.ProcessTableMessage(msg)
260            grid.EndBatch()
261        if self.onChange:
262            self.onChange()
263
264# def MakeDrawAtom(data,atom):
265#     'Convert atom to format needed to draw it'
266#     generalData = data['General']
267#     deftype = G2obj.validateAtomDrawType(
268#         GSASIIpath.GetConfigValue('DrawAtoms_default'),generalData)
269#     if generalData['Type'] in ['nuclear','faulted',]:
270#         atomInfo = [atom[:2]+atom[3:6]+['1']+[deftype]+
271#                     ['']+[[255,255,255]]+atom[9:]+[[],[]]][0]
272#     ct,cs = [1,8]         #type & color
273#     atNum = generalData['AtomTypes'].index(atom[ct])
274#     atomInfo[cs] = list(generalData['Color'][atNum])
275#     return atomInfo
276
277class ConstraintDialog(wx.Dialog):
278    '''Window to edit Constraint values
279    '''
280    def __init__(self,parent,title,text,data,separator='*',varname="",varyflag=False):
281        wx.Dialog.__init__(self,parent,-1,'Edit '+title, 
282            pos=wx.DefaultPosition,style=wx.DEFAULT_DIALOG_STYLE)
283        self.data = data[:]
284        self.newvar = [varname,varyflag]
285        panel = wx.Panel(self)
286        mainSizer = wx.BoxSizer(wx.VERTICAL)
287        topLabl = wx.StaticText(panel,-1,text)
288        mainSizer.Add((10,10),1)
289        mainSizer.Add(topLabl,0,wx.LEFT,10)
290        mainSizer.Add((10,10),1)
291        dataGridSizer = wx.FlexGridSizer(cols=3,hgap=2,vgap=2)
292        self.OkBtn = wx.Button(panel,wx.ID_OK)
293        for id in range(len(self.data)):
294            lbl1 = lbl = str(self.data[id][1])
295            if lbl[-1] != '=': lbl1 = lbl + ' ' + separator + ' '
296            name = wx.StaticText(panel,wx.ID_ANY,lbl1,style=wx.ALIGN_RIGHT)
297            scale = G2G.ValidatedTxtCtrl(panel,self.data[id],0,OKcontrol=self.DisableOK)
298            dataGridSizer.Add(name,0,wx.LEFT|wx.RIGHT|WACV,5)
299            dataGridSizer.Add(scale,0,wx.RIGHT,3)
300            if ':' in lbl:
301                dataGridSizer.Add(
302                    wx.StaticText(panel,-1,G2obj.fmtVarDescr(lbl)),
303                    0,wx.RIGHT|WACV,3)
304            else:
305                dataGridSizer.Add((-1,-1))
306        if title == 'New Variable':
307            name = wx.StaticText(panel,wx.ID_ANY,"New variable's\nname (optional)",
308                style=wx.ALIGN_CENTER)
309            scale = G2G.ValidatedTxtCtrl(panel,self.newvar,0,notBlank=False)
310            dataGridSizer.Add(name,0,wx.LEFT|wx.RIGHT|WACV,5)
311            dataGridSizer.Add(scale,0,wx.RIGHT|WACV,3)
312            self.refine = wx.CheckBox(panel,label='Refine?')
313            self.refine.SetValue(self.newvar[1]==True)
314            self.refine.Bind(wx.EVT_CHECKBOX, self.OnCheckBox)
315            dataGridSizer.Add(self.refine,0,wx.RIGHT|WACV,3)
316           
317        mainSizer.Add(dataGridSizer,0,wx.EXPAND)
318        self.OkBtn.Bind(wx.EVT_BUTTON, self.OnOk)
319        self.OkBtn.SetDefault()
320        cancelBtn = wx.Button(panel,wx.ID_CANCEL)
321        cancelBtn.Bind(wx.EVT_BUTTON, self.OnCancel)
322        btnSizer = wx.BoxSizer(wx.HORIZONTAL)
323        btnSizer.Add((20,20),1)
324        btnSizer.Add(self.OkBtn)
325        btnSizer.Add((20,20),1)
326        btnSizer.Add(cancelBtn)
327        btnSizer.Add((20,20),1)
328
329        mainSizer.Add(btnSizer,0,wx.EXPAND, 10)
330        panel.SetSizer(mainSizer)
331        panel.Fit()
332        self.Fit()
333        self.CenterOnParent()
334       
335    def DisableOK(self,setting):
336        for id in range(len(self.data)):  # coefficient cannot be zero
337            if abs(self.data[id][0]) < 1.e-20:
338                setting  = False
339                break   
340        if setting:
341            self.OkBtn.Enable()
342        else:
343            self.OkBtn.Disable()
344
345    def OnCheckBox(self,event):
346        self.newvar[1] = self.refine.GetValue()
347
348    def OnOk(self,event):
349        parent = self.GetParent()
350        parent.Raise()
351        self.EndModal(wx.ID_OK)             
352
353    def OnCancel(self,event):
354        parent = self.GetParent()
355        parent.Raise()
356        self.EndModal(wx.ID_CANCEL)             
357
358    def GetData(self):
359        return self.data
360       
361#####  Constraints ################################################################################           
362def UpdateConstraints(G2frame,data):
363    '''Called when Constraints tree item is selected.
364    Displays the constraints in the data window
365    '''
366       
367    def FindEquivVarb(name,nameList):
368        'Creates a list of variables appropriate to constrain with name'
369        outList = []
370        #phlist = []
371        items = name.split(':')
372        namelist = [items[2],]
373        if 'dA' in name:
374            namelist = ['dAx','dAy','dAz']
375        elif 'AU' in name:
376            namelist = ['AUiso','AU11','AU22','AU33','AU12','AU13','AU23']
377        elif 'AM' in name:
378            namelist = ['AMx','AMy','AMz']
379        elif items[-1] in ['A0','A1','A2','A3','A4','A5']:
380            namelist = ['A0','A1','A2','A3','A4','A5']
381        elif items[-1] in ['D11','D22','D33','D12','D13','D23']:
382            namelist = ['D11','D22','D33','D12','D13','D23']
383        elif 'Tm' in name:
384            namelist = ['Tmin','Tmax']
385        elif 'MX' in name or 'MY' in name or 'MZ' in name:
386            namelist = ['MXcos','MYcos','MZcos','MXsin','MYsin','MZsin']
387        elif 'mV' in name:
388            namelist = ['mV0','mV1','mV2']
389        elif 'RB' in name:
390            rbfx = 'RB'+items[2][2]
391            if 'T' in name and 'Tr' not in name:
392                namelist = [rbfx+'T11',rbfx+'T22',rbfx+'T33',rbfx+'T12',rbfx+'T13',rbfx+'T23']
393            if 'L' in name:
394                namelist = [rbfx+'L11',rbfx+'L22',rbfx+'L33',rbfx+'L12',rbfx+'L13',rbfx+'L23']
395            if 'S' in name:
396                namelist = [rbfx+'S12',rbfx+'S13',rbfx+'S21',rbfx+'S23',rbfx+'S31',rbfx+'S32',rbfx+'SAA',rbfx+'SBB']
397            if 'U' in name:
398                namelist = [rbfx+'U',]
399
400        for item in nameList:
401            keys = item.split(':')
402            #if keys[0] not in phlist:
403            #    phlist.append(keys[0])
404            if items[1] == '*' and keys[2] in namelist: # wildcard -- select only sequential options
405                keys[1] = '*'
406                mitem = ':'.join(keys)
407                if mitem == name: continue
408                if mitem not in outList: outList.append(mitem)
409            elif keys[2] in namelist and item != name:
410                outList.append(item)
411        return outList
412       
413    def SelectVarbs(page,FrstVarb,varList,legend,constType):
414        '''Select variables used in constraints after one variable has
415        been selected. This routine determines the appropriate variables to be
416        used based on the one that has been selected and asks for more to be added.
417
418        It then creates the constraint and adds it to the constraints list.
419       
420        Called from OnAddEquivalence, OnAddFunction & OnAddConstraint (all but
421        OnAddHold)
422
423        :param list page: defines calling page -- type of variables to be used
424        :parm GSASIIobj.G2VarObj FrstVarb: reference to first selected variable
425        :param list varList: list of other appropriate variables to select from
426        :param str legend: header for selection dialog
427        :param str constType: type of constraint to be generated
428        :returns: a constraint, as defined in
429          :ref:`GSASIIobj <Constraint_definitions_table>`
430        '''
431        choices = [[i]+list(G2obj.VarDescr(i)) for i in varList]
432        meaning = G2obj.getDescr(FrstVarb.name)
433        if not meaning:
434            meaning = "(no definition found!)"
435        l = str(FrstVarb).split(':')
436        # make lists of phases & histograms to iterate over
437        phaselist = [l[0]]
438        if l[0]:
439            phaselbl = ['phase #'+l[0]]
440            if len(Phases) > 1:
441                phaselist += ['all'] 
442                phaselbl += ['all phases']
443        else:
444            phaselbl = ['']
445        histlist = [l[1]]
446        if l[1] == '*':
447            pass
448        elif l[1]:
449            histlbl = ['histogram #'+l[1]]
450            if len(Histograms) > 1:
451                histlist += ['all']
452                histlbl += ['all histograms']
453                typ = Histograms[G2obj.LookupHistName(l[1])[0]]['Instrument Parameters'][0]['Type'][1]
454                i = 0
455                for hist in Histograms:
456                    if Histograms[hist]['Instrument Parameters'][0]['Type'][1] == typ: i += 1
457                if i > 1:
458                    histlist += ['all='+typ]
459                    histlbl += ['all '+typ+' histograms']
460        else:
461            histlbl = ['']
462        # make a list of equivalent parameter names
463        nameList = [FrstVarb.name]
464        for var in varList:
465            nam = var.split(":")[2]
466            if nam not in nameList: nameList += [nam]
467        # add "wild-card" names to the list of variables
468        if l[1] == '*':
469            pass
470        elif page[1] == 'phs':
471            if 'RB' in FrstVarb.name:
472                pass
473            elif FrstVarb.atom is None:
474                for nam in nameList:
475                    for ph,plbl in zip(phaselist,phaselbl):
476                        if plbl: plbl = 'For ' + plbl
477                        var = ph+"::"+nam
478                        if var == str(FrstVarb) or var in varList: continue
479                        varList += [var]
480                        choices.append([var,plbl,meaning])
481            else:
482                for nam in nameList:
483                    for ph,plbl in zip(phaselist,phaselbl):
484                        if plbl: plbl = ' in ' + plbl
485                        for atype in ['']+TypeList:
486                            if atype:
487                                albl = "For "+atype+" atoms"
488                                akey = "all="+atype                       
489                            else:
490                                albl = "For all atoms"
491                                akey = "all"
492                            var = ph+"::"+nam+":"+akey
493                            if var == str(FrstVarb) or var in varList: continue
494                            varList += [var]
495                            choices.append([var,albl+plbl,meaning])
496        elif page[1] == 'hap':
497            if FrstVarb.name == "Scale":
498                meaning = "Phase fraction"
499            for nam in nameList:
500                for ph,plbl in zip(phaselist,phaselbl):
501                    if plbl: plbl = 'For ' + plbl
502                    for hst,hlbl in zip(histlist,histlbl):
503                        if hlbl:
504                            if plbl:
505                                hlbl = ' in ' + hlbl
506                            else:
507                                hlbl = 'For ' + hlbl                               
508                            var = ph+":"+hst+":"+nam
509                            if var == str(FrstVarb) or var in varList: continue
510                            varList += [var]
511                            choices.append([var,plbl+hlbl,meaning])
512        elif page[1] == 'hst':
513            if FrstVarb.name == "Scale":
514                meaning = "Scale factor"
515            for nam in nameList:
516                for hst,hlbl in zip(histlist,histlbl):
517                    if hlbl:
518                        hlbl = 'For ' + hlbl                               
519                        var = ":"+hst+":"+nam
520                        if var == str(FrstVarb) or var in varList: continue
521                        varList += [var]
522                        choices.append([var,hlbl,meaning])
523        elif page[1] == 'glb' or page[1] == 'sym':
524            pass
525        else:
526            raise Exception('Unknown constraint page '+ page[1])                   
527        if len(choices):
528            l1 = l2 = 1
529            for i1,i2,i3 in choices:
530                l1 = max(l1,len(i1))
531                l2 = max(l2,len(i2))
532            fmt = "{:"+str(l1)+"s} {:"+str(l2)+"s} {:s}"
533            atchoices = [fmt.format(*i1) for i1 in choices] # reformat list as str with columns
534            dlg = G2G.G2MultiChoiceDialog(
535                G2frame,legend,
536                'Constrain '+str(FrstVarb)+' with...',atchoices,
537                toggle=False,size=(625,400),monoFont=True)
538            dlg.CenterOnParent()
539            res = dlg.ShowModal()
540            Selections = dlg.GetSelections()[:]
541            dlg.Destroy()
542            if res != wx.ID_OK: return []
543            if len(Selections) == 0:
544                dlg = wx.MessageDialog(
545                    G2frame,
546                    'No variables were selected to include with '+str(FrstVarb),
547                    'No variables')
548                dlg.CenterOnParent()
549                dlg.ShowModal()
550                dlg.Destroy()
551                return []
552        else:
553            dlg = wx.MessageDialog(
554                G2frame,
555                'There are no appropriate variables to include with '+str(FrstVarb),
556                'No variables')
557            dlg.CenterOnParent()
558            dlg.ShowModal()
559            dlg.Destroy()
560            return []
561        # now process the variables provided by the user
562        varbs = [str(FrstVarb),] # list of selected variables
563        for sel in Selections:
564            var = varList[sel]
565            # phase(s) included
566            l = var.split(':')
567            if l[0] == "all":
568                phlist = [str(Phases[phase]['pId']) for phase in Phases]
569            else:
570                phlist = [l[0]]
571            # histogram(s) included
572            if l[1] == "all":
573                hstlist = [str(Histograms[hist]['hId']) for hist in Histograms]
574            elif '=' in l[1]:
575                htyp = l[1].split('=')[1]
576                hstlist = [str(Histograms[hist]['hId']) for hist in Histograms if 
577                           Histograms[hist]['Instrument Parameters'][0]['Type'][1] == htyp]
578            else:
579                hstlist = [l[1]]
580            if len(l) == 3:
581                for ph in phlist:
582                    for hst in hstlist:
583                        var = ph + ":" + hst + ":" + l[2]
584                        if var in varbs: continue
585                        varbs.append(var)
586            else: # constraints with atoms or rigid bodies
587                if len(l) == 5: # rigid body parameter
588                    var = ':'.join(l)
589                    if var in varbs: continue
590                    varbs.append(var)
591                elif l[3] == "all":
592                    for ph in phlist:
593                        key = G2obj.LookupPhaseName(ph)[0]
594                        for hst in hstlist: # should be blank
595                            for iatm,at in enumerate(Phases[key]['Atoms']):
596                                var = ph + ":" + hst + ":" + l[2] + ":" + str(iatm)
597                                if var in varbs: continue
598                                varbs.append(var)
599                elif '=' in l[3]:
600                    for ph in phlist:
601                        key = G2obj.LookupPhaseName(ph)[0]
602                        cx,ct,cs,cia = Phases[key]['General']['AtomPtrs']
603                        for hst in hstlist: # should be blank
604                            atyp = l[3].split('=')[1]
605                            for iatm,at in enumerate(Phases[key]['Atoms']):
606                                if at[ct] != atyp: continue
607                                var = ph + ":" + hst + ":" + l[2] + ":" + str(iatm)
608                                if var in varbs: continue
609                                varbs.append(var)
610                else:
611                    for ph in phlist:
612                        key = G2obj.LookupPhaseName(ph)[0]
613                        for hst in hstlist: # should be blank
614                            var = ph + ":" + hst + ":" + l[2] + ":" + l[3]
615                            if var in varbs: continue
616                            varbs.append(var)
617        if len(varbs) >= 1 or 'constraint' in constType:
618            constr = [[1.0,FrstVarb]]
619            for item in varbs[1:]:
620                constr += [[1.0,G2obj.G2VarObj(item)]]
621            if 'equivalence' in constType:
622                return [constr+[None,None,'e']]
623            elif 'function' in constType:
624                return [constr+[None,False,'f']]
625            elif 'constraint' in constType:
626                return [constr+[1.0,None,'c']]
627            else:
628                raise Exception('Unknown constraint type: '+str(constType))
629        else:
630            dlg = wx.MessageDialog(
631                G2frame,
632                'There are no selected variables to include with '+str(FrstVarb),
633                'No variables')
634            dlg.CenterOnParent()
635            dlg.ShowModal()
636            dlg.Destroy()
637        return []
638   
639    def ConstraintsLoad(data,newcons=[]):
640        '''Load all constraints. Constraints based on symmetry (etc.)
641        are generated by running :func:`GSASIIstrIO.GetPhaseData`.
642        '''
643        G2mv.InitVars()
644        #Find all constraints
645        constraintSet = []
646        for key in data:
647            if key.startswith('_'): continue
648            constraintSet += data[key]
649        if newcons:
650            constraintSet = constraintSet + newcons
651        constDictList,fixedList,ignored = G2stIO.ProcessConstraints(constraintSet)
652        # generate symmetry constraints to check for conflicts
653        rigidbodyDict = G2frame.GPXtree.GetItemPyData(   
654            G2gd.GetGPXtreeItemId(G2frame,G2frame.root,'Rigid bodies'))
655        rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
656        rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,Print=False)
657        (Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,
658             BLtables,MFtables,maxSSwave) = G2stIO.GetPhaseData(
659                 Phases,RestraintDict=None,rbIds=rbIds,Print=False) # generates atom symmetry constraints
660        return constDictList,phaseDict,fixedList
661           
662    def ConstraintsCheck(data,newcons=[],reqVaryList=None):
663        '''Load constraints & check them for errors. Since error checking
664        can cause changes in constraints in case of repairable conflicts
665        between equivalences, reload the constraints again after the check.
666        This could probably be done more effectively, but only reloading when
667        needed, but the reload should be quick.
668       
669        When reqVaryList is included (see WarnConstraintLimit) then
670        parameters with limits are checked against constraints and a
671        warning is shown.
672        '''
673        constDictList,phaseDict,fixedList = ConstraintsLoad(data,newcons)
674        msg = G2mv.EvaluateMultipliers(constDictList,phaseDict)
675        if msg:
676            return 'Unable to interpret multiplier(s): '+msg,''
677        res = G2mv.CheckConstraints('',constDictList,fixedList)
678        # reload constraints in case any were merged in MoveConfEquiv
679        ConstraintsLoad(data,newcons)
680        impossible = []
681        if reqVaryList:
682            Controls = G2frame.GPXtree.GetItemPyData(G2gd.GetGPXtreeItemId(G2frame,G2frame.root, 'Controls'))
683            for key in ('parmMinDict','parmMaxDict','parmFrozen'):
684                if key not in Controls: Controls[key] = {}
685            varyList = reqVaryList[:]
686            try:
687                G2mv.GenerateConstraints(varyList,constDictList,fixedList,phaseDict)
688                G2mv.Map2Dict(phaseDict,varyList)
689                # check for limits on dependent vars
690                consVars = [i for i in reqVaryList if i not in varyList]
691                impossible = set(
692                    [str(i) for i in Controls['parmMinDict'] if i in consVars] + 
693                    [str(i) for i in Controls['parmMaxDict'] if i in consVars])
694            except G2mv.ConstraintException:
695                pass
696        if impossible:
697            msg = ''
698            for i in sorted(impossible):
699                if msg: msg += ', '
700                msg += i
701            msg =  ' &'.join(msg.rsplit(',',1))
702            msg = ('Note: limits on variable(s) '+msg+
703            ' will be ignored because they are constrained.')
704            G2G.G2MessageBox(G2frame,msg,'Limits ignored for constrained vars')
705        return res
706
707    def CheckAddedConstraint(newcons):
708        '''Check a new constraint that has just been input.
709        If there is an error display a message and discard the last entry
710
711        Since the varylist is not available, no warning messages
712        should be generated here
713
714        :returns: True if constraint should be added
715        '''
716       
717        errmsg,warnmsg = ConstraintsCheck(data,newcons)
718        if errmsg:
719            G2frame.ErrorDialog('Constraint Error',
720                'Error with newly added constraint:\n'+errmsg+
721                '\nIgnoring newly added constraint',parent=G2frame)
722            # reset error status
723            errmsg,warnmsg = ConstraintsCheck(data)
724            if errmsg:
725                print (errmsg)
726                print (G2mv.VarRemapShow([],True))
727            return False
728        elif warnmsg:
729            print ('Unexpected contraint warning:\n'+warnmsg)
730        return True
731
732    def WarnConstraintLimit():
733        '''Check if constraints reference variables with limits.
734        Displays a warning message, but does nothing
735        '''
736       
737        try:
738            parmDict,reqVaryList = G2frame.MakeLSParmDict()
739            errmsg,warnmsg = ConstraintsCheck(data,[],reqVaryList)
740        except:
741            print('Error retrieving parameters')           
742
743
744    def CheckChangedConstraint():
745        '''Check all constraints after an edit has been made.
746        If there is an error display a message and reject the change.
747
748        Since the varylist is not available, no warning messages
749        should be generated.
750       
751        :returns: True if the edit should be retained
752        '''
753        errmsg,warnmsg = ConstraintsCheck(data)
754        if errmsg:
755            G2frame.ErrorDialog('Constraint Error',
756                'Error after editing constraint:\n'+errmsg+
757                '\nDiscarding last constraint edit',parent=G2frame)
758            # reset error status
759            errmsg,warnmsg = ConstraintsCheck(data)
760            if errmsg:
761                print (errmsg)
762                print (G2mv.VarRemapShow([],True))
763            return False
764        elif warnmsg:
765            print ('Unexpected contraint warning:\n'+warnmsg)
766        return True
767             
768    def PageSelection(page):
769        'Decode page reference'
770        if page[1] == "phs":
771            vartype = "phase"
772            varList = G2obj.removeNonRefined(phaseList)  # remove any non-refinable prms from list
773            constrDictEnt = 'Phase'
774        elif page[1] == "hap":
775            vartype = "Histogram*Phase"
776            varList = G2obj.removeNonRefined(hapList)  # remove any non-refinable prms from list
777            constrDictEnt = 'HAP'
778        elif page[1] == "hst":
779            vartype = "Histogram"
780            varList = G2obj.removeNonRefined(histList)  # remove any non-refinable prms from list
781            constrDictEnt = 'Hist'
782        elif page[1] == "glb":
783            vartype = "Global"
784            varList = G2obj.removeNonRefined(globalList)   # remove any non-refinable prms from list
785
786            constrDictEnt = 'Global'
787        elif page[1] == "sym":
788            return None,None,None
789        else:
790            raise Exception('Should not happen!')
791        return vartype,varList,constrDictEnt
792
793    def OnAddHold(event):
794        '''Create a new Hold constraint
795
796        Hold constraints allows the user to select one variable (the list of available
797        variables depends on which tab is currently active).
798        '''
799        page = G2frame.Page
800        vartype,varList,constrDictEnt = PageSelection(page)
801        if vartype is None: return
802        varList = G2obj.SortVariables(varList)
803        title1 = "Hold "+vartype+" variable"
804        if not varList:
805            G2frame.ErrorDialog('No variables','There are no variables of type '+vartype,
806                parent=G2frame)
807            return
808        l2 = l1 = 1
809        for i in varList:
810            l1 = max(l1,len(i))
811            loc,desc = G2obj.VarDescr(i)
812            l2 = max(l2,len(loc))
813        fmt = "{:"+str(l1)+"s} {:"+str(l2)+"s} {:s}"
814        varListlbl = [fmt.format(i,*G2obj.VarDescr(i)) for i in varList]
815        #varListlbl = ["("+i+") "+G2obj.fmtVarDescr(i) for i in varList]
816        legend = "Select variables to hold (Will not be varied, even if vary flag is set)"
817        dlg = G2G.G2MultiChoiceDialog(G2frame,
818            legend,title1,varListlbl,toggle=False,size=(625,400),monoFont=True)
819        dlg.CenterOnParent()
820        if dlg.ShowModal() == wx.ID_OK:
821            for sel in dlg.GetSelections():
822                Varb = varList[sel]
823                VarObj = G2obj.G2VarObj(Varb)
824                newcons = [[[0.0,VarObj],None,None,'h']]
825                if CheckAddedConstraint(newcons):
826                    data[constrDictEnt] += newcons
827        dlg.Destroy()
828        wx.CallAfter(OnPageChanged,None)
829       
830    def OnAddEquivalence(event):
831        '''add an Equivalence constraint'''
832        page = G2frame.Page
833        vartype,varList,constrDictEnt = PageSelection(page)
834        if vartype is None: return
835        title1 = "Create equivalence constraint between "+vartype+" variables"
836        title2 = "Select additional "+vartype+" variable(s) to be equivalent with "
837        if not varList:
838            G2frame.ErrorDialog('No variables','There are no variables of type '+vartype,
839                parent=G2frame)
840            return
841#        legend = "Select variables to make equivalent (only one of the variables will be varied when all are set to be varied)"
842        GetAddVars(page,title1,title2,varList,constrDictEnt,'equivalence')
843       
844    def OnAddAtomEquiv(event):
845        ''' Add equivalences between all parameters on atoms '''
846        page = G2frame.Page
847        vartype,varList,constrDictEnt = PageSelection(page)
848        if vartype is None: return
849        title1 = "Setup equivalent atom variables"
850        title2 = "Select additional atoms(s) to be equivalent with "
851        if not varList:
852            G2frame.ErrorDialog('No variables','There are no variables of type '+vartype,
853                parent=G2frame)
854            return
855#        legend = "Select atoms to make equivalent (only one of the atom variables will be varied when all are set to be varied)"
856        GetAddAtomVars(page,title1,title2,varList,constrDictEnt,'equivalence')
857       
858    def OnAddRiding(event):
859        ''' Add riding equivalences between all parameters on atoms '''
860        page = G2frame.Page
861        vartype,varList,constrDictEnt = PageSelection(page)
862        if vartype is None: return
863        title1 = "Setup riding atoms "
864        title2 = "Select additional atoms(s) to ride on "
865        if not varList:
866            G2frame.ErrorDialog('No variables','There are no variables of type '+vartype,
867                parent=G2frame)
868            return
869#        legend = "Select atoms to ride (only one of the atom variables will be varied when all are set to be varied)"
870        GetAddAtomVars(page,title1,title2,varList,constrDictEnt,'riding')
871   
872    def OnAddFunction(event):
873        '''add a Function (new variable) constraint'''
874        page = G2frame.Page
875        vartype,varList,constrDictEnt = PageSelection(page)
876        if vartype is None: return
877        title1 = "Setup new variable based on "+vartype+" variables"
878        title2 = "Include additional "+vartype+" variable(s) to be included with "
879        if not varList:
880            G2frame.ErrorDialog('No variables','There are no variables of type '+vartype,
881                parent=G2frame)
882            return
883#        legend = "Select variables to include in a new variable (the new variable will be varied when all included variables are varied)"
884        GetAddVars(page,title1,title2,varList,constrDictEnt,'function')
885                       
886    def OnAddConstraint(event):
887        '''add a constraint equation to the constraints list'''
888        page = G2frame.Page
889        vartype,varList,constrDictEnt = PageSelection(page)
890        if vartype is None: return
891        title1 = "Creating constraint on "+vartype+" variables"
892        title2 = "Select additional "+vartype+" variable(s) to include in constraint with "
893        if not varList:
894            G2frame.ErrorDialog('No variables','There are no variables of type '+vartype,
895                parent=G2frame)
896            return
897#        legend = "Select variables to include in a constraint equation (the values will be constrainted to equal a specified constant)"
898        GetAddVars(page,title1,title2,varList,constrDictEnt,'constraint')
899
900    def GetAddVars(page,title1,title2,varList,constrDictEnt,constType):
901        '''Get the variables to be added for OnAddEquivalence, OnAddFunction,
902        and OnAddConstraint. Then create and check the constraint.
903        '''
904        #varListlbl = ["("+i+") "+G2obj.fmtVarDescr(i) for i in varList]
905        if constType == 'equivalence':
906            omitVars = G2mv.GetDependentVars()
907        else:
908            omitVars = []
909        varList = G2obj.SortVariables([i for i in varList if i not in omitVars])
910        l2 = l1 = 1
911        for i in varList:
912            l1 = max(l1,len(i))
913            loc,desc = G2obj.VarDescr(i)
914            l2 = max(l2,len(loc))
915        fmt = "{:"+str(l1)+"s} {:"+str(l2)+"s} {:s}"
916        varListlbl = [fmt.format(i,*G2obj.VarDescr(i)) for i in varList]
917        dlg = G2G.G2SingleChoiceDialog(G2frame,'Select 1st variable:',
918            title1,varListlbl,monoFont=True,size=(625,400))
919        dlg.CenterOnParent()
920        if dlg.ShowModal() == wx.ID_OK:
921            if constType == 'equivalence':
922                omitVars = G2mv.GetDependentVars() + G2mv.GetIndependentVars()
923            sel = dlg.GetSelection()
924            FrstVarb = varList[sel]
925            VarObj = G2obj.G2VarObj(FrstVarb)
926            moreVarb = G2obj.SortVariables(FindEquivVarb(FrstVarb,[i for i in varList if i not in omitVars]))
927            newcons = SelectVarbs(page,VarObj,moreVarb,title2+FrstVarb,constType)
928            if len(newcons) > 0:
929                if CheckAddedConstraint(newcons):
930                    data[constrDictEnt] += newcons
931        dlg.Destroy()
932        WarnConstraintLimit()
933        wx.CallAfter(OnPageChanged,None)
934                       
935    def FindNeighbors(phase,FrstName,AtNames):
936        General = phase['General']
937        cx,ct,cs,cia = General['AtomPtrs']
938        Atoms = phase['Atoms']
939        atNames = [atom[ct-1] for atom in Atoms]
940        Cell = General['Cell'][1:7]
941        Amat,Bmat = G2lat.cell2AB(Cell)
942        atTypes = General['AtomTypes']
943        Radii = np.array(General['BondRadii'])
944        AtInfo = dict(zip(atTypes,Radii)) #or General['BondRadii']
945        Orig = atNames.index(FrstName.split()[1])
946        OType = Atoms[Orig][ct]
947        XYZ = G2mth.getAtomXYZ(Atoms,cx)       
948        Neigh = []
949        Dx = np.inner(Amat,XYZ-XYZ[Orig]).T
950        dist = np.sqrt(np.sum(Dx**2,axis=1))
951        sumR = AtInfo[OType]+0.5    #H-atoms only!
952        IndB = ma.nonzero(ma.masked_greater(dist-0.85*sumR,0.))
953        for j in IndB[0]:
954            if j != Orig:
955                Neigh.append(AtNames[j])
956        return Neigh
957       
958    def GetAddAtomVars(page,title1,title2,varList,constrDictEnt,constType):
959        '''Get the atom variables to be added for OnAddAtomEquiv. Then create and
960        check the constraints. Riding for H atoms only.
961        '''
962        Atoms = {G2obj.VarDescr(i)[0]:[] for i in varList if 'Atom' in G2obj.VarDescr(i)[0]}
963        for item in varList:
964            atName = G2obj.VarDescr(item)[0]
965            if atName in Atoms:
966                Atoms[atName].append(item)
967        AtNames = list(Atoms.keys())
968        AtNames.sort()
969        dlg = G2G.G2SingleChoiceDialog(G2frame,'Select 1st atom:',
970            title1,AtNames,monoFont=True,size=(625,400))
971        dlg.CenterOnParent()
972        FrstAtom = ''
973        if dlg.ShowModal() == wx.ID_OK:
974            sel = dlg.GetSelection()
975            FrstAtom = AtNames[sel]
976            if 'riding' in constType:
977                phaseName = (FrstAtom.split(' in ')[1]).strip()
978                phase = Phases[phaseName]
979                AtNames = FindNeighbors(phase,FrstAtom,AtNames)
980            else:
981                AtNames.remove(FrstAtom)
982        dlg.Destroy()
983        if FrstAtom == '':
984            print ('no atom selected')
985            return
986        dlg = G2G.G2MultiChoiceDialog(
987            G2frame,title2+FrstAtom,
988            'Constrain '+str(FrstAtom)+' with...',AtNames,
989            toggle=False,size=(625,400),monoFont=True)
990        if dlg.ShowModal() == wx.ID_OK:
991            Selections = dlg.GetSelections()[:]
992        else:
993            print ('no target atom selected')
994            dlg.Destroy()
995            return
996        dlg.Destroy()
997        for name in Atoms[FrstAtom]:
998            newcons = []
999            constr = []
1000            if 'riding' in constType:
1001                if 'AUiso' in name:
1002                    constr = [[1.0,G2obj.G2VarObj(name)]]
1003                elif 'AU11' in name:
1004                    pass
1005                elif 'AU' not in name:
1006                    constr = [[1.0,G2obj.G2VarObj(name)]]
1007            else:
1008                constr = [[1.0,G2obj.G2VarObj(name)]]
1009            pref = ':'+name.rsplit(':',1)[0].split(':',1)[1]    #get stuff between phase id & atom id
1010            for sel in Selections:
1011                name2 = Atoms[AtNames[sel]][0]
1012                pid = name2.split(':',1)[0]                     #get phase id for 2nd atom
1013                id = name2.rsplit(':',1)[-1]                    #get atom no. for 2nd atom
1014                if 'riding' in constType:
1015                    pref = pid+pref
1016                    if 'AUiso' in pref:
1017                        parts = pref.split('AUiso')
1018                        constr += [[1.2,G2obj.G2VarObj('%s:%s'%(parts[0]+'AUiso',id))]]
1019                    elif 'AU' not in pref:
1020                        constr += [[1.0,G2obj.G2VarObj('%s:%s'%(pref,id))]]
1021                else:
1022                    constr += [[1.0,G2obj.G2VarObj('%s:%s'%(pid+pref,id))]]
1023            if not constr:
1024                continue
1025            if 'frac' in pref and 'riding' not in constType:
1026                newcons = [constr+[1.0,None,'c']]
1027            else:
1028                newcons = [constr+[None,None,'e']]
1029            if len(newcons) > 0:
1030                if CheckAddedConstraint(newcons):
1031                    data[constrDictEnt] += newcons
1032        WarnConstraintLimit()
1033        wx.CallAfter(OnPageChanged,None)
1034                       
1035    def MakeConstraintsSizer(name,pageDisplay):
1036        '''Creates a sizer displaying all of the constraints entered of
1037        the specified type.
1038
1039        :param str name: the type of constraints to be displayed ('HAP',
1040          'Hist', 'Phase', 'Global', 'Sym-Generated')
1041        :param wx.Panel pageDisplay: parent panel for sizer
1042        :returns: wx.Sizer created by method
1043        '''
1044        if name == 'Sym-Generated':         #show symmetry generated constraints
1045            Sizer1 =  wx.BoxSizer(wx.VERTICAL)
1046            if symHolds:
1047                Sizer1.Add(wx.StaticText(pageDisplay,wx.ID_ANY,
1048                    'Position variables fixed by space group symmetry'))
1049                Sizer1.Add((-1,5))
1050                Sizer = wx.FlexGridSizer(0,2,0,0)
1051                Sizer1.Add(Sizer)
1052                for var in symHolds:
1053                    Sizer.Add(wx.StaticText(pageDisplay,wx.ID_ANY,'  FIXED'),
1054                              0,WACV|wx.ALIGN_CENTER|wx.RIGHT|wx.LEFT,3)
1055                    Sizer.Add(wx.StaticText(pageDisplay,wx.ID_ANY,var))
1056                    Sizer.Add((-1,-1))
1057                    Sizer.Add((-1,2))
1058            else:
1059                Sizer1.Add(wx.StaticText(pageDisplay,wx.ID_ANY,
1060                    'No holds generated'))
1061            Sizer1.Add((-1,10))
1062            symGen = G2mv.GetSymEquiv()
1063            if len(symGen) == 0:
1064                Sizer1.Add(wx.StaticText(pageDisplay,wx.ID_ANY,
1065                    'No equvalences generated'))
1066                return Sizer1
1067            Sizer1.Add(wx.StaticText(pageDisplay,wx.ID_ANY,
1068                'Equivalences generated based on cell/space group input'))
1069            Sizer1.Add((-1,5))
1070            Sizer = wx.FlexGridSizer(0,2,0,0)
1071            Sizer1.Add(Sizer)
1072            for sym in symGen:
1073                Sizer.Add(wx.StaticText(pageDisplay,wx.ID_ANY,'  EQUIV'),
1074                    0,WACV|wx.ALIGN_CENTER|wx.RIGHT|wx.LEFT,3)
1075                Sizer.Add(wx.StaticText(pageDisplay,wx.ID_ANY,sym))
1076                Sizer.Add((-1,-1))
1077                Sizer.Add((-1,2))
1078            return Sizer1
1079        constSizer = wx.FlexGridSizer(0,6,0,0)
1080        maxlen = 50 # characters before wrapping a constraint
1081        for Id,item in enumerate(data[name]):
1082            refineflag = False
1083            helptext = ""
1084            eqString = ['',]
1085            problemItem = False
1086            badVar = False
1087            for term in item[:-3]:
1088                if str(term[1]) in G2mv.problemVars:
1089                    problemItem = True
1090            if item[-1] == 'h': # Hold on variable
1091                constSizer.Add((-1,-1),0)              # blank space for edit button
1092                typeString = 'FIXED'
1093                var = str(item[0][1])
1094                if '?' in var: badVar = True
1095                varMean = G2obj.fmtVarDescr(var)
1096                eqString[-1] =  var +'   '
1097                helptext = "Prevents variable:\n"+ var + " ("+ varMean + ")\nfrom being changed"
1098            elif item[-1] == 'f' or item[-1] == 'e' or item[-1] == 'c': # not true on original-style (2011?) constraints
1099                constEdit = wx.Button(pageDisplay,wx.ID_ANY,'Edit',style=wx.BU_EXACTFIT)
1100                constEdit.Bind(wx.EVT_BUTTON,OnConstEdit)
1101                constSizer.Add(constEdit,0,wx.LEFT|wx.RIGHT|WACV|wx.ALIGN_CENTER,1)            # edit button
1102                Indx[constEdit.GetId()] = [Id,name]
1103                if item[-1] == 'f':
1104                    helptext = "A new variable"
1105                    if item[-3]:
1106                        helptext += " named "+str(item[-3])
1107                    helptext += " is created from a linear combination of the following variables:\n"
1108                    for term in item[:-3]:
1109                        var = str(term[1])
1110                        if '?' in var: badVar = True
1111                        if len(eqString[-1]) > maxlen:
1112                            eqString.append(' ')
1113                        m = term[0]
1114                        if eqString[-1] != '':
1115                            if m >= 0:
1116                                eqString[-1] += ' + '
1117                            else:
1118                                eqString[-1] += ' - '
1119                                m = abs(m)
1120                        if m == 1:
1121                            eqString[-1] += '{:} '.format(var)
1122                        else:
1123                            eqString[-1] += '{:g}*{:} '.format(m,var)
1124                        varMean = G2obj.fmtVarDescr(var)
1125                        helptext += '\n{:3g} * {:} '.format(m,var) + " ("+ varMean + ")"
1126                    # Add ISODISTORT help items
1127                    if '_Explain' in data:
1128                        # this ignores the phase number. TODO: refactor that in
1129                        hlptxt = None
1130                        try:
1131                            hlptxt = data['_Explain'].get(item[-3])
1132                        except TypeError:
1133                            # note that phase RanId is item[-3].phase
1134                            hlptxt = data['_Explain'].get(str(item[-3].phase)+item[-3].name)
1135                        if hlptxt:
1136                            helptext += '\n\n'
1137                            helptext += hlptxt
1138                    if item[-3]:
1139                        typeString = '(NEWVAR) ' + str(item[-3]) + ' = '
1140                    else:
1141                        typeString = 'New Variable = '
1142                    #print 'refine',item[-2]
1143                    refineflag = True
1144                elif item[-1] == 'c':
1145                    helptext = "The following variables constrained to equal a constant:"
1146                    for term in item[:-3]:
1147                        var = str(term[1])
1148                        if '?' in var: badVar = True
1149                        if len(eqString[-1]) > maxlen:
1150                            eqString.append(' ')
1151                        m = term[0]
1152                        if eqString[-1] != '':
1153                            if term[0] > 0:
1154                                eqString[-1] += ' + '
1155                            else:
1156                                eqString[-1] += ' - '
1157                                m = -term[0]
1158                        if m == 1:
1159                            eqString[-1] += '{:} '.format(var)
1160                        else:
1161                            eqString[-1] += '{:g}*{:} '.format(m,var)
1162                        varMean = G2obj.fmtVarDescr(var)
1163                        helptext += '\n{:3g} * {:} '.format(m,var) + " ("+ varMean + ")"
1164                    typeString = 'CONST'
1165                    eqString[-1] += ' = '+str(item[-3])
1166                elif item[-1] == 'e' and len(item[:-3]) == 2:
1167                    if item[0][0] == 0: item[0][0] = 1.0
1168                    if item[1][0] == 0: item[1][0] = 1.0
1169                    var = str(item[0][1])
1170                    if '?' in var: badVar = True
1171                    helptext = 'Variable {:} '.format(var) + " ("+ G2obj.fmtVarDescr(var) + ")"
1172                    helptext += "\n\nis equivalent to "
1173                    m = item[0][0]/item[1][0]
1174                    var1 = str(item[1][1])
1175                    helptext += '\n{:3g} * {:} '.format(m,var1) + " ("+ G2obj.fmtVarDescr(var1) + ")"
1176                    eqString[-1] += '{:} = {:}'.format(var1,var)
1177                    if m != 1:
1178                        eqString[-1] += ' / ' + str(m)
1179                    typeString = 'EQUIV'
1180                elif item[-1] == 'e':
1181                    helptext = "The following variable:"
1182                    normval = item[0][0]
1183                    indepterm = item[0][1]
1184                    for i,term in enumerate(item[:-3]):
1185                        var = str(term[1])
1186                        if '?' in var: badVar = True
1187                        if term[0] == 0: term[0] = 1.0
1188                        if len(eqString[-1]) > maxlen:
1189                            eqString.append(' ')
1190                        varMean = G2obj.fmtVarDescr(var)
1191                        if i == 0: # move independent variable to end, as requested by Bob
1192                            helptext += '\n{:} '.format(var) + " ("+ varMean + ")"
1193                            helptext += "\n\nis equivalent to the following, noting multipliers:"
1194                            continue
1195                        elif eqString[-1] != '':
1196                            eqString[-1] += ' = '
1197                        #m = normval/term[0]
1198                        m = term[0]/normval
1199                        if m == 1:
1200                            eqString[-1] += '{:}'.format(var)
1201                        else:
1202                            eqString[-1] += '{:g}*{:} '.format(m,var)
1203                        helptext += '\n{:3g} * {:} '.format(m,var) + " ("+ varMean + ")"
1204                    eqString[-1] += ' = {:} '.format(indepterm)
1205                    typeString = 'EQUIV'
1206                else:
1207                    print ('Unexpected constraint'+item)
1208               
1209            else:
1210                print ('Removing old-style constraints')
1211                data[name] = []
1212                return constSizer
1213            constDel = wx.Button(pageDisplay,wx.ID_ANY,'Delete',style=wx.BU_EXACTFIT)
1214            constDel.Bind(wx.EVT_BUTTON,OnConstDel)
1215            Indx[constDel.GetId()] = [Id,name]
1216            constSizer.Add(constDel,0,wx.LEFT|wx.RIGHT|WACV|wx.ALIGN_CENTER,1)             # delete button
1217            if helptext:
1218                ch = G2G.HelpButton(pageDisplay,helptext)
1219                constSizer.Add(ch,0,wx.LEFT|wx.RIGHT|WACV|wx.ALIGN_CENTER,1)
1220            else:
1221                constSizer.Add((-1,-1))
1222            if refineflag:
1223                ch = G2G.G2CheckBox(pageDisplay,'',item,-2)
1224                constSizer.Add(ch,0,wx.LEFT|wx.RIGHT|WACV|wx.ALIGN_CENTER,1)
1225            else:
1226                constSizer.Add((-1,-1))               
1227            constSizer.Add(wx.StaticText(pageDisplay,wx.ID_ANY,typeString),
1228                           0,WACV|wx.ALIGN_CENTER|wx.RIGHT|wx.LEFT,3)
1229            if badVar: eqString[-1] += ' -- Error: variable removed'
1230            if problemItem: eqString[-1] += ' -- Conflict: see console'
1231            if len(eqString) > 1:
1232                Eq = wx.BoxSizer(wx.VERTICAL)
1233                for s in eqString:
1234                    line = wx.StaticText(pageDisplay,wx.ID_ANY,s)
1235                    if problemItem: line.SetBackgroundColour(wx.YELLOW)
1236                    Eq.Add(line,0)
1237                Eq.Add((-1,4))
1238            else:
1239                Eq = wx.StaticText(pageDisplay,wx.ID_ANY,eqString[0])
1240                if problemItem or badVar: Eq.SetBackgroundColour(wx.YELLOW)
1241            constSizer.Add(Eq,1,WACV)
1242        return constSizer
1243               
1244    def OnConstDel(event):
1245        'Delete a constraint'
1246        Obj = event.GetEventObject()
1247        Id,name = Indx[Obj.GetId()]
1248        del data[name][Id]
1249        ConstraintsLoad(data)
1250        wx.CallAfter(OnPageChanged,None)
1251
1252    def OnConstEdit(event):
1253        '''Called to edit an individual contraint in response to a
1254        click on its Edit button
1255        '''
1256        Obj = event.GetEventObject()
1257        Id,name = Indx[Obj.GetId()]
1258        if data[name][Id][-1] == 'f':
1259            items = data[name][Id][:-3]
1260            constType = 'New Variable'
1261            if data[name][Id][-3]:
1262                varname = str(data[name][Id][-3])
1263            else:
1264                varname = ""
1265            lbl = 'Enter value for each term in constraint; sum = new variable'
1266            dlg = ConstraintDialog(G2frame,constType,lbl,items,
1267                                   varname=varname,varyflag=data[name][Id][-2])
1268        elif data[name][Id][-1] == 'c':
1269            items = data[name][Id][:-3]+[
1270                [str(data[name][Id][-3]),'fixed value =']]
1271            constType = 'Constraint'
1272            lbl = 'Edit value for each term in constant constraint sum'
1273            dlg = ConstraintDialog(G2frame,constType,lbl,items)
1274        elif data[name][Id][-1] == 'e':
1275            items = data[name][Id][:-3]
1276            constType = 'Equivalence'
1277            lbl = 'The following terms are set to be equal:'
1278            dlg = ConstraintDialog(G2frame,constType,lbl,items,'*')
1279        else:
1280            return
1281        try:
1282            prev = copy.deepcopy(data[name][Id])
1283            if dlg.ShowModal() == wx.ID_OK:
1284                result = dlg.GetData()
1285                for i in range(len(data[name][Id][:-3])):
1286                    if type(data[name][Id][i]) is tuple: # fix non-mutable construct
1287                        data[name][Id][i] = list(data[name][Id][i])
1288                    data[name][Id][i][0] = result[i][0]
1289                if data[name][Id][-1] == 'c':
1290                    data[name][Id][-3] = str(result[-1][0])
1291                elif data[name][Id][-1] == 'f':
1292                    data[name][Id][-2] = dlg.newvar[1]
1293                    if dlg.newvar[0]:
1294                        # process the variable name to put in global form (::var)
1295                        varname = str(dlg.newvar[0]).strip().replace(' ','_')
1296                        if varname.startswith('::'):
1297                            varname = varname[2:]
1298                        varname = varname.replace(':',';')
1299                        if varname:
1300                            data[name][Id][-3] = varname
1301                        else:
1302                            data[name][Id][-3] = ''
1303                if not CheckChangedConstraint():
1304                    data[name][Id] = prev
1305            else:
1306                data[name][Id] = prev
1307        except:
1308            import traceback
1309            print (traceback.format_exc())
1310        finally:
1311            dlg.Destroy()
1312        wx.CallAfter(OnPageChanged,None)
1313   
1314    def UpdateConstraintPanel(panel,typ):
1315        '''Update the contents of the selected Constraint
1316        notebook tab. Called in :func:`OnPageChanged`
1317        '''
1318        if panel.GetSizer(): panel.GetSizer().Clear(True)
1319        Siz = wx.BoxSizer(wx.VERTICAL)
1320        Siz.Add((5,5),0)
1321        Siz.Add(MakeConstraintsSizer(typ,panel),1,wx.EXPAND)
1322        panel.SetSizer(Siz,True)
1323        Size = Siz.GetMinSize()
1324        Size[0] += 40
1325        Size[1] = max(Size[1],450) + 20
1326        panel.SetSize(Size)
1327        panel.SetScrollbars(10,10,Size[0]/10-4,Size[1]/10-1)
1328        panel.Show()
1329
1330    def OnPageChanged(event):
1331        '''Called when a tab is pressed or when a "select tab" menu button is
1332        used (see RaisePage), or to refresh the current tab contents (event=None)
1333        '''
1334        if event:       #page change event!
1335            page = event.GetSelection()
1336        else: # called directly, get current page
1337            try:
1338                page = G2frame.constr.GetSelection()
1339            except:
1340                if GSASIIpath.GetConfigValue('debug'): print('DBG_gpx open error:C++ Run time error - skipped')
1341                return
1342        G2frame.constr.ChangeSelection(page)
1343        text = G2frame.constr.GetPageText(page)
1344        G2frame.dataWindow.ConstraintEdit.Enable(G2G.wxID_EQUIVALANCEATOMS,False)
1345#        G2frame.dataWindow.ConstraintEdit.Enable(G2G.wxID_ADDRIDING,False)
1346        if text == 'Histogram/Phase':
1347            enableEditCons = [False]+4*[True]
1348            G2frame.Page = [page,'hap']
1349            UpdateConstraintPanel(HAPConstr,'HAP')
1350        elif text == 'Histogram':
1351            enableEditCons = [False]+4*[True]
1352            G2frame.Page = [page,'hst']
1353            UpdateConstraintPanel(HistConstr,'Hist')
1354        elif text == 'Phase':
1355            enableEditCons = 5*[True]
1356            G2frame.Page = [page,'phs']
1357            G2frame.dataWindow.ConstraintEdit.Enable(G2G.wxID_EQUIVALANCEATOMS,True)
1358#            G2frame.dataWindow.ConstraintEdit.Enable(G2G.wxID_ADDRIDING,True)
1359            if 'DELETED' in str(PhaseConstr):   #seems to be no other way to do this (wx bug)
1360                if GSASIIpath.GetConfigValue('debug'):
1361                    print ('DBG_wx error: PhaseConstr not cleanly deleted after Refine')
1362                return
1363            UpdateConstraintPanel(PhaseConstr,'Phase')
1364        elif text == 'Global':
1365            enableEditCons = [False]+4*[True]
1366            G2frame.Page = [page,'glb']
1367            UpdateConstraintPanel(GlobalConstr,'Global')
1368        else:
1369            enableEditCons = 5*[False]
1370            G2frame.Page = [page,'sym']
1371            UpdateConstraintPanel(SymConstr,'Sym-Generated')
1372        # remove menu items when not allowed
1373        for obj,flag in zip(G2frame.dataWindow.ConstraintEdit.GetMenuItems(),enableEditCons): 
1374            obj.Enable(flag)
1375        G2frame.dataWindow.SetDataSize()
1376
1377    def RaisePage(event):
1378        'Respond to a "select tab" menu button'
1379        try:
1380            i = (G2G.wxID_CONSPHASE,
1381                 G2G.wxID_CONSHAP,
1382                 G2G.wxID_CONSHIST,
1383                 G2G.wxID_CONSGLOBAL,
1384                 G2G.wxID_CONSSYM,
1385                ).index(event.GetId())
1386            G2frame.constr.SetSelection(i)
1387            wx.CallAfter(OnPageChanged,None)
1388        except ValueError:
1389            print('Unexpected event in RaisePage')
1390
1391    def SetStatusLine(text):
1392        G2frame.GetStatusBar().SetStatusText(text,1)
1393       
1394    def OnShowISODISTORT(event):
1395        ShowIsoDistortCalc(G2frame)
1396
1397    #### UpdateConstraints execution starts here ##############################
1398    if not data:
1399        data.update({'Hist':[],'HAP':[],'Phase':[],'Global':[]})       #empty dict - fill it
1400    if 'Global' not in data:                                            #patch
1401        data['Global'] = []
1402    # DEBUG code #=====================================
1403    #import GSASIIconstrGUI
1404    #reload(GSASIIconstrGUI)
1405    #reload(G2obj)
1406    #reload(G2stIO)
1407    #import GSASIIstrMain
1408    #reload(GSASIIstrMain)   
1409    #reload(G2mv)
1410    #reload(G2gd)
1411    #===================================================
1412    Histograms,Phases = G2frame.GetUsedHistogramsAndPhasesfromTree()
1413    if not len(Phases) or not len(Histograms):
1414        dlg = wx.MessageDialog(G2frame,'You need both phases and histograms to see Constraints',
1415            'No phases or histograms')
1416        dlg.CenterOnParent()
1417        dlg.ShowModal()
1418        dlg.Destroy()
1419        return
1420    G2obj.IndexAllIds(Histograms,Phases)
1421    for p in Phases:
1422        if 'ISODISTORT' in Phases[p]:
1423            G2frame.dataWindow.ConstraintEdit.Enable(G2G.wxID_SHOWISO,True)
1424            break
1425    ###### patch: convert old-style (str) variables in constraints to G2VarObj objects #####
1426    for key,value in data.items():
1427        if key.startswith('_'): continue
1428        j = 0
1429        for cons in value:
1430            #print cons             # DEBUG
1431            for i in range(len(cons[:-3])):
1432                if type(cons[i][1]) is str:
1433                    cons[i][1] = G2obj.G2VarObj(cons[i][1])
1434                    j += 1
1435        if j:
1436            print (str(key) + ': '+str(j)+' variable(s) as strings converted to objects')
1437    ##### end patch #############################################################################
1438    rigidbodyDict = G2frame.GPXtree.GetItemPyData(
1439        G2gd.GetGPXtreeItemId(G2frame,G2frame.root,'Rigid bodies'))
1440    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
1441    rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,Print=False)
1442    badPhaseParms = ['Ax','Ay','Az','Amul','AI/A','Atype','SHorder','AwaveType','FwaveType','PwaveType','MwaveType','Vol','isMag',]
1443    globalList = list(rbDict.keys())
1444    globalList.sort()
1445    try:
1446        AtomDict = dict([Phases[phase]['pId'],Phases[phase]['Atoms']] for phase in Phases)
1447    except KeyError:
1448        G2frame.ErrorDialog('Constraint Error','Constraints cannot be set until a cycle of least squares'+
1449                            ' has been run.\nWe suggest you refine a scale factor.')
1450        return
1451
1452    # create a list of the phase variables
1453    symHolds = []
1454    (Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtable,BLtable,
1455         MFtable,maxSSwave) = G2stIO.GetPhaseData(
1456             Phases,rbIds=rbIds,Print=False,symHold=symHolds)
1457    phaseList = []
1458    for item in phaseDict:
1459        if item.split(':')[2] not in badPhaseParms:
1460            phaseList.append(item)
1461    phaseList.sort()
1462    phaseAtNames = {}
1463    phaseAtTypes = {}
1464    TypeList = []
1465    for item in phaseList:
1466        Split = item.split(':')
1467        if Split[2][:2] in ['AU','Af','dA','AM']:
1468            Id = int(Split[0])
1469            phaseAtNames[item] = AtomDict[Id][int(Split[3])][0]
1470            phaseAtTypes[item] = AtomDict[Id][int(Split[3])][1]
1471            if phaseAtTypes[item] not in TypeList:
1472                TypeList.append(phaseAtTypes[item])
1473        else:
1474            phaseAtNames[item] = ''
1475            phaseAtTypes[item] = ''
1476             
1477    # create a list of the hist*phase variables
1478    seqList = G2frame.testSeqRefineMode()
1479    if seqList: # for sequential refinement, only process 1st histgram in list
1480        histDict = {seqList[0]:Histograms[seqList[0]]}
1481    else:
1482        histDict = Histograms
1483    hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,histDict,Print=False,resetRefList=False)
1484    hapList = sorted([i for i in hapDict.keys() if i.split(':')[2] not in ('Type',)])
1485    if seqList: # convert histogram # to wildcard
1486        wildList = [] # list of variables with "*" for histogram number
1487        for i in hapList:
1488            s = i.split(':')
1489            if s[1] == "": continue
1490            s[1] = '*'
1491            sj = ':'.join(s)
1492            if sj not in wildList: wildList.append(sj)
1493        hapList = wildList
1494    histVary,histDict,controlDict = G2stIO.GetHistogramData(histDict,Print=False)
1495    histList = list(histDict.keys())
1496    histList.sort()
1497    if seqList: # convert histogram # to wildcard
1498        wildList = [] # list of variables with "*" for histogram number
1499        for i in histList:
1500            s = i.split(':')
1501            if s[1] == "": continue
1502            s[1] = '*'
1503            sj = ':'.join(s)
1504            if sj not in wildList: wildList.append(sj)
1505        histList = wildList       
1506    Indx = {}
1507    G2frame.Page = [0,'phs']
1508   
1509    G2gd.SetDataMenuBar(G2frame,G2frame.dataWindow.ConstraintMenu)
1510    SetStatusLine('')
1511   
1512    G2frame.Bind(wx.EVT_MENU, OnAddConstraint, id=G2G.wxID_CONSTRAINTADD)
1513    G2frame.Bind(wx.EVT_MENU, OnAddFunction, id=G2G.wxID_FUNCTADD)
1514    G2frame.Bind(wx.EVT_MENU, OnAddEquivalence, id=G2G.wxID_EQUIVADD)
1515    G2frame.Bind(wx.EVT_MENU, OnAddHold, id=G2G.wxID_HOLDADD)
1516    G2frame.Bind(wx.EVT_MENU, OnAddAtomEquiv, id=G2G.wxID_EQUIVALANCEATOMS)
1517#    G2frame.Bind(wx.EVT_MENU, OnAddRiding, id=G2G.wxID_ADDRIDING)
1518    G2frame.Bind(wx.EVT_MENU, OnShowISODISTORT, id=G2G.wxID_SHOWISO)
1519    # tab commands
1520    for id in (G2G.wxID_CONSPHASE,
1521               G2G.wxID_CONSHAP,
1522               G2G.wxID_CONSHIST,
1523               G2G.wxID_CONSGLOBAL,
1524               G2G.wxID_CONSSYM,
1525               ):
1526        G2frame.Bind(wx.EVT_MENU, RaisePage,id=id)
1527
1528    #G2frame.constr = G2G.GSNoteBook(parent=G2frame.dataWindow,size=G2frame.dataWindow.GetClientSize())
1529    G2frame.constr = G2G.GSNoteBook(parent=G2frame.dataWindow)
1530    G2frame.dataWindow.GetSizer().Add(G2frame.constr,1,wx.ALL|wx.EXPAND)
1531    # note that order of pages is hard-coded in RaisePage
1532    PhaseConstr = wx.ScrolledWindow(G2frame.constr)
1533    G2frame.constr.AddPage(PhaseConstr,'Phase')
1534    HAPConstr = wx.ScrolledWindow(G2frame.constr)
1535    G2frame.constr.AddPage(HAPConstr,'Histogram/Phase')
1536    HistConstr = wx.ScrolledWindow(G2frame.constr)
1537    G2frame.constr.AddPage(HistConstr,'Histogram')
1538    GlobalConstr = wx.ScrolledWindow(G2frame.constr)
1539    G2frame.constr.AddPage(GlobalConstr,'Global')
1540    SymConstr = wx.ScrolledWindow(G2frame.constr)
1541    G2frame.constr.AddPage(SymConstr,'Sym-Generated')
1542    wx.CallAfter(OnPageChanged,None)
1543    G2frame.constr.Bind(wx.aui.EVT_AUINOTEBOOK_PAGE_CHANGED, OnPageChanged)
1544    # validate all the constrants -- should not see any errors here normally
1545    errmsg,warnmsg = ConstraintsCheck(data)
1546    if errmsg:
1547        G2frame.ErrorDialog('Constraint Error',
1548                            'Error in constraints:\n'+errmsg+'\nCheck console output for more information',
1549                            parent=G2frame)
1550        print (errmsg)
1551        print (G2mv.VarRemapShow([],True))
1552        return
1553    elif warnmsg:
1554        print ('Unexpected contraint warning:\n'+warnmsg)
1555    WarnConstraintLimit()
1556       
1557###### check scale & phase fractions, create constraint if needed #############
1558def CheckAllScalePhaseFractions(G2frame):
1559    '''Check if scale factor and all phase fractions are refined without a constraint
1560    for all used histograms, if so, offer the user a chance to create a constraint
1561    on the sum of phase fractions
1562    '''
1563    histograms, phases = G2frame.GetUsedHistogramsAndPhasesfromTree()
1564    for i,hist in enumerate(histograms):
1565        CheckScalePhaseFractions(G2frame,hist,histograms,phases)
1566       
1567def CheckScalePhaseFractions(G2frame,hist,histograms,phases):
1568    '''Check if scale factor and all phase fractions are refined without a constraint
1569    for histogram hist, if so, offer the user a chance to create a constraint
1570    on the sum of phase fractions
1571    '''
1572    if G2frame.testSeqRefineMode():
1573        histStr = '*'
1574    else:
1575        histStr = str(histograms[hist]['hId'])
1576    # Is this powder?
1577    if not hist.startswith('PWDR '): return
1578    # do this only if the scale factor is varied
1579    if not histograms[hist]['Sample Parameters']['Scale'][1]: return
1580    # are all phase fractions varied in all used histograms?
1581    phaseCount = 0
1582    for p in phases:
1583        if hist not in phases[p]['Histograms']: continue
1584        if phases[p]['Histograms'][hist]['Use'] and not phases[p]['Histograms'][hist]['Scale'][1]:
1585            return
1586        else:
1587            phaseCount += 1
1588   
1589    # all phase fractions and scale factor varied, now scan through constraints
1590    sub = G2gd.GetGPXtreeItemId(G2frame,G2frame.root,'Constraints') 
1591    Constraints = G2frame.GPXtree.GetItemPyData(sub)
1592    for c in Constraints.get('HAP',[]):
1593        if c[-1] != 'c': continue
1594        if not c[-3]: continue
1595        if len(c[:-3]) != phaseCount: continue
1596        # got a constraint equation with right number of terms, is it on phase fractions for
1597        # the correct histogram?
1598        if all([(i[1].name == 'Scale' and i[1].varname().split(':')[1] == histStr) for i in c[:-3]]):
1599            # got a constraint, this is OK
1600            return
1601    dlg = wx.MessageDialog(G2frame,'You are refining the scale factor and all phase fractions for histogram #'+
1602        histStr+'. This will produce an unstable refinement. '+
1603        'Do you want to constrain the sum of phase fractions?','Create constraint?',wx.OK|wx.CANCEL)
1604    if dlg.ShowModal() != wx.ID_OK:
1605        dlg.Destroy()
1606        return
1607    dlg.Destroy()
1608
1609    constr = []
1610    for p in phases:
1611        if hist not in phases[p]['Histograms']: continue
1612        if not phases[p]['Histograms'][hist]['Use']: continue
1613        constr += [[1.0,G2obj.G2VarObj(':'.join((str(phases[p]['pId']),histStr,'Scale')))]]
1614    constr += [1.0,None,'c']
1615    Constraints['HAP'] += [constr]
1616       
1617#### Make nuclear/magnetic phase transition constraints - called by OnTransform in G2phsGUI ##########
1618def TransConstraints(G2frame,oldPhase,newPhase,Trans,Vec,atCodes):
1619    '''Add constraints for new magnetic phase created via transformation of old
1620    nuclear one
1621    NB: A = [G11,G22,G33,2*G12,2*G13,2*G23]
1622    '''
1623   
1624    def SetUniqAj(pId,iA,SGData):
1625        SGLaue = SGData['SGLaue']
1626        SGUniq = SGData['SGUniq']
1627        if SGLaue in ['m3','m3m']:
1628            if iA in [0,1,2]:
1629                parm = '%d::%s'%(pId,'A0')
1630            else:
1631                parm = None
1632        elif SGLaue in ['4/m','4/mmm']:
1633            if iA in [0,1]:
1634                parm = '%d::%s'%(pId,'A0')
1635            elif iA == 2:
1636                parm = '%d::%s'%(pId,'A2')
1637            else:
1638                parm = None
1639        elif SGLaue in ['6/m','6/mmm','3m1', '31m', '3']:
1640            if iA in [0,1,3]:
1641                parm = '%d::%s'%(pId,'A0')
1642            elif iA == 2:
1643                parm = '%d::%s'%(pId,'A2')
1644            else:
1645                parm = None
1646        elif SGLaue in ['3R', '3mR']:
1647            if ia in [0,1,2]:
1648                parm = '%d::%s'%(pId,'A0')
1649            else:
1650                parm = '%d::%s'%(pId,'A3')
1651        elif SGLaue in ['mmm',]:
1652            if iA in [0,1,2]:
1653                parm = '%d::A%s'%(pId,iA)
1654            else:
1655                parm = None
1656        elif SGLaue == '2/m':
1657            if iA in [0,1,2]:
1658                parm = '%d::A%s'%(pId,iA)
1659            elif iA == 3 and SGUniq == 'c':
1660                parm = '%d::A%s'%(pId,iA)
1661            elif iA == 4 and SGUniq == 'b':
1662                parm = '%d::A%s'%(pId,iA)
1663            elif iA == 5 and SGUniq == 'a':
1664                parm = '%d::A%s'%(pId,iA)
1665            else:
1666                parm = None           
1667        else:
1668            parm = '%d::A%s'%(pId,iA)
1669        return parm
1670   
1671    Histograms,Phases = G2frame.GetUsedHistogramsAndPhasesfromTree()
1672    UseList = newPhase['Histograms']
1673    detTrans = np.abs(nl.det(Trans))
1674    opId = oldPhase['pId']
1675    npId = newPhase['pId']
1676    cx,ct,cs,cia = newPhase['General']['AtomPtrs']
1677    nAtoms = newPhase['Atoms']
1678    nSGData = newPhase['General']['SGData']
1679    #oAcof = G2lat.cell2A(oldPhase['General']['Cell'][1:7])
1680    #nAcof = G2lat.cell2A(newPhase['General']['Cell'][1:7])
1681    item = G2gd.GetGPXtreeItemId(G2frame,G2frame.root,'Constraints')
1682    if not item:
1683        print('Error: no constraints in Data Tree')
1684        return
1685    constraints = G2frame.GPXtree.GetItemPyData(item)
1686    xnames = ['dAx','dAy','dAz']
1687    # constraints on matching atom params between phases
1688    for ia,code in enumerate(atCodes):
1689        atom = nAtoms[ia]
1690        if not ia and atom[cia] == 'A':
1691            wx.MessageDialog(G2frame,
1692                'Anisotropic thermal motion constraints are not developed at the present time',
1693                'Anisotropic thermal constraint?',style=wx.ICON_INFORMATION).ShowModal()
1694        siteSym = G2spc.SytSym(atom[cx:cx+3],nSGData)[0]
1695        CSX = G2spc.GetCSxinel(siteSym)
1696#        CSU = G2spc.GetCSuinel(siteSym)
1697        item = code.split('+')[0]
1698        iat,opr = item.split(':')
1699        Nop = abs(int(opr))%100-1
1700        if '-' in opr:
1701            Nop *= -1
1702        Opr = oldPhase['General']['SGData']['SGOps'][abs(Nop)][0]
1703        if Nop < 0:         #inversion
1704            Opr *= -1
1705        XOpr = np.inner(Opr,Trans)
1706        for i,ix in enumerate(list(CSX[0])):
1707            if not ix:
1708                continue
1709            name = xnames[i]
1710            IndpCon = [1.0,G2obj.G2VarObj('%d::%s:%d'%(npId,name,ia))]
1711            DepCons = []
1712            for iop,opval in enumerate(XOpr[i]):
1713                if opval:
1714                    DepCons.append([opval,G2obj.G2VarObj('%d::%s:%s'%(opId,xnames[iop],iat))])
1715            if len(DepCons) == 1:
1716                constraints['Phase'].append([DepCons[0],IndpCon,None,None,'e'])
1717            elif len(DepCons) > 1:
1718                for Dep in DepCons:
1719                    Dep[0] *= -1
1720                constraints['Phase'].append([IndpCon]+DepCons+[0.0,None,'c'])
1721        for name in ['Afrac','AUiso']:
1722            IndpCon = [1.0,G2obj.G2VarObj('%d::%s:%d'%(npId,name,ia))]
1723            DepCons = [1.0,G2obj.G2VarObj('%d::%s:%s'%(opId,name,iat))]
1724            constraints['Phase'].append([DepCons,IndpCon,None,None,'e'])
1725           
1726        # unfinished Anisotropic constraint generation
1727#        Uids = [[0,0,'AU11'],[1,1,'AU22'],[2,2,'AU33'],[0,1,'AU12'],[0,2,'AU13'],[1,2,'AU23']]
1728#        DepConsDict = dict(zip(Us,[[],[],[],[],[],[]]))
1729#        for iu,Uid in enumerate(Uids):
1730#            UMT = np.zeros((3,3))
1731#            UMT[Uid[0],Uid[1]] = 1
1732#            nUMT = G2lat.prodMGMT(UMT,invTrans)
1733#            nUT = G2lat.UijtoU6(nUMT)
1734#            for iu,nU in enumerate(nUT):
1735#                if abs(nU) > 1.e-8:
1736#                    parm = '%d::%s;%s'%(opId,Us[iu],iat)
1737#                    DepConsDict[Uid[2]].append([abs(nU%1.),G2obj.G2VarObj(parm)])
1738#        nUcof = atom[iu:iu+6]
1739#        conStrings = []
1740#        for iU,Usi in enumerate(Us):
1741#            parm = '%d::%s;%d'%(npId,Usi,ia)
1742#            parmDict[parm] = nUcof[iU]
1743#            varyList.append(parm)
1744#            IndpCon = [1.0,G2obj.G2VarObj(parm)]
1745#            conStr = str([IndpCon,DepConsDict[Usi]])
1746#            if conStr in conStrings:
1747#                continue
1748#            conStrings.append(conStr)
1749#            if len(DepConsDict[Usi]) == 1:
1750#                if DepConsDict[Usi][0]:
1751#                    constraints['Phase'].append([IndpCon,DepConsDict[Usi][0],None,None,'e'])
1752#            elif len(DepConsDict[Usi]) > 1:       
1753#                for Dep in DepConsDict[Usi]:
1754#                    Dep[0] *= -1
1755#                constraints['Phase'].append([IndpCon]+DepConsDict[Usi]+[0.0,None,'c'])
1756           
1757        #how do I do Uij's for most Trans?
1758
1759    # constraints on lattice parameters between phases
1760#    T = nl.inv(Trans).T
1761    # T = Trans.T
1762    # conMat = [
1763    #     [T[0,0]**2,T[0,1]**2,T[0,2]**2,T[0,0]*T[0,1],T[0,0]*T[0,2],T[0,1]*T[0,2]],
1764    #     [T[1,0]**2,T[1,1]**2,T[1,2]**2,T[1,0]*T[1,1],T[1,0]*T[1,2],T[1,1]*T[1,2]],
1765    #     [T[2,0]**2,T[2,1]**2,T[2,2]**2,T[2,0]*T[2,1],T[2,0]*T[2,2],T[2,1]*T[2,2]],
1766    #     [2.*T[0,0]*T[1,0],2.*T[0,1]*T[1,1],2.*T[0,2]*T[1,2],T[0,0]*T[1,1]+T[0,1]*T[1,0],T[0,0]*T[1,2]+T[0,2]*T[1,0],T[0,1]*T[1,2]+T[0,2]*T[1,1]],
1767    #     [2.*T[0,0]*T[2,0],2.*T[0,1]*T[2,1],2.*T[0,2]*T[2,2],T[0,0]*T[2,1]+T[0,1]*T[2,0],T[0,0]*T[2,2]+T[0,2]*T[2,0],T[0,1]*T[2,2]+T[0,2]*T[2,1]],
1768    #     [2.*T[1,0]*T[2,0],2.*T[1,1]*T[2,1],2.*T[1,2]*T[2,2],T[1,0]*T[2,1]+T[1,1]*T[2,0],T[1,0]*T[2,2]+T[1,2]*T[2,0],T[1,1]*T[2,2]+T[1,2]*T[2,1]]
1769    #     ]
1770    # Gnew = conMat * A:
1771#         T00**2*a0  T01**2*a1 T02**2*a2 T00*T01*a3    T00*T02*a4    T01*T02*a5
1772#         T10**2*a0  T11**2*a1 T12**2*a2 T10*T11*a3    T10*T12*a4    T11*T12*a5
1773#         T20**2*a0  T21**2*a1 T22**2*a2 T20*T21*a3    T20*T22*a4    T21*T22*a5
1774#         2*T00*T10*a0      2*T01*T11*a1     2*T02*T12*a2     (T00*T11 + T01*T10)*a3      (T00*T12 + T02*T10)*a4      (T01*T12 + T02*T11)*a5
1775#         2*T00*T20*a0      2*T01*T21*a1     2*T02*T22*a2     (T00*T21 + T01*T20)*a3      (T00*T22 + T02*T20)*a4      (T01*T22 + T02*T21)*a5
1776#         2*T10*T20*a0      2*T11*T21*a1     2*T12*T22*a2     (T10*T21 + T11*T20)*a3      (T10*T22 + T12*T20)*a4      (T11*T22 + T12*T21)*a5
1777    # Generated as symbolic code using:
1778    # import sym
1779    # A0, A1, A2, A3, A4, A5 = sym.symbols('A0, A1, A2, A3, A4, A5')
1780    # G = sym.Matrix([ [A0,    A3/2.,  A4/2.], [A3/2.,  A1,    A5/2.], [A4/2., A5/2.,    A2]])
1781    # transformation matrix
1782    # T00, T10, T20, T01, T11, T21, T02, T12, T22 = sym.symbols('T00, T10, T20, T01, T11, T21, T02, T12, T22')
1783    # Tr = sym.Matrix([ [T00, T10, T20], [T01, T11, T21], [T02, T12, T22],])
1784    # Gnew = (Tr.T*G)*Tr
1785   
1786    #print('old A',G2lat.cell2A(oldPhase['General']['Cell'][1:7]))
1787    #print('new A',G2lat.cell2A(newPhase['General']['Cell'][1:7]))
1788   
1789#this is still incorrect for hex/trig/ortho/tetragonal --> monoclinic
1790   
1791#    for iAnew,Asi in enumerate(['A0','A1','A2','A3','A4','A5']): # loop through A[i] for new cell
1792#        Nparm = str(npId) + '::' + Asi
1793#        if Nparm != SetUniqAj(npId,iAnew,nSGData):
1794#            continue # skip: Ai constrained from Aj or must be zero
1795#        multDict = {}
1796#        for iAorg in range(6):
1797#            cA = conMat[iAnew][iAorg] # coeff for A[i] in constraint matrix
1798#            if abs(cA) < 1.e-8: continue
1799#            parm = SetUniqAj(opId,iAorg,oSGData) # translate to unique A[i] in original cell
1800#            if not parm: continue # must be zero
1801#            # sum coeff
1802#            if parm in multDict:
1803#                multDict[parm] += cA
1804#            else:
1805#                multDict[parm] = cA
1806#        # any non-zero multipliers?
1807#        maxMult = 0
1808#        for i in multDict:
1809#            maxMult = max(maxMult,abs(multDict[i]))
1810#        if maxMult <= 0:  # Nparm computes as zero; Fix this parameter
1811#            constraints['Phase'] += [[
1812#                [0.0,G2obj.G2VarObj(Nparm)],
1813#                None,None,'h']]
1814#        elif len(multDict) == 1:        # create equivalence
1815#            key = list(multDict.keys())[0]
1816#            constraints['Phase'] += [[
1817#                [1.0,G2obj.G2VarObj(key)],
1818#                [multDict[key],G2obj.G2VarObj(Nparm)],
1819#                None,None,'e']]
1820#        else:                           # create constraint
1821#            constr = [[-1.0,G2obj.G2VarObj(Nparm)]]
1822#            for key in multDict:
1823#                constr += [[multDict[key],G2obj.G2VarObj(key)]]
1824#            constr += [0.0,None,'c']
1825#            constraints['Phase'] += [constr]
1826   
1827    # constraints on HAP Scale, etc.
1828    for hId,hist in enumerate(UseList):    #HAP - seems OK
1829        ohapkey = '%d:%d:'%(opId,hId)
1830        nhapkey = '%d:%d:'%(npId,hId)
1831        IndpCon = [1.0,G2obj.G2VarObj(ohapkey+'Scale')]
1832        DepCons = [detTrans,G2obj.G2VarObj(nhapkey+'Scale')]
1833        constraints['HAP'].append([DepCons,IndpCon,None,None,'e'])
1834        for name in ['Size;i','Mustrain;i']:
1835            IndpCon = [1.0,G2obj.G2VarObj(ohapkey+name)]
1836            DepCons = [1.0,G2obj.G2VarObj(nhapkey+name)]
1837            constraints['HAP'].append([IndpCon,DepCons,None,None,'e'])
1838       
1839#### Rigid bodies #############################################################
1840resRBsel = None
1841def UpdateRigidBodies(G2frame,data):
1842    '''Called when Rigid bodies tree item is selected.
1843    Displays the rigid bodies in the data window
1844    ''' 
1845    def OnPageChanged(event):
1846        global resList
1847        resList = []
1848        if event:       #page change event!
1849            page = event.GetSelection()
1850        else:
1851            try:
1852                page = G2frame.rbBook.GetSelection()
1853            except:
1854                if GSASIIpath.GetConfigValue('debug'): print('DBG_gpx open error:C++ Run time error - skipped')
1855                return
1856        G2frame.rbBook.ChangeSelection(page)
1857        text = G2frame.rbBook.GetPageText(page)
1858        if text == 'Vector rigid bodies':
1859            G2gd.SetDataMenuBar(G2frame,G2frame.dataWindow.VectorBodyMenu)
1860            G2frame.Bind(wx.EVT_MENU, AddVectorRB, id=G2G.wxID_VECTORBODYADD)
1861            G2frame.Bind(wx.EVT_MENU, ExtractPhaseRB, id=G2G.wxID_VECTORBODYIMP)
1862            G2frame.Bind(wx.EVT_MENU, AddVectTrans, id=G2G.wxID_VECTORBODYEXTD)
1863            G2frame.Bind(wx.EVT_MENU, SaveVectorRB, id=G2G.wxID_VECTORBODYSAV)
1864            G2frame.Bind(wx.EVT_MENU, ReadVectorRB, id=G2G.wxID_VECTORBODYRD)
1865            G2frame.Page = [page,'vrb']
1866            UpdateVectorRB()
1867        elif text == 'Residue rigid bodies':
1868            G2gd.SetDataMenuBar(G2frame,G2frame.dataWindow.RigidBodyMenu)
1869            G2frame.Bind(wx.EVT_MENU, AddResidueRB, id=G2G.wxID_RIGIDBODYADD)
1870            G2frame.Bind(wx.EVT_MENU, ExtractPhaseRB, id=G2G.wxID_RIGIDBODYIMP)
1871            G2frame.Bind(wx.EVT_MENU, OnImportRigidBody, id=G2G.wxID_RIGIDBODYIMPORT)
1872            G2frame.Bind(wx.EVT_MENU, OnSaveRigidBody, id=G2G.wxID_RIGIDBODYSAVE)
1873            G2frame.Bind(wx.EVT_MENU, OnDefineTorsSeq, id=G2G.wxID_RESIDUETORSSEQ) #enable only if residue RBs exist?
1874            G2frame.Bind(wx.EVT_MENU, DumpVectorRB, id=G2G.wxID_RESBODYSAV)
1875            G2frame.Bind(wx.EVT_MENU, LoadVectorRB, id=G2G.wxID_RESBODYRD)
1876            G2frame.Page = [page,'rrb']
1877            UpdateResidueRB()
1878        else:
1879            G2gd.SetDataMenuBar(G2frame)
1880            #G2frame.Page = [page,'rrb']
1881           
1882    def getMacroFile(macName):
1883        defDir = os.path.join(os.path.split(__file__)[0],'GSASIImacros')
1884        dlg = wx.FileDialog(G2frame,message='Choose '+macName+' rigid body macro file',
1885            defaultDir=defDir,defaultFile="",wildcard="GSAS-II macro file (*.mac)|*.mac",
1886            style=wx.FD_OPEN | wx.FD_CHANGE_DIR)
1887        try:
1888            if dlg.ShowModal() == wx.ID_OK:
1889                macfile = dlg.GetPath()
1890                macro = open(macfile,'r')
1891                head = macro.readline()
1892                if macName not in head:
1893                    print (head)
1894                    print ('**** ERROR - wrong restraint macro file selected, try again ****')
1895                    macro = []
1896            else: # cancel was pressed
1897                macro = []
1898        finally:
1899            dlg.Destroy()
1900        return macro        #advanced past 1st line
1901       
1902    def getTextFile():
1903        dlg = wx.FileDialog(G2frame,'Choose rigid body text file', G2frame.LastGPXdir, '',
1904            "GSAS-II text file (*.txt)|*.txt|XYZ file (*.xyz)|*.xyz|"
1905            "Sybyl mol2 file (*.mol2)|*.mol2|PDB file (*.pdb;*.ent)|*.pdb;*.ent",
1906            wx.FD_OPEN | wx.FD_CHANGE_DIR)
1907        try:
1908            if dlg.ShowModal() == wx.ID_OK:
1909                txtfile = dlg.GetPath()
1910                ext = os.path.splitext(txtfile)[1]
1911                text = open(txtfile,'r')
1912            else: # cancel was pressed
1913                ext = ''
1914                text = []
1915        finally:
1916            dlg.Destroy()
1917        if 'ent' in ext:
1918            ext = '.pdb'
1919        return text,ext.lower()
1920       
1921    def OnImportRigidBody(event):
1922        page = G2frame.rbBook.GetSelection()
1923        if 'Vector' in G2frame.rbBook.GetPageText(page):
1924            pass
1925        elif 'Residue' in G2frame.rbBook.GetPageText(page):
1926            ImportResidueRB()
1927           
1928    def OnSaveRigidBody(event):
1929        page = G2frame.rbBook.GetSelection()
1930        if 'Vector' in G2frame.rbBook.GetPageText(page):
1931            pass
1932        elif 'Residue' in G2frame.rbBook.GetPageText(page):
1933            SaveResidueRB()
1934           
1935    def DumpVectorRB(event):
1936        global resRBsel
1937        if resRBsel not in data['Residue']:
1938            return
1939        rbData = data['Residue'][resRBsel]
1940        pth = G2G.GetExportPath(G2frame)
1941        dlg = wx.FileDialog(G2frame, 'Choose file to save residue rigid body',
1942            pth, '', 'RRB files (*.resbody)|*.resbody',
1943            wx.FD_SAVE|wx.FD_OVERWRITE_PROMPT)
1944        try:
1945            if dlg.ShowModal() == wx.ID_OK:
1946                filename = dlg.GetPath()
1947                filename = os.path.splitext(filename)[0]+'.resbody'  # set extension
1948                fp = open(filename,'w')
1949                fp.write('Name: '+rbData['RBname']+'\n')
1950                fp.write('atNames: ')
1951                for i in rbData['atNames']:
1952                    fp.write(str(i)+" ") 
1953                fp.write('\n')
1954                for item in rbData['rbSeq']:
1955                    fp.write('rbSeq: ') 
1956                    fp.write('{:d} {:d} {:.1f}: '.format(*item[:3]))
1957                    for num in item[3]:
1958                        fp.write('{:d} '.format(num))
1959                    fp.write('\n')
1960                for i,sym in enumerate(rbData['rbTypes']):
1961                    fp.write("{:3s}".format(sym))
1962                    fp.write('{:8.5f}{:9.5f}{:9.5f}   '
1963                            .format(*rbData['rbXYZ'][i]))
1964                    fp.write('\n')
1965                fp.close()
1966                print ('Vector rigid body saved to: '+filename)
1967               
1968        finally:
1969            dlg.Destroy()
1970       
1971    def LoadVectorRB(event):
1972        AtInfo = data['Residue']['AtInfo']
1973        pth = G2G.GetExportPath(G2frame)
1974        dlg = wx.FileDialog(G2frame, 'Choose file to read vector rigid body',
1975            pth, '', 'RRB files (*.resbody)|*.resbody',
1976            wx.FD_OPEN)
1977        try:
1978            if dlg.ShowModal() == wx.ID_OK:
1979                filename = dlg.GetPath()
1980                filename = os.path.splitext(filename)[0]+'.resbody'  # set extension
1981                fp = open(filename,'r')
1982                l = fp.readline().strip()
1983                if 'Name' not in l:
1984                    fp.close()
1985                    G2frame.ErrorDialog('Read Error',
1986                        'File '+filename+' does not start with Name\nFirst line ='
1987                        +l+'\ninvalid file',parent=G2frame)
1988                    return
1989                name = l.split(':')[1].strip()
1990                line = fp.readline().strip().split(':')[1].split()
1991                atNames = [i for i in line]
1992                types = []
1993                coords = []
1994                l = fp.readline().strip()
1995                rbSeq = []
1996                while 'rbSeq' in l:
1997                    tag,vals,lst = l.split(':')
1998                    seq = []
1999                    for t,v in zip((int,int,float),vals.split()): 
2000                        seq.append(t(v))
2001                    seq.append([])
2002                    for num in lst.split():
2003                        seq[-1].append(int(num))
2004                    rbSeq.append(seq)
2005                    l = fp.readline().strip()                   
2006                while l:
2007                    nums = l.strip().split()
2008                    types.append(nums.pop(0))
2009                    t = types[-1]
2010                    if t not in AtInfo:
2011                        Info = G2elem.GetAtomInfo(t)
2012                        AtInfo[t] = [Info['Drad'],Info['Color']]
2013                    coords.append([float(nums.pop(0)) for j in range(3)])
2014                    l = fp.readline().strip()
2015                fp.close()
2016            else:
2017                return       
2018        finally:
2019            dlg.Destroy()
2020        coords = np.array(coords)
2021        rbid = ran.randint(0,sys.maxsize)
2022        namelist = [data['Residue'][key]['RBname'] for key in data['Residue']
2023                        if 'RBname' in data['Residue'][key]]
2024        name = G2obj.MakeUniqueLabel(name,namelist)
2025        data['Residue'][rbid] = {'RBname':name,
2026                'rbXYZ': coords,
2027                'rbRef':[0,1,2,False],
2028                'rbTypes':types, 'atNames':atNames,
2029                'useCount':0,
2030                'rbSeq':rbSeq, 'SelSeq':[0,0],}
2031        data['RBIds']['Residue'].append(rbid)
2032        UpdateResidueRB()
2033   
2034    def AddVectorRB(event):
2035        'Create a new vector rigid body'
2036        AtInfo = data['Vector']['AtInfo']
2037        dlg = G2G.MultiIntegerDialog(G2frame,'New Rigid Body',['No. atoms','No. translations'],[3,1])
2038        if dlg.ShowModal() == wx.ID_OK:
2039            nAtoms,nTrans = dlg.GetValues()
2040            if nAtoms < 3:
2041                dlg.Destroy()
2042                G2G.G2MessageBox(G2frame,'A vector rigid body must have 3 or more atoms')
2043                return
2044            rbid = ran.randint(0,sys.maxsize)
2045            vecMag = [1.0 for i in range(nTrans)]
2046            vecRef = [False for i in range(nTrans)]
2047            vecVal = [np.zeros((nAtoms,3)) for j in range(nTrans)]
2048            rbTypes = ['C' for i in range(nAtoms)]
2049            Info = G2elem.GetAtomInfo('C')
2050            AtInfo['C'] = [Info['Drad'],Info['Color']]
2051            name = 'UNKRB'
2052            namelist = [data['Vector'][key]['RBname'] for key in data['Vector']
2053                        if 'RBname' in data['Vector'][key]]
2054            name = G2obj.MakeUniqueLabel(name,namelist)
2055            data['Vector'][rbid] = {'RBname':name,'VectMag':vecMag,'rbXYZ':np.zeros((nAtoms,3)),
2056                'rbRef':[0,1,2,False],'VectRef':vecRef,'rbTypes':rbTypes,'rbVect':vecVal,'useCount':0}
2057            data['RBIds']['Vector'].append(rbid)
2058        dlg.Destroy()
2059        UpdateVectorRB()
2060
2061    def ExtractPhaseRB(event):
2062        'Extract a rigid body from a file with a phase'
2063        def SetupDrawing(atmData):
2064            '''Add the dicts needed for G2plt.PlotStructure to work to the
2065            reader .Phase object
2066            '''
2067            generalData = atmData['General']
2068            generalData['BondRadii'] = []
2069
2070            G2phG.SetDrawingDefaults(atmData['Drawing'])
2071            atmData['Drawing'].update(
2072                {'oldxy':[0.,0.],'Quaternion':[0.,0.,0.,1.],'cameraPos':150.,
2073                     'viewDir':[0,0,1],'atomPtrs': [2, 1, 6, 17],
2074                     })
2075            atmData['Drawing']['showRigidBodies'] = False
2076            generalData['Map'] = {'MapType':False, 'rho':[]}
2077            generalData['AtomTypes'] = []
2078            generalData['BondRadii'] = []
2079            generalData['AngleRadii'] = []
2080            generalData['vdWRadii'] = []
2081            generalData['Color'] = []
2082            generalData['Isotopes'] = {}
2083            generalData['Isotope'] = {}
2084            cx,ct,cs,cia = generalData['AtomPtrs']
2085            generalData['Mydir'] = G2frame.dirname
2086            for iat,atom in enumerate(atmData['Atoms']):
2087                atom[ct] = atom[ct].lower().capitalize()      #force elem symbol to standard form
2088                if atom[ct] not in generalData['AtomTypes'] and atom[ct] != 'UNK':
2089                    Info = G2elem.GetAtomInfo(atom[ct])
2090                    if not Info:
2091                        atom[ct] = 'UNK'
2092                        continue
2093                    atom[ct] = Info['Symbol'] # N.B. symbol might be changed by GetAtomInfo
2094                    generalData['AtomTypes'].append(atom[ct])
2095                    generalData['Z'] = Info['Z']
2096                    generalData['Isotopes'][atom[ct]] = Info['Isotopes']
2097                    generalData['BondRadii'].append(Info['Drad'])
2098                    generalData['AngleRadii'].append(Info['Arad'])
2099                    generalData['vdWRadii'].append(Info['Vdrad'])
2100                    if atom[ct] in generalData['Isotope']:
2101                        if generalData['Isotope'][atom[ct]] not in generalData['Isotopes'][atom[ct]]:
2102                            isotope = list(generalData['Isotopes'][atom[ct]].keys())[-1]
2103                            generalData['Isotope'][atom[ct]] = isotope
2104                    else:
2105                        generalData['Isotope'][atom[ct]] = 'Nat. Abund.'
2106                        if 'Nat. Abund.' not in generalData['Isotopes'][atom[ct]]:
2107                            isotope = list(generalData['Isotopes'][atom[ct]].keys())[-1]
2108                            generalData['Isotope'][atom[ct]] = isotope
2109                    generalData['Color'].append(Info['Color'])
2110                    # if generalData['Type'] == 'magnetic':
2111                    #     if len(landeg) < len(generalData['AtomTypes']):
2112                    #         landeg.append(2.0)
2113            atmData['Drawing']['Atoms'] = []
2114            for atom in atmData['Atoms']:
2115                atmData['Drawing']['Atoms'].append(G2mth.MakeDrawAtom(atmData,atom))
2116
2117        def onCancel(event,page=0):
2118            'complete or bail out from RB define, cleaning up'
2119            G2frame.rbBook.DeletePage(G2frame.rbBook.FindPage(pagename))
2120            G2frame.rbBook.SetSelection(page)
2121
2122        def Page1():
2123            '''Show the GUI for first stage of the rigid body with all atoms in
2124            phase in crystal coordinates. Select the atoms to go onto the
2125            next stage
2126            '''
2127            def ShowSelection(selections):
2128                'respond to change in atom selections'
2129                ct,cs = [1,8]
2130                generalData = rd.Phase['General']
2131                for i,atom in enumerate(rd.Phase['Drawing']['Atoms']):
2132                    if i in selections:
2133                        factor = 1
2134                    else:
2135                        factor = 2.5
2136                    atNum = generalData['AtomTypes'].index(atom[ct]) 
2137                    atom[cs] = list(np.array(generalData['Color'][atNum])//factor) 
2138                draw(*drawArgs)
2139            def onPage1OK(event):
2140                '1st section has been completed, move onto next'
2141                G2frame.G2plotNB.Delete(rd.Phase['General']['Name'])
2142                GetCoords(atmsel)
2143                Page2()
2144
2145            if 'macromolecular' == rd.Phase['General']['Type']:
2146                # for PDB imports, lets see if a quick reformat of atoms list will work
2147                rd.Phase['Atoms'] = [a[3:] for a in rd.Phase['Atoms']]
2148                rd.Phase['General']['AtomPtrs']  = [i-3 for i in rd.Phase['General']['AtomPtrs']]
2149                rd.Phase['General']['Type'] = 'nuclear'
2150            SetupDrawing(rd.Phase) # add information to reader object to allow plotting
2151            atomlist = [atom[0] for atom in rd.Phase['Atoms']]
2152            atmsel = list(range(len(rd.Phase['Atoms'])))
2153            # broken -- # why no bonds?
2154            #for atm in rd.Phase['Drawing']['Atoms']:
2155            #    atm[6] = 'balls & sticks'
2156
2157            draw,drawArgs = G2plt.PlotStructure(G2frame,rd.Phase,True)
2158            ShowSelection(atmsel)
2159
2160            if G2frame.rbBook.FindPage(pagename) is not None:
2161                G2frame.rbBook.DeletePage(G2frame.rbBook.FindPage(pagename))
2162
2163            RBImp = wx.ScrolledWindow(G2frame.rbBook)
2164            RBImpPnl = wx.Panel(RBImp)
2165            G2frame.rbBook.AddPage(RBImp,pagename)
2166            G2frame.rbBook.SetSelection(G2frame.rbBook.FindPage(pagename))
2167
2168            HelpInfo = '''
2169This window shows all the atoms that were read from the
2170selected phase file. Select the atoms that will be used in the
2171rigid body processing (this may include atoms needed to
2172define an axis or origin that will not be included in the
2173eventual rigid body.) Note that in the plot window,
2174unselected atoms appear much darker than selected atoms.
2175'''
2176            mainSizer = G2G.G2MultiChoiceWindow(RBImpPnl,
2177                            'Select atoms to import',
2178                            atomlist,atmsel,OnChange=ShowSelection,
2179                            helpText=HelpInfo)
2180
2181            # OK/Cancel buttons       
2182            btnsizer = wx.StdDialogButtonSizer()
2183            OKbtn = wx.Button(RBImpPnl, wx.ID_OK, 'Continue')
2184            OKbtn.SetDefault()
2185            btnsizer.AddButton(OKbtn)
2186            OKbtn.Bind(wx.EVT_BUTTON,onPage1OK)
2187            btn = wx.Button(RBImpPnl, wx.ID_CANCEL)
2188            btn.Bind(wx.EVT_BUTTON,onCancel)
2189            btnsizer.AddButton(btn)
2190            btnsizer.Realize()
2191            mainSizer.Add(btnsizer,0,wx.ALIGN_CENTER,50)
2192
2193            RBImpPnl.SetSizer(mainSizer,True)
2194
2195            mainSizer.Layout()   
2196            Size = mainSizer.GetMinSize()
2197            Size[0] += 40
2198            Size[1] = max(Size[1],G2frame.GetSize()[1]-200) + 20
2199            RBImpPnl.SetSize(Size)
2200            RBImp.SetScrollbars(10,10,Size[0]/10-4,Size[1]/10-1)
2201            RBImp.Scroll(0,0)
2202
2203        def Page2():
2204            '''Show the GUI for the second stage, where selected atoms are
2205            now in Cartesian space, manipulate the axes and export selected
2206            atoms to a vector or residue rigid body.
2207            '''
2208            def UpdateDraw(event=None):
2209                'Called when info changes in grid, replots'
2210                UpdateVectorBody(rbData)
2211                DrawCallback()
2212               
2213            def onSetAll(event):
2214                'Set all atoms as selected'
2215                grid.completeEdits()
2216                for i in range(len(rd.Phase['RBselection'])):
2217                    rd.Phase['RBselection'][i] = 1 # table needs 0/1 for T/F
2218                grid.ForceRefresh()
2219                UpdateDraw()
2220               
2221            def onToggle(event):
2222                'Toggles selection state for all atoms'
2223                grid.completeEdits()
2224                for i in range(len(rd.Phase['RBselection'])):
2225                    rd.Phase['RBselection'][i] = int(not rd.Phase['RBselection'][i])
2226                grid.ForceRefresh()
2227                UpdateDraw()
2228               
2229            def onSetOrigin(event):
2230                'Resets origin to midpoint between all selected atoms'
2231                grid.completeEdits()
2232                center = np.array([0.,0.,0.])
2233                count = 0
2234                for i in range(len(rd.Phase['RBselection'])):
2235                    if rd.Phase['RBselection'][i]:
2236                        count += 1
2237                        center += rd.Phase['RBcoords'][i]
2238                if count:
2239                    rd.Phase['RBcoords'] -= center/count
2240                grid.ForceRefresh()
2241                UpdateDraw()
2242               
2243            def onSetX(event):
2244                grid.completeEdits()
2245                center = np.array([0.,0.,0.])
2246                count = 0
2247                for i in range(len(rd.Phase['RBselection'])):
2248                    if rd.Phase['RBselection'][i]:
2249                        count += 1
2250                        center += rd.Phase['RBcoords'][i]
2251                if not count:
2252                    G2G.G2MessageBox(G2frame,'No atoms selected',
2253                                    'Selection required')
2254                    return
2255                XYZP = center/count
2256                if np.sqrt(sum(XYZP**2)) < 0.1:
2257                    G2G.G2MessageBox(G2frame,
2258                            'The selected atom(s) are too close to the origin',
2259                            'near origin')
2260                    return
2261                if bntOpts['direction'] == 'y':
2262                    YP = XYZP / np.sqrt(np.sum(XYZP**2))
2263                    ZP = np.cross((1,0,0),YP)
2264                    if sum(ZP*ZP) < .1: # pathological condition: Y' along X
2265                        ZP = np.cross((0,0,1),YP)
2266                    XP = np.cross(YP,ZP)
2267                elif bntOpts['direction'] == 'z':
2268                    ZP = XYZP / np.sqrt(np.sum(XYZP**2))
2269                    XP = np.cross((0,1,0),ZP)
2270                    if sum(XP*XP) < .1: # pathological condition: X' along Y
2271                        XP = np.cross((0,0,1),ZP)
2272                    YP = np.cross(ZP,XP)
2273                else:
2274                    XP = XYZP / np.sqrt(np.sum(XYZP**2))
2275                    YP = np.cross((0,0,1),XP)
2276                    if sum(YP*YP) < .1: # pathological condition: X' along Z
2277                        YP = np.cross((0,1,0),XP)
2278                    ZP = np.cross(XP,YP)
2279                trans = np.array((XP,YP,ZP))
2280                # update atoms in place
2281                rd.Phase['RBcoords'][:] = np.inner(trans,rd.Phase['RBcoords']).T
2282                grid.ForceRefresh()
2283                UpdateDraw()
2284
2285            def onSetPlane(event): 
2286                '''Finds eigen vector/matrix for best "ellipsoid" about atoms;
2287                rotate atoms so that smallest axis is along choice.
2288                '''
2289                grid.completeEdits()
2290                selList = [i==1 for i in rd.Phase['RBselection']]
2291                XYZ = rd.Phase['RBcoords'][selList]
2292                Natoms = len(XYZ)
2293                if Natoms < 3: 
2294                    G2G.G2MessageBox(G2frame,'A plane requires three or more atoms','Need more atoms')
2295                    return
2296                Zmat = np.zeros((3,3))
2297                for xyz in XYZ:
2298                    Zmat += np.outer(xyz.T,xyz)
2299                Evec,Emat = nl.eig(Zmat)
2300                Order = np.argsort(np.nan_to_num(Evec))     #short-long order
2301                if bntOpts['plane'] == 'xy':        #short along z
2302                    trans = np.array([Emat[Order[2]],Emat[Order[1]],Emat[Order[0]]])
2303                elif bntOpts['plane'] == 'yz':      #short along x
2304                    trans = np.array([Emat[Order[0]],Emat[Order[2]],Emat[Order[1]]])
2305                elif bntOpts['plane'] == 'xz':      #short along y
2306                    trans = np.array([Emat[Order[1]],Emat[Order[0]],Emat[Order[2]]])
2307                else:
2308                    print('unexpected plane',bntOpts['plane'])
2309                    return
2310                # update atoms in place
2311                rd.Phase['RBcoords'][:] = np.inner(trans,rd.Phase['RBcoords']).T
2312                grid.ForceRefresh()
2313                UpdateDraw()
2314
2315            def onWriteXYZ(event):
2316                '''Writes selected atoms in a .xyz file for use in Avogadro, etc.
2317                '''
2318                grid.completeEdits()
2319                center = np.array([0.,0.,0.])
2320                count = 0
2321                for i in range(len(rd.Phase['RBselection'])):
2322                    if rd.Phase['RBselection'][i]:
2323                        count += 1
2324                        center += rd.Phase['RBcoords'][i]
2325                if count:
2326                    center /= count
2327                else:
2328                    print('nothing selected')
2329                    return
2330                obj = G2IO.ExportBaseclass(G2frame,'XYZ','.xyz')
2331                #obj.InitExport(None)
2332                if obj.ExportSelect():    # set export parameters; ask for file name
2333                    return
2334                obj.OpenFile()
2335                obj.Write(str(count))
2336                obj.Write('')
2337                for i in range(len(rd.Phase['RBselection'])):
2338                    if rd.Phase['RBselection'][i]:
2339                        line = ' ' + rd.Phase['RBtypes'][i]
2340                        for xyz in rd.Phase['RBcoords'][i]:
2341                            line += ' ' + str(xyz)
2342                        obj.Write(line)
2343                obj.CloseFile()
2344                #GSASIIpath.IPyBreak()
2345               
2346            def onAddVector(event):
2347                '''Adds selected atoms as a new vector rigid body.
2348                Closes out the importer tab when done.
2349                '''
2350                grid.completeEdits()
2351                name = os.path.splitext(os.path.split(filename)[1])[0]
2352                namelist = [data['Vector'][key]['RBname'] for key in
2353                        data['Vector'] if 'RBname' in data['Vector'][key]]
2354                name = G2obj.MakeUniqueLabel(name,namelist)
2355                rb = MakeVectorBody(name)
2356                UpdateVectorBody(rb,True)
2357                if len(rb['rbTypes']) < 3: return # must have at least 3 atoms
2358                rbid = ran.randint(0,sys.maxsize)
2359                data['Vector'][rbid] = rb
2360                data['RBIds']['Vector'].append(rbid)
2361                for t in rb['rbTypes']:
2362                    if t in data['Vector']['AtInfo']: continue
2363                    Info = G2elem.GetAtomInfo(t)
2364                    data['Vector']['AtInfo'][t] = [Info['Drad'],Info['Color']]
2365                G2frame.G2plotNB.Delete('Rigid body')
2366                onCancel(event,0)
2367               
2368            def onAddResidue(event):
2369                '''Adds selected atoms as a new residue rigid body.
2370                Closes out the importer tab when done.
2371                '''
2372                grid.completeEdits()
2373                name = os.path.split(filename)[1]
2374                rbXYZ = []
2375                rbTypes = []
2376                atNames = []
2377                for i in rd.Phase['RBindex']:
2378                    if rd.Phase['RBselection'][i]:
2379                        rbXYZ.append(rd.Phase['RBcoords'][i])
2380                        rbTypes.append(rd.Phase['RBtypes'][i])
2381                        atNames.append(rd.Phase['RBlbls'][i])
2382                if len(rbTypes) < 3: return # must have at least 3 atoms
2383                rbXYZ = np.array(rbXYZ)
2384                rbid = ran.randint(0,sys.maxsize)
2385                namelist = [data['Residue'][key]['RBname'] for key in
2386                        data['Residue'] if 'RBname' in data['Residue'][key]]
2387                name = G2obj.MakeUniqueLabel(name,namelist)
2388                data['Residue'][rbid] = {'RBname':name,'rbXYZ':rbXYZ,
2389                    'rbTypes':rbTypes,'atNames':atNames,'rbRef':[0,1,2,False],
2390                    'rbSeq':[],'SelSeq':[0,0],'useCount':0}
2391                data['RBIds']['Residue'].append(rbid)
2392                for t in rbTypes:
2393                    if t in data['Residue']['AtInfo']: continue
2394                    Info = G2elem.GetAtomInfo(t)
2395                    data['Residue']['AtInfo'][t] = [Info['Drad'],Info['Color']]
2396
2397                print ('Rigid body added')
2398                G2frame.G2plotNB.Delete('Rigid body')
2399                onCancel(event,1)
2400
2401            if G2frame.rbBook.FindPage(pagename) is not None:
2402                G2frame.rbBook.DeletePage(G2frame.rbBook.FindPage(pagename))
2403            RBImp = wx.ScrolledWindow(G2frame.rbBook)
2404            RBImpPnl = wx.Panel(RBImp)
2405            G2frame.rbBook.AddPage(RBImp,pagename)
2406            G2frame.rbBook.SetSelection(G2frame.rbBook.FindPage(pagename))
2407            AtInfo = {}
2408            for t in rd.Phase['RBtypes']:
2409                if t in AtInfo: continue
2410                Info = G2elem.GetAtomInfo(t)
2411                AtInfo[t] = [Info['Drad'],Info['Color']]
2412            plotDefaults = {'oldxy':[0.,0.],'Quaternion':[0.,0.,0.,1.],'cameraPos':30.,'viewDir':[0,0,1],}
2413
2414            rd.Phase['RBindex'] = list(range(len(rd.Phase['RBtypes'])))
2415            rd.Phase['RBselection'] = len(rd.Phase['RBtypes']) * [1]
2416            name = 'UNKRB'
2417            namelist = [data['Vector'][key]['RBname'] for key in
2418                        data['Vector'] if 'RBname' in data['Vector'][key]]
2419            name = G2obj.MakeUniqueLabel(name,namelist)
2420            rbData = MakeVectorBody()
2421            DrawCallback = G2plt.PlotRigidBody(G2frame,'Vector',AtInfo,rbData,plotDefaults)
2422
2423            mainSizer = wx.BoxSizer(wx.HORIZONTAL)
2424            btnSizer = wx.BoxSizer(wx.VERTICAL)
2425            helpText = '''
2426In this window, if wanted,
2427one can select one or more atoms and use them
2428to define an origin, a specified axis or place the selected atoms into
2429a selected plane. (Different sets of atoms can be used for each
2430operation.)
2431%%Once that is done, atoms can be selected and can be exported in a
2432"XYZ" file for use in a program such as Avogadro or can be used to
2433create a Vector or Residue rigid body.
2434'''
2435            btnSizer.Add(G2G.HelpButton(RBImpPnl,helpText,wrap=400),
2436                             0,wx.ALIGN_RIGHT)
2437            btnSizer.Add(wx.StaticText(RBImpPnl,wx.ID_ANY,'Reorder atoms by dragging'),0,wx.ALL)
2438            btnSizer.Add((-1,15))
2439            btn = wx.Button(RBImpPnl, wx.ID_ANY, 'Set All')
2440            btn.Bind(wx.EVT_BUTTON,onSetAll)
2441            btnSizer.Add(btn,0,wx.ALIGN_CENTER)
2442            btn = wx.Button(RBImpPnl, wx.ID_ANY, 'Toggle')
2443            btn.Bind(wx.EVT_BUTTON,onToggle)
2444            btnSizer.Add(btn,0,wx.ALIGN_CENTER)
2445            btnSizer.Add((-1,15))
2446            btnSizer.Add(wx.StaticText(RBImpPnl,wx.ID_ANY,'Reorient using selected\natoms...'),0,wx.ALL)
2447            btnSizer.Add((-1,5))
2448            btn = wx.Button(RBImpPnl, wx.ID_ANY, 'Set origin')
2449            btn.Bind(wx.EVT_BUTTON,onSetOrigin)
2450            btnSizer.Add(btn,0,wx.ALIGN_CENTER)
2451
2452            bntOpts = {'plane':'xy','direction':'x'}
2453            inSizer = wx.BoxSizer(wx.HORIZONTAL)
2454            btn = wx.Button(RBImpPnl, wx.ID_ANY, 'Place in plane')
2455            btn.Bind(wx.EVT_BUTTON,onSetPlane)
2456            inSizer.Add(btn)
2457            inSizer.Add(G2G.G2ChoiceButton(RBImpPnl,('xy','yz','xz'),None,None,bntOpts,'plane'))
2458            btnSizer.Add(inSizer,0,wx.ALIGN_CENTER)
2459           
2460            inSizer = wx.BoxSizer(wx.HORIZONTAL)
2461            btn = wx.Button(RBImpPnl, wx.ID_ANY, 'Define as')
2462            btn.Bind(wx.EVT_BUTTON,onSetX)
2463            inSizer.Add(btn)
2464            inSizer.Add(G2G.G2ChoiceButton(RBImpPnl,('x','y','z'),None,None,bntOpts,'direction'))
2465            btnSizer.Add(inSizer,0,wx.ALIGN_CENTER)
2466           
2467            btnSizer.Add((-1,15))
2468            btnSizer.Add(wx.StaticText(RBImpPnl,wx.ID_ANY,'Use selected atoms to\ncreate...'),0,wx.ALL)
2469            btnSizer.Add((-1,5))
2470            btn = wx.Button(RBImpPnl, wx.ID_ANY, 'export as xyz')
2471            btn.Bind(wx.EVT_BUTTON,onWriteXYZ)
2472            btnSizer.Add(btn,0,wx.ALIGN_CENTER)
2473            btnSizer.Add((-1,10))
2474            btn = wx.Button(RBImpPnl, wx.ID_ANY, 'a Vector Body')
2475            btn.Bind(wx.EVT_BUTTON,onAddVector)
2476            btnSizer.Add(btn,0,wx.ALIGN_CENTER)
2477            btn = wx.Button(RBImpPnl, wx.ID_ANY, 'a Residue Body')
2478            btn.Bind(wx.EVT_BUTTON,onAddResidue)
2479            btnSizer.Add(btn,0,wx.ALIGN_CENTER)
2480            btn = wx.Button(RBImpPnl, wx.ID_CANCEL)
2481            btn.Bind(wx.EVT_BUTTON,onCancel)
2482            btnSizer.Add((-1,10))
2483            btnSizer.Add(btn,0,wx.ALIGN_CENTER)
2484
2485            mainSizer.Add(btnSizer)
2486            mainSizer.Add((5,5))
2487            grid = DragableRBGrid(RBImpPnl,rd.Phase,UpdateDraw)
2488            mainSizer.Add(grid)
2489            RBImpPnl.SetSizer(mainSizer,True)
2490            mainSizer.Layout()   
2491            Size = mainSizer.GetMinSize()
2492            Size[0] += 40
2493            Size[1] = max(Size[1],G2frame.GetSize()[1]-200) + 20
2494            RBImpPnl.SetSize(Size)
2495            RBImp.SetScrollbars(10,10,Size[0]/10-4,Size[1]/10-1)
2496            RBImp.Scroll(0,0)
2497
2498        def GetCoords(atmsel):
2499            '''Create orthogonal coordinates for selected atoms.
2500            Place the origin at the center of the body
2501            '''
2502            atms = rd.Phase['Atoms']
2503            cell = rd.Phase['General']['Cell'][1:7]
2504            Amat,Bmat = G2lat.cell2AB(cell)
2505            rd.Phase['RBcoords'] = np.array([np.inner(Amat,atms[i][3:6]) for i in atmsel])
2506            rd.Phase['RBcoords'] -= rd.Phase['RBcoords'].mean(axis=0)  # origin to middle
2507            rd.Phase['RBtypes'] = [atms[i][1] for i in atmsel]
2508            rd.Phase['RBlbls'] = [atms[i][0] for i in atmsel]
2509
2510        def UpdateVectorBody(rb,useSelection=False):
2511            '''Put the atoms in order to pass for plotting or for storage as
2512            a vector rigid body.
2513
2514            :param dict rb: rigid body contents created in :func:`MakeVectorBody`
2515            :param bool useSelection: True if the rd.Phase['RBselection']
2516              values will be used to select which atoms are included in the
2517              rigid body. If False (default) they are included in rb
2518              and are used for plotting.         
2519            '''
2520            coordlist = []
2521            typeslist = []
2522            sellist = []
2523            for i in rd.Phase['RBindex']:
2524                use = True
2525                if useSelection and not rd.Phase['RBselection'][i]: use = False
2526                if use:
2527                    coordlist.append(rd.Phase['RBcoords'][i])
2528                    typeslist.append(rd.Phase['RBtypes'][i])
2529                    sellist.append(rd.Phase['RBselection'][i])
2530            coordlist = np.array(coordlist)
2531            rb['rbXYZ'] = coordlist
2532            rb['rbVect'] = [coordlist]
2533            rb['rbTypes'] = typeslist
2534            if not useSelection:
2535                rb['Selection'] = sellist
2536            elif 'Selection' in rb:
2537                del rb['Selection']
2538
2539        def MakeVectorBody(name=''):
2540            '''Make the basic vector rigid body dict (w/o coordinates) used for
2541            export and for plotting
2542            '''
2543            vecMag = [1.0]
2544            vecRef = [False]
2545            rb = {'RBname':name,'VectMag':vecMag,
2546                    'rbRef':[0,1,2,False],'VectRef':vecRef,
2547                    'useCount':0}
2548            UpdateVectorBody(rb)
2549            return rb
2550
2551        # too lazy to figure out why wx crashes
2552        if wx.__version__.split('.')[0] != '4':
2553            wx.MessageBox('Sorry, wxPython 4.x is required to run this command',
2554                                  caption='Update Python',
2555                                  style=wx.ICON_EXCLAMATION)
2556            return
2557        if platform.python_version()[:1] == '2':
2558            wx.MessageBox('Sorry, Python >=3.x is required to run this command',
2559                                  caption='Update Python',
2560                                  style=wx.ICON_EXCLAMATION)
2561            return
2562
2563        # get importer type and a phase file of that type
2564        G2sc.LoadG2fil()
2565        choices = [rd.formatName for  rd in G2sc.Readers['Phase']] 
2566        dlg = G2G.G2SingleChoiceDialog(G2frame,'Select the format of the file',
2567                                     'select format',choices)
2568        dlg.CenterOnParent()
2569        try:
2570            if dlg.ShowModal() == wx.ID_OK:
2571                col = dlg.GetSelection()
2572            else:
2573                col = None
2574                return
2575        finally:
2576            dlg.Destroy()
2577        reader = G2sc.Readers['Phase'][col]
2578
2579        choices = reader.formatName + " file ("
2580        w = ""
2581        for extn in reader.extensionlist:
2582            if w != "": w += ";"
2583            w += "*" + extn
2584        choices += w + ")|" + w
2585        #choices += "|zip archive (.zip)|*.zip"
2586        if not reader.strictExtension:
2587            choices += "|any file (*.*)|*.*"
2588        typ = '( type '+reader.formatName+')'
2589        filelist = G2G.GetImportFile(G2frame,
2590                        message="Choose phase input file"+typ,
2591                        defaultFile="",wildcard=choices,style=wx.FD_OPEN)
2592        if len(filelist) != 1: return
2593
2594        # read in the phase file
2595        filename = filelist[0]
2596        rd = reader
2597        with open(filename, 'r'):
2598            rd.ReInitialize()
2599            rd.errors = ""
2600            if not rd.ContentsValidator(filename):   # Report error
2601                G2fl.G2Print("Warning: File {} has a validation error".format(filename))
2602                return
2603            if len(rd.selections) > 1:
2604                print("File {} has {} phases. This is unexpected."
2605                                    .format(filename,len(rd.selections)))
2606                return
2607
2608            rd.objname = os.path.basename(filename)
2609            try:
2610                rd.Reader(filename)
2611            except Exception as msg:
2612                G2fl.G2Print("Warning: read of file {} failed\n{}".format(
2613                    filename,rd.errors))
2614                if GSASIIpath.GetConfigValue('debug'):
2615                    print(msg)
2616                    import traceback
2617                    print (traceback.format_exc())
2618                    GSASIIpath.IPyBreak()
2619                return
2620
2621        pagename = 'Rigid body importer'
2622        Page1()
2623        return
2624
2625    def AddVectTrans(event):
2626        'Add a translation to an existing vector rigid body'
2627        choices = []
2628        rbIdlist = []
2629        for rbid in data['RBIds']['Vector']:
2630            if rbid != 'AtInfo':
2631                rbIdlist.append(rbid)
2632                choices.append(data['Vector'][rbid]['RBname'])
2633        if len(choices) == 0:
2634            G2G.G2MessageBox(G2frame,'No Vector Rigid Bodies found',
2635                                 'No VR Bodies')
2636            return
2637        elif len(choices) == 1:
2638            rbid = rbIdlist[0]
2639        else:
2640            dlg = G2G.G2SingleChoiceDialog(G2frame,'Select the rigid body to save',
2641                                  'select format',choices)
2642            try:
2643                if dlg.ShowModal() == wx.ID_OK:
2644                    rbid = rbIdlist[dlg.GetSelection()]
2645                else:
2646                    return
2647            finally:
2648                dlg.Destroy()
2649        data['Vector'][rbid]['VectMag'] += [1.0]
2650        data['Vector'][rbid]['VectRef'] += [False]
2651        nAtoms = len(data['Vector'][rbid]['rbXYZ'])
2652        data['Vector'][rbid]['rbVect'] += [np.zeros((nAtoms,3))]
2653        UpdateVectorRB()
2654       
2655    def SaveVectorRB(event):
2656        choices = []
2657        rbIdlist = []
2658        for rbid in data['RBIds']['Vector']:
2659            if rbid != 'AtInfo':
2660                rbIdlist.append(rbid)
2661                choices.append(data['Vector'][rbid]['RBname'])
2662        if len(choices) == 0:
2663            G2G.G2MessageBox(G2frame,'No Vector Rigid Bodies found',
2664                                 'No VR Bodies')
2665            return
2666        elif len(choices) == 1:
2667            rbid = rbIdlist[0]
2668        else:
2669            dlg = G2G.G2SingleChoiceDialog(G2frame,'Select the rigid body to save',
2670                                  'select format',choices)
2671            try:
2672                if dlg.ShowModal() == wx.ID_OK:
2673                    rbid = rbIdlist[dlg.GetSelection()]
2674                else:
2675                    return
2676            finally:
2677                dlg.Destroy()
2678             
2679        pth = G2G.GetExportPath(G2frame)
2680        dlg = wx.FileDialog(G2frame, 'Choose file to save vector rigid body',
2681            pth, '', 'VRB files (*.vecbody)|*.vecbody',
2682            wx.FD_SAVE|wx.FD_OVERWRITE_PROMPT)
2683        try:
2684            if dlg.ShowModal() == wx.ID_OK:
2685                filename = dlg.GetPath()
2686                filename = os.path.splitext(filename)[0]+'.vecbody'  # set extension
2687                fp = open(filename,'w')
2688                fp.write('Name: '+data['Vector'][rbid]['RBname']+'\n')
2689                fp.write('Trans: ')
2690                for i in data['Vector'][rbid]['VectMag']:
2691                    fp.write(str(i)+" ") 
2692                fp.write('\n')
2693                ntrans = len(data['Vector'][rbid]['VectMag'])
2694                for i,sym in enumerate(data['Vector'][rbid]['rbTypes']):
2695                    fp.write("{:3s}".format(sym))
2696                    for j in range(ntrans):
2697                        fp.write('{:8.5f}{:9.5f}{:9.5f}   '
2698                            .format(*data['Vector'][rbid]['rbVect'][j][i]))
2699                    fp.write('\n')
2700                fp.close()
2701                print ('Vector rigid body saved to: '+filename)
2702        finally:
2703            dlg.Destroy()
2704           
2705    def ReadVectorRB(event):
2706        AtInfo = data['Vector']['AtInfo']
2707        pth = G2G.GetExportPath(G2frame)
2708        dlg = wx.FileDialog(G2frame, 'Choose file to read vector rigid body',
2709            pth, '', 'VRB files (*.vecbody)|*.vecbody',
2710            wx.FD_OPEN)
2711        try:
2712            if dlg.ShowModal() == wx.ID_OK:
2713                filename = dlg.GetPath()
2714                filename = os.path.splitext(filename)[0]+'.vecbody'  # set extension
2715                fp = open(filename,'r')
2716                l = fp.readline().strip()
2717                if 'Name' not in l:
2718                    fp.close()
2719                    G2frame.ErrorDialog('Read Error',
2720                        'File '+filename+' does not start with Name\nFirst line ='
2721                        +l+'\ninvalid file',parent=G2frame)
2722                    return
2723                name = l.split(':')[1].strip()
2724                trans = fp.readline().strip().split(':')[1].split()
2725                vecMag = [float(i) for i in trans]
2726                ntrans = len(trans)
2727                vecs = [[] for i in range(ntrans)]
2728                types = []
2729                l = fp.readline().strip()
2730                while l:
2731                    nums = l.strip().split()
2732                    types.append(nums.pop(0))
2733                    t = types[-1]
2734                    if t not in AtInfo:
2735                        Info = G2elem.GetAtomInfo(t)
2736                        AtInfo[t] = [Info['Drad'],Info['Color']]
2737                    for i in range(ntrans):
2738                        vecs[i].append([float(nums.pop(0)) for j in range(3)])
2739                    l = fp.readline().strip()
2740                fp.close()
2741            else:
2742                return       
2743        finally:
2744            dlg.Destroy()
2745        natoms = len(types)
2746        vecs = [np.array(vecs[i]) for i in range(ntrans)]
2747        rbid = ran.randint(0,sys.maxsize)
2748        namelist = [data['Vector'][key]['RBname'] for key in data['Vector']
2749                        if 'RBname' in data['Vector'][key]]
2750        name = G2obj.MakeUniqueLabel(name,namelist)
2751        data['Vector'][rbid] = {'RBname':name,'VectMag':vecMag,
2752                'rbXYZ':np.zeros((natoms,3)),
2753                'rbRef':[0,1,2,False],'VectRef':ntrans*[False],
2754                'rbTypes':types,
2755                'rbVect':vecs,'useCount':0}
2756        data['RBIds']['Vector'].append(rbid)
2757        UpdateVectorRB()
2758       
2759    def AddResidueRB(event):
2760        global resRBsel
2761        AtInfo = data['Residue']['AtInfo']
2762        macro = getMacroFile('rigid body')
2763        if not macro:
2764            return
2765        macStr = macro.readline()
2766        while macStr:
2767            items = macStr.split()
2768            if 'I' == items[0]:
2769                resRBsel = ran.randint(0,sys.maxsize)
2770                rbName = items[1]
2771                rbTypes = []
2772                rbXYZ = []
2773                rbSeq = []
2774                atNames = []
2775                nAtms,nSeq,nOrig,mRef,nRef = [int(items[i]) for i in [2,3,4,5,6]]
2776                for iAtm in range(nAtms):
2777                    macStr = macro.readline().split()
2778                    atName = macStr[0]
2779                    atType = macStr[1]
2780                    atNames.append(atName)
2781                    rbXYZ.append([float(macStr[i]) for i in [2,3,4]])
2782                    rbTypes.append(atType)
2783                    if atType not in AtInfo:
2784                        Info = G2elem.GetAtomInfo(atType)
2785                        AtInfo[atType] = [Info['Drad'],Info['Color']]
2786                rbXYZ = np.array(rbXYZ)-np.array(rbXYZ[nOrig-1])
2787                for iSeq in range(nSeq):
2788                    macStr = macro.readline().split()
2789                    mSeq = int(macStr[0])
2790                    for jSeq in range(mSeq):
2791                        macStr = macro.readline().split()
2792                        iBeg = int(macStr[0])-1
2793                        iFin = int(macStr[1])-1
2794                        angle = 0.0
2795                        nMove = int(macStr[2])
2796                        iMove = [int(macStr[i])-1 for i in range(3,nMove+3)]
2797                        rbSeq.append([iBeg,iFin,angle,iMove])
2798                namelist = [data['Residue'][key]['RBname'] for key in
2799                           data['Residue'] if 'RBname' in data['Residue'][key]]
2800                rbName = G2obj.MakeUniqueLabel(rbName,namelist)
2801                data['Residue'][resRBsel] = {'RBname':rbName,'rbXYZ':rbXYZ,'rbTypes':rbTypes,
2802                    'atNames':atNames,'rbRef':[nOrig-1,mRef-1,nRef-1,True],'rbSeq':rbSeq,
2803                    'SelSeq':[0,0],'useCount':0,'molCent':None}
2804                data['RBIds']['Residue'].append(resRBsel)
2805                print ('Rigid body '+rbName+' added')
2806            macStr = macro.readline()
2807        macro.close()
2808        UpdateResidueRB()
2809       
2810    def ImportResidueRB():
2811        global resRBsel
2812        AtInfo = data['Residue']['AtInfo']
2813        text,ext = getTextFile()
2814        if not text:
2815            return
2816        resRBsel = ran.randint(0,sys.maxsize)
2817        rbTypes = []
2818        rbXYZ = []
2819        atNames = []
2820        txtStr = text.readline()
2821        if 'xyz' in ext:
2822            txtStr = text.readline()
2823            txtStr = text.readline()
2824        elif 'mol2' in ext:
2825            while 'ATOM' not in txtStr:
2826                txtStr = text.readline()
2827            txtStr = text.readline()
2828        elif 'pdb' in ext:
2829            while 'ATOM' not in txtStr[:6] and 'HETATM' not in txtStr[:6]:
2830                txtStr = text.readline()
2831        items = txtStr.split()
2832        nat = 1
2833        while len(items):
2834            if 'txt' in ext:
2835                atName = items[0]
2836                atType = items[1]
2837                rbXYZ.append([float(items[i]) for i in [2,3,4]])
2838            elif 'xyz' in ext:
2839                atType = items[0]
2840                rbXYZ.append([float(items[i]) for i in [1,2,3]])
2841                atName = '%s%d'%(atType,nat)
2842            elif 'mol2' in ext:
2843                atType = items[1]
2844                atName = items[1]+items[0]
2845                rbXYZ.append([float(items[i]) for i in [2,3,4]])
2846            elif 'pdb' in ext:
2847                atType = items[-1]
2848                if not items[2][-1].isnumeric():
2849                    atName = '%s%d'%(items[2],nat)
2850                else:
2851                    atName = '5s'%items[2]
2852                xyz = txtStr[30:55].split()                   
2853                rbXYZ.append([float(x) for x in xyz])
2854            atNames.append(atName)
2855            rbTypes.append(atType)
2856            if atType not in AtInfo:
2857                Info = G2elem.GetAtomInfo(atType)
2858                AtInfo[atType] = [Info['Drad'],Info['Color']]
2859            txtStr = text.readline()
2860            if 'mol2' in ext and 'BOND' in txtStr:
2861                break
2862            if 'pdb' in ext and ('ATOM' not in txtStr[:6] and 'HETATM' not in txtStr[:6]):
2863                break
2864            items = txtStr.split()
2865            nat += 1
2866        if len(atNames) < 3:
2867            G2G.G2MessageBox(G2frame,'Not enough atoms in rigid body; must be 3 or more')
2868        else:
2869            rbXYZ = np.array(rbXYZ)-np.array(rbXYZ[0])
2870            Xxyz = rbXYZ[1]
2871            X = Xxyz/np.sqrt(np.sum(Xxyz**2))
2872            Yxyz = rbXYZ[2]
2873            Y = Yxyz/np.sqrt(np.sum(Yxyz**2))
2874            Mat = G2mth.getRBTransMat(X,Y)
2875            rbXYZ = np.inner(Mat,rbXYZ).T
2876            name = 'UNKRB'
2877            namelist = [data['Residue'][key]['RBname'] for key in data['Residue']
2878                        if 'RBname' in data['Residue'][key]]
2879            name = G2obj.MakeUniqueLabel(name,namelist)
2880            data['Residue'][resRBsel] = {'RBname':name,'rbXYZ':rbXYZ,'rbTypes':rbTypes,
2881                'atNames':atNames,'rbRef':[0,1,2,False],'rbSeq':[],'SelSeq':[0,0],'useCount':0,'molCent':False}
2882            data['RBIds']['Residue'].append(resRBsel)
2883            print ('Rigid body UNKRB added')
2884        text.close()
2885        UpdateResidueRB()
2886       
2887    def SaveResidueRB():
2888        global resRBsel
2889        pth = G2G.GetExportPath(G2frame)
2890        dlg = wx.FileDialog(G2frame, 'Choose PDB file for Atom XYZ', pth, '', 
2891            'PDB files (*.pdb)|*.pdb',wx.FD_SAVE|wx.FD_OVERWRITE_PROMPT)
2892        try:
2893            if dlg.ShowModal() == wx.ID_OK:
2894                filename = dlg.GetPath()
2895                filename = os.path.splitext(filename)[0]+'.pdb'  # make extension .pdb
2896                File = open(filename,'w')       
2897                rbData =  data['Residue'][resRBsel]
2898                for iat,xyz in enumerate(rbData['rbXYZ']):
2899                    File.write('ATOM %6d  %-4s%3s     1    %8.3f%8.3f%8.3f  1.00  0.00          %2s\n'%(
2900                        iat,rbData['atNames'][iat],rbData['RBname'][:3],xyz[0],xyz[1],xyz[2],rbData['rbTypes'][iat]))
2901                File.close()
2902                print ('Atom XYZ saved to: '+filename)
2903        finally:
2904            dlg.Destroy()
2905       
2906       
2907    def FindNeighbors(Orig,XYZ,atTypes,atNames,AtInfo):
2908        Radii = []
2909        for Atype in atTypes:
2910            Radii.append(AtInfo[Atype][0])
2911        Radii = np.array(Radii)
2912        Neigh = []
2913        Dx = XYZ-XYZ[Orig]
2914        dist = np.sqrt(np.sum(Dx**2,axis=1))
2915        sumR = Radii[Orig]+Radii
2916        IndB = ma.nonzero(ma.masked_greater(dist-0.85*sumR,0.))
2917        for j in IndB[0]:
2918            if j != Orig and atTypes[j] != 'H':
2919                Neigh.append(atNames[j])
2920        return Neigh
2921       
2922    def FindAllNeighbors(XYZ,atTypes,atNames,AtInfo):
2923        NeighDict = {}
2924        for iat,xyz in enumerate(atNames):
2925            NeighDict[atNames[iat]] = FindNeighbors(iat,XYZ,atTypes,atNames,AtInfo)
2926        return NeighDict
2927       
2928    def FindRiding(Orig,Pivot,NeighDict):
2929        riding = [Orig,Pivot]
2930        iAdd = 1
2931        new = True
2932        while new:
2933            newAtms = NeighDict[riding[iAdd]]
2934            for At in newAtms:
2935                new = False
2936                if At not in riding:
2937                    riding.append(At)
2938                    new = True
2939            iAdd += 1
2940            if iAdd < len(riding):
2941                new = True
2942        return riding[2:]
2943                       
2944    def OnDefineTorsSeq(event):
2945        global resRBsel
2946        rbData = data['Residue'][resRBsel]
2947        if not len(rbData):
2948            return
2949        atNames = rbData['atNames']
2950        AtInfo = data['Residue']['AtInfo']
2951        atTypes = rbData['rbTypes']
2952        XYZ = rbData['rbXYZ']
2953        neighDict = FindAllNeighbors(XYZ,atTypes,atNames,AtInfo)
2954        TargList = []           
2955        dlg = wx.SingleChoiceDialog(G2frame,'Select origin atom for torsion sequence','Origin atom',rbData['atNames'])
2956        if dlg.ShowModal() == wx.ID_OK:
2957            Orig = dlg.GetSelection()
2958            TargList = neighDict[atNames[Orig]]
2959        dlg.Destroy()
2960        if not len(TargList):
2961            return
2962        dlg = wx.SingleChoiceDialog(G2frame,'Select pivot atom for torsion sequence','Pivot atom',TargList)
2963        if dlg.ShowModal() == wx.ID_OK:
2964            Piv = atNames.index(TargList[dlg.GetSelection()])
2965            riding = FindRiding(atNames[Orig],atNames[Piv],neighDict)
2966            Riding = []
2967            for atm in riding:
2968                Riding.append(atNames.index(atm))
2969            rbData['rbSeq'].append([Orig,Piv,0.0,Riding])           
2970        dlg.Destroy()
2971        UpdateResidueRB()
2972
2973    def UpdateVectorRB(Scroll=0):
2974        '''Display & edit a selected Vector RB
2975        '''
2976        global resRBsel
2977        def rbNameSizer(rbid,rbData):
2978
2979            def OnRBName(event):
2980                event.Skip()
2981                Obj = event.GetEventObject()
2982                name = Obj.GetValue()
2983                if name == rbData['RBname']: return # no change
2984                namelist = [data['Vector'][key]['RBname'] for key in
2985                        data['Vector'] if 'RBname' in data['Vector'][key]]
2986                name = G2obj.MakeUniqueLabel(name,namelist)
2987                rbData['RBname'] = name
2988                wx.CallAfter(UpdateVectorRB)
2989               
2990            def OnDelRB(event):
2991                Obj = event.GetEventObject()
2992                rbid = Indx[Obj.GetId()]
2993                if rbid in data['Vector']:
2994                    del data['Vector'][rbid]
2995                    data['RBIds']['Vector'].remove(rbid)
2996                    rbData['useCount'] -= 1
2997                wx.CallAfter(UpdateVectorRB)
2998               
2999            def OnPlotRB(event):
3000                G2plt.PlotRigidBody(G2frame,'Vector',AtInfo,rbData,plotDefaults)
3001           
3002            # start of rbNameSizer
3003            nameSizer = wx.BoxSizer(wx.HORIZONTAL)
3004            nameSizer.Add(wx.StaticText(VectorRBDisplay,-1,'Rigid body name: '),0,WACV)
3005            RBname = wx.TextCtrl(VectorRBDisplay,-1,rbData['RBname'],
3006                                     style=wx.TE_PROCESS_ENTER)
3007            RBname.Bind(wx.EVT_LEAVE_WINDOW, OnRBName)
3008            RBname.Bind(wx.EVT_TEXT_ENTER,OnRBName)
3009            RBname.Bind(wx.EVT_KILL_FOCUS,OnRBName)
3010            nameSizer.Add(RBname,0,WACV)
3011            nameSizer.Add((5,0),)
3012            plotRB =  wx.Button(VectorRBDisplay,wx.ID_ANY,'Plot',style=wx.BU_EXACTFIT)
3013            plotRB.Bind(wx.EVT_BUTTON, OnPlotRB)
3014            Indx[plotRB.GetId()] = rbid
3015            nameSizer.Add(plotRB,0,WACV)
3016            nameSizer.Add((5,0),)
3017            if not rbData['useCount']:
3018                delRB = wx.Button(VectorRBDisplay,wx.ID_ANY,"Delete",style=wx.BU_EXACTFIT)
3019                delRB.Bind(wx.EVT_BUTTON, OnDelRB)
3020                Indx[delRB.GetId()] = rbid
3021                nameSizer.Add(delRB,0,WACV)
3022            return nameSizer
3023           
3024        def rbRefAtmSizer(rbid,rbData):
3025           
3026            def OnRefSel(event):
3027                Obj = event.GetEventObject()
3028                iref = Indx[Obj.GetId()]
3029                sel = Obj.GetValue()
3030                rbData['rbRef'][iref] = atNames.index(sel)
3031                FillRefChoice(rbid,rbData)
3032           
3033            refAtmSizer = wx.BoxSizer(wx.HORIZONTAL)
3034            atNames = [name+str(i) for i,name in enumerate(rbData['rbTypes'])]
3035            rbRef = rbData.get('rbRef',[0,1,2,False])
3036            rbData['rbRef'] = rbRef
3037            if rbData['useCount']:
3038                refAtmSizer.Add(wx.StaticText(VectorRBDisplay,-1,
3039                    'Orientation reference atoms A-B-C: %s, %s, %s'%(atNames[rbRef[0]], \
3040                     atNames[rbRef[1]],atNames[rbRef[2]])),0)
3041            else:
3042                refAtmSizer.Add(wx.StaticText(VectorRBDisplay,-1,
3043                    'Orientation reference atoms A-B-C: '),0,WACV)
3044                for i in range(3):
3045                    choices = [atNames[j] for j in refChoice[rbid][i]]
3046                    refSel = wx.ComboBox(VectorRBDisplay,-1,value='',
3047                        choices=choices,style=wx.CB_READONLY|wx.CB_DROPDOWN)
3048                    refSel.SetValue(atNames[rbRef[i]])
3049                    refSel.Bind(wx.EVT_COMBOBOX, OnRefSel)
3050                    Indx[refSel.GetId()] = i
3051                    refAtmSizer.Add(refSel,0,WACV)
3052                refHelpInfo = '''
3053* The "Orientation Reference" control defines the Cartesian
3054axes for rigid bodies with the three atoms, A, B and C.
3055The vector from B to A defines the x-axis and the y axis is placed
3056in the plane defined by B to A and C to A. A,B,C must not be collinear.
3057'''
3058                hlp = G2G.HelpButton(VectorRBDisplay,refHelpInfo,wrap=400)
3059                refAtmSizer.Add(hlp,0,wx.LEFT|wx.RIGHT|wx.ALIGN_CENTER_VERTICAL,2)
3060            return refAtmSizer
3061                       
3062        def rbVectMag(rbid,imag,rbData):
3063           
3064            def OnRBVectorMag(event):
3065                event.Skip()
3066                Obj = event.GetEventObject()
3067                rbid,imag = Indx[Obj.GetId()]
3068                try:
3069                    val = float(Obj.GetValue())
3070                    if val <= 0.:
3071                        raise ValueError
3072                    rbData['VectMag'][imag] = val
3073                except ValueError:
3074                    pass
3075                Obj.SetValue('%8.4f'%(val))
3076                wx.CallAfter(UpdateVectorRB,VectorRB.GetScrollPos(wx.VERTICAL))
3077                G2plt.PlotRigidBody(G2frame,'Vector',AtInfo,data['Vector'][rbid],plotDefaults)
3078               
3079            def OnRBVectorRef(event):
3080                Obj = event.GetEventObject()
3081                rbid,imag = Indx[Obj.GetId()]
3082                rbData['VectRef'][imag] = Obj.GetValue()
3083                       
3084            magSizer = wx.BoxSizer(wx.HORIZONTAL)
3085            magSizer.Add(wx.StaticText(VectorRBDisplay,-1,'Translation magnitude: '),0,WACV)
3086            magValue = wx.TextCtrl(VectorRBDisplay,-1,'%8.4f'%(rbData['VectMag'][imag]))
3087            Indx[magValue.GetId()] = [rbid,imag]
3088            magValue.Bind(wx.EVT_TEXT_ENTER,OnRBVectorMag)
3089            magValue.Bind(wx.EVT_KILL_FOCUS,OnRBVectorMag)
3090            magSizer.Add(magValue,0,WACV)
3091            magSizer.Add((5,0),)
3092            magref = wx.CheckBox(VectorRBDisplay,label=' Refine?') 
3093            magref.SetValue(rbData['VectRef'][imag])
3094            magref.Bind(wx.EVT_CHECKBOX,OnRBVectorRef)
3095            Indx[magref.GetId()] = [rbid,imag]
3096            magSizer.Add(magref,0,WACV)
3097            return magSizer
3098           
3099        def rbVectors(rbid,imag,mag,XYZ,rbData):
3100
3101            def TypeSelect(event):
3102                AtInfo = data['Vector']['AtInfo']
3103                r,c = event.GetRow(),event.GetCol()
3104                if vecGrid.GetColLabelValue(c) == 'Type':
3105                    PE = G2elemGUI.PickElement(G2frame,oneOnly=True)
3106                    if PE.ShowModal() == wx.ID_OK:
3107                        if PE.Elem != 'None':
3108                            El = PE.Elem.strip().lower().capitalize()
3109                            if El not in AtInfo:
3110                                Info = G2elem.GetAtomInfo(El)
3111                                AtInfo[El] = [Info['Drad'],Info['Color']]
3112                            rbData['rbTypes'][r] = El
3113                            vecGrid.SetCellValue(r,c,El)
3114                    PE.Destroy()
3115                wx.CallAfter(UpdateVectorRB,VectorRB.GetScrollPos(wx.VERTICAL))
3116
3117            def ChangeCell(event):
3118                r,c =  event.GetRow(),event.GetCol()
3119                if r >= 0 and (0 <= c < 3):
3120                    try:
3121                        val = float(vecGrid.GetCellValue(r,c))
3122                        rbData['rbVect'][imag][r][c] = val
3123                    except ValueError:
3124                        pass
3125                G2plt.PlotRigidBody(G2frame,'Vector',AtInfo,data['Vector'][rbid],plotDefaults)
3126                wx.CallAfter(UpdateVectorRB,VectorRB.GetScrollPos(wx.VERTICAL))
3127
3128            vecSizer = wx.BoxSizer()
3129            Types = 3*[wg.GRID_VALUE_FLOAT+':10,5',]+[wg.GRID_VALUE_STRING,]+3*[wg.GRID_VALUE_FLOAT+':10,5',]
3130            colLabels = ['Vector x','Vector y','Vector z','Type','Cart x','Cart y','Cart z']
3131            table = []
3132            rowLabels = []
3133            atNames = []
3134            for ivec,xyz in enumerate(rbData['rbVect'][imag]):
3135                table.append(list(xyz)+[rbData['rbTypes'][ivec],]+list(XYZ[ivec]))
3136                rowLabels.append(str(ivec))
3137                atNames.append(rbData['rbTypes'][ivec]+str(ivec))
3138            rbData['atNames'] = atNames
3139            vecTable = G2G.Table(table,rowLabels=rowLabels,colLabels=colLabels,types=Types)
3140            vecGrid = G2G.GSGrid(VectorRBDisplay)
3141            vecGrid.SetTable(vecTable, True)
3142            if 'phoenix' in wx.version():
3143                vecGrid.Bind(wg.EVT_GRID_CELL_CHANGED, ChangeCell)
3144            else:
3145                vecGrid.Bind(wg.EVT_GRID_CELL_CHANGE, ChangeCell)
3146            if not imag:
3147                vecGrid.Bind(wg.EVT_GRID_CELL_LEFT_DCLICK, TypeSelect)
3148            attr = wx.grid.GridCellAttr()
3149            attr.IncRef()
3150            attr.SetEditor(G2G.GridFractionEditor(vecGrid))
3151            for c in range(3):
3152                attr.IncRef()
3153                vecGrid.SetColAttr(c, attr)
3154            for row in range(vecTable.GetNumberRows()):
3155                if imag:
3156                    vecGrid.SetCellStyle(row,3,VERY_LIGHT_GREY,True)                   
3157                for col in [4,5,6]:
3158                    vecGrid.SetCellStyle(row,col,VERY_LIGHT_GREY,True)
3159#            vecGrid.SetScrollRate(0,0)
3160            vecGrid.AutoSizeColumns(False)
3161            vecSizer.Add(vecGrid)
3162            return vecSizer
3163       
3164        def FillRefChoice(rbid,rbData):
3165            choiceIds = [i for i in range(len(rbData['rbTypes']))]
3166           
3167            rbRef = rbData.get('rbRef',[-1,-1,-1,False])
3168            for i in range(3):
3169                if rbRef[i] in choiceIds: choiceIds.remove(rbRef[i])
3170            refChoice[rbid] = [choiceIds[:],choiceIds[:],choiceIds[:]]
3171            for i in range(3):
3172                refChoice[rbid][i].append(rbRef[i])
3173                refChoice[rbid][i].sort()
3174               
3175        def OnRBSelect(event):
3176            global resRBsel
3177            sel = rbSelect.GetSelection()
3178            if sel == 0: return # 1st entry is blank
3179            rbname = rbchoice[sel-1]
3180            resRBsel = RBnames[rbname]
3181            wx.CallLater(100,UpdateVectorRB)
3182           
3183        # beginning of UpdateVectorRB
3184        AtInfo = data['Vector']['AtInfo']
3185        refChoice = {}
3186        RefObjs = []
3187
3188        GS = VectorRBDisplay.GetSizer()
3189        if GS: 
3190            try:        #get around a c++ error in wx 4.0; doing is again seems to be OK
3191                GS.Clear(True)
3192            except:
3193                GS.Clear(True)
3194       
3195        RBnames = {}
3196        for rbid in data['RBIds']['Vector']:
3197            RBnames.update({data['Vector'][rbid]['RBname']:rbid,})
3198        if not RBnames:
3199            return
3200        rbchoice = list(RBnames.keys())
3201        rbchoice.sort()
3202        if GS: 
3203            VectorRBSizer = GS
3204        else:
3205            VectorRBSizer = wx.BoxSizer(wx.VERTICAL)
3206        if resRBsel not in data['RBIds']['Vector']:
3207            resRBsel = RBnames[rbchoice[0]]
3208        if len(RBnames) > 1:
3209            selSizer = wx.BoxSizer(wx.HORIZONTAL)
3210            selSizer.Add(wx.StaticText(VectorRBDisplay,label=' Select rigid body to view:'),0)
3211            rbSelect = wx.ComboBox(VectorRBDisplay,choices=['']+rbchoice)
3212            name = data['Vector'][resRBsel]['RBname']
3213            rbSelect.SetSelection(1+rbchoice.index(name))
3214            rbSelect.Bind(wx.EVT_COMBOBOX,OnRBSelect)
3215            selSizer.Add(rbSelect,0)
3216            VectorRBSizer.Add(selSizer,0)
3217        rbData = data['Vector'][resRBsel]
3218        if 'DELETED' in str(G2frame.GetStatusBar()):   #seems to be no other way to do this (wx bug)
3219            if GSASIIpath.GetConfigValue('debug'):
3220                print ('DBG_wx error: Rigid Body/Status not cleanly deleted after Refine')
3221            return
3222        SetStatusLine(' You may use e.g. "c60" or "s60" for a vector entry')
3223        FillRefChoice(resRBsel,rbData)
3224        VectorRBSizer.Add(rbNameSizer(resRBsel,rbData),0)
3225        VectorRBSizer.Add(rbRefAtmSizer(resRBsel,rbData),0)
3226        XYZ = np.array([[0.,0.,0.] for Ty in rbData['rbTypes']])
3227        for imag,mag in enumerate(rbData['VectMag']):
3228            XYZ += mag*rbData['rbVect'][imag]
3229            VectorRBSizer.Add(rbVectMag(rbid,imag,rbData),0)
3230            VectorRBSizer.Add(rbVectors(rbid,imag,mag,XYZ,rbData),0)
3231        VectorRBSizer.Add((5,5),0)
3232        data['Vector'][rbid]['rbXYZ'] = XYZ       
3233       
3234        VectorRBSizer.Add((5,25),)
3235        VectorRBSizer.Layout()   
3236        VectorRBDisplay.SetSizer(VectorRBSizer,True)
3237        VectorRBDisplay.SetAutoLayout(True)
3238        Size = VectorRBSizer.GetMinSize()
3239        VectorRBDisplay.SetSize(Size)
3240       
3241        Size[0] += 40
3242        Size[1] = max(Size[1],450) + 20
3243        VectorRB.SetSize(Size)
3244        VectorRB.SetScrollbars(10,10,Size[0]/10-4,Size[1]/10-1)
3245        G2frame.dataWindow.SendSizeEvent()
3246       
3247        VectorRBDisplay.Show()
3248       
3249       
3250    def UpdateResidueRB():
3251        '''Draw the contents of the Residue Rigid Body tab for Rigid Bodies tree entry
3252        '''
3253        global resRBsel
3254        def rbNameSizer(rbid,rbData):
3255           
3256            def OnRBName(event):
3257                event.Skip()
3258                Obj = event.GetEventObject()
3259                name = Obj.GetValue()
3260                if name == rbData['RBname']: return # no change
3261                namelist = [data['Residue'][key]['RBname'] for key in
3262                        data['Residue'] if 'RBname' in data['Residue'][key]]
3263                name = G2obj.MakeUniqueLabel(name,namelist)
3264                rbData['RBname'] = name
3265                wx.CallAfter(UpdateResidueRB)
3266               
3267            def OnDelRB(event):
3268                Obj = event.GetEventObject()
3269                rbid = Indx[Obj.GetId()]
3270                if rbid in data['Residue']: 
3271                    del data['Residue'][rbid]
3272                    data['RBIds']['Residue'].remove(rbid)
3273                wx.CallAfter(UpdateResidueRB)
3274               
3275            def OnStripH(event):
3276                Obj = event.GetEventObject()
3277                rbid = Indx[Obj.GetId()]
3278                if rbid in data['Residue']:
3279                    newNames = []
3280                    newTypes = []
3281                    newXYZ = []
3282                    for i,atype in enumerate(rbData['rbTypes']):
3283                        if atype != 'H':
3284                            newNames.append(rbData['atNames'][i])
3285                            newTypes.append(rbData['rbTypes'][i])
3286                            newXYZ.append(rbData['rbXYZ'][i])
3287                    rbData['atNames'] = newNames
3288                    rbData['rbTypes'] = newTypes
3289                    rbData['rbXYZ'] = newXYZ
3290                G2plt.PlotRigidBody(G2frame,'Residue',AtInfo,rbData,plotDefaults)
3291                wx.CallAfter(UpdateResidueRB)
3292
3293            def OnPlotRB(event):
3294                G2plt.PlotRigidBody(G2frame,'Residue',AtInfo,rbData,plotDefaults)
3295               
3296            # start of rbNameSizer
3297            nameSizer = wx.BoxSizer(wx.HORIZONTAL)
3298            nameSizer.Add(wx.StaticText(ResidueRBDisplay,-1,'Residue name: '),0,WACV)
3299            RBname = wx.TextCtrl(ResidueRBDisplay,-1,rbData['RBname'],
3300                                     style=wx.TE_PROCESS_ENTER)
3301            RBname.Bind(wx.EVT_LEAVE_WINDOW, OnRBName)
3302            RBname.Bind(wx.EVT_TEXT_ENTER,OnRBName)
3303            RBname.Bind(wx.EVT_KILL_FOCUS,OnRBName)
3304            nameSizer.Add(RBname,0,WACV)
3305            nameSizer.Add((5,0),)
3306            plotRB =  wx.Button(ResidueRBDisplay,wx.ID_ANY,'Plot',style=wx.BU_EXACTFIT)
3307            plotRB.Bind(wx.EVT_BUTTON, OnPlotRB)
3308            Indx[plotRB.GetId()] = rbid
3309            nameSizer.Add(plotRB,0,WACV)
3310            nameSizer.Add((5,0),)
3311            if not rbData['useCount']:
3312                delRB = wx.Button(ResidueRBDisplay,wx.ID_ANY,"Delete",style=wx.BU_EXACTFIT)
3313                delRB.Bind(wx.EVT_BUTTON, OnDelRB)
3314                Indx[delRB.GetId()] = rbid
3315                nameSizer.Add(delRB,0,WACV)
3316                if 'H'  in rbData['rbTypes']:
3317                    stripH = wx.Button(ResidueRBDisplay,wx.ID_ANY,"Strip H-atoms",style=wx.BU_EXACTFIT)
3318                    stripH.Bind(wx.EVT_BUTTON, OnStripH)
3319                    Indx[stripH.GetId()] = rbid
3320                    nameSizer.Add(stripH,0,WACV)
3321            nameSizer.Add(wx.StaticText(ResidueRBDisplay,-1,'  body type #'+
3322                                        str(data['RBIds']['Residue'].index(rbid))),0,WACV)
3323            return nameSizer
3324           
3325        def rbResidues(rbid,rbData):
3326           
3327            def TypeSelect(event):
3328                AtInfo = data['Residue']['AtInfo']
3329                r,c = event.GetRow(),event.GetCol()
3330                if resGrid.GetColLabelValue(c) == 'Type':
3331                    PE = G2elemGUI.PickElement(G2frame,oneOnly=True)
3332                    if PE.ShowModal() == wx.ID_OK:
3333                        if PE.Elem != 'None':
3334                            El = PE.Elem.strip().lower().capitalize()
3335                            if El not in AtInfo:
3336                                Info = G2elem.GetAtomInfo(El)
3337                                AtInfo[El] = [Info['Drad'],Info['Color']]
3338                            rbData['rbTypes'][r] = El
3339                            resGrid.SetCellValue(r,c,El)
3340                    PE.Destroy()
3341
3342            def ChangeCell(event):
3343                r,c =  event.GetRow(),event.GetCol()
3344                if c == 0:
3345                    rbData['atNames'][r] = resGrid.GetCellValue(r,c)
3346                if r >= 0 and (2 <= c <= 4):
3347                    try:
3348                        val = float(resGrid.GetCellValue(r,c))
3349                        rbData['rbXYZ'][r][c-2] = val
3350                    except ValueError:
3351                        pass
3352                G2plt.PlotRigidBody(G2frame,'Residue',AtInfo,rbData,plotDefaults)
3353                       
3354            def OnRefSel(event):
3355                Obj = event.GetEventObject()
3356                iref,res,jref = Indx[Obj.GetId()]
3357                sel = Obj.GetValue()
3358                ind = atNames.index(sel)
3359                if rbData['rbTypes'][ind] == 'H':
3360                    G2G.G2MessageBox(G2frame,'You should not select an H-atom for rigid body orientation')
3361                rbData['rbRef'][iref] = ind
3362                FillRefChoice(rbid,rbData)
3363                for i,ref in enumerate(RefObjs[jref]):
3364                    ref.SetItems([atNames[j] for j in refChoice[rbid][i]])
3365                    ref.SetValue(atNames[rbData['rbRef'][i]])                   
3366                rbXYZ = rbData['rbXYZ']
3367                if not iref:     #origin change
3368                    rbXYZ -= rbXYZ[ind]
3369                Xxyz = rbXYZ[rbData['rbRef'][1]]
3370                X = Xxyz/np.sqrt(np.sum(Xxyz**2))
3371                Yxyz = rbXYZ[rbData['rbRef'][2]]
3372                Y = Yxyz/np.sqrt(np.sum(Yxyz**2))
3373                Mat = G2mth.getRBTransMat(X,Y)
3374                rbXYZ = np.inner(Mat,rbXYZ).T
3375                rbData['rbXYZ'] = rbXYZ
3376                rbData['molCent'] = False
3377                res.ClearSelection()
3378                resTable = res.GetTable()
3379                for r in range(res.GetNumberRows()):
3380                    row = resTable.GetRowValues(r)
3381                    row[2:4] = rbXYZ[r]
3382                    resTable.SetRowValues(r,row)
3383                res.ForceRefresh()
3384                G2plt.PlotRigidBody(G2frame,'Residue',AtInfo,rbData,plotDefaults)
3385               
3386            def OnMolCent(event):
3387                rbData['molCent'] = not rbData['molCent']
3388                if rbData['molCent']:
3389                    Obj = event.GetEventObject()
3390                    res = Indx[Obj.GetId()]
3391                    rbXYZ = rbData['rbXYZ']
3392                    rbCent = np.array([np.sum(rbXYZ[:,0]),np.sum(rbXYZ[:,1]),np.sum(rbXYZ[:,2])])/rbXYZ.shape[0]
3393                    rbXYZ -= rbCent
3394                    rbData['rbXYZ'] = rbXYZ
3395                    res.ClearSelection()
3396                    resTable = res.GetTable()
3397                    for r in range(res.GetNumberRows()):
3398                        row = resTable.GetRowValues(r)
3399                        row[2:4] = rbXYZ[r]
3400                        resTable.SetRowValues(r,row)
3401                    res.ForceRefresh()
3402                    G2plt.PlotRigidBody(G2frame,'Residue',AtInfo,rbData,plotDefaults)
3403                   
3404            def OnCycleXYZ(event):
3405                Obj = event.GetEventObject()
3406                res = Indx[Obj.GetId()]
3407                rbXYZ = rbData['rbXYZ']
3408                resTable = res.GetTable()
3409                for r in range(res.GetNumberRows()):
3410                    rbXYZ[r] = np.roll(rbXYZ[r],1)
3411                    row = resTable.GetRowValues(r)
3412                    row[2:4] = rbXYZ[r]
3413                    resTable.SetRowValues(r,row)
3414                res.ForceRefresh()
3415                G2plt.PlotRigidBody(G2frame,'Residue',AtInfo,rbData,plotDefaults)
3416               
3417            Types = 2*[wg.GRID_VALUE_STRING,]+3*[wg.GRID_VALUE_FLOAT+':10,5',]
3418            colLabels = ['Name','Type','Cart x','Cart y','Cart z']
3419            table = []
3420            rowLabels = []
3421            for ivec,xyz in enumerate(rbData['rbXYZ']):
3422                table.append([rbData['atNames'][ivec],]+[rbData['rbTypes'][ivec],]+list(xyz))
3423                rowLabels.append(str(ivec))
3424            vecTable = G2G.Table(table,rowLabels=rowLabels,colLabels=colLabels,types=Types)
3425            resGrid = G2G.GSGrid(ResidueRBDisplay)
3426            Indx[resGrid.GetId()] = rbid
3427            resList.append(resGrid)
3428            resGrid.SetTable(vecTable, True)
3429            if 'phoenix' in wx.version():
3430                resGrid.Bind(wg.EVT_GRID_CELL_CHANGED, ChangeCell)
3431            else:
3432                resGrid.Bind(wg.EVT_GRID_CELL_CHANGE, ChangeCell)
3433            resGrid.Bind(wg.EVT_GRID_CELL_LEFT_DCLICK, TypeSelect)
3434            for c in range(2,5):
3435                attr = wx.grid.GridCellAttr()
3436                attr.IncRef()
3437                attr.SetEditor(G2G.GridFractionEditor(resGrid))
3438                resGrid.SetColAttr(c, attr)
3439            for row in range(vecTable.GetNumberRows()):
3440                resGrid.SetReadOnly(row,1,True)
3441            resGrid.AutoSizeColumns(False)
3442            vecSizer = wx.BoxSizer()
3443            vecSizer.Add(resGrid)
3444           
3445            refAtmSizer = wx.BoxSizer(wx.HORIZONTAL)
3446            atNames = rbData['atNames']
3447            rbRef = rbData['rbRef']
3448            if rbData['rbRef'][3] or rbData['useCount']:
3449                refAtmSizer.Add(wx.StaticText(ResidueRBDisplay,-1,
3450                    'Orientation reference non-H atoms A-B-C: %s, %s, %s'%(atNames[rbRef[0]], \
3451                     atNames[rbRef[1]],atNames[rbRef[2]])),0)
3452            else:
3453                refAtmSizer.Add(wx.StaticText(ResidueRBDisplay,-1,
3454                    'Orientation reference non-H atoms A-B-C: '),0,WACV)
3455                refObj = [0,0,0]
3456                for i in range(3):
3457                    choices = [atNames[j] for j in refChoice[rbid][i]]
3458                    refSel = wx.ComboBox(ResidueRBDisplay,-1,value='',
3459                        choices=choices,style=wx.CB_READONLY|wx.CB_DROPDOWN)
3460                    refSel.SetValue(atNames[rbRef[i]])
3461                    refSel.Bind(wx.EVT_COMBOBOX, OnRefSel)
3462                    Indx[refSel.GetId()] = [i,resGrid,len(RefObjs)]
3463                    refObj[i] = refSel
3464                    refAtmSizer.Add(refSel,0,WACV)
3465                RefObjs.append(refObj)
3466                cycleXYZ = wx.Button(ResidueRBDisplay,label=' Cycle XYZ?')
3467                cycleXYZ.Bind(wx.EVT_BUTTON,OnCycleXYZ)
3468                Indx[cycleXYZ.GetId()] = resGrid
3469                refAtmSizer.Add(cycleXYZ,0,WACV)
3470                if 'molCent' not in rbData: rbData['molCent'] = False           #patch
3471                molcent = wx.Button(ResidueRBDisplay,label=' Center RB?')
3472                molcent.Bind(wx.EVT_BUTTON,OnMolCent)
3473                Indx[molcent.GetId()] = resGrid
3474                refAtmSizer.Add(molcent,0,WACV)
3475                refHelpInfo = '''
3476* The "Orientation Reference" control defines the Cartesian
3477axes for rigid bodies with the three atoms, A, B and C.
3478The vector from B to A defines the x-axis and the y axis is placed
3479in the plane defined by B to A and C to A. A,B,C must not be collinear.
3480 
3481%%* The origin is at A unless the "Center RB?" button is pressed.
3482
3483%%* The 'Cycle XYZ' button will permute the rigid body XYZ coordinates so
3484XYZ --> ZXY. Repeat if needed.
3485
3486%%* The "Center RB?" button will shift the origin of the
3487rigid body to be the midpoint of all atoms in the body (not mass weighted).
3488'''
3489                hlp = G2G.HelpButton(ResidueRBDisplay,refHelpInfo,wrap=400)
3490                refAtmSizer.Add(hlp,0,wx.LEFT|wx.RIGHT|wx.ALIGN_CENTER_VERTICAL,2)
3491           
3492            mainSizer = wx.BoxSizer(wx.VERTICAL)
3493            mainSizer.Add(refAtmSizer)
3494            mainSizer.Add(vecSizer)
3495            return mainSizer
3496           
3497        def Add2SeqSizer(seqSizer,angSlide,rbid,iSeq,Seq,atNames):
3498           
3499            def ChangeAngle(event):
3500                event.Skip()
3501                Obj = event.GetEventObject()
3502                rbid,Seq = Indx[Obj.GetId()][:2]
3503                val = Seq[2]
3504                try:
3505                    val = float(Obj.GetValue())
3506                    Seq[2] = val
3507                except ValueError:
3508                    pass
3509                Obj.SetValue('%8.2f'%(val))
3510                G2plt.PlotRigidBody(G2frame,'Residue',AtInfo,data['Residue'][rbid],plotDefaults)
3511               
3512            def OnRadBtn(event):
3513                Obj = event.GetEventObject()
3514                Seq,iSeq,angId = Indx[Obj.GetId()]
3515                data['Residue'][rbid]['SelSeq'] = [iSeq,angId]
3516                angSlide.SetValue(int(100*Seq[2]))
3517               
3518            def OnDelBtn(event):
3519                Obj = event.GetEventObject()
3520                rbid,Seq = Indx[Obj.GetId()]
3521                data['Residue'][rbid]['rbSeq'].remove(Seq)       
3522                wx.CallAfter(UpdateResidueRB)
3523           
3524            iBeg,iFin,angle,iMove = Seq
3525            ang = wx.TextCtrl(ResidueRBDisplay,wx.ID_ANY,
3526                    '%8.2f'%(angle),size=(70,-1),style=wx.TE_PROCESS_ENTER)
3527            if not iSeq:
3528                radBt = wx.RadioButton(ResidueRBDisplay,wx.ID_ANY,
3529                                           '',style=wx.RB_GROUP)
3530                data['Residue'][rbid]['SelSeq'] = [iSeq,ang.GetId()]
3531                radBt.SetValue(True)
3532            else:
3533                radBt = wx.RadioButton(ResidueRBDisplay,wx.ID_ANY,'')
3534            radBt.Bind(wx.EVT_RADIOBUTTON,OnRadBtn)                   
3535            seqSizer.Add(radBt)
3536            delBt =  wx.Button(ResidueRBDisplay,wx.ID_ANY,'Del',
3537                                style=wx.BU_EXACTFIT)
3538            delBt.Bind(wx.EVT_BUTTON,OnDelBtn)
3539            seqSizer.Add(delBt)
3540            bond = wx.StaticText(ResidueRBDisplay,wx.ID_ANY,
3541                        '%s %s'%(atNames[iBeg],atNames[iFin]),size=(50,20))
3542            seqSizer.Add(bond,0,WACV)
3543            Indx[radBt.GetId()] = [Seq,iSeq,ang.GetId()]
3544            Indx[delBt.GetId()] = [rbid,Seq]
3545            Indx[ang.GetId()] = [rbid,Seq,ang]
3546            ang.Bind(wx.EVT_TEXT_ENTER,ChangeAngle)
3547            ang.Bind(wx.EVT_KILL_FOCUS,ChangeAngle)
3548            seqSizer.Add(ang,0,WACV)
3549            atms = ''
3550            for i in iMove:   
3551                atms += ' %s,'%(atNames[i])
3552            moves = wx.StaticText(ResidueRBDisplay,wx.ID_ANY,
3553                            atms[:-1],size=(200,20))
3554            seqSizer.Add(moves,1,wx.EXPAND|wx.RIGHT)
3555            return seqSizer
3556           
3557        def SlideSizer():
3558           
3559            def OnSlider(event):
3560                Obj = event.GetEventObject()
3561                rbData = Indx[Obj.GetId()]
3562                iSeq,angId = rbData['SelSeq']
3563                val = float(Obj.GetValue())/100.
3564                rbData['rbSeq'][iSeq][2] = val
3565                Indx[angId][2].SetValue('%8.2f'%(val))
3566                G2plt.PlotRigidBody(G2frame,'Residue',AtInfo,rbData,plotDefaults)
3567           
3568            slideSizer = wx.BoxSizer(wx.HORIZONTAL)
3569            slideSizer.Add(wx.StaticText(ResidueRBDisplay,-1,'Selected torsion angle:'),0)
3570            iSeq,angId = rbData['SelSeq']
3571            angSlide = wx.Slider(ResidueRBDisplay,-1,
3572                int(100*rbData['rbSeq'][iSeq][2]),0,36000,size=(200,20),
3573                style=wx.SL_HORIZONTAL)
3574            angSlide.Bind(wx.EVT_SLIDER, OnSlider)
3575            Indx[angSlide.GetId()] = rbData
3576            slideSizer.Add(angSlide,0)           
3577            return slideSizer,angSlide
3578           
3579        def FillRefChoice(rbid,rbData):
3580            choiceIds = [i for i in range(len(rbData['atNames']))]
3581            for seq in rbData['rbSeq']:
3582                for i in seq[3]:
3583                    try:
3584                        choiceIds.remove(i)
3585                    except ValueError:
3586                        pass
3587            rbRef = rbData['rbRef']
3588            for i in range(3):
3589                try:
3590                    choiceIds.remove(rbRef[i])
3591                except ValueError:
3592                    pass
3593            refChoice[rbid] = [choiceIds[:],choiceIds[:],choiceIds[:]]
3594            for i in range(3):
3595                refChoice[rbid][i].append(rbRef[i])
3596                refChoice[rbid][i].sort()
3597               
3598        def OnRBSelect(event):
3599            global resRBsel
3600            sel = rbSelect.GetSelection()
3601            if sel == 0: return # 1st entry is blank
3602            rbname = rbchoice[sel-1]
3603            resRBsel = RBnames[rbname]
3604            G2plt.PlotRigidBody(G2frame,'Residue',AtInfo,rbData,plotDefaults)
3605            wx.CallLater(100,UpdateResidueRB)
3606           
3607        #----- beginning of UpdateResidueRB -----------------------------------------------
3608        AtInfo = data['Residue']['AtInfo']
3609        refChoice = {}
3610        RefObjs = []
3611
3612        GS = ResidueRBDisplay.GetSizer()
3613        if GS: 
3614            try:        #get around a c++ error in wx 4.0; doing is again seems to be OK
3615                GS.Clear(True)
3616            except:
3617                GS.Clear(True)
3618       
3619        RBnames = {}
3620        for rbid in data['RBIds']['Residue']:
3621            RBnames.update({data['Residue'][rbid]['RBname']:rbid,})
3622        if not RBnames:
3623            return
3624        rbchoice = list(RBnames.keys())
3625        rbchoice.sort()
3626        if GS: 
3627            ResidueRBSizer = GS
3628        else:
3629            ResidueRBSizer = wx.BoxSizer(wx.VERTICAL)
3630        if resRBsel not in data['RBIds']['Residue']:
3631            resRBsel = RBnames[rbchoice[0]]
3632        if len(RBnames) > 1:
3633            selSizer = wx.BoxSizer(wx.HORIZONTAL)
3634            selSizer.Add(wx.StaticText(ResidueRBDisplay,
3635                                label=' Select residue to view:'),0)
3636            rbSelect = wx.ComboBox(ResidueRBDisplay,choices=['']+rbchoice)
3637            name = data['Residue'][resRBsel]['RBname']
3638            rbSelect.SetSelection(1+rbchoice.index(name))
3639            rbSelect.Bind(wx.EVT_COMBOBOX,OnRBSelect)
3640            selSizer.Add(rbSelect,0)
3641            ResidueRBSizer.Add(selSizer,0)
3642        rbData = data['Residue'][resRBsel]
3643        FillRefChoice(resRBsel,rbData)
3644        ResidueRBSizer.Add(rbNameSizer(resRBsel,rbData),0)
3645        ResidueRBSizer.Add(rbResidues(resRBsel,rbData),0)
3646        if len(rbData['rbSeq']):
3647            ResidueRBSizer.Add((-1,15),0)
3648            slideSizer,angSlide = SlideSizer()
3649            seqSizer = wx.FlexGridSizer(0,5,4,8)
3650            for lbl in 'Sel','','Bond','Angle','Riding atoms':
3651                seqSizer.Add(wx.StaticText(ResidueRBDisplay,wx.ID_ANY,lbl))
3652            ResidueRBSizer.Add(seqSizer)
3653#            for iSeq,Seq in enumerate(rbData['rbSeq']):
3654#                ResidueRBSizer.Add(SeqSizer(angSlide,resRBsel,iSeq,Seq,rbData['atNames']))
3655            for iSeq,Seq in enumerate(rbData['rbSeq']):
3656                Add2SeqSizer(seqSizer,angSlide,resRBsel,iSeq,Seq,rbData['atNames'])
3657            ResidueRBSizer.Add(slideSizer,)
3658
3659        ResidueRBSizer.Add((5,25),)
3660        ResidueRBSizer.Layout()   
3661        ResidueRBDisplay.SetSizer(ResidueRBSizer,True)
3662        ResidueRBDisplay.SetAutoLayout(True)
3663        Size = ResidueRBSizer.GetMinSize()
3664        ResidueRBDisplay.SetSize(Size)
3665       
3666        Size[0] += 40
3667        Size[1] = max(Size[1],450) + 20
3668        ResidueRB.SetSize(Size)
3669        ResidueRB.SetScrollbars(10,10,Size[0]/10-4,Size[1]/10-1)
3670        G2frame.dataWindow.SendSizeEvent()
3671       
3672        ResidueRBDisplay.Show()
3673       
3674    def SetStatusLine(text):
3675        G2frame.GetStatusBar().SetStatusText(text,1)                                     
3676
3677    #### UpdateRigidBodies starts here =========
3678    global resList,resRBsel           
3679    if not data.get('RBIds') or not data:
3680        data.update({'Vector':{'AtInfo':{}},'Residue':{'AtInfo':{}},
3681            'RBIds':{'Vector':[],'Residue':[]}})       #empty/bad dict - fill it
3682    Indx = {}
3683    resList = []
3684    plotDefaults = {'oldxy':[0.,0.],'Quaternion':[0.,0.,0.,1.],'cameraPos':30.,'viewDir':[0,0,1],}
3685    G2frame.rbBook = G2G.GSNoteBook(parent=G2frame.dataWindow)
3686    G2frame.dataWindow.GetSizer().Add(G2frame.rbBook,1,wx.ALL|wx.EXPAND)
3687    VectorRB = wx.ScrolledWindow(G2frame.rbBook)
3688    VectorRBDisplay = wx.Panel(VectorRB)
3689    G2frame.rbBook.AddPage(VectorRB,'Vector rigid bodies')
3690    ResidueRB = wx.ScrolledWindow(G2frame.rbBook)
3691    ResidueRBDisplay = wx.Panel(ResidueRB)
3692    G2frame.rbBook.AddPage(ResidueRB,'Residue rigid bodies')
3693    # vector RBs are not too common, so select Residue as the default when one is present
3694    if len(data['RBIds']['Residue']) > 0 and len(data['RBIds']['Vector']) == 0:
3695        G2frame.rbBook.ChangeSelection(1)
3696        OnPageChanged(None)
3697    G2gd.SetDataMenuBar(G2frame,G2frame.dataWindow.RigidBodyMenu)
3698    SetStatusLine('')
3699    UpdateVectorRB()
3700    G2frame.rbBook.Bind(wx.aui.EVT_AUINOTEBOOK_PAGE_CHANGED, OnPageChanged)
3701    wx.CallAfter(OnPageChanged,None)
3702
3703def ShowIsoDistortCalc(G2frame,phase=None):
3704    '''Compute the ISODISTORT mode values from the current coordinates.
3705    Called in response to the (Phase/Atoms tab) AtomCompute or
3706    Constraints/Edit Constr. "Show ISODISTORT modes" menu item, which
3707    should be enabled only when Phase['ISODISTORT'] is defined.
3708    '''
3709    def _onClose(event):
3710        dlg.EndModal(wx.ID_CANCEL)
3711    def fmtHelp(item,fullname):
3712        helptext = "A new variable"
3713        if item[-3]:
3714            helptext += " named "+str(item[-3])
3715        helptext += " is a linear combination of the following parameters:\n"
3716        first = True
3717        for term in item[:-3]:
3718            line = ''
3719            var = str(term[1])
3720            m = term[0]
3721            if first:
3722                first = False
3723                line += ' = '
3724            else:
3725                if m >= 0:
3726                    line += ' + '
3727                else:
3728                    line += ' - '
3729                m = abs(m)
3730            line += '%.3f*%s '%(m,var)
3731            varMean = G2obj.fmtVarDescr(var)
3732            helptext += "\n" + line + " ("+ varMean + ")"
3733        helptext += '\n\nISODISTORT full name: '+str(fullname)
3734        return helptext
3735
3736    Histograms,Phases = G2frame.GetUsedHistogramsAndPhasesfromTree() # init for constraint
3737    isophases = [p for p in Phases if 'ISODISTORT' in Phases[p]]
3738   
3739    if not isophases:
3740        G2frame.ErrorDialog('no ISODISTORT phases',
3741                            'Unexpected error: no ISODISTORT phases')
3742        return
3743    if phase and phase not in isophases:
3744        G2frame.ErrorDialog('not ISODISTORT phase',
3745                            'Unexpected error: selected phase is not an ISODISTORT phase')
3746        print('Unexpected error: selected phase is not an ISODISTORT phase',
3747                  phase,isophases)
3748    elif not phase and len(isophases) == 1:
3749        phase = isophases[0]
3750    elif not phase:
3751        dlg = wx.SingleChoiceDialog(G2frame,'Select phase from ISODISTORT phases',
3752                                        'Select Phase',isophases)
3753        if dlg.ShowModal() == wx.ID_OK:
3754            sel = dlg.GetSelection()
3755            phase = isophases[sel]
3756        else:
3757            return
3758    # if len(data.get('Histograms',[])) == 0:
3759    #     G2frame.ErrorDialog(
3760    #         'No data',
3761    #         'Sorry, this computation requires that a histogram first be added to the phase'
3762    #         )
3763    #     return
3764   
3765    covdata = G2frame.GPXtree.GetItemPyData(
3766        G2gd.GetGPXtreeItemId(G2frame,G2frame.root,'Covariance'))
3767    # make a lookup table for named NewVar Phase constraints
3768    sub = G2gd.GetGPXtreeItemId(G2frame,G2frame.root,'Constraints') 
3769    Constraints = G2frame.GPXtree.GetItemPyData(sub)
3770    constDict = {}
3771    for c in Constraints['Phase']:
3772        if c[-1] != 'f' or not c[-3]: continue
3773        constDict[str(c[-3])] = c
3774
3775    parmDict,varyList = G2frame.MakeLSParmDict()
3776
3777    dlg = wx.Dialog(G2frame,wx.ID_ANY,'ISODISTORT mode values',#size=(630,400),
3778                       style=wx.DEFAULT_DIALOG_STYLE|wx.RESIZE_BORDER)
3779    mainSizer = wx.BoxSizer(wx.VERTICAL)
3780    if 'ISODISTORT' not in Phases[phase]:
3781        G2frame.ErrorDialog('not ISODISTORT phase',
3782                            'Unexpected error: selected phase is not an ISODISTORT phase')
3783        return
3784    data = Phases[phase]
3785    ISO = data['ISODISTORT']
3786    mainSizer.Add(wx.StaticText(dlg,wx.ID_ANY,
3787                                'ISODISTORT mode computation for cordinates in phase '+
3788                                str(data['General'].get('Name'))))
3789    aSizer = wx.BoxSizer(wx.HORIZONTAL)
3790    panel1 = wxscroll.ScrolledPanel(
3791        dlg, wx.ID_ANY,#size=(100,200),
3792        style = wx.TAB_TRAVERSAL|wx.SUNKEN_BORDER)
3793    subSizer1 = wx.FlexGridSizer(cols=2,hgap=5,vgap=2)
3794    panel2 = wxscroll.ScrolledPanel(
3795        dlg, wx.ID_ANY,#size=(100,200),
3796        style = wx.TAB_TRAVERSAL|wx.SUNKEN_BORDER)
3797    subSizer2 = wx.FlexGridSizer(cols=3,hgap=5,vgap=2)
3798    subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,'Parameter name  '))
3799    subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,' value'),0,wx.ALIGN_RIGHT)
3800    subSizer2.Add((-1,-1))
3801    subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,'Mode name  '))
3802    subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,' value'),0,wx.ALIGN_RIGHT)
3803    # ISODISTORT displacive modes
3804    if 'G2VarList' in ISO:
3805        deltaList = []
3806        for gv,Ilbl in zip(ISO['G2VarList'],ISO['IsoVarList']):
3807            dvar = gv.varname()
3808            var = dvar.replace('::dA','::A')
3809            albl = Ilbl[:Ilbl.rfind('_')]
3810            v = Ilbl[Ilbl.rfind('_')+1:]
3811            pval = ISO['ParentStructure'][albl][['dx','dy','dz'].index(v)]
3812            if var in parmDict:
3813                cval = parmDict[var][0]
3814            else:
3815                dlg.EndModal(wx.ID_CANCEL)
3816                G2frame.ErrorDialog('Atom not found',"No value found for parameter "+str(var))
3817                return
3818            deltaList.append(cval-pval)
3819        modeVals = np.inner(ISO['Var2ModeMatrix'],deltaList)
3820        for lbl,xyz,var,val,norm,G2mode in zip(
3821                ISO['IsoVarList'],deltaList,
3822                ISO['IsoModeList'],modeVals,ISO['NormList'],ISO['G2ModeList'] ):
3823            #GSASIIpath.IPyBreak()
3824            if str(G2mode) in constDict:
3825                ch = G2G.HelpButton(panel2,fmtHelp(constDict[str(G2mode)],var))
3826                subSizer2.Add(ch,0,wx.LEFT|wx.RIGHT|WACV|wx.ALIGN_CENTER,1)
3827            else:
3828                subSizer2.Add((-1,-1))
3829            subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,str(lbl)))
3830            try:
3831                value = G2py3.FormatSigFigs(xyz)
3832            except TypeError:
3833                value = str(xyz)           
3834            subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,value),0,wx.ALIGN_RIGHT)
3835            subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,str(var)))
3836            try:
3837                value = G2py3.FormatSigFigs(val/norm)
3838                if 'varyList' in covdata:
3839                    if str(G2mode) in covdata['varyList']:
3840                        sig = covdata['sig'][covdata['varyList'].index(str(G2mode))]
3841                        value = G2mth.ValEsd(val/norm,sig/norm)
3842            except TypeError:
3843                value = '?'
3844            subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,value),0,wx.ALIGN_RIGHT)
3845            #GSASIIpath.IPyBreak()
3846    # ISODISTORT occupancy modes
3847    if 'G2OccVarList' in ISO:
3848        deltaList = []
3849        for gv,Ilbl in zip(ISO['G2OccVarList'],ISO['OccVarList']):
3850            var = gv.varname()
3851            albl = Ilbl[:Ilbl.rfind('_')]
3852            pval = ISO['BaseOcc'][albl]
3853            if var in parmDict:
3854                cval = parmDict[var][0]
3855            else:
3856                dlg.EndModal(wx.ID_CANCEL)
3857                G2frame.ErrorDialog('Atom not found',"No value found for parameter "+str(var))
3858                return
3859            deltaList.append(cval-pval)
3860        modeVals = np.inner(ISO['Var2OccMatrix'],deltaList)
3861        for lbl,delocc,var,val,norm,G2mode in zip(
3862                ISO['OccVarList'],deltaList,
3863                ISO['OccModeList'],modeVals,ISO['OccNormList'],ISO['G2OccModeList']):
3864            if str(G2mode) in constDict:
3865                ch = G2G.HelpButton(panel2,fmtHelp(constDict[str(G2mode)],var))
3866                subSizer2.Add(ch,0,wx.LEFT|wx.RIGHT|WACV|wx.ALIGN_CENTER,1)
3867            else:
3868                subSizer2.Add((-1,-1))
3869            subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,str(lbl)))
3870            try:
3871                value = G2py3.FormatSigFigs(delocc)
3872            except TypeError:
3873                value = str(delocc)
3874            subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,value),0,wx.ALIGN_RIGHT)
3875            #subSizer.Add((10,-1))
3876            subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,str(var)))
3877            try:
3878                value = G2py3.FormatSigFigs(val/norm)
3879                if 'varyList' in covdata:
3880                    if str(G2mode) in covdata['varyList']:
3881                        sig = covdata['sig'][covdata['varyList'].index(str(G2mode))]
3882                        value = G2mth.ValEsd(val/norm,sig/norm)
3883            except TypeError:
3884                value = '?'
3885            subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,value),0,wx.ALIGN_RIGHT)
3886
3887    # finish up ScrolledPanel
3888    panel1.SetSizer(subSizer1)
3889    panel2.SetSizer(subSizer2)
3890    panel1.SetAutoLayout(1)
3891    panel1.SetupScrolling()
3892    panel2.SetAutoLayout(1)
3893    panel2.SetupScrolling()
3894    # Allow window to be enlarged but not made smaller
3895    dlg.SetSizer(mainSizer)
3896    w1,l1 = subSizer1.GetSize()
3897    w2,l2 = subSizer2.GetSize()
3898    panel1.SetMinSize((w1+10,200))
3899    panel2.SetMinSize((w2+20,200))
3900    aSizer.Add(panel1,1, wx.ALL|wx.EXPAND,1)
3901    aSizer.Add(panel2,2, wx.ALL|wx.EXPAND,1)
3902    mainSizer.Add(aSizer,1, wx.ALL|wx.EXPAND,1)
3903
3904    # make OK button
3905    btnsizer = wx.BoxSizer(wx.HORIZONTAL)
3906    btn = wx.Button(dlg, wx.ID_CLOSE) 
3907    btn.Bind(wx.EVT_BUTTON,_onClose)
3908    btnsizer.Add(btn)
3909    mainSizer.Add(btnsizer, 0, wx.ALIGN_CENTER|wx.ALL, 5)
3910
3911    mainSizer.Fit(dlg)
3912    dlg.SetMinSize(dlg.GetSize())
3913    dlg.ShowModal()
3914    dlg.Destroy()
Note: See TracBrowser for help on using the repository browser.