source: trunk/GSASIIconstrGUI.py @ 4906

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