source: trunk/GSASIIconstrGUI.py @ 4822

Last change on this file since 4822 was 4822, checked in by toby, 9 months ago

Warn if param limits on constrained vars

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