source: trunk/GSASIIconstrGUI.py @ 4907

Last change on this file since 4907 was 4907, checked in by vondreele, 5 months ago

use eigenvalue/vector math to find molecule orientation in onSetPlane

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