source: trunk/GSASIIconstrGUI.py @ 4809

Last change on this file since 4809 was 4809, checked in by toby, 10 months ago

display more Residue rigid body info; better docs

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