source: trunk/GSASIIconstrGUI.py @ 4851

Last change on this file since 4851 was 4851, checked in by toby, 7 months ago

reformat equivalence pairs as Bob prefers

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