source: trunk/GSASIIconstrGUI.py @ 4905

Last change on this file since 4905 was 4905, checked in by toby, 5 months ago

fix delete for general; flag constraints where vars are removed

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