source: trunk/GSASIIconstrGUI.py @ 4983

Last change on this file since 4983 was 4983, checked in by toby, 21 months ago

fix initialization problem when EQUIV are changed to constraints; improve warnings

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