source: trunk/GSASIIconstrGUI.py @ 4849

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

apply symmetry constraints on rigid body origin; show position vars that are fixed on sym-gen tab for constraints; minor bug fixes for missing keys

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 170.9 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASIIconstrGUI - constraint GUI routines
3########### SVN repository information ###################
4# $Date: 2021-03-11 03:07:37 +0000 (Thu, 11 Mar 2021) $
5# $Author: toby $
6# $Revision: 4849 $
7# $URL: trunk/GSASIIconstrGUI.py $
8# $Id: GSASIIconstrGUI.py 4849 2021-03-11 03:07:37Z 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: 4849 $")
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':
1157                    helptext = "The following variable:"
1158                    normval = item[:-3][0][0]
1159                    indepterm = item[:-3][0][1]
1160                    for i,term in enumerate(item[:-3]):
1161                        var = str(term[1])
1162                        if term[0] == 0: term[0] = 1.0
1163                        if len(eqString[-1]) > maxlen:
1164                            eqString.append(' ')
1165                        varMean = G2obj.fmtVarDescr(var)
1166                        if i == 0: # move independent variable to end, as requested by Bob
1167                            helptext += '\n{:} '.format(var) + " ("+ varMean + ")"
1168                            helptext += "\n\nis equivalent to the following, noting multipliers:"
1169                            continue
1170                        elif eqString[-1] != '':
1171                            eqString[-1] += ' = '
1172                        m = normval/term[0]
1173                        if m == 1:
1174                            eqString[-1] += '{:}'.format(var)
1175                        else:
1176                            eqString[-1] += '{:g}*{:} '.format(m,var)
1177                        helptext += '\n{:3g} * {:} '.format(m,var) + " ("+ varMean + ")"
1178                    eqString[-1] += ' = {:} '.format(indepterm)
1179                    typeString = 'EQUIV'
1180                else:
1181                    print ('Unexpected constraint'+item)
1182               
1183            else:
1184                print ('Removing old-style constraints')
1185                data[name] = []
1186                return constSizer
1187            constDel = wx.Button(pageDisplay,wx.ID_ANY,'Delete',style=wx.BU_EXACTFIT)
1188            constDel.Bind(wx.EVT_BUTTON,OnConstDel)
1189            Indx[constDel.GetId()] = [Id,name]
1190            constSizer.Add(constDel,0,wx.LEFT|wx.RIGHT|WACV|wx.ALIGN_CENTER,1)             # delete button
1191            if helptext:
1192                ch = G2G.HelpButton(pageDisplay,helptext)
1193                constSizer.Add(ch,0,wx.LEFT|wx.RIGHT|WACV|wx.ALIGN_CENTER,1)
1194            else:
1195                constSizer.Add((-1,-1))
1196            if refineflag:
1197                ch = G2G.G2CheckBox(pageDisplay,'',item,-2)
1198                constSizer.Add(ch,0,wx.LEFT|wx.RIGHT|WACV|wx.ALIGN_CENTER,1)
1199            else:
1200                constSizer.Add((-1,-1))               
1201            constSizer.Add(wx.StaticText(pageDisplay,wx.ID_ANY,typeString),
1202                           0,WACV|wx.ALIGN_CENTER|wx.RIGHT|wx.LEFT,3)
1203            if problemItem: eqString[-1] += ' -- Conflict: see console'
1204            if len(eqString) > 1:
1205                Eq = wx.BoxSizer(wx.VERTICAL)
1206                for s in eqString:
1207                    line = wx.StaticText(pageDisplay,wx.ID_ANY,s)
1208                    if problemItem: line.SetBackgroundColour(wx.YELLOW)
1209                    Eq.Add(line,0)
1210                Eq.Add((-1,4))
1211            else:
1212                Eq = wx.StaticText(pageDisplay,wx.ID_ANY,eqString[0])
1213                if problemItem: Eq.SetBackgroundColour(wx.YELLOW)
1214            constSizer.Add(Eq,1,WACV)
1215        return constSizer
1216               
1217    def OnConstDel(event):
1218        'Delete a constraint'
1219        Obj = event.GetEventObject()
1220        Id,name = Indx[Obj.GetId()]
1221        del data[name][Id]
1222        ConstraintsLoad(data)
1223        wx.CallAfter(OnPageChanged,None)
1224
1225    def OnConstEdit(event):
1226        '''Called to edit an individual contraint in response to a
1227        click on its Edit button
1228        '''
1229        Obj = event.GetEventObject()
1230        Id,name = Indx[Obj.GetId()]
1231        if data[name][Id][-1] == 'f':
1232            items = data[name][Id][:-3]
1233            constType = 'New Variable'
1234            if data[name][Id][-3]:
1235                varname = str(data[name][Id][-3])
1236            else:
1237                varname = ""
1238            lbl = 'Enter value for each term in constraint; sum = new variable'
1239            dlg = ConstraintDialog(G2frame,constType,lbl,items,
1240                                   varname=varname,varyflag=data[name][Id][-2])
1241        elif data[name][Id][-1] == 'c':
1242            items = data[name][Id][:-3]+[
1243                [str(data[name][Id][-3]),'fixed value =']]
1244            constType = 'Constraint'
1245            lbl = 'Edit value for each term in constant constraint sum'
1246            dlg = ConstraintDialog(G2frame,constType,lbl,items)
1247        elif data[name][Id][-1] == 'e':
1248            items = data[name][Id][:-3]
1249            constType = 'Equivalence'
1250            lbl = 'The following terms are set to be equal:'
1251            dlg = ConstraintDialog(G2frame,constType,lbl,items,'/')
1252        else:
1253            return
1254        try:
1255            prev = copy.deepcopy(data[name][Id])
1256            if dlg.ShowModal() == wx.ID_OK:
1257                result = dlg.GetData()
1258                for i in range(len(data[name][Id][:-3])):
1259                    if type(data[name][Id][i]) is tuple: # fix non-mutable construct
1260                        data[name][Id][i] = list(data[name][Id][i])
1261                    data[name][Id][i][0] = result[i][0]
1262                if data[name][Id][-1] == 'c':
1263                    data[name][Id][-3] = str(result[-1][0])
1264                elif data[name][Id][-1] == 'f':
1265                    data[name][Id][-2] = dlg.newvar[1]
1266                    if dlg.newvar[0]:
1267                        # process the variable name to put in global form (::var)
1268                        varname = str(dlg.newvar[0]).strip().replace(' ','_')
1269                        if varname.startswith('::'):
1270                            varname = varname[2:]
1271                        varname = varname.replace(':',';')
1272                        if varname:
1273                            data[name][Id][-3] = varname
1274                        else:
1275                            data[name][Id][-3] = ''
1276                if not CheckChangedConstraint():
1277                    data[name][Id] = prev
1278            else:
1279                data[name][Id] = prev
1280        except:
1281            import traceback
1282            print (traceback.format_exc())
1283        finally:
1284            dlg.Destroy()
1285        wx.CallAfter(OnPageChanged,None)
1286   
1287    def UpdateConstraintPanel(panel,typ):
1288        '''Update the contents of the selected Constraint
1289        notebook tab. Called in :func:`OnPageChanged`
1290        '''
1291        if panel.GetSizer(): panel.GetSizer().Clear(True)
1292        Siz = wx.BoxSizer(wx.VERTICAL)
1293        Siz.Add((5,5),0)
1294        Siz.Add(MakeConstraintsSizer(typ,panel),1,wx.EXPAND)
1295        panel.SetSizer(Siz,True)
1296        Size = Siz.GetMinSize()
1297        Size[0] += 40
1298        Size[1] = max(Size[1],450) + 20
1299        panel.SetSize(Size)
1300        panel.SetScrollbars(10,10,Size[0]/10-4,Size[1]/10-1)
1301        panel.Show()
1302
1303    def OnPageChanged(event):
1304        '''Called when a tab is pressed or when a "select tab" menu button is
1305        used (see RaisePage), or to refresh the current tab contents (event=None)
1306        '''
1307        if event:       #page change event!
1308            page = event.GetSelection()
1309        else: # called directly, get current page
1310            try:
1311                page = G2frame.constr.GetSelection()
1312            except:
1313                if GSASIIpath.GetConfigValue('debug'): print('DBG_gpx open error:C++ Run time error - skipped')
1314                return
1315        G2frame.constr.ChangeSelection(page)
1316        text = G2frame.constr.GetPageText(page)
1317        G2frame.dataWindow.ConstraintEdit.Enable(G2G.wxID_EQUIVALANCEATOMS,False)
1318#        G2frame.dataWindow.ConstraintEdit.Enable(G2G.wxID_ADDRIDING,False)
1319        if text == 'Histogram/Phase':
1320            enableEditCons = [False]+4*[True]
1321            G2frame.Page = [page,'hap']
1322            UpdateConstraintPanel(HAPConstr,'HAP')
1323        elif text == 'Histogram':
1324            enableEditCons = [False]+4*[True]
1325            G2frame.Page = [page,'hst']
1326            UpdateConstraintPanel(HistConstr,'Hist')
1327        elif text == 'Phase':
1328            enableEditCons = 5*[True]
1329            G2frame.Page = [page,'phs']
1330            G2frame.dataWindow.ConstraintEdit.Enable(G2G.wxID_EQUIVALANCEATOMS,True)
1331#            G2frame.dataWindow.ConstraintEdit.Enable(G2G.wxID_ADDRIDING,True)
1332            if 'DELETED' in str(PhaseConstr):   #seems to be no other way to do this (wx bug)
1333                if GSASIIpath.GetConfigValue('debug'):
1334                    print ('DBG_wx error: PhaseConstr not cleanly deleted after Refine')
1335                return
1336            UpdateConstraintPanel(PhaseConstr,'Phase')
1337        elif text == 'Global':
1338            enableEditCons = [False]+4*[True]
1339            G2frame.Page = [page,'glb']
1340            UpdateConstraintPanel(GlobalConstr,'Global')
1341        else:
1342            enableEditCons = 5*[False]
1343            G2frame.Page = [page,'sym']
1344            UpdateConstraintPanel(SymConstr,'Sym-Generated')
1345        # remove menu items when not allowed
1346        for obj,flag in zip(G2frame.dataWindow.ConstraintEdit.GetMenuItems(),enableEditCons): 
1347            obj.Enable(flag)
1348        G2frame.dataWindow.SetDataSize()
1349
1350    def RaisePage(event):
1351        'Respond to a "select tab" menu button'
1352        try:
1353            i = (G2G.wxID_CONSPHASE,
1354                 G2G.wxID_CONSHAP,
1355                 G2G.wxID_CONSHIST,
1356                 G2G.wxID_CONSGLOBAL,
1357                 G2G.wxID_CONSSYM,
1358                ).index(event.GetId())
1359            G2frame.constr.SetSelection(i)
1360            wx.CallAfter(OnPageChanged,None)
1361        except ValueError:
1362            print('Unexpected event in RaisePage')
1363
1364    def SetStatusLine(text):
1365        G2frame.GetStatusBar().SetStatusText(text,1)
1366       
1367    def OnShowISODISTORT(event):
1368        ShowIsoDistortCalc(G2frame)
1369
1370    #### UpdateConstraints execution starts here ##############################
1371    if not data:
1372        data.update({'Hist':[],'HAP':[],'Phase':[],'Global':[]})       #empty dict - fill it
1373    if 'Global' not in data:                                            #patch
1374        data['Global'] = []
1375    # DEBUG code #=====================================
1376    #import GSASIIconstrGUI
1377    #reload(GSASIIconstrGUI)
1378    #reload(G2obj)
1379    #reload(G2stIO)
1380    #import GSASIIstrMain
1381    #reload(GSASIIstrMain)   
1382    #reload(G2mv)
1383    #reload(G2gd)
1384    #===================================================
1385    Histograms,Phases = G2frame.GetUsedHistogramsAndPhasesfromTree()
1386    if not len(Phases) or not len(Histograms):
1387        dlg = wx.MessageDialog(G2frame,'You need both phases and histograms to see Constraints',
1388            'No phases or histograms')
1389        dlg.CenterOnParent()
1390        dlg.ShowModal()
1391        dlg.Destroy()
1392        return
1393    G2obj.IndexAllIds(Histograms,Phases)
1394    for p in Phases:
1395        if 'ISODISTORT' in Phases[p]:
1396            G2frame.dataWindow.ConstraintEdit.Enable(G2G.wxID_SHOWISO,True)
1397            break
1398    ###### patch: convert old-style (str) variables in constraints to G2VarObj objects #####
1399    for key,value in data.items():
1400        if key.startswith('_'): continue
1401        j = 0
1402        for cons in value:
1403            #print cons             # DEBUG
1404            for i in range(len(cons[:-3])):
1405                if type(cons[i][1]) is str:
1406                    cons[i][1] = G2obj.G2VarObj(cons[i][1])
1407                    j += 1
1408        if j:
1409            print (str(key) + ': '+str(j)+' variable(s) as strings converted to objects')
1410    ##### end patch #############################################################################
1411    rigidbodyDict = G2frame.GPXtree.GetItemPyData(
1412        G2gd.GetGPXtreeItemId(G2frame,G2frame.root,'Rigid bodies'))
1413    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
1414    rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,Print=False)
1415    badPhaseParms = ['Ax','Ay','Az','Amul','AI/A','Atype','SHorder','AwaveType','FwaveType','PwaveType','MwaveType','Vol','isMag',]
1416    globalList = list(rbDict.keys())
1417    globalList.sort()
1418    try:
1419        AtomDict = dict([Phases[phase]['pId'],Phases[phase]['Atoms']] for phase in Phases)
1420    except KeyError:
1421        G2frame.ErrorDialog('Constraint Error','Constraints cannot be set until a cycle of least squares'+
1422                            ' has been run.\nWe suggest you refine a scale factor.')
1423        return
1424
1425    # create a list of the phase variables
1426    symHolds = []
1427    (Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtable,BLtable,
1428         MFtable,maxSSwave) = G2stIO.GetPhaseData(
1429             Phases,rbIds=rbIds,Print=False,symHold=symHolds)
1430    phaseList = []
1431    for item in phaseDict:
1432        if item.split(':')[2] not in badPhaseParms:
1433            phaseList.append(item)
1434    phaseList.sort()
1435    phaseAtNames = {}
1436    phaseAtTypes = {}
1437    TypeList = []
1438    for item in phaseList:
1439        Split = item.split(':')
1440        if Split[2][:2] in ['AU','Af','dA','AM']:
1441            Id = int(Split[0])
1442            phaseAtNames[item] = AtomDict[Id][int(Split[3])][0]
1443            phaseAtTypes[item] = AtomDict[Id][int(Split[3])][1]
1444            if phaseAtTypes[item] not in TypeList:
1445                TypeList.append(phaseAtTypes[item])
1446        else:
1447            phaseAtNames[item] = ''
1448            phaseAtTypes[item] = ''
1449             
1450    # create a list of the hist*phase variables
1451    seqList = G2frame.testSeqRefineMode()
1452    if seqList: # for sequential refinement, only process 1st histgram in list
1453        histDict = {seqList[0]:Histograms[seqList[0]]}
1454    else:
1455        histDict = Histograms
1456    hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,histDict,Print=False,resetRefList=False)
1457    hapList = sorted([i for i in hapDict.keys() if i.split(':')[2] not in ('Type',)])
1458    if seqList: # convert histogram # to wildcard
1459        wildList = [] # list of variables with "*" for histogram number
1460        for i in hapList:
1461            s = i.split(':')
1462            if s[1] == "": continue
1463            s[1] = '*'
1464            sj = ':'.join(s)
1465            if sj not in wildList: wildList.append(sj)
1466        hapList = wildList
1467    histVary,histDict,controlDict = G2stIO.GetHistogramData(histDict,Print=False)
1468    histList = list(histDict.keys())
1469    histList.sort()
1470    if seqList: # convert histogram # to wildcard
1471        wildList = [] # list of variables with "*" for histogram number
1472        for i in histList:
1473            s = i.split(':')
1474            if s[1] == "": continue
1475            s[1] = '*'
1476            sj = ':'.join(s)
1477            if sj not in wildList: wildList.append(sj)
1478        histList = wildList       
1479    Indx = {}
1480    G2frame.Page = [0,'phs']
1481   
1482    G2gd.SetDataMenuBar(G2frame,G2frame.dataWindow.ConstraintMenu)
1483    SetStatusLine('')
1484   
1485    G2frame.Bind(wx.EVT_MENU, OnAddConstraint, id=G2G.wxID_CONSTRAINTADD)
1486    G2frame.Bind(wx.EVT_MENU, OnAddFunction, id=G2G.wxID_FUNCTADD)
1487    G2frame.Bind(wx.EVT_MENU, OnAddEquivalence, id=G2G.wxID_EQUIVADD)
1488    G2frame.Bind(wx.EVT_MENU, OnAddHold, id=G2G.wxID_HOLDADD)
1489    G2frame.Bind(wx.EVT_MENU, OnAddAtomEquiv, id=G2G.wxID_EQUIVALANCEATOMS)
1490#    G2frame.Bind(wx.EVT_MENU, OnAddRiding, id=G2G.wxID_ADDRIDING)
1491    G2frame.Bind(wx.EVT_MENU, OnShowISODISTORT, id=G2G.wxID_SHOWISO)
1492    # tab commands
1493    for id in (G2G.wxID_CONSPHASE,
1494               G2G.wxID_CONSHAP,
1495               G2G.wxID_CONSHIST,
1496               G2G.wxID_CONSGLOBAL,
1497               G2G.wxID_CONSSYM,
1498               ):
1499        G2frame.Bind(wx.EVT_MENU, RaisePage,id=id)
1500
1501    #G2frame.constr = G2G.GSNoteBook(parent=G2frame.dataWindow,size=G2frame.dataWindow.GetClientSize())
1502    G2frame.constr = G2G.GSNoteBook(parent=G2frame.dataWindow)
1503    G2frame.dataWindow.GetSizer().Add(G2frame.constr,1,wx.ALL|wx.EXPAND)
1504    # note that order of pages is hard-coded in RaisePage
1505    PhaseConstr = wx.ScrolledWindow(G2frame.constr)
1506    G2frame.constr.AddPage(PhaseConstr,'Phase')
1507    HAPConstr = wx.ScrolledWindow(G2frame.constr)
1508    G2frame.constr.AddPage(HAPConstr,'Histogram/Phase')
1509    HistConstr = wx.ScrolledWindow(G2frame.constr)
1510    G2frame.constr.AddPage(HistConstr,'Histogram')
1511    GlobalConstr = wx.ScrolledWindow(G2frame.constr)
1512    G2frame.constr.AddPage(GlobalConstr,'Global')
1513    SymConstr = wx.ScrolledWindow(G2frame.constr)
1514    G2frame.constr.AddPage(SymConstr,'Sym-Generated')
1515    wx.CallAfter(OnPageChanged,None)
1516    G2frame.constr.Bind(wx.aui.EVT_AUINOTEBOOK_PAGE_CHANGED, OnPageChanged)
1517    # validate all the constrants -- should not see any errors here normally
1518    errmsg,warnmsg = ConstraintsCheck(data)
1519    if errmsg:
1520        G2frame.ErrorDialog('Constraint Error',
1521                            'Error in constraints:\n'+errmsg+'\nCheck console output for more information',
1522                            parent=G2frame)
1523        print (errmsg)
1524        print (G2mv.VarRemapShow([],True))
1525        return
1526    elif warnmsg:
1527        print ('Unexpected contraint warning:\n'+warnmsg)
1528    WarnConstraintLimit()
1529       
1530###### check scale & phase fractions, create constraint if needed #############
1531def CheckAllScalePhaseFractions(G2frame):
1532    '''Check if scale factor and all phase fractions are refined without a constraint
1533    for all used histograms, if so, offer the user a chance to create a constraint
1534    on the sum of phase fractions
1535    '''
1536    histograms, phases = G2frame.GetUsedHistogramsAndPhasesfromTree()
1537    for i,hist in enumerate(histograms):
1538        CheckScalePhaseFractions(G2frame,hist,histograms,phases)
1539       
1540def CheckScalePhaseFractions(G2frame,hist,histograms,phases):
1541    '''Check if scale factor and all phase fractions are refined without a constraint
1542    for histogram hist, if so, offer the user a chance to create a constraint
1543    on the sum of phase fractions
1544    '''
1545    if G2frame.testSeqRefineMode():
1546        histStr = '*'
1547    else:
1548        histStr = str(histograms[hist]['hId'])
1549    # Is this powder?
1550    if not hist.startswith('PWDR '): return
1551    # do this only if the scale factor is varied
1552    if not histograms[hist]['Sample Parameters']['Scale'][1]: return
1553    # are all phase fractions varied in all used histograms?
1554    phaseCount = 0
1555    for p in phases:
1556        if hist not in phases[p]['Histograms']: continue
1557        if phases[p]['Histograms'][hist]['Use'] and not phases[p]['Histograms'][hist]['Scale'][1]:
1558            return
1559        else:
1560            phaseCount += 1
1561   
1562    # all phase fractions and scale factor varied, now scan through constraints
1563    sub = G2gd.GetGPXtreeItemId(G2frame,G2frame.root,'Constraints') 
1564    Constraints = G2frame.GPXtree.GetItemPyData(sub)
1565    for c in Constraints.get('HAP',[]):
1566        if c[-1] != 'c': continue
1567        if not c[-3]: continue
1568        if len(c[:-3]) != phaseCount: continue
1569        # got a constraint equation with right number of terms, is it on phase fractions for
1570        # the correct histogram?
1571        if all([(i[1].name == 'Scale' and i[1].varname().split(':')[1] == histStr) for i in c[:-3]]):
1572            # got a constraint, this is OK
1573            return
1574    dlg = wx.MessageDialog(G2frame,'You are refining the scale factor and all phase fractions for histogram #'+
1575        histStr+'. This will produce an unstable refinement. '+
1576        'Do you want to constrain the sum of phase fractions?','Create constraint?',wx.OK|wx.CANCEL)
1577    if dlg.ShowModal() != wx.ID_OK:
1578        dlg.Destroy()
1579        return
1580    dlg.Destroy()
1581
1582    constr = []
1583    for p in phases:
1584        if hist not in phases[p]['Histograms']: continue
1585        if not phases[p]['Histograms'][hist]['Use']: continue
1586        constr += [[1.0,G2obj.G2VarObj(':'.join((str(phases[p]['pId']),histStr,'Scale')))]]
1587    constr += [1.0,None,'c']
1588    Constraints['HAP'] += [constr]
1589       
1590#### Make nuclear/magnetic phase transition constraints - called by OnTransform in G2phsGUI ##########
1591def TransConstraints(G2frame,oldPhase,newPhase,Trans,Vec,atCodes):
1592    '''Add constraints for new magnetic phase created via transformation of old
1593    nuclear one
1594    NB: A = [G11,G22,G33,2*G12,2*G13,2*G23]
1595    '''
1596   
1597    def SetUniqAj(pId,iA,SGData):
1598        SGLaue = SGData['SGLaue']
1599        SGUniq = SGData['SGUniq']
1600        if SGLaue in ['m3','m3m']:
1601            if iA in [0,1,2]:
1602                parm = '%d::%s'%(pId,'A0')
1603            else:
1604                parm = None
1605        elif SGLaue in ['4/m','4/mmm']:
1606            if iA in [0,1]:
1607                parm = '%d::%s'%(pId,'A0')
1608            elif iA == 2:
1609                parm = '%d::%s'%(pId,'A2')
1610            else:
1611                parm = None
1612        elif SGLaue in ['6/m','6/mmm','3m1', '31m', '3']:
1613            if iA in [0,1,3]:
1614                parm = '%d::%s'%(pId,'A0')
1615            elif iA == 2:
1616                parm = '%d::%s'%(pId,'A2')
1617            else:
1618                parm = None
1619        elif SGLaue in ['3R', '3mR']:
1620            if ia in [0,1,2]:
1621                parm = '%d::%s'%(pId,'A0')
1622            else:
1623                parm = '%d::%s'%(pId,'A3')
1624        elif SGLaue in ['mmm',]:
1625            if iA in [0,1,2]:
1626                parm = '%d::A%s'%(pId,iA)
1627            else:
1628                parm = None
1629        elif SGLaue == '2/m':
1630            if iA in [0,1,2]:
1631                parm = '%d::A%s'%(pId,iA)
1632            elif iA == 3 and SGUniq == 'c':
1633                parm = '%d::A%s'%(pId,iA)
1634            elif iA == 4 and SGUniq == 'b':
1635                parm = '%d::A%s'%(pId,iA)
1636            elif iA == 5 and SGUniq == 'a':
1637                parm = '%d::A%s'%(pId,iA)
1638            else:
1639                parm = None           
1640        else:
1641            parm = '%d::A%s'%(pId,iA)
1642        return parm
1643   
1644    Histograms,Phases = G2frame.GetUsedHistogramsAndPhasesfromTree()
1645    UseList = newPhase['Histograms']
1646    detTrans = np.abs(nl.det(Trans))
1647    opId = oldPhase['pId']
1648    npId = newPhase['pId']
1649    cx,ct,cs,cia = newPhase['General']['AtomPtrs']
1650    nAtoms = newPhase['Atoms']
1651    nSGData = newPhase['General']['SGData']
1652    #oAcof = G2lat.cell2A(oldPhase['General']['Cell'][1:7])
1653    #nAcof = G2lat.cell2A(newPhase['General']['Cell'][1:7])
1654    item = G2gd.GetGPXtreeItemId(G2frame,G2frame.root,'Constraints')
1655    if not item:
1656        print('Error: no constraints in Data Tree')
1657        return
1658    constraints = G2frame.GPXtree.GetItemPyData(item)
1659    xnames = ['dAx','dAy','dAz']
1660    # constraints on matching atom params between phases
1661    for ia,code in enumerate(atCodes):
1662        atom = nAtoms[ia]
1663        if not ia and atom[cia] == 'A':
1664            wx.MessageDialog(G2frame,
1665                'Anisotropic thermal motion constraints are not developed at the present time',
1666                'Anisotropic thermal constraint?',style=wx.ICON_INFORMATION).ShowModal()
1667        siteSym = G2spc.SytSym(atom[cx:cx+3],nSGData)[0]
1668        CSX = G2spc.GetCSxinel(siteSym)
1669#        CSU = G2spc.GetCSuinel(siteSym)
1670        item = code.split('+')[0]
1671        iat,opr = item.split(':')
1672        Nop = abs(int(opr))%100-1
1673        if '-' in opr:
1674            Nop *= -1
1675        Opr = oldPhase['General']['SGData']['SGOps'][abs(Nop)][0]
1676        if Nop < 0:         #inversion
1677            Opr *= -1
1678        XOpr = np.inner(Opr,Trans)
1679        for i,ix in enumerate(list(CSX[0])):
1680            if not ix:
1681                continue
1682            name = xnames[i]
1683            IndpCon = [1.0,G2obj.G2VarObj('%d::%s:%d'%(npId,name,ia))]
1684            DepCons = []
1685            for iop,opval in enumerate(XOpr[i]):
1686                if opval:
1687                    DepCons.append([opval,G2obj.G2VarObj('%d::%s:%s'%(opId,xnames[iop],iat))])
1688            if len(DepCons) == 1:
1689                constraints['Phase'].append([DepCons[0],IndpCon,None,None,'e'])
1690            elif len(DepCons) > 1:
1691                for Dep in DepCons:
1692                    Dep[0] *= -1
1693                constraints['Phase'].append([IndpCon]+DepCons+[0.0,None,'c'])
1694        for name in ['Afrac','AUiso']:
1695            IndpCon = [1.0,G2obj.G2VarObj('%d::%s:%d'%(npId,name,ia))]
1696            DepCons = [1.0,G2obj.G2VarObj('%d::%s:%s'%(opId,name,iat))]
1697            constraints['Phase'].append([DepCons,IndpCon,None,None,'e'])
1698           
1699        # unfinished Anisotropic constraint generation
1700#        Uids = [[0,0,'AU11'],[1,1,'AU22'],[2,2,'AU33'],[0,1,'AU12'],[0,2,'AU13'],[1,2,'AU23']]
1701#        DepConsDict = dict(zip(Us,[[],[],[],[],[],[]]))
1702#        for iu,Uid in enumerate(Uids):
1703#            UMT = np.zeros((3,3))
1704#            UMT[Uid[0],Uid[1]] = 1
1705#            nUMT = G2lat.prodMGMT(UMT,invTrans)
1706#            nUT = G2lat.UijtoU6(nUMT)
1707#            for iu,nU in enumerate(nUT):
1708#                if abs(nU) > 1.e-8:
1709#                    parm = '%d::%s;%s'%(opId,Us[iu],iat)
1710#                    DepConsDict[Uid[2]].append([abs(nU%1.),G2obj.G2VarObj(parm)])
1711#        nUcof = atom[iu:iu+6]
1712#        conStrings = []
1713#        for iU,Usi in enumerate(Us):
1714#            parm = '%d::%s;%d'%(npId,Usi,ia)
1715#            parmDict[parm] = nUcof[iU]
1716#            varyList.append(parm)
1717#            IndpCon = [1.0,G2obj.G2VarObj(parm)]
1718#            conStr = str([IndpCon,DepConsDict[Usi]])
1719#            if conStr in conStrings:
1720#                continue
1721#            conStrings.append(conStr)
1722#            if len(DepConsDict[Usi]) == 1:
1723#                if DepConsDict[Usi][0]:
1724#                    constraints['Phase'].append([IndpCon,DepConsDict[Usi][0],None,None,'e'])
1725#            elif len(DepConsDict[Usi]) > 1:       
1726#                for Dep in DepConsDict[Usi]:
1727#                    Dep[0] *= -1
1728#                constraints['Phase'].append([IndpCon]+DepConsDict[Usi]+[0.0,None,'c'])
1729           
1730        #how do I do Uij's for most Trans?
1731
1732    # constraints on lattice parameters between phases
1733#    T = nl.inv(Trans).T
1734    # T = Trans.T
1735    # conMat = [
1736    #     [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]],
1737    #     [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]],
1738    #     [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]],
1739    #     [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]],
1740    #     [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]],
1741    #     [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]]
1742    #     ]
1743    # Gnew = conMat * A:
1744#         T00**2*a0  T01**2*a1 T02**2*a2 T00*T01*a3    T00*T02*a4    T01*T02*a5
1745#         T10**2*a0  T11**2*a1 T12**2*a2 T10*T11*a3    T10*T12*a4    T11*T12*a5
1746#         T20**2*a0  T21**2*a1 T22**2*a2 T20*T21*a3    T20*T22*a4    T21*T22*a5
1747#         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
1748#         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
1749#         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
1750    # Generated as symbolic code using:
1751    # import sym
1752    # A0, A1, A2, A3, A4, A5 = sym.symbols('A0, A1, A2, A3, A4, A5')
1753    # G = sym.Matrix([ [A0,    A3/2.,  A4/2.], [A3/2.,  A1,    A5/2.], [A4/2., A5/2.,    A2]])
1754    # transformation matrix
1755    # T00, T10, T20, T01, T11, T21, T02, T12, T22 = sym.symbols('T00, T10, T20, T01, T11, T21, T02, T12, T22')
1756    # Tr = sym.Matrix([ [T00, T10, T20], [T01, T11, T21], [T02, T12, T22],])
1757    # Gnew = (Tr.T*G)*Tr
1758   
1759    #print('old A',G2lat.cell2A(oldPhase['General']['Cell'][1:7]))
1760    #print('new A',G2lat.cell2A(newPhase['General']['Cell'][1:7]))
1761   
1762#this is still incorrect for hex/trig/ortho/tetragonal --> monoclinic
1763   
1764#    for iAnew,Asi in enumerate(['A0','A1','A2','A3','A4','A5']): # loop through A[i] for new cell
1765#        Nparm = str(npId) + '::' + Asi
1766#        if Nparm != SetUniqAj(npId,iAnew,nSGData):
1767#            continue # skip: Ai constrained from Aj or must be zero
1768#        multDict = {}
1769#        for iAorg in range(6):
1770#            cA = conMat[iAnew][iAorg] # coeff for A[i] in constraint matrix
1771#            if abs(cA) < 1.e-8: continue
1772#            parm = SetUniqAj(opId,iAorg,oSGData) # translate to unique A[i] in original cell
1773#            if not parm: continue # must be zero
1774#            # sum coeff
1775#            if parm in multDict:
1776#                multDict[parm] += cA
1777#            else:
1778#                multDict[parm] = cA
1779#        # any non-zero multipliers?
1780#        maxMult = 0
1781#        for i in multDict:
1782#            maxMult = max(maxMult,abs(multDict[i]))
1783#        if maxMult <= 0:  # Nparm computes as zero; Fix this parameter
1784#            constraints['Phase'] += [[
1785#                [0.0,G2obj.G2VarObj(Nparm)],
1786#                None,None,'h']]
1787#        elif len(multDict) == 1:        # create equivalence
1788#            key = list(multDict.keys())[0]
1789#            constraints['Phase'] += [[
1790#                [1.0,G2obj.G2VarObj(key)],
1791#                [multDict[key],G2obj.G2VarObj(Nparm)],
1792#                None,None,'e']]
1793#        else:                           # create constraint
1794#            constr = [[-1.0,G2obj.G2VarObj(Nparm)]]
1795#            for key in multDict:
1796#                constr += [[multDict[key],G2obj.G2VarObj(key)]]
1797#            constr += [0.0,None,'c']
1798#            constraints['Phase'] += [constr]
1799   
1800    # constraints on HAP Scale, etc.
1801    for hId,hist in enumerate(UseList):    #HAP - seems OK
1802        ohapkey = '%d:%d:'%(opId,hId)
1803        nhapkey = '%d:%d:'%(npId,hId)
1804        IndpCon = [1.0,G2obj.G2VarObj(ohapkey+'Scale')]
1805        DepCons = [detTrans,G2obj.G2VarObj(nhapkey+'Scale')]
1806        constraints['HAP'].append([DepCons,IndpCon,None,None,'e'])
1807        for name in ['Size;i','Mustrain;i']:
1808            IndpCon = [1.0,G2obj.G2VarObj(ohapkey+name)]
1809            DepCons = [1.0,G2obj.G2VarObj(nhapkey+name)]
1810            constraints['HAP'].append([IndpCon,DepCons,None,None,'e'])
1811       
1812#### Rigid bodies #############################################################
1813resRBsel = None
1814def UpdateRigidBodies(G2frame,data):
1815    '''Called when Rigid bodies tree item is selected.
1816    Displays the rigid bodies in the data window
1817    ''' 
1818    def OnPageChanged(event):
1819        global resList
1820        resList = []
1821        if event:       #page change event!
1822            page = event.GetSelection()
1823        else:
1824            try:
1825                page = G2frame.rbBook.GetSelection()
1826            except:
1827                if GSASIIpath.GetConfigValue('debug'): print('DBG_gpx open error:C++ Run time error - skipped')
1828                return
1829        G2frame.rbBook.ChangeSelection(page)
1830        text = G2frame.rbBook.GetPageText(page)
1831        if text == 'Vector rigid bodies':
1832            G2gd.SetDataMenuBar(G2frame,G2frame.dataWindow.VectorBodyMenu)
1833            G2frame.Bind(wx.EVT_MENU, AddVectorRB, id=G2G.wxID_VECTORBODYADD)
1834            G2frame.Bind(wx.EVT_MENU, ExtractPhaseRB, id=G2G.wxID_VECTORBODYIMP)
1835            G2frame.Bind(wx.EVT_MENU, AddVectTrans, id=G2G.wxID_VECTORBODYEXTD)
1836            G2frame.Bind(wx.EVT_MENU, SaveVectorRB, id=G2G.wxID_VECTORBODYSAV)
1837            G2frame.Bind(wx.EVT_MENU, ReadVectorRB, id=G2G.wxID_VECTORBODYRD)
1838            G2frame.Page = [page,'vrb']
1839            UpdateVectorRB()
1840        elif text == 'Residue rigid bodies':
1841            G2gd.SetDataMenuBar(G2frame,G2frame.dataWindow.RigidBodyMenu)
1842            G2frame.Bind(wx.EVT_MENU, AddResidueRB, id=G2G.wxID_RIGIDBODYADD)
1843            G2frame.Bind(wx.EVT_MENU, ExtractPhaseRB, id=G2G.wxID_RIGIDBODYIMP)
1844            G2frame.Bind(wx.EVT_MENU, OnImportRigidBody, id=G2G.wxID_RIGIDBODYIMPORT)
1845            G2frame.Bind(wx.EVT_MENU, OnSaveRigidBody, id=G2G.wxID_RIGIDBODYSAVE)
1846            G2frame.Bind(wx.EVT_MENU, OnDefineTorsSeq, id=G2G.wxID_RESIDUETORSSEQ) #enable only if residue RBs exist?
1847            G2frame.Bind(wx.EVT_MENU, DumpVectorRB, id=G2G.wxID_RESBODYSAV)
1848            G2frame.Bind(wx.EVT_MENU, LoadVectorRB, id=G2G.wxID_RESBODYRD)
1849            G2frame.Page = [page,'rrb']
1850            UpdateResidueRB()
1851        else:
1852            G2gd.SetDataMenuBar(G2frame)
1853            #G2frame.Page = [page,'rrb']
1854           
1855    def getMacroFile(macName):
1856        defDir = os.path.join(os.path.split(__file__)[0],'GSASIImacros')
1857        dlg = wx.FileDialog(G2frame,message='Choose '+macName+' rigid body macro file',
1858            defaultDir=defDir,defaultFile="",wildcard="GSAS-II macro file (*.mac)|*.mac",
1859            style=wx.FD_OPEN | wx.FD_CHANGE_DIR)
1860        try:
1861            if dlg.ShowModal() == wx.ID_OK:
1862                macfile = dlg.GetPath()
1863                macro = open(macfile,'r')
1864                head = macro.readline()
1865                if macName not in head:
1866                    print (head)
1867                    print ('**** ERROR - wrong restraint macro file selected, try again ****')
1868                    macro = []
1869            else: # cancel was pressed
1870                macro = []
1871        finally:
1872            dlg.Destroy()
1873        return macro        #advanced past 1st line
1874       
1875    def getTextFile():
1876        dlg = wx.FileDialog(G2frame,'Choose rigid body text file', G2frame.LastGPXdir, '',
1877            "GSAS-II text file (*.txt)|*.txt|XYZ file (*.xyz)|*.xyz|"
1878            "Sybyl mol2 file (*.mol2)|*.mol2|PDB file (*.pdb;*.ent)|*.pdb;*.ent",
1879            wx.FD_OPEN | wx.FD_CHANGE_DIR)
1880        try:
1881            if dlg.ShowModal() == wx.ID_OK:
1882                txtfile = dlg.GetPath()
1883                ext = os.path.splitext(txtfile)[1]
1884                text = open(txtfile,'r')
1885            else: # cancel was pressed
1886                ext = ''
1887                text = []
1888        finally:
1889            dlg.Destroy()
1890        if 'ent' in ext:
1891            ext = '.pdb'
1892        return text,ext.lower()
1893       
1894    def OnImportRigidBody(event):
1895        page = G2frame.rbBook.GetSelection()
1896        if 'Vector' in G2frame.rbBook.GetPageText(page):
1897            pass
1898        elif 'Residue' in G2frame.rbBook.GetPageText(page):
1899            ImportResidueRB()
1900           
1901    def OnSaveRigidBody(event):
1902        page = G2frame.rbBook.GetSelection()
1903        if 'Vector' in G2frame.rbBook.GetPageText(page):
1904            pass
1905        elif 'Residue' in G2frame.rbBook.GetPageText(page):
1906            SaveResidueRB()
1907           
1908    def DumpVectorRB(event):
1909        global resRBsel
1910        if resRBsel not in data['Residue']:
1911            return
1912        rbData = data['Residue'][resRBsel]
1913        pth = G2G.GetExportPath(G2frame)
1914        dlg = wx.FileDialog(G2frame, 'Choose file to save residue rigid body',
1915            pth, '', 'RRB files (*.resbody)|*.resbody',
1916            wx.FD_SAVE|wx.FD_OVERWRITE_PROMPT)
1917        try:
1918            if dlg.ShowModal() == wx.ID_OK:
1919                filename = dlg.GetPath()
1920                filename = os.path.splitext(filename)[0]+'.resbody'  # set extension
1921                fp = open(filename,'w')
1922                fp.write('Name: '+rbData['RBname']+'\n')
1923                fp.write('atNames: ')
1924                for i in rbData['atNames']:
1925                    fp.write(str(i)+" ") 
1926                fp.write('\n')
1927                for item in rbData['rbSeq']:
1928                    fp.write('rbSeq: ') 
1929                    fp.write('{:d} {:d} {:.1f}: '.format(*item[:3]))
1930                    for num in item[3]:
1931                        fp.write('{:d} '.format(num))
1932                    fp.write('\n')
1933                for i,sym in enumerate(rbData['rbTypes']):
1934                    fp.write("{:3s}".format(sym))
1935                    fp.write('{:8.5f}{:9.5f}{:9.5f}   '
1936                            .format(*rbData['rbXYZ'][i]))
1937                    fp.write('\n')
1938                fp.close()
1939                print ('Vector rigid body saved to: '+filename)
1940               
1941        finally:
1942            dlg.Destroy()
1943       
1944    def LoadVectorRB(event):
1945        AtInfo = data['Residue']['AtInfo']
1946        pth = G2G.GetExportPath(G2frame)
1947        dlg = wx.FileDialog(G2frame, 'Choose file to read vector rigid body',
1948            pth, '', 'RRB files (*.resbody)|*.resbody',
1949            wx.FD_OPEN)
1950        try:
1951            if dlg.ShowModal() == wx.ID_OK:
1952                filename = dlg.GetPath()
1953                filename = os.path.splitext(filename)[0]+'.resbody'  # set extension
1954                fp = open(filename,'r')
1955                l = fp.readline().strip()
1956                if 'Name' not in l:
1957                    fp.close()
1958                    G2frame.ErrorDialog('Read Error',
1959                        'File '+filename+' does not start with Name\nFirst line ='
1960                        +l+'\ninvalid file',parent=G2frame)
1961                    return
1962                name = l.split(':')[1].strip()
1963                line = fp.readline().strip().split(':')[1].split()
1964                atNames = [i for i in line]
1965                types = []
1966                coords = []
1967                l = fp.readline().strip()
1968                rbSeq = []
1969                while 'rbSeq' in l:
1970                    tag,vals,lst = l.split(':')
1971                    seq = []
1972                    for t,v in zip((int,int,float),vals.split()): 
1973                        seq.append(t(v))
1974                    seq.append([])
1975                    for num in lst.split():
1976                        seq[-1].append(int(num))
1977                    rbSeq.append(seq)
1978                    l = fp.readline().strip()                   
1979                while l:
1980                    nums = l.strip().split()
1981                    types.append(nums.pop(0))
1982                    t = types[-1]
1983                    if t not in AtInfo:
1984                        Info = G2elem.GetAtomInfo(t)
1985                        AtInfo[t] = [Info['Drad'],Info['Color']]
1986                    coords.append([float(nums.pop(0)) for j in range(3)])
1987                    l = fp.readline().strip()
1988                fp.close()
1989            else:
1990                return       
1991        finally:
1992            dlg.Destroy()
1993        coords = np.array(coords)
1994        rbid = ran.randint(0,sys.maxsize)
1995        namelist = [data['Residue'][key]['RBname'] for key in data['Residue']
1996                        if 'RBname' in data['Residue'][key]]
1997        name = G2obj.MakeUniqueLabel(name,namelist)
1998        data['Residue'][rbid] = {'RBname':name,
1999                'rbXYZ': coords,
2000                'rbRef':[0,1,2,False],
2001                'rbTypes':types, 'atNames':atNames,
2002                'useCount':0,
2003                'rbSeq':rbSeq, 'SelSeq':[0,0],}
2004        data['RBIds']['Residue'].append(rbid)
2005        UpdateResidueRB()
2006   
2007    def AddVectorRB(event):
2008        'Create a new vector rigid body'
2009        AtInfo = data['Vector']['AtInfo']
2010        dlg = G2G.MultiIntegerDialog(G2frame,'New Rigid Body',['No. atoms','No. translations'],[3,1])
2011        if dlg.ShowModal() == wx.ID_OK:
2012            nAtoms,nTrans = dlg.GetValues()
2013            if nAtoms < 3:
2014                dlg.Destroy()
2015                G2G.G2MessageBox(G2frame,'A vector rigid body must have 3 or more atoms')
2016                return
2017            rbid = ran.randint(0,sys.maxsize)
2018            vecMag = [1.0 for i in range(nTrans)]
2019            vecRef = [False for i in range(nTrans)]
2020            vecVal = [np.zeros((nAtoms,3)) for j in range(nTrans)]
2021            rbTypes = ['C' for i in range(nAtoms)]
2022            Info = G2elem.GetAtomInfo('C')
2023            AtInfo['C'] = [Info['Drad'],Info['Color']]
2024            name = 'UNKRB'
2025            namelist = [data['Vector'][key]['RBname'] for key in data['Vector']
2026                        if 'RBname' in data['Vector'][key]]
2027            name = G2obj.MakeUniqueLabel(name,namelist)
2028            data['Vector'][rbid] = {'RBname':name,'VectMag':vecMag,'rbXYZ':np.zeros((nAtoms,3)),
2029                'rbRef':[0,1,2,False],'VectRef':vecRef,'rbTypes':rbTypes,'rbVect':vecVal,'useCount':0}
2030            data['RBIds']['Vector'].append(rbid)
2031        dlg.Destroy()
2032        UpdateVectorRB()
2033
2034    def ExtractPhaseRB(event):
2035        'Extract a rigid body from a file with a phase'
2036        def SetupDrawing(atmData):
2037            '''Add the dicts needed for G2plt.PlotStructure to work to the
2038            reader .Phase object
2039            '''
2040            generalData = atmData['General']
2041            generalData['BondRadii'] = []
2042
2043            G2phG.SetDrawingDefaults(atmData['Drawing'])
2044            atmData['Drawing'].update(
2045                {'oldxy':[0.,0.],'Quaternion':[0.,0.,0.,1.],'cameraPos':150.,
2046                     'viewDir':[0,0,1],'atomPtrs': [2, 1, 6, 17],
2047                     })
2048            atmData['Drawing']['showRigidBodies'] = False
2049            generalData['Map'] = {'MapType':False, 'rho':[]}
2050            generalData['AtomTypes'] = []
2051            generalData['BondRadii'] = []
2052            generalData['AngleRadii'] = []
2053            generalData['vdWRadii'] = []
2054            generalData['Color'] = []
2055            generalData['Isotopes'] = {}
2056            generalData['Isotope'] = {}
2057            cx,ct,cs,cia = generalData['AtomPtrs']
2058            generalData['Mydir'] = G2frame.dirname
2059            for iat,atom in enumerate(atmData['Atoms']):
2060                atom[ct] = atom[ct].lower().capitalize()      #force elem symbol to standard form
2061                if atom[ct] not in generalData['AtomTypes'] and atom[ct] != 'UNK':
2062                    Info = G2elem.GetAtomInfo(atom[ct])
2063                    if not Info:
2064                        atom[ct] = 'UNK'
2065                        continue
2066                    atom[ct] = Info['Symbol'] # N.B. symbol might be changed by GetAtomInfo
2067                    generalData['AtomTypes'].append(atom[ct])
2068                    generalData['Z'] = Info['Z']
2069                    generalData['Isotopes'][atom[ct]] = Info['Isotopes']
2070                    generalData['BondRadii'].append(Info['Drad'])
2071                    generalData['AngleRadii'].append(Info['Arad'])
2072                    generalData['vdWRadii'].append(Info['Vdrad'])
2073                    if atom[ct] in generalData['Isotope']:
2074                        if generalData['Isotope'][atom[ct]] not in generalData['Isotopes'][atom[ct]]:
2075                            isotope = list(generalData['Isotopes'][atom[ct]].keys())[-1]
2076                            generalData['Isotope'][atom[ct]] = isotope
2077                    else:
2078                        generalData['Isotope'][atom[ct]] = 'Nat. Abund.'
2079                        if 'Nat. Abund.' not in generalData['Isotopes'][atom[ct]]:
2080                            isotope = list(generalData['Isotopes'][atom[ct]].keys())[-1]
2081                            generalData['Isotope'][atom[ct]] = isotope
2082                    generalData['Color'].append(Info['Color'])
2083                    # if generalData['Type'] == 'magnetic':
2084                    #     if len(landeg) < len(generalData['AtomTypes']):
2085                    #         landeg.append(2.0)
2086            atmData['Drawing']['Atoms'] = []
2087            for atom in atmData['Atoms']:
2088                atmData['Drawing']['Atoms'].append(G2mth.MakeDrawAtom(atmData,atom))
2089
2090        def onCancel(event,page=0):
2091            'complete or bail out from RB define, cleaning up'
2092            G2frame.rbBook.DeletePage(G2frame.rbBook.FindPage(pagename))
2093            G2frame.rbBook.SetSelection(page)
2094
2095        def Page1():
2096            '''Show the GUI for first stage of the rigid body with all atoms in
2097            phase in crystal coordinates. Select the atoms to go onto the
2098            next stage
2099            '''
2100            def ShowSelection(selections):
2101                'respond to change in atom selections'
2102                ct,cs = [1,8]
2103                generalData = rd.Phase['General']
2104                for i,atom in enumerate(rd.Phase['Drawing']['Atoms']):
2105                    if i in selections:
2106                        factor = 1
2107                    else:
2108                        factor = 2.5
2109                    atNum = generalData['AtomTypes'].index(atom[ct]) 
2110                    atom[cs] = list(np.array(generalData['Color'][atNum])//factor) 
2111                draw(*drawArgs)
2112            def onPage1OK(event):
2113                '1st section has been completed, move onto next'
2114                G2frame.G2plotNB.Delete(rd.Phase['General']['Name'])
2115                GetCoords(atmsel)
2116                Page2()
2117
2118            if 'macromolecular' == rd.Phase['General']['Type']:
2119                # for PDB imports, lets see if a quick reformat of atoms list will work
2120                rd.Phase['Atoms'] = [a[3:] for a in rd.Phase['Atoms']]
2121                rd.Phase['General']['AtomPtrs']  = [i-3 for i in rd.Phase['General']['AtomPtrs']]
2122                rd.Phase['General']['Type'] = 'nuclear'
2123            SetupDrawing(rd.Phase) # add information to reader object to allow plotting
2124            atomlist = [atom[0] for atom in rd.Phase['Atoms']]
2125            atmsel = list(range(len(rd.Phase['Atoms'])))
2126            # broken -- # why no bonds?
2127            #for atm in rd.Phase['Drawing']['Atoms']:
2128            #    atm[6] = 'balls & sticks'
2129
2130            draw,drawArgs = G2plt.PlotStructure(G2frame,rd.Phase,True)
2131            ShowSelection(atmsel)
2132
2133            if G2frame.rbBook.FindPage(pagename) is not None:
2134                G2frame.rbBook.DeletePage(G2frame.rbBook.FindPage(pagename))
2135
2136            RBImp = wx.ScrolledWindow(G2frame.rbBook)
2137            RBImpPnl = wx.Panel(RBImp)
2138            G2frame.rbBook.AddPage(RBImp,pagename)
2139            G2frame.rbBook.SetSelection(G2frame.rbBook.FindPage(pagename))
2140
2141            HelpInfo = '''
2142This window shows all the atoms that were read from the
2143selected phase file. Select the atoms that will be used in the
2144rigid body processing (this may include atoms needed to
2145define an axis or origin that will not be included in the
2146eventual rigid body.) Note that in the plot window,
2147unselected atoms appear much darker than selected atoms.
2148'''
2149            mainSizer = G2G.G2MultiChoiceWindow(RBImpPnl,
2150                            'Select atoms to import',
2151                            atomlist,atmsel,OnChange=ShowSelection,
2152                            helpText=HelpInfo)
2153
2154            # OK/Cancel buttons       
2155            btnsizer = wx.StdDialogButtonSizer()
2156            OKbtn = wx.Button(RBImpPnl, wx.ID_OK, 'Continue')
2157            OKbtn.SetDefault()
2158            btnsizer.AddButton(OKbtn)
2159            OKbtn.Bind(wx.EVT_BUTTON,onPage1OK)
2160            btn = wx.Button(RBImpPnl, wx.ID_CANCEL)
2161            btn.Bind(wx.EVT_BUTTON,onCancel)
2162            btnsizer.AddButton(btn)
2163            btnsizer.Realize()
2164            mainSizer.Add(btnsizer,0,wx.ALIGN_CENTER,50)
2165
2166            RBImpPnl.SetSizer(mainSizer,True)
2167
2168            mainSizer.Layout()   
2169            Size = mainSizer.GetMinSize()
2170            Size[0] += 40
2171            Size[1] = max(Size[1],G2frame.GetSize()[1]-200) + 20
2172            RBImpPnl.SetSize(Size)
2173            RBImp.SetScrollbars(10,10,Size[0]/10-4,Size[1]/10-1)
2174            RBImp.Scroll(0,0)
2175
2176        def Page2():
2177            '''Show the GUI for the second stage, where selected atoms are
2178            now in Cartesian space, manipulate the axes and export selected
2179            atoms to a vector or residue rigid body.
2180            '''
2181            def UpdateDraw(event=None):
2182                'Called when info changes in grid, replots'
2183                UpdateVectorBody(rbData)
2184                DrawCallback()
2185               
2186            def onSetAll(event):
2187                'Set all atoms as selected'
2188                grid.completeEdits()
2189                for i in range(len(rd.Phase['RBselection'])):
2190                    rd.Phase['RBselection'][i] = 1 # table needs 0/1 for T/F
2191                grid.ForceRefresh()
2192                UpdateDraw()
2193               
2194            def onToggle(event):
2195                'Toggles selection state for all atoms'
2196                grid.completeEdits()
2197                for i in range(len(rd.Phase['RBselection'])):
2198                    rd.Phase['RBselection'][i] = int(not rd.Phase['RBselection'][i])
2199                grid.ForceRefresh()
2200                UpdateDraw()
2201               
2202            def onSetOrigin(event):
2203                'Resets origin to midpoint between all selected atoms'
2204                grid.completeEdits()
2205                center = np.array([0.,0.,0.])
2206                count = 0
2207                for i in range(len(rd.Phase['RBselection'])):
2208                    if rd.Phase['RBselection'][i]:
2209                        count += 1
2210                        center += rd.Phase['RBcoords'][i]
2211                if count:
2212                    rd.Phase['RBcoords'] -= center/count
2213                grid.ForceRefresh()
2214                UpdateDraw()
2215               
2216            def onSetX(event):
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 not count:
2225                    G2G.G2MessageBox(G2frame,'No atoms selected',
2226                                    'Selection required')
2227                    return
2228                XYZP = center/count
2229                if np.sqrt(sum(XYZP**2)) < 0.1:
2230                    G2G.G2MessageBox(G2frame,
2231                            'The selected atom(s) are too close to the origin',
2232                            'near origin')
2233                    return
2234                if bntOpts['direction'] == 'y':
2235                    YP = XYZP / np.sqrt(np.sum(XYZP**2))
2236                    ZP = np.cross((1,0,0),YP)
2237                    if sum(ZP*ZP) < .1: # pathological condition: Y' along X
2238                        ZP = np.cross((0,0,1),YP)
2239                    XP = np.cross(YP,ZP)
2240                elif bntOpts['direction'] == 'z':
2241                    ZP = XYZP / np.sqrt(np.sum(XYZP**2))
2242                    XP = np.cross((0,1,0),ZP)
2243                    if sum(XP*XP) < .1: # pathological condition: X' along Y
2244                        XP = np.cross((0,0,1),ZP)
2245                    YP = np.cross(ZP,XP)
2246                else:
2247                    XP = XYZP / np.sqrt(np.sum(XYZP**2))
2248                    YP = np.cross((0,0,1),XP)
2249                    if sum(YP*YP) < .1: # pathological condition: X' along Z
2250                        YP = np.cross((0,1,0),XP)
2251                    ZP = np.cross(XP,YP)
2252                trans = np.array((XP,YP,ZP))
2253                # update atoms in place
2254                rd.Phase['RBcoords'][:] = np.inner(trans,rd.Phase['RBcoords']).T
2255                grid.ForceRefresh()
2256                UpdateDraw()
2257
2258            def onSetPlane(event): 
2259                '''Compute least-squares plane for selected atoms;
2260                move atoms so that LS plane aligned with x-y plane,
2261                with minimum required change to x
2262                '''
2263                grid.completeEdits()
2264                selList = [i==1 for i in rd.Phase['RBselection']]
2265                XYZ = rd.Phase['RBcoords'][selList]
2266                if len(XYZ) < 3: 
2267                    G2G.G2MessageBox(G2frame,'A plane requires three or more atoms',
2268                                     'Need more atoms')
2269                    return
2270                # fit 3 ways (in case of singularity) and take result with lowest residual
2271                X,Y,Z = [XYZ[:,i] for i in (0,1,2)]
2272                XZ = copy.copy(XYZ)
2273                XZ[:,1] = 1
2274                (a,d,b), resd, rank, sing = nl.lstsq(XZ, -Y)
2275                resid_min = resd
2276                normal = a,1,b
2277                YZ = copy.copy(XYZ)
2278                YZ[:,0] = 1
2279                (d,a,b), resd, rank, sing = nl.lstsq(YZ, -X)
2280                if resid_min > resd:
2281                    resid_min = resd
2282                    normal = 1,a,b
2283                XY = copy.copy(XYZ)
2284                XY[:,2] = 1
2285                (a,b,d), resd, rank, sing = nl.lstsq(XY, -Z)
2286                if resid_min > resd:
2287                    resid_min = resd
2288                    normal = a,b,1
2289                # solve for  ax + bx + z + c = 0 or equivalently ax + bx + c = -z
2290                # try:
2291                # except:
2292                #     G2G.G2MessageBox(G2frame,
2293                #             'Error computing plane; are atoms in a line?',
2294                #             'Computation error')
2295                #     return
2296                if bntOpts['plane'] == 'xy':
2297                    # new coordinate system is
2298                    #   ZP, z' normal to plane
2299                    #   YP, y' = z' cross x (= [0,1,-b])
2300                    #   XP, x' = (z' cross x) cross z'
2301                    # this puts XP as close as possible to X with XP & YP in plane
2302                    ZP = np.array(normal)
2303                    ZP /= np.sqrt(np.sum(ZP**2))
2304                    YP = np.cross(ZP,[1,0,0])
2305                    if sum(YP*YP) < .1: # pathological condition: z' along x
2306                        YP = np.cross(ZP,[0,1,0])
2307                    YP /= np.sqrt(np.sum(YP**2))
2308                    XP = np.cross(YP,ZP)
2309                elif bntOpts['plane'] == 'yz':
2310                    # new coordinate system is
2311                    #   XP, x' normal to plane
2312                    #   ZP, z' = x' cross y
2313                    #   YP, y' = (x' cross y) cross x'
2314                    # this puts y' as close as possible to y with z' & y' in plane
2315                    XP = np.array(normal)
2316                    XP /= np.sqrt(np.sum(XP**2))
2317                    ZP = np.cross(XP,[0,1,0])
2318                    if sum(ZP*ZP) < .1: # pathological condition: x' along y
2319                        ZP = np.cross(XP,(0,0,1))
2320                    ZP /= np.sqrt(np.sum(ZP**2))
2321                    YP = np.cross(ZP,XP)
2322                elif bntOpts['plane'] == 'xz':
2323                    # new coordinate system is
2324                    #   YP, y' normal to plane
2325                    #   ZP, z' = x cross y'
2326                    #   XP, y' = (x cross y') cross z'
2327                    # this puts XP as close as possible to X with XP & YP in plane
2328                    YP = np.array(normal)
2329                    YP /= np.sqrt(np.sum(YP**2))
2330                    ZP = np.cross([1,0,0],YP)
2331                    if sum(ZP*ZP) < .1: # pathological condition: y' along x
2332                        ZP = np.cross([0,1,0],YP)
2333                    ZP /= np.sqrt(np.sum(ZP**2))
2334                    XP = np.cross(YP,ZP)
2335                else:
2336                    print('unexpected plane',bntOpts['plane'])
2337                    return
2338                trans = np.array((XP,YP,ZP))
2339                # update atoms in place
2340                rd.Phase['RBcoords'][:] = np.inner(trans,rd.Phase['RBcoords']).T
2341                grid.ForceRefresh()
2342                UpdateDraw()
2343
2344            def onWriteXYZ(event):
2345                '''Writes selected atoms in a .xyz file for use in Avogadro, etc.
2346                '''
2347                grid.completeEdits()
2348                center = np.array([0.,0.,0.])
2349                count = 0
2350                for i in range(len(rd.Phase['RBselection'])):
2351                    if rd.Phase['RBselection'][i]:
2352                        count += 1
2353                        center += rd.Phase['RBcoords'][i]
2354                if count:
2355                    center /= count
2356                else:
2357                    print('nothing selected')
2358                    return
2359                obj = G2IO.ExportBaseclass(G2frame,'XYZ','.xyz')
2360                #obj.InitExport(None)
2361                if obj.ExportSelect():    # set export parameters; ask for file name
2362                    return
2363                obj.OpenFile()
2364                obj.Write(str(count))
2365                obj.Write('')
2366                for i in range(len(rd.Phase['RBselection'])):
2367                    if rd.Phase['RBselection'][i]:
2368                        line = ' ' + rd.Phase['RBtypes'][i]
2369                        for xyz in rd.Phase['RBcoords'][i]:
2370                            line += ' ' + str(xyz)
2371                        obj.Write(line)
2372                obj.CloseFile()
2373                #GSASIIpath.IPyBreak()
2374               
2375            def onAddVector(event):
2376                '''Adds selected atoms as a new vector rigid body.
2377                Closes out the importer tab when done.
2378                '''
2379                grid.completeEdits()
2380                name = os.path.splitext(os.path.split(filename)[1])[0]
2381                namelist = [data['Vector'][key]['RBname'] for key in
2382                        data['Vector'] if 'RBname' in data['Vector'][key]]
2383                name = G2obj.MakeUniqueLabel(name,namelist)
2384                rb = MakeVectorBody(name)
2385                UpdateVectorBody(rb,True)
2386                if len(rb['rbTypes']) < 3: return # must have at least 3 atoms
2387                rbid = ran.randint(0,sys.maxsize)
2388                data['Vector'][rbid] = rb
2389                data['RBIds']['Vector'].append(rbid)
2390                for t in rb['rbTypes']:
2391                    if t in data['Vector']['AtInfo']: continue
2392                    Info = G2elem.GetAtomInfo(t)
2393                    data['Vector']['AtInfo'][t] = [Info['Drad'],Info['Color']]
2394                G2frame.G2plotNB.Delete('Rigid body')
2395                onCancel(event,0)
2396               
2397            def onAddResidue(event):
2398                '''Adds selected atoms as a new residue rigid body.
2399                Closes out the importer tab when done.
2400                '''
2401                grid.completeEdits()
2402                name = os.path.split(filename)[1]
2403                rbXYZ = []
2404                rbTypes = []
2405                atNames = []
2406                for i in rd.Phase['RBindex']:
2407                    if rd.Phase['RBselection'][i]:
2408                        rbXYZ.append(rd.Phase['RBcoords'][i])
2409                        rbTypes.append(rd.Phase['RBtypes'][i])
2410                        atNames.append(rd.Phase['RBlbls'][i])
2411                if len(rbTypes) < 3: return # must have at least 3 atoms
2412                rbXYZ = np.array(rbXYZ)
2413                rbid = ran.randint(0,sys.maxsize)
2414                namelist = [data['Residue'][key]['RBname'] for key in
2415                        data['Residue'] if 'RBname' in data['Residue'][key]]
2416                name = G2obj.MakeUniqueLabel(name,namelist)
2417                data['Residue'][rbid] = {'RBname':name,'rbXYZ':rbXYZ,
2418                    'rbTypes':rbTypes,'atNames':atNames,'rbRef':[0,1,2,False],
2419                    'rbSeq':[],'SelSeq':[0,0],'useCount':0}
2420                data['RBIds']['Residue'].append(rbid)
2421                for t in rbTypes:
2422                    if t in data['Residue']['AtInfo']: continue
2423                    Info = G2elem.GetAtomInfo(t)
2424                    data['Residue']['AtInfo'][t] = [Info['Drad'],Info['Color']]
2425
2426                print ('Rigid body added')
2427                G2frame.G2plotNB.Delete('Rigid body')
2428                onCancel(event,1)
2429
2430            if G2frame.rbBook.FindPage(pagename) is not None:
2431                G2frame.rbBook.DeletePage(G2frame.rbBook.FindPage(pagename))
2432            RBImp = wx.ScrolledWindow(G2frame.rbBook)
2433            RBImpPnl = wx.Panel(RBImp)
2434            G2frame.rbBook.AddPage(RBImp,pagename)
2435            G2frame.rbBook.SetSelection(G2frame.rbBook.FindPage(pagename))
2436            AtInfo = {}
2437            for t in rd.Phase['RBtypes']:
2438                if t in AtInfo: continue
2439                Info = G2elem.GetAtomInfo(t)
2440                AtInfo[t] = [Info['Drad'],Info['Color']]
2441            plotDefaults = {'oldxy':[0.,0.],'Quaternion':[0.,0.,0.,1.],'cameraPos':30.,'viewDir':[0,0,1],}
2442
2443            rd.Phase['RBindex'] = list(range(len(rd.Phase['RBtypes'])))
2444            rd.Phase['RBselection'] = len(rd.Phase['RBtypes']) * [1]
2445            name = 'UNKRB'
2446            namelist = [data['Vector'][key]['RBname'] for key in
2447                        data['Vector'] if 'RBname' in data['Vector'][key]]
2448            name = G2obj.MakeUniqueLabel(name,namelist)
2449            rbData = MakeVectorBody()
2450            DrawCallback = G2plt.PlotRigidBody(G2frame,'Vector',AtInfo,rbData,plotDefaults)
2451
2452            mainSizer = wx.BoxSizer(wx.HORIZONTAL)
2453            btnSizer = wx.BoxSizer(wx.VERTICAL)
2454            helpText = '''
2455In this window, if wanted,
2456one can select one or more atoms and use them
2457to define an origin, a specified axis or place the selected atoms into
2458a selected plane. (Different sets of atoms can be used for each
2459operation.)
2460%%Once that is done, atoms can be selected and can be exported in a
2461"XYZ" file for use in a program such as Avogadro or can be used to
2462create a Vector or Residue rigid body.
2463'''
2464            btnSizer.Add(G2G.HelpButton(RBImpPnl,helpText,wrap=400),
2465                             0,wx.ALIGN_RIGHT)
2466            btnSizer.Add(wx.StaticText(RBImpPnl,wx.ID_ANY,'Reorder atoms by dragging'),0,wx.ALL)
2467            btnSizer.Add((-1,15))
2468            btn = wx.Button(RBImpPnl, wx.ID_ANY, 'Set All')
2469            btn.Bind(wx.EVT_BUTTON,onSetAll)
2470            btnSizer.Add(btn,0,wx.ALIGN_CENTER)
2471            btn = wx.Button(RBImpPnl, wx.ID_ANY, 'Toggle')
2472            btn.Bind(wx.EVT_BUTTON,onToggle)
2473            btnSizer.Add(btn,0,wx.ALIGN_CENTER)
2474            btnSizer.Add((-1,15))
2475            btnSizer.Add(wx.StaticText(RBImpPnl,wx.ID_ANY,'Reorient using selected\natoms...'),0,wx.ALL)
2476            btnSizer.Add((-1,5))
2477            btn = wx.Button(RBImpPnl, wx.ID_ANY, 'Set origin')
2478            btn.Bind(wx.EVT_BUTTON,onSetOrigin)
2479            btnSizer.Add(btn,0,wx.ALIGN_CENTER)
2480
2481            bntOpts = {'plane':'xy','direction':'x'}
2482            inSizer = wx.BoxSizer(wx.HORIZONTAL)
2483            btn = wx.Button(RBImpPnl, wx.ID_ANY, 'Place in plane')
2484            btn.Bind(wx.EVT_BUTTON,onSetPlane)
2485            inSizer.Add(btn)
2486            inSizer.Add(G2G.G2ChoiceButton(RBImpPnl,('xy','yz','xz'),None,None,bntOpts,'plane'))
2487            btnSizer.Add(inSizer,0,wx.ALIGN_CENTER)
2488           
2489            inSizer = wx.BoxSizer(wx.HORIZONTAL)
2490            btn = wx.Button(RBImpPnl, wx.ID_ANY, 'Define as')
2491            btn.Bind(wx.EVT_BUTTON,onSetX)
2492            inSizer.Add(btn)
2493            inSizer.Add(G2G.G2ChoiceButton(RBImpPnl,('x','y','z'),None,None,bntOpts,'direction'))
2494            btnSizer.Add(inSizer,0,wx.ALIGN_CENTER)
2495           
2496            btnSizer.Add((-1,15))
2497            btnSizer.Add(wx.StaticText(RBImpPnl,wx.ID_ANY,'Use selected atoms to\ncreate...'),0,wx.ALL)
2498            btnSizer.Add((-1,5))
2499            btn = wx.Button(RBImpPnl, wx.ID_ANY, 'export as xyz')
2500            btn.Bind(wx.EVT_BUTTON,onWriteXYZ)
2501            btnSizer.Add(btn,0,wx.ALIGN_CENTER)
2502            btnSizer.Add((-1,10))
2503            btn = wx.Button(RBImpPnl, wx.ID_ANY, 'a Vector Body')
2504            btn.Bind(wx.EVT_BUTTON,onAddVector)
2505            btnSizer.Add(btn,0,wx.ALIGN_CENTER)
2506            btn = wx.Button(RBImpPnl, wx.ID_ANY, 'a Residue Body')
2507            btn.Bind(wx.EVT_BUTTON,onAddResidue)
2508            btnSizer.Add(btn,0,wx.ALIGN_CENTER)
2509            btn = wx.Button(RBImpPnl, wx.ID_CANCEL)
2510            btn.Bind(wx.EVT_BUTTON,onCancel)
2511            btnSizer.Add((-1,10))
2512            btnSizer.Add(btn,0,wx.ALIGN_CENTER)
2513
2514            mainSizer.Add(btnSizer)
2515            mainSizer.Add((5,5))
2516            grid = DragableRBGrid(RBImpPnl,rd.Phase,UpdateDraw)
2517            mainSizer.Add(grid)
2518            RBImpPnl.SetSizer(mainSizer,True)
2519            mainSizer.Layout()   
2520            Size = mainSizer.GetMinSize()
2521            Size[0] += 40
2522            Size[1] = max(Size[1],G2frame.GetSize()[1]-200) + 20
2523            RBImpPnl.SetSize(Size)
2524            RBImp.SetScrollbars(10,10,Size[0]/10-4,Size[1]/10-1)
2525            RBImp.Scroll(0,0)
2526
2527        def GetCoords(atmsel):
2528            '''Create orthogonal coordinates for selected atoms.
2529            Place the origin at the center of the body
2530            '''
2531            atms = rd.Phase['Atoms']
2532            cell = rd.Phase['General']['Cell'][1:7]
2533            Amat,Bmat = G2lat.cell2AB(cell)
2534            rd.Phase['RBcoords'] = np.array([np.inner(Amat,atms[i][3:6]) for i in atmsel])
2535            rd.Phase['RBcoords'] -= rd.Phase['RBcoords'].mean(axis=0)  # origin to middle
2536            rd.Phase['RBtypes'] = [atms[i][1] for i in atmsel]
2537            rd.Phase['RBlbls'] = [atms[i][0] for i in atmsel]
2538
2539        def UpdateVectorBody(rb,useSelection=False):
2540            '''Put the atoms in order to pass for plotting or for storage as
2541            a vector rigid body.
2542
2543            :param dict rb: rigid body contents created in :func:`MakeVectorBody`
2544            :param bool useSelection: True if the rd.Phase['RBselection']
2545              values will be used to select which atoms are included in the
2546              rigid body. If False (default) they are included in rb
2547              and are used for plotting.         
2548            '''
2549            coordlist = []
2550            typeslist = []
2551            sellist = []
2552            for i in rd.Phase['RBindex']:
2553                use = True
2554                if useSelection and not rd.Phase['RBselection'][i]: use = False
2555                if use:
2556                    coordlist.append(rd.Phase['RBcoords'][i])
2557                    typeslist.append(rd.Phase['RBtypes'][i])
2558                    sellist.append(rd.Phase['RBselection'][i])
2559            coordlist = np.array(coordlist)
2560            rb['rbXYZ'] = coordlist
2561            rb['rbVect'] = [coordlist]
2562            rb['rbTypes'] = typeslist
2563            if not useSelection:
2564                rb['Selection'] = sellist
2565            elif 'Selection' in rb:
2566                del rb['Selection']
2567
2568        def MakeVectorBody(name=''):
2569            '''Make the basic vector rigid body dict (w/o coordinates) used for
2570            export and for plotting
2571            '''
2572            vecMag = [1.0]
2573            vecRef = [False]
2574            rb = {'RBname':name,'VectMag':vecMag,
2575                    'rbRef':[0,1,2,False],'VectRef':vecRef,
2576                    'useCount':0}
2577            UpdateVectorBody(rb)
2578            return rb
2579
2580        # too lazy to figure out why wx crashes
2581        if wx.__version__.split('.')[0] != '4':
2582            wx.MessageBox('Sorry, wxPython 4.x is required to run this command',
2583                                  caption='Update Python',
2584                                  style=wx.ICON_EXCLAMATION)
2585            return
2586        if platform.python_version()[:1] == '2':
2587            wx.MessageBox('Sorry, Python >=3.x is required to run this command',
2588                                  caption='Update Python',
2589                                  style=wx.ICON_EXCLAMATION)
2590            return
2591
2592        # get importer type and a phase file of that type
2593        G2sc.LoadG2fil()
2594        choices = [rd.formatName for  rd in G2sc.Readers['Phase']] 
2595        dlg = G2G.G2SingleChoiceDialog(G2frame,'Select the format of the file',
2596                                     'select format',choices)
2597        dlg.CenterOnParent()
2598        try:
2599            if dlg.ShowModal() == wx.ID_OK:
2600                col = dlg.GetSelection()
2601            else:
2602                col = None
2603                return
2604        finally:
2605            dlg.Destroy()
2606        reader = G2sc.Readers['Phase'][col]
2607
2608        choices = reader.formatName + " file ("
2609        w = ""
2610        for extn in reader.extensionlist:
2611            if w != "": w += ";"
2612            w += "*" + extn
2613        choices += w + ")|" + w
2614        #choices += "|zip archive (.zip)|*.zip"
2615        if not reader.strictExtension:
2616            choices += "|any file (*.*)|*.*"
2617        typ = '( type '+reader.formatName+')'
2618        filelist = G2G.GetImportFile(G2frame,
2619                        message="Choose phase input file"+typ,
2620                        defaultFile="",wildcard=choices,style=wx.FD_OPEN)
2621        if len(filelist) != 1: return
2622
2623        # read in the phase file
2624        filename = filelist[0]
2625        rd = reader
2626        with open(filename, 'r'):
2627            rd.ReInitialize()
2628            rd.errors = ""
2629            if not rd.ContentsValidator(filename):   # Report error
2630                G2fl.G2Print("Warning: File {} has a validation error".format(filename))
2631                return
2632            if len(rd.selections) > 1:
2633                print("File {} has {} phases. This is unexpected."
2634                                    .format(filename,len(rd.selections)))
2635                return
2636
2637            rd.objname = os.path.basename(filename)
2638            try:
2639                rd.Reader(filename)
2640            except Exception as msg:
2641                G2fl.G2Print("Warning: read of file {} failed\n{}".format(
2642                    filename,rd.errors))
2643                if GSASIIpath.GetConfigValue('debug'):
2644                    print(msg)
2645                    import traceback
2646                    print (traceback.format_exc())
2647                    GSASIIpath.IPyBreak()
2648                return
2649
2650        pagename = 'Rigid body importer'
2651        Page1()
2652        return
2653
2654    def AddVectTrans(event):
2655        'Add a translation to an existing vector rigid body'
2656        choices = []
2657        rbIdlist = []
2658        for rbid in data['RBIds']['Vector']:
2659            if rbid != 'AtInfo':
2660                rbIdlist.append(rbid)
2661                choices.append(data['Vector'][rbid]['RBname'])
2662        if len(choices) == 0:
2663            G2G.G2MessageBox(G2frame,'No Vector Rigid Bodies found',
2664                                 'No VR Bodies')
2665            return
2666        elif len(choices) == 1:
2667            rbid = rbIdlist[0]
2668        else:
2669            dlg = G2G.G2SingleChoiceDialog(G2frame,'Select the rigid body to save',
2670                                  'select format',choices)
2671            try:
2672                if dlg.ShowModal() == wx.ID_OK:
2673                    rbid = rbIdlist[dlg.GetSelection()]
2674                else:
2675                    return
2676            finally:
2677                dlg.Destroy()
2678        data['Vector'][rbid]['VectMag'] += [1.0]
2679        data['Vector'][rbid]['VectRef'] += [False]
2680        nAtoms = len(data['Vector'][rbid]['rbXYZ'])
2681        data['Vector'][rbid]['rbVect'] += [np.zeros((nAtoms,3))]
2682        UpdateVectorRB()
2683       
2684    def SaveVectorRB(event):
2685        choices = []
2686        rbIdlist = []
2687        for rbid in data['RBIds']['Vector']:
2688            if rbid != 'AtInfo':
2689                rbIdlist.append(rbid)
2690                choices.append(data['Vector'][rbid]['RBname'])
2691        if len(choices) == 0:
2692            G2G.G2MessageBox(G2frame,'No Vector Rigid Bodies found',
2693                                 'No VR Bodies')
2694            return
2695        elif len(choices) == 1:
2696            rbid = rbIdlist[0]
2697        else:
2698            dlg = G2G.G2SingleChoiceDialog(G2frame,'Select the rigid body to save',
2699                                  'select format',choices)
2700            try:
2701                if dlg.ShowModal() == wx.ID_OK:
2702                    rbid = rbIdlist[dlg.GetSelection()]
2703                else:
2704                    return
2705            finally:
2706                dlg.Destroy()
2707             
2708        pth = G2G.GetExportPath(G2frame)
2709        dlg = wx.FileDialog(G2frame, 'Choose file to save vector rigid body',
2710            pth, '', 'VRB files (*.vecbody)|*.vecbody',
2711            wx.FD_SAVE|wx.FD_OVERWRITE_PROMPT)
2712        try:
2713            if dlg.ShowModal() == wx.ID_OK:
2714                filename = dlg.GetPath()
2715                filename = os.path.splitext(filename)[0]+'.vecbody'  # set extension
2716                fp = open(filename,'w')
2717                fp.write('Name: '+data['Vector'][rbid]['RBname']+'\n')
2718                fp.write('Trans: ')
2719                for i in data['Vector'][rbid]['VectMag']:
2720                    fp.write(str(i)+" ") 
2721                fp.write('\n')
2722                ntrans = len(data['Vector'][rbid]['VectMag'])
2723                for i,sym in enumerate(data['Vector'][rbid]['rbTypes']):
2724                    fp.write("{:3s}".format(sym))
2725                    for j in range(ntrans):
2726                        fp.write('{:8.5f}{:9.5f}{:9.5f}   '
2727                            .format(*data['Vector'][rbid]['rbVect'][j][i]))
2728                    fp.write('\n')
2729                fp.close()
2730                print ('Vector rigid body saved to: '+filename)
2731        finally:
2732            dlg.Destroy()
2733           
2734    def ReadVectorRB(event):
2735        AtInfo = data['Vector']['AtInfo']
2736        pth = G2G.GetExportPath(G2frame)
2737        dlg = wx.FileDialog(G2frame, 'Choose file to read vector rigid body',
2738            pth, '', 'VRB files (*.vecbody)|*.vecbody',
2739            wx.FD_OPEN)
2740        try:
2741            if dlg.ShowModal() == wx.ID_OK:
2742                filename = dlg.GetPath()
2743                filename = os.path.splitext(filename)[0]+'.vecbody'  # set extension
2744                fp = open(filename,'r')
2745                l = fp.readline().strip()
2746                if 'Name' not in l:
2747                    fp.close()
2748                    G2frame.ErrorDialog('Read Error',
2749                        'File '+filename+' does not start with Name\nFirst line ='
2750                        +l+'\ninvalid file',parent=G2frame)
2751                    return
2752                name = l.split(':')[1].strip()
2753                trans = fp.readline().strip().split(':')[1].split()
2754                vecMag = [float(i) for i in trans]
2755                ntrans = len(trans)
2756                vecs = [[] for i in range(ntrans)]
2757                types = []
2758                l = fp.readline().strip()
2759                while l:
2760                    nums = l.strip().split()
2761                    types.append(nums.pop(0))
2762                    t = types[-1]
2763                    if t not in AtInfo:
2764                        Info = G2elem.GetAtomInfo(t)
2765                        AtInfo[t] = [Info['Drad'],Info['Color']]
2766                    for i in range(ntrans):
2767                        vecs[i].append([float(nums.pop(0)) for j in range(3)])
2768                    l = fp.readline().strip()
2769                fp.close()
2770            else:
2771                return       
2772        finally:
2773            dlg.Destroy()
2774        natoms = len(types)
2775        vecs = [np.array(vecs[i]) for i in range(ntrans)]
2776        rbid = ran.randint(0,sys.maxsize)
2777        namelist = [data['Vector'][key]['RBname'] for key in data['Vector']
2778                        if 'RBname' in data['Vector'][key]]
2779        name = G2obj.MakeUniqueLabel(name,namelist)
2780        data['Vector'][rbid] = {'RBname':name,'VectMag':vecMag,
2781                'rbXYZ':np.zeros((natoms,3)),
2782                'rbRef':[0,1,2,False],'VectRef':ntrans*[False],
2783                'rbTypes':types,
2784                'rbVect':vecs,'useCount':0}
2785        data['RBIds']['Vector'].append(rbid)
2786        UpdateVectorRB()
2787       
2788    def AddResidueRB(event):
2789        global resRBsel
2790        AtInfo = data['Residue']['AtInfo']
2791        macro = getMacroFile('rigid body')
2792        if not macro:
2793            return
2794        macStr = macro.readline()
2795        while macStr:
2796            items = macStr.split()
2797            if 'I' == items[0]:
2798                resRBsel = ran.randint(0,sys.maxsize)
2799                rbName = items[1]
2800                rbTypes = []
2801                rbXYZ = []
2802                rbSeq = []
2803                atNames = []
2804                nAtms,nSeq,nOrig,mRef,nRef = [int(items[i]) for i in [2,3,4,5,6]]
2805                for iAtm in range(nAtms):
2806                    macStr = macro.readline().split()
2807                    atName = macStr[0]
2808                    atType = macStr[1]
2809                    atNames.append(atName)
2810                    rbXYZ.append([float(macStr[i]) for i in [2,3,4]])
2811                    rbTypes.append(atType)
2812                    if atType not in AtInfo:
2813                        Info = G2elem.GetAtomInfo(atType)
2814                        AtInfo[atType] = [Info['Drad'],Info['Color']]
2815                rbXYZ = np.array(rbXYZ)-np.array(rbXYZ[nOrig-1])
2816                for iSeq in range(nSeq):
2817                    macStr = macro.readline().split()
2818                    mSeq = int(macStr[0])
2819                    for jSeq in range(mSeq):
2820                        macStr = macro.readline().split()
2821                        iBeg = int(macStr[0])-1
2822                        iFin = int(macStr[1])-1
2823                        angle = 0.0
2824                        nMove = int(macStr[2])
2825                        iMove = [int(macStr[i])-1 for i in range(3,nMove+3)]
2826                        rbSeq.append([iBeg,iFin,angle,iMove])
2827                namelist = [data['Residue'][key]['RBname'] for key in
2828                           data['Residue'] if 'RBname' in data['Residue'][key]]
2829                rbName = G2obj.MakeUniqueLabel(rbName,namelist)
2830                data['Residue'][resRBsel] = {'RBname':rbName,'rbXYZ':rbXYZ,'rbTypes':rbTypes,
2831                    'atNames':atNames,'rbRef':[nOrig-1,mRef-1,nRef-1,True],'rbSeq':rbSeq,
2832                    'SelSeq':[0,0],'useCount':0,'molCent':None}
2833                data['RBIds']['Residue'].append(resRBsel)
2834                print ('Rigid body '+rbName+' added')
2835            macStr = macro.readline()
2836        macro.close()
2837        UpdateResidueRB()
2838       
2839    def ImportResidueRB():
2840        global resRBsel
2841        AtInfo = data['Residue']['AtInfo']
2842        text,ext = getTextFile()
2843        if not text:
2844            return
2845        resRBsel = ran.randint(0,sys.maxsize)
2846        rbTypes = []
2847        rbXYZ = []
2848        atNames = []
2849        txtStr = text.readline()
2850        if 'xyz' in ext:
2851            txtStr = text.readline()
2852            txtStr = text.readline()
2853        elif 'mol2' in ext:
2854            while 'ATOM' not in txtStr:
2855                txtStr = text.readline()
2856            txtStr = text.readline()
2857        elif 'pdb' in ext:
2858            while 'ATOM' not in txtStr[:6] and 'HETATM' not in txtStr[:6]:
2859                txtStr = text.readline()
2860        items = txtStr.split()
2861        nat = 1
2862        while len(items):
2863            if 'txt' in ext:
2864                atName = items[0]
2865                atType = items[1]
2866                rbXYZ.append([float(items[i]) for i in [2,3,4]])
2867            elif 'xyz' in ext:
2868                atType = items[0]
2869                rbXYZ.append([float(items[i]) for i in [1,2,3]])
2870                atName = '%s%d'%(atType,nat)
2871            elif 'mol2' in ext:
2872                atType = items[1]
2873                atName = items[1]+items[0]
2874                rbXYZ.append([float(items[i]) for i in [2,3,4]])
2875            elif 'pdb' in ext:
2876                atType = items[-1]
2877                if not items[2][-1].isnumeric():
2878                    atName = '%s%d'%(items[2],nat)
2879                else:
2880                    atName = '5s'%items[2]
2881                xyz = txtStr[30:55].split()                   
2882                rbXYZ.append([float(x) for x in xyz])
2883            atNames.append(atName)
2884            rbTypes.append(atType)
2885            if atType not in AtInfo:
2886                Info = G2elem.GetAtomInfo(atType)
2887                AtInfo[atType] = [Info['Drad'],Info['Color']]
2888            txtStr = text.readline()
2889            if 'mol2' in ext and 'BOND' in txtStr:
2890                break
2891            if 'pdb' in ext and ('ATOM' not in txtStr[:6] and 'HETATM' not in txtStr[:6]):
2892                break
2893            items = txtStr.split()
2894            nat += 1
2895        if len(atNames) < 3:
2896            G2G.G2MessageBox(G2frame,'Not enough atoms in rigid body; must be 3 or more')
2897        else:
2898            rbXYZ = np.array(rbXYZ)-np.array(rbXYZ[0])
2899            Xxyz = rbXYZ[1]
2900            X = Xxyz/np.sqrt(np.sum(Xxyz**2))
2901            Yxyz = rbXYZ[2]
2902            Y = Yxyz/np.sqrt(np.sum(Yxyz**2))
2903            Mat = G2mth.getRBTransMat(X,Y)
2904            rbXYZ = np.inner(Mat,rbXYZ).T
2905            name = 'UNKRB'
2906            namelist = [data['Residue'][key]['RBname'] for key in data['Residue']
2907                        if 'RBname' in data['Residue'][key]]
2908            name = G2obj.MakeUniqueLabel(name,namelist)
2909            data['Residue'][resRBsel] = {'RBname':name,'rbXYZ':rbXYZ,'rbTypes':rbTypes,
2910                'atNames':atNames,'rbRef':[0,1,2,False],'rbSeq':[],'SelSeq':[0,0],'useCount':0,'molCent':False}
2911            data['RBIds']['Residue'].append(resRBsel)
2912            print ('Rigid body UNKRB added')
2913        text.close()
2914        UpdateResidueRB()
2915       
2916    def SaveResidueRB():
2917        global resRBsel
2918        pth = G2G.GetExportPath(G2frame)
2919        dlg = wx.FileDialog(G2frame, 'Choose PDB file for Atom XYZ', pth, '', 
2920            'PDB files (*.pdb)|*.pdb',wx.FD_SAVE|wx.FD_OVERWRITE_PROMPT)
2921        try:
2922            if dlg.ShowModal() == wx.ID_OK:
2923                filename = dlg.GetPath()
2924                filename = os.path.splitext(filename)[0]+'.pdb'  # make extension .pdb
2925                File = open(filename,'w')       
2926                rbData =  data['Residue'][resRBsel]
2927                for iat,xyz in enumerate(rbData['rbXYZ']):
2928                    File.write('ATOM %6d  %-4s%3s     1    %8.3f%8.3f%8.3f  1.00  0.00          %2s\n'%(
2929                        iat,rbData['atNames'][iat],rbData['RBname'][:3],xyz[0],xyz[1],xyz[2],rbData['rbTypes'][iat]))
2930                File.close()
2931                print ('Atom XYZ saved to: '+filename)
2932        finally:
2933            dlg.Destroy()
2934       
2935       
2936    def FindNeighbors(Orig,XYZ,atTypes,atNames,AtInfo):
2937        Radii = []
2938        for Atype in atTypes:
2939            Radii.append(AtInfo[Atype][0])
2940        Radii = np.array(Radii)
2941        Neigh = []
2942        Dx = XYZ-XYZ[Orig]
2943        dist = np.sqrt(np.sum(Dx**2,axis=1))
2944        sumR = Radii[Orig]+Radii
2945        IndB = ma.nonzero(ma.masked_greater(dist-0.85*sumR,0.))
2946        for j in IndB[0]:
2947            if j != Orig and atTypes[j] != 'H':
2948                Neigh.append(atNames[j])
2949        return Neigh
2950       
2951    def FindAllNeighbors(XYZ,atTypes,atNames,AtInfo):
2952        NeighDict = {}
2953        for iat,xyz in enumerate(atNames):
2954            NeighDict[atNames[iat]] = FindNeighbors(iat,XYZ,atTypes,atNames,AtInfo)
2955        return NeighDict
2956       
2957    def FindRiding(Orig,Pivot,NeighDict):
2958        riding = [Orig,Pivot]
2959        iAdd = 1
2960        new = True
2961        while new:
2962            newAtms = NeighDict[riding[iAdd]]
2963            for At in newAtms:
2964                new = False
2965                if At not in riding:
2966                    riding.append(At)
2967                    new = True
2968            iAdd += 1
2969            if iAdd < len(riding):
2970                new = True
2971        return riding[2:]
2972                       
2973    def OnDefineTorsSeq(event):
2974        global resRBsel
2975        rbData = data['Residue'][resRBsel]
2976        if not len(rbData):
2977            return
2978        atNames = rbData['atNames']
2979        AtInfo = data['Residue']['AtInfo']
2980        atTypes = rbData['rbTypes']
2981        XYZ = rbData['rbXYZ']
2982        neighDict = FindAllNeighbors(XYZ,atTypes,atNames,AtInfo)
2983        TargList = []           
2984        dlg = wx.SingleChoiceDialog(G2frame,'Select origin atom for torsion sequence','Origin atom',rbData['atNames'])
2985        if dlg.ShowModal() == wx.ID_OK:
2986            Orig = dlg.GetSelection()
2987            TargList = neighDict[atNames[Orig]]
2988        dlg.Destroy()
2989        if not len(TargList):
2990            return
2991        dlg = wx.SingleChoiceDialog(G2frame,'Select pivot atom for torsion sequence','Pivot atom',TargList)
2992        if dlg.ShowModal() == wx.ID_OK:
2993            Piv = atNames.index(TargList[dlg.GetSelection()])
2994            riding = FindRiding(atNames[Orig],atNames[Piv],neighDict)
2995            Riding = []
2996            for atm in riding:
2997                Riding.append(atNames.index(atm))
2998            rbData['rbSeq'].append([Orig,Piv,0.0,Riding])           
2999        dlg.Destroy()
3000        UpdateResidueRB()
3001
3002    def UpdateVectorRB(Scroll=0):
3003        '''Display & edit a selected Vector RB
3004        '''
3005        global resRBsel
3006        def rbNameSizer(rbid,rbData):
3007
3008            def OnRBName(event):
3009                event.Skip()
3010                Obj = event.GetEventObject()
3011                name = Obj.GetValue()
3012                if name == rbData['RBname']: return # no change
3013                namelist = [data['Vector'][key]['RBname'] for key in
3014                        data['Vector'] if 'RBname' in data['Vector'][key]]
3015                name = G2obj.MakeUniqueLabel(name,namelist)
3016                rbData['RBname'] = name
3017                wx.CallAfter(UpdateVectorRB)
3018               
3019            def OnDelRB(event):
3020                Obj = event.GetEventObject()
3021                rbid = Indx[Obj.GetId()]
3022                if rbid in data['Vector']:
3023                    del data['Vector'][rbid]
3024                    data['RBIds']['Vector'].remove(rbid)
3025                    rbData['useCount'] -= 1
3026                wx.CallAfter(UpdateVectorRB)
3027               
3028            def OnPlotRB(event):
3029                G2plt.PlotRigidBody(G2frame,'Vector',AtInfo,rbData,plotDefaults)
3030           
3031            # start of rbNameSizer
3032            nameSizer = wx.BoxSizer(wx.HORIZONTAL)
3033            nameSizer.Add(wx.StaticText(VectorRBDisplay,-1,'Rigid body name: '),0,WACV)
3034            RBname = wx.TextCtrl(VectorRBDisplay,-1,rbData['RBname'],
3035                                     style=wx.TE_PROCESS_ENTER)
3036            RBname.Bind(wx.EVT_LEAVE_WINDOW, OnRBName)
3037            RBname.Bind(wx.EVT_TEXT_ENTER,OnRBName)
3038            RBname.Bind(wx.EVT_KILL_FOCUS,OnRBName)
3039            nameSizer.Add(RBname,0,WACV)
3040            nameSizer.Add((5,0),)
3041            plotRB =  wx.Button(VectorRBDisplay,wx.ID_ANY,'Plot',style=wx.BU_EXACTFIT)
3042            plotRB.Bind(wx.EVT_BUTTON, OnPlotRB)
3043            Indx[plotRB.GetId()] = rbid
3044            nameSizer.Add(plotRB,0,WACV)
3045            nameSizer.Add((5,0),)
3046            if not rbData['useCount']:
3047                delRB = wx.Button(VectorRBDisplay,wx.ID_ANY,"Delete",style=wx.BU_EXACTFIT)
3048                delRB.Bind(wx.EVT_BUTTON, OnDelRB)
3049                Indx[delRB.GetId()] = rbid
3050                nameSizer.Add(delRB,0,WACV)
3051            return nameSizer
3052           
3053        def rbRefAtmSizer(rbid,rbData):
3054           
3055            def OnRefSel(event):
3056                Obj = event.GetEventObject()
3057                iref = Indx[Obj.GetId()]
3058                sel = Obj.GetValue()
3059                rbData['rbRef'][iref] = atNames.index(sel)
3060                FillRefChoice(rbid,rbData)
3061           
3062            refAtmSizer = wx.BoxSizer(wx.HORIZONTAL)
3063            atNames = [name+str(i) for i,name in enumerate(rbData['rbTypes'])]
3064            rbRef = rbData.get('rbRef',[0,1,2,False])
3065            rbData['rbRef'] = rbRef
3066            if rbData['useCount']:
3067                refAtmSizer.Add(wx.StaticText(VectorRBDisplay,-1,
3068                    'Orientation reference atoms A-B-C: %s, %s, %s'%(atNames[rbRef[0]], \
3069                     atNames[rbRef[1]],atNames[rbRef[2]])),0)
3070            else:
3071                refAtmSizer.Add(wx.StaticText(VectorRBDisplay,-1,
3072                    'Orientation reference atoms A-B-C: '),0,WACV)
3073                for i in range(3):
3074                    choices = [atNames[j] for j in refChoice[rbid][i]]
3075                    refSel = wx.ComboBox(VectorRBDisplay,-1,value='',
3076                        choices=choices,style=wx.CB_READONLY|wx.CB_DROPDOWN)
3077                    refSel.SetValue(atNames[rbRef[i]])
3078                    refSel.Bind(wx.EVT_COMBOBOX, OnRefSel)
3079                    Indx[refSel.GetId()] = i
3080                    refAtmSizer.Add(refSel,0,WACV)
3081                refHelpInfo = '''
3082* The "Orientation Reference" control defines the Cartesian
3083axes for rigid bodies with the three atoms, A, B and C.
3084The vector from B to A defines the x-axis and the y axis is placed
3085in the plane defined by B to A and C to A. A,B,C must not be collinear.
3086'''
3087                hlp = G2G.HelpButton(VectorRBDisplay,refHelpInfo,wrap=400)
3088                refAtmSizer.Add(hlp,0,wx.LEFT|wx.RIGHT|wx.ALIGN_CENTER_VERTICAL,2)
3089            return refAtmSizer
3090                       
3091        def rbVectMag(rbid,imag,rbData):
3092           
3093            def OnRBVectorMag(event):
3094                event.Skip()
3095                Obj = event.GetEventObject()
3096                rbid,imag = Indx[Obj.GetId()]
3097                try:
3098                    val = float(Obj.GetValue())
3099                    if val <= 0.:
3100                        raise ValueError
3101                    rbData['VectMag'][imag] = val
3102                except ValueError:
3103                    pass
3104                Obj.SetValue('%8.4f'%(val))
3105                wx.CallAfter(UpdateVectorRB,VectorRB.GetScrollPos(wx.VERTICAL))
3106                G2plt.PlotRigidBody(G2frame,'Vector',AtInfo,data['Vector'][rbid],plotDefaults)
3107               
3108            def OnRBVectorRef(event):
3109                Obj = event.GetEventObject()
3110                rbid,imag = Indx[Obj.GetId()]
3111                rbData['VectRef'][imag] = Obj.GetValue()
3112                       
3113            magSizer = wx.BoxSizer(wx.HORIZONTAL)
3114            magSizer.Add(wx.StaticText(VectorRBDisplay,-1,'Translation magnitude: '),0,WACV)
3115            magValue = wx.TextCtrl(VectorRBDisplay,-1,'%8.4f'%(rbData['VectMag'][imag]))
3116            Indx[magValue.GetId()] = [rbid,imag]
3117            magValue.Bind(wx.EVT_TEXT_ENTER,OnRBVectorMag)
3118            magValue.Bind(wx.EVT_KILL_FOCUS,OnRBVectorMag)
3119            magSizer.Add(magValue,0,WACV)
3120            magSizer.Add((5,0),)
3121            magref = wx.CheckBox(VectorRBDisplay,label=' Refine?') 
3122            magref.SetValue(rbData['VectRef'][imag])
3123            magref.Bind(wx.EVT_CHECKBOX,OnRBVectorRef)
3124            Indx[magref.GetId()] = [rbid,imag]
3125            magSizer.Add(magref,0,WACV)
3126            return magSizer
3127           
3128        def rbVectors(rbid,imag,mag,XYZ,rbData):
3129
3130            def TypeSelect(event):
3131                AtInfo = data['Vector']['AtInfo']
3132                r,c = event.GetRow(),event.GetCol()
3133                if vecGrid.GetColLabelValue(c) == 'Type':
3134                    PE = G2elemGUI.PickElement(G2frame,oneOnly=True)
3135                    if PE.ShowModal() == wx.ID_OK:
3136                        if PE.Elem != 'None':
3137                            El = PE.Elem.strip().lower().capitalize()
3138                            if El not in AtInfo:
3139                                Info = G2elem.GetAtomInfo(El)
3140                                AtInfo[El] = [Info['Drad'],Info['Color']]
3141                            rbData['rbTypes'][r] = El
3142                            vecGrid.SetCellValue(r,c,El)
3143                    PE.Destroy()
3144                wx.CallAfter(UpdateVectorRB,VectorRB.GetScrollPos(wx.VERTICAL))
3145
3146            def ChangeCell(event):
3147                r,c =  event.GetRow(),event.GetCol()
3148                if r >= 0 and (0 <= c < 3):
3149                    try:
3150                        val = float(vecGrid.GetCellValue(r,c))
3151                        rbData['rbVect'][imag][r][c] = val
3152                    except ValueError:
3153                        pass
3154                G2plt.PlotRigidBody(G2frame,'Vector',AtInfo,data['Vector'][rbid],plotDefaults)
3155                wx.CallAfter(UpdateVectorRB,VectorRB.GetScrollPos(wx.VERTICAL))
3156
3157            vecSizer = wx.BoxSizer()
3158            Types = 3*[wg.GRID_VALUE_FLOAT+':10,5',]+[wg.GRID_VALUE_STRING,]+3*[wg.GRID_VALUE_FLOAT+':10,5',]
3159            colLabels = ['Vector x','Vector y','Vector z','Type','Cart x','Cart y','Cart z']
3160            table = []
3161            rowLabels = []
3162            atNames = []
3163            for ivec,xyz in enumerate(rbData['rbVect'][imag]):
3164                table.append(list(xyz)+[rbData['rbTypes'][ivec],]+list(XYZ[ivec]))
3165                rowLabels.append(str(ivec))
3166                atNames.append(rbData['rbTypes'][ivec]+str(ivec))
3167            rbData['atNames'] = atNames
3168            vecTable = G2G.Table(table,rowLabels=rowLabels,colLabels=colLabels,types=Types)
3169            vecGrid = G2G.GSGrid(VectorRBDisplay)
3170            vecGrid.SetTable(vecTable, True)
3171            if 'phoenix' in wx.version():
3172                vecGrid.Bind(wg.EVT_GRID_CELL_CHANGED, ChangeCell)
3173            else:
3174                vecGrid.Bind(wg.EVT_GRID_CELL_CHANGE, ChangeCell)
3175            if not imag:
3176                vecGrid.Bind(wg.EVT_GRID_CELL_LEFT_DCLICK, TypeSelect)
3177            attr = wx.grid.GridCellAttr()
3178            attr.IncRef()
3179            attr.SetEditor(G2G.GridFractionEditor(vecGrid))
3180            for c in range(3):
3181                attr.IncRef()
3182                vecGrid.SetColAttr(c, attr)
3183            for row in range(vecTable.GetNumberRows()):
3184                if imag:
3185                    vecGrid.SetCellStyle(row,3,VERY_LIGHT_GREY,True)                   
3186                for col in [4,5,6]:
3187                    vecGrid.SetCellStyle(row,col,VERY_LIGHT_GREY,True)
3188#            vecGrid.SetScrollRate(0,0)
3189            vecGrid.AutoSizeColumns(False)
3190            vecSizer.Add(vecGrid)
3191            return vecSizer
3192       
3193        def FillRefChoice(rbid,rbData):
3194            choiceIds = [i for i in range(len(rbData['rbTypes']))]
3195           
3196            rbRef = rbData.get('rbRef',[-1,-1,-1,False])
3197            for i in range(3):
3198                if rbRef[i] in choiceIds: choiceIds.remove(rbRef[i])
3199            refChoice[rbid] = [choiceIds[:],choiceIds[:],choiceIds[:]]
3200            for i in range(3):
3201                refChoice[rbid][i].append(rbRef[i])
3202                refChoice[rbid][i].sort()
3203               
3204        def OnRBSelect(event):
3205            global resRBsel
3206            sel = rbSelect.GetSelection()
3207            if sel == 0: return # 1st entry is blank
3208            rbname = rbchoice[sel-1]
3209            resRBsel = RBnames[rbname]
3210            wx.CallLater(100,UpdateVectorRB)
3211           
3212        # beginning of UpdateVectorRB
3213        AtInfo = data['Vector']['AtInfo']
3214        refChoice = {}
3215        RefObjs = []
3216
3217        GS = VectorRBDisplay.GetSizer()
3218        if GS: 
3219            try:        #get around a c++ error in wx 4.0; doing is again seems to be OK
3220                GS.Clear(True)
3221            except:
3222                GS.Clear(True)
3223       
3224        RBnames = {}
3225        for rbid in data['RBIds']['Vector']:
3226            RBnames.update({data['Vector'][rbid]['RBname']:rbid,})
3227        if not RBnames:
3228            return
3229        rbchoice = list(RBnames.keys())
3230        rbchoice.sort()
3231        if GS: 
3232            VectorRBSizer = GS
3233        else:
3234            VectorRBSizer = wx.BoxSizer(wx.VERTICAL)
3235        if resRBsel not in data['RBIds']['Vector']:
3236            resRBsel = RBnames[rbchoice[0]]
3237        if len(RBnames) > 1:
3238            selSizer = wx.BoxSizer(wx.HORIZONTAL)
3239            selSizer.Add(wx.StaticText(VectorRBDisplay,label=' Select rigid body to view:'),0)
3240            rbSelect = wx.ComboBox(VectorRBDisplay,choices=['']+rbchoice)
3241            name = data['Vector'][resRBsel]['RBname']
3242            rbSelect.SetSelection(1+rbchoice.index(name))
3243            rbSelect.Bind(wx.EVT_COMBOBOX,OnRBSelect)
3244            selSizer.Add(rbSelect,0)
3245            VectorRBSizer.Add(selSizer,0)
3246        rbData = data['Vector'][resRBsel]
3247        if 'DELETED' in str(G2frame.GetStatusBar()):   #seems to be no other way to do this (wx bug)
3248            if GSASIIpath.GetConfigValue('debug'):
3249                print ('DBG_wx error: Rigid Body/Status not cleanly deleted after Refine')
3250            return
3251        SetStatusLine(' You may use e.g. "c60" or "s60" for a vector entry')
3252        FillRefChoice(resRBsel,rbData)
3253        VectorRBSizer.Add(rbNameSizer(resRBsel,rbData),0)
3254        VectorRBSizer.Add(rbRefAtmSizer(resRBsel,rbData),0)
3255        XYZ = np.array([[0.,0.,0.] for Ty in rbData['rbTypes']])
3256        for imag,mag in enumerate(rbData['VectMag']):
3257            XYZ += mag*rbData['rbVect'][imag]
3258            VectorRBSizer.Add(rbVectMag(rbid,imag,rbData),0)
3259            VectorRBSizer.Add(rbVectors(rbid,imag,mag,XYZ,rbData),0)
3260        VectorRBSizer.Add((5,5),0)
3261        data['Vector'][rbid]['rbXYZ'] = XYZ       
3262       
3263        VectorRBSizer.Add((5,25),)
3264        VectorRBSizer.Layout()   
3265        VectorRBDisplay.SetSizer(VectorRBSizer,True)
3266        VectorRBDisplay.SetAutoLayout(True)
3267        Size = VectorRBSizer.GetMinSize()
3268        VectorRBDisplay.SetSize(Size)
3269       
3270        Size[0] += 40
3271        Size[1] = max(Size[1],450) + 20
3272        VectorRB.SetSize(Size)
3273        VectorRB.SetScrollbars(10,10,Size[0]/10-4,Size[1]/10-1)
3274        G2frame.dataWindow.SendSizeEvent()
3275       
3276        VectorRBDisplay.Show()
3277       
3278       
3279    def UpdateResidueRB():
3280        '''Draw the contents of the Residue Rigid Body tab for Rigid Bodies tree entry
3281        '''
3282        global resRBsel
3283        def rbNameSizer(rbid,rbData):
3284           
3285            def OnRBName(event):
3286                event.Skip()
3287                Obj = event.GetEventObject()
3288                name = Obj.GetValue()
3289                if name == rbData['RBname']: return # no change
3290                namelist = [data['Residue'][key]['RBname'] for key in
3291                        data['Residue'] if 'RBname' in data['Residue'][key]]
3292                name = G2obj.MakeUniqueLabel(name,namelist)
3293                rbData['RBname'] = name
3294                wx.CallAfter(UpdateResidueRB)
3295               
3296            def OnDelRB(event):
3297                Obj = event.GetEventObject()
3298                rbid = Indx[Obj.GetId()]
3299                if rbid in data['Residue']: 
3300                    del data['Residue'][rbid]
3301                    data['RBIds']['Residue'].remove(rbid)
3302                wx.CallAfter(UpdateResidueRB)
3303               
3304            def OnStripH(event):
3305                Obj = event.GetEventObject()
3306                rbid = Indx[Obj.GetId()]
3307                if rbid in data['Residue']:
3308                    newNames = []
3309                    newTypes = []
3310                    newXYZ = []
3311                    for i,atype in enumerate(rbData['rbTypes']):
3312                        if atype != 'H':
3313                            newNames.append(rbData['atNames'][i])
3314                            newTypes.append(rbData['rbTypes'][i])
3315                            newXYZ.append(rbData['rbXYZ'][i])
3316                    rbData['atNames'] = newNames
3317                    rbData['rbTypes'] = newTypes
3318                    rbData['rbXYZ'] = newXYZ
3319                G2plt.PlotRigidBody(G2frame,'Residue',AtInfo,rbData,plotDefaults)
3320                wx.CallAfter(UpdateResidueRB)
3321
3322            def OnPlotRB(event):
3323                G2plt.PlotRigidBody(G2frame,'Residue',AtInfo,rbData,plotDefaults)
3324               
3325            # start of rbNameSizer
3326            nameSizer = wx.BoxSizer(wx.HORIZONTAL)
3327            nameSizer.Add(wx.StaticText(ResidueRBDisplay,-1,'Residue name: '),0,WACV)
3328            RBname = wx.TextCtrl(ResidueRBDisplay,-1,rbData['RBname'],
3329                                     style=wx.TE_PROCESS_ENTER)
3330            RBname.Bind(wx.EVT_LEAVE_WINDOW, OnRBName)
3331            RBname.Bind(wx.EVT_TEXT_ENTER,OnRBName)
3332            RBname.Bind(wx.EVT_KILL_FOCUS,OnRBName)
3333            nameSizer.Add(RBname,0,WACV)
3334            nameSizer.Add((5,0),)
3335            plotRB =  wx.Button(ResidueRBDisplay,wx.ID_ANY,'Plot',style=wx.BU_EXACTFIT)
3336            plotRB.Bind(wx.EVT_BUTTON, OnPlotRB)
3337            Indx[plotRB.GetId()] = rbid
3338            nameSizer.Add(plotRB,0,WACV)
3339            nameSizer.Add((5,0),)
3340            if not rbData['useCount']:
3341                delRB = wx.Button(ResidueRBDisplay,wx.ID_ANY,"Delete",style=wx.BU_EXACTFIT)
3342                delRB.Bind(wx.EVT_BUTTON, OnDelRB)
3343                Indx[delRB.GetId()] = rbid
3344                nameSizer.Add(delRB,0,WACV)
3345                if 'H'  in rbData['rbTypes']:
3346                    stripH = wx.Button(ResidueRBDisplay,wx.ID_ANY,"Strip H-atoms",style=wx.BU_EXACTFIT)
3347                    stripH.Bind(wx.EVT_BUTTON, OnStripH)
3348                    Indx[stripH.GetId()] = rbid
3349                    nameSizer.Add(stripH,0,WACV)
3350            nameSizer.Add(wx.StaticText(ResidueRBDisplay,-1,'  body type #'+
3351                                        str(data['RBIds']['Residue'].index(rbid))),0,WACV)
3352            return nameSizer
3353           
3354        def rbResidues(rbid,rbData):
3355           
3356            def TypeSelect(event):
3357                AtInfo = data['Residue']['AtInfo']
3358                r,c = event.GetRow(),event.GetCol()
3359                if resGrid.GetColLabelValue(c) == 'Type':
3360                    PE = G2elemGUI.PickElement(G2frame,oneOnly=True)
3361                    if PE.ShowModal() == wx.ID_OK:
3362                        if PE.Elem != 'None':
3363                            El = PE.Elem.strip().lower().capitalize()
3364                            if El not in AtInfo:
3365                                Info = G2elem.GetAtomInfo(El)
3366                                AtInfo[El] = [Info['Drad'],Info['Color']]
3367                            rbData['rbTypes'][r] = El
3368                            resGrid.SetCellValue(r,c,El)
3369                    PE.Destroy()
3370
3371            def ChangeCell(event):
3372                r,c =  event.GetRow(),event.GetCol()
3373                if c == 0:
3374                    rbData['atNames'][r] = resGrid.GetCellValue(r,c)
3375                if r >= 0 and (2 <= c <= 4):
3376                    try:
3377                        val = float(resGrid.GetCellValue(r,c))
3378                        rbData['rbXYZ'][r][c-2] = val
3379                    except ValueError:
3380                        pass
3381                G2plt.PlotRigidBody(G2frame,'Residue',AtInfo,rbData,plotDefaults)
3382                       
3383            def OnRefSel(event):
3384                Obj = event.GetEventObject()
3385                iref,res,jref = Indx[Obj.GetId()]
3386                sel = Obj.GetValue()
3387                ind = atNames.index(sel)
3388                if rbData['rbTypes'][ind] == 'H':
3389                    G2G.G2MessageBox(G2frame,'You should not select an H-atom for rigid body orientation')
3390                rbData['rbRef'][iref] = ind
3391                FillRefChoice(rbid,rbData)
3392                for i,ref in enumerate(RefObjs[jref]):
3393                    ref.SetItems([atNames[j] for j in refChoice[rbid][i]])
3394                    ref.SetValue(atNames[rbData['rbRef'][i]])                   
3395                rbXYZ = rbData['rbXYZ']
3396                if not iref:     #origin change
3397                    rbXYZ -= rbXYZ[ind]
3398                Xxyz = rbXYZ[rbData['rbRef'][1]]
3399                X = Xxyz/np.sqrt(np.sum(Xxyz**2))
3400                Yxyz = rbXYZ[rbData['rbRef'][2]]
3401                Y = Yxyz/np.sqrt(np.sum(Yxyz**2))
3402                Mat = G2mth.getRBTransMat(X,Y)
3403                rbXYZ = np.inner(Mat,rbXYZ).T
3404                rbData['rbXYZ'] = rbXYZ
3405                rbData['molCent'] = False
3406                res.ClearSelection()
3407                resTable = res.GetTable()
3408                for r in range(res.GetNumberRows()):
3409                    row = resTable.GetRowValues(r)
3410                    row[2:4] = rbXYZ[r]
3411                    resTable.SetRowValues(r,row)
3412                res.ForceRefresh()
3413                G2plt.PlotRigidBody(G2frame,'Residue',AtInfo,rbData,plotDefaults)
3414               
3415            def OnMolCent(event):
3416                rbData['molCent'] = not rbData['molCent']
3417                if rbData['molCent']:
3418                    Obj = event.GetEventObject()
3419                    res = Indx[Obj.GetId()]
3420                    rbXYZ = rbData['rbXYZ']
3421                    rbCent = np.array([np.sum(rbXYZ[:,0]),np.sum(rbXYZ[:,1]),np.sum(rbXYZ[:,2])])/rbXYZ.shape[0]
3422                    rbXYZ -= rbCent
3423                    rbData['rbXYZ'] = rbXYZ
3424                    res.ClearSelection()
3425                    resTable = res.GetTable()
3426                    for r in range(res.GetNumberRows()):
3427                        row = resTable.GetRowValues(r)
3428                        row[2:4] = rbXYZ[r]
3429                        resTable.SetRowValues(r,row)
3430                    res.ForceRefresh()
3431                    G2plt.PlotRigidBody(G2frame,'Residue',AtInfo,rbData,plotDefaults)
3432                   
3433            def OnCycleXYZ(event):
3434                Obj = event.GetEventObject()
3435                res = Indx[Obj.GetId()]
3436                rbXYZ = rbData['rbXYZ']
3437                resTable = res.GetTable()
3438                for r in range(res.GetNumberRows()):
3439                    rbXYZ[r] = np.roll(rbXYZ[r],1)
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            Types = 2*[wg.GRID_VALUE_STRING,]+3*[wg.GRID_VALUE_FLOAT+':10,5',]
3447            colLabels = ['Name','Type','Cart x','Cart y','Cart z']
3448            table = []
3449            rowLabels = []
3450            for ivec,xyz in enumerate(rbData['rbXYZ']):
3451                table.append([rbData['atNames'][ivec],]+[rbData['rbTypes'][ivec],]+list(xyz))
3452                rowLabels.append(str(ivec))
3453            vecTable = G2G.Table(table,rowLabels=rowLabels,colLabels=colLabels,types=Types)
3454            resGrid = G2G.GSGrid(ResidueRBDisplay)
3455            Indx[resGrid.GetId()] = rbid
3456            resList.append(resGrid)
3457            resGrid.SetTable(vecTable, True)
3458            if 'phoenix' in wx.version():
3459                resGrid.Bind(wg.EVT_GRID_CELL_CHANGED, ChangeCell)
3460            else:
3461                resGrid.Bind(wg.EVT_GRID_CELL_CHANGE, ChangeCell)
3462            resGrid.Bind(wg.EVT_GRID_CELL_LEFT_DCLICK, TypeSelect)
3463            for c in range(2,5):
3464                attr = wx.grid.GridCellAttr()
3465                attr.IncRef()
3466                attr.SetEditor(G2G.GridFractionEditor(resGrid))
3467                resGrid.SetColAttr(c, attr)
3468            for row in range(vecTable.GetNumberRows()):
3469                resGrid.SetReadOnly(row,1,True)
3470            resGrid.AutoSizeColumns(False)
3471            vecSizer = wx.BoxSizer()
3472            vecSizer.Add(resGrid)
3473           
3474            refAtmSizer = wx.BoxSizer(wx.HORIZONTAL)
3475            atNames = rbData['atNames']
3476            rbRef = rbData['rbRef']
3477            if rbData['rbRef'][3] or rbData['useCount']:
3478                refAtmSizer.Add(wx.StaticText(ResidueRBDisplay,-1,
3479                    'Orientation reference non-H atoms A-B-C: %s, %s, %s'%(atNames[rbRef[0]], \
3480                     atNames[rbRef[1]],atNames[rbRef[2]])),0)
3481            else:
3482                refAtmSizer.Add(wx.StaticText(ResidueRBDisplay,-1,
3483                    'Orientation reference non-H atoms A-B-C: '),0,WACV)
3484                refObj = [0,0,0]
3485                for i in range(3):
3486                    choices = [atNames[j] for j in refChoice[rbid][i]]
3487                    refSel = wx.ComboBox(ResidueRBDisplay,-1,value='',
3488                        choices=choices,style=wx.CB_READONLY|wx.CB_DROPDOWN)
3489                    refSel.SetValue(atNames[rbRef[i]])
3490                    refSel.Bind(wx.EVT_COMBOBOX, OnRefSel)
3491                    Indx[refSel.GetId()] = [i,resGrid,len(RefObjs)]
3492                    refObj[i] = refSel
3493                    refAtmSizer.Add(refSel,0,WACV)
3494                RefObjs.append(refObj)
3495                cycleXYZ = wx.Button(ResidueRBDisplay,label=' Cycle XYZ?')
3496                cycleXYZ.Bind(wx.EVT_BUTTON,OnCycleXYZ)
3497                Indx[cycleXYZ.GetId()] = resGrid
3498                refAtmSizer.Add(cycleXYZ,0,WACV)
3499                if 'molCent' not in rbData: rbData['molCent'] = False           #patch
3500                molcent = wx.Button(ResidueRBDisplay,label=' Center RB?')
3501                molcent.Bind(wx.EVT_BUTTON,OnMolCent)
3502                Indx[molcent.GetId()] = resGrid
3503                refAtmSizer.Add(molcent,0,WACV)
3504                refHelpInfo = '''
3505* The "Orientation Reference" control defines the Cartesian
3506axes for rigid bodies with the three atoms, A, B and C.
3507The vector from B to A defines the x-axis and the y axis is placed
3508in the plane defined by B to A and C to A. A,B,C must not be collinear.
3509 
3510%%* The origin is at A unless the "Center RB?" button is pressed.
3511
3512%%* The 'Cycle XYZ' button will permute the rigid body XYZ coordinates so
3513XYZ --> ZXY. Repeat if needed.
3514
3515%%* The "Center RB?" button will shift the origin of the
3516rigid body to be the midpoint of all atoms in the body (not mass weighted).
3517'''
3518                hlp = G2G.HelpButton(ResidueRBDisplay,refHelpInfo,wrap=400)
3519                refAtmSizer.Add(hlp,0,wx.LEFT|wx.RIGHT|wx.ALIGN_CENTER_VERTICAL,2)
3520           
3521            mainSizer = wx.BoxSizer(wx.VERTICAL)
3522            mainSizer.Add(refAtmSizer)
3523            mainSizer.Add(vecSizer)
3524            return mainSizer
3525           
3526        def Add2SeqSizer(seqSizer,angSlide,rbid,iSeq,Seq,atNames):
3527           
3528            def ChangeAngle(event):
3529                event.Skip()
3530                Obj = event.GetEventObject()
3531                rbid,Seq = Indx[Obj.GetId()][:2]
3532                val = Seq[2]
3533                try:
3534                    val = float(Obj.GetValue())
3535                    Seq[2] = val
3536                except ValueError:
3537                    pass
3538                Obj.SetValue('%8.2f'%(val))
3539                G2plt.PlotRigidBody(G2frame,'Residue',AtInfo,data['Residue'][rbid],plotDefaults)
3540               
3541            def OnRadBtn(event):
3542                Obj = event.GetEventObject()
3543                Seq,iSeq,angId = Indx[Obj.GetId()]
3544                data['Residue'][rbid]['SelSeq'] = [iSeq,angId]
3545                angSlide.SetValue(int(100*Seq[2]))
3546               
3547            def OnDelBtn(event):
3548                Obj = event.GetEventObject()
3549                rbid,Seq = Indx[Obj.GetId()]
3550                data['Residue'][rbid]['rbSeq'].remove(Seq)       
3551                wx.CallAfter(UpdateResidueRB)
3552           
3553            iBeg,iFin,angle,iMove = Seq
3554            ang = wx.TextCtrl(ResidueRBDisplay,wx.ID_ANY,
3555                    '%8.2f'%(angle),size=(70,-1),style=wx.TE_PROCESS_ENTER)
3556            if not iSeq:
3557                radBt = wx.RadioButton(ResidueRBDisplay,wx.ID_ANY,
3558                                           '',style=wx.RB_GROUP)
3559                data['Residue'][rbid]['SelSeq'] = [iSeq,ang.GetId()]
3560                radBt.SetValue(True)
3561            else:
3562                radBt = wx.RadioButton(ResidueRBDisplay,wx.ID_ANY,'')
3563            radBt.Bind(wx.EVT_RADIOBUTTON,OnRadBtn)                   
3564            seqSizer.Add(radBt)
3565            delBt =  wx.Button(ResidueRBDisplay,wx.ID_ANY,'Del',
3566                                style=wx.BU_EXACTFIT)
3567            delBt.Bind(wx.EVT_BUTTON,OnDelBtn)
3568            seqSizer.Add(delBt)
3569            bond = wx.StaticText(ResidueRBDisplay,wx.ID_ANY,
3570                        '%s %s'%(atNames[iBeg],atNames[iFin]),size=(50,20))
3571            seqSizer.Add(bond,0,WACV)
3572            Indx[radBt.GetId()] = [Seq,iSeq,ang.GetId()]
3573            Indx[delBt.GetId()] = [rbid,Seq]
3574            Indx[ang.GetId()] = [rbid,Seq,ang]
3575            ang.Bind(wx.EVT_TEXT_ENTER,ChangeAngle)
3576            ang.Bind(wx.EVT_KILL_FOCUS,ChangeAngle)
3577            seqSizer.Add(ang,0,WACV)
3578            atms = ''
3579            for i in iMove:   
3580                atms += ' %s,'%(atNames[i])
3581            moves = wx.StaticText(ResidueRBDisplay,wx.ID_ANY,
3582                            atms[:-1],size=(200,20))
3583            seqSizer.Add(moves,1,wx.EXPAND|wx.RIGHT)
3584            return seqSizer
3585           
3586        def SlideSizer():
3587           
3588            def OnSlider(event):
3589                Obj = event.GetEventObject()
3590                rbData = Indx[Obj.GetId()]
3591                iSeq,angId = rbData['SelSeq']
3592                val = float(Obj.GetValue())/100.
3593                rbData['rbSeq'][iSeq][2] = val
3594                Indx[angId][2].SetValue('%8.2f'%(val))
3595                G2plt.PlotRigidBody(G2frame,'Residue',AtInfo,rbData,plotDefaults)
3596           
3597            slideSizer = wx.BoxSizer(wx.HORIZONTAL)
3598            slideSizer.Add(wx.StaticText(ResidueRBDisplay,-1,'Selected torsion angle:'),0)
3599            iSeq,angId = rbData['SelSeq']
3600            angSlide = wx.Slider(ResidueRBDisplay,-1,
3601                int(100*rbData['rbSeq'][iSeq][2]),0,36000,size=(200,20),
3602                style=wx.SL_HORIZONTAL)
3603            angSlide.Bind(wx.EVT_SLIDER, OnSlider)
3604            Indx[angSlide.GetId()] = rbData
3605            slideSizer.Add(angSlide,0)           
3606            return slideSizer,angSlide
3607           
3608        def FillRefChoice(rbid,rbData):
3609            choiceIds = [i for i in range(len(rbData['atNames']))]
3610            for seq in rbData['rbSeq']:
3611                for i in seq[3]:
3612                    try:
3613                        choiceIds.remove(i)
3614                    except ValueError:
3615                        pass
3616            rbRef = rbData['rbRef']
3617            for i in range(3):
3618                try:
3619                    choiceIds.remove(rbRef[i])
3620                except ValueError:
3621                    pass
3622            refChoice[rbid] = [choiceIds[:],choiceIds[:],choiceIds[:]]
3623            for i in range(3):
3624                refChoice[rbid][i].append(rbRef[i])
3625                refChoice[rbid][i].sort()
3626               
3627        def OnRBSelect(event):
3628            global resRBsel
3629            sel = rbSelect.GetSelection()
3630            if sel == 0: return # 1st entry is blank
3631            rbname = rbchoice[sel-1]
3632            resRBsel = RBnames[rbname]
3633            G2plt.PlotRigidBody(G2frame,'Residue',AtInfo,rbData,plotDefaults)
3634            wx.CallLater(100,UpdateResidueRB)
3635           
3636        #----- beginning of UpdateResidueRB -----------------------------------------------
3637        AtInfo = data['Residue']['AtInfo']
3638        refChoice = {}
3639        RefObjs = []
3640
3641        GS = ResidueRBDisplay.GetSizer()
3642        if GS: 
3643            try:        #get around a c++ error in wx 4.0; doing is again seems to be OK
3644                GS.Clear(True)
3645            except:
3646                GS.Clear(True)
3647       
3648        RBnames = {}
3649        for rbid in data['RBIds']['Residue']:
3650            RBnames.update({data['Residue'][rbid]['RBname']:rbid,})
3651        if not RBnames:
3652            return
3653        rbchoice = list(RBnames.keys())
3654        rbchoice.sort()
3655        if GS: 
3656            ResidueRBSizer = GS
3657        else:
3658            ResidueRBSizer = wx.BoxSizer(wx.VERTICAL)
3659        if resRBsel not in data['RBIds']['Residue']:
3660            resRBsel = RBnames[rbchoice[0]]
3661        if len(RBnames) > 1:
3662            selSizer = wx.BoxSizer(wx.HORIZONTAL)
3663            selSizer.Add(wx.StaticText(ResidueRBDisplay,
3664                                label=' Select residue to view:'),0)
3665            rbSelect = wx.ComboBox(ResidueRBDisplay,choices=['']+rbchoice)
3666            name = data['Residue'][resRBsel]['RBname']
3667            rbSelect.SetSelection(1+rbchoice.index(name))
3668            rbSelect.Bind(wx.EVT_COMBOBOX,OnRBSelect)
3669            selSizer.Add(rbSelect,0)
3670            ResidueRBSizer.Add(selSizer,0)
3671        rbData = data['Residue'][resRBsel]
3672        FillRefChoice(resRBsel,rbData)
3673        ResidueRBSizer.Add(rbNameSizer(resRBsel,rbData),0)
3674        ResidueRBSizer.Add(rbResidues(resRBsel,rbData),0)
3675        if len(rbData['rbSeq']):
3676            ResidueRBSizer.Add((-1,15),0)
3677            slideSizer,angSlide = SlideSizer()
3678            seqSizer = wx.FlexGridSizer(0,5,4,8)
3679            for lbl in 'Sel','','Bond','Angle','Riding atoms':
3680                seqSizer.Add(wx.StaticText(ResidueRBDisplay,wx.ID_ANY,lbl))
3681            ResidueRBSizer.Add(seqSizer)
3682#            for iSeq,Seq in enumerate(rbData['rbSeq']):
3683#                ResidueRBSizer.Add(SeqSizer(angSlide,resRBsel,iSeq,Seq,rbData['atNames']))
3684            for iSeq,Seq in enumerate(rbData['rbSeq']):
3685                Add2SeqSizer(seqSizer,angSlide,resRBsel,iSeq,Seq,rbData['atNames'])
3686            ResidueRBSizer.Add(slideSizer,)
3687
3688        ResidueRBSizer.Add((5,25),)
3689        ResidueRBSizer.Layout()   
3690        ResidueRBDisplay.SetSizer(ResidueRBSizer,True)
3691        ResidueRBDisplay.SetAutoLayout(True)
3692        Size = ResidueRBSizer.GetMinSize()
3693        ResidueRBDisplay.SetSize(Size)
3694       
3695        Size[0] += 40
3696        Size[1] = max(Size[1],450) + 20
3697        ResidueRB.SetSize(Size)
3698        ResidueRB.SetScrollbars(10,10,Size[0]/10-4,Size[1]/10-1)
3699        G2frame.dataWindow.SendSizeEvent()
3700       
3701        ResidueRBDisplay.Show()
3702       
3703    def SetStatusLine(text):
3704        G2frame.GetStatusBar().SetStatusText(text,1)                                     
3705
3706    #### UpdateRigidBodies starts here =========
3707    global resList,resRBsel           
3708    if not data.get('RBIds') or not data:
3709        data.update({'Vector':{'AtInfo':{}},'Residue':{'AtInfo':{}},
3710            'RBIds':{'Vector':[],'Residue':[]}})       #empty/bad dict - fill it
3711    Indx = {}
3712    resList = []
3713    plotDefaults = {'oldxy':[0.,0.],'Quaternion':[0.,0.,0.,1.],'cameraPos':30.,'viewDir':[0,0,1],}
3714    G2frame.rbBook = G2G.GSNoteBook(parent=G2frame.dataWindow)
3715    G2frame.dataWindow.GetSizer().Add(G2frame.rbBook,1,wx.ALL|wx.EXPAND)
3716    VectorRB = wx.ScrolledWindow(G2frame.rbBook)
3717    VectorRBDisplay = wx.Panel(VectorRB)
3718    G2frame.rbBook.AddPage(VectorRB,'Vector rigid bodies')
3719    ResidueRB = wx.ScrolledWindow(G2frame.rbBook)
3720    ResidueRBDisplay = wx.Panel(ResidueRB)
3721    G2frame.rbBook.AddPage(ResidueRB,'Residue rigid bodies')
3722    # vector RBs are not too common, so select Residue as the default when one is present
3723    if len(data['RBIds']['Residue']) > 0 and len(data['RBIds']['Vector']) == 0:
3724        G2frame.rbBook.ChangeSelection(1)
3725        OnPageChanged(None)
3726    G2gd.SetDataMenuBar(G2frame,G2frame.dataWindow.RigidBodyMenu)
3727    SetStatusLine('')
3728    UpdateVectorRB()
3729    G2frame.rbBook.Bind(wx.aui.EVT_AUINOTEBOOK_PAGE_CHANGED, OnPageChanged)
3730    wx.CallAfter(OnPageChanged,None)
3731
3732def ShowIsoDistortCalc(G2frame,phase=None):
3733    '''Compute the ISODISTORT mode values from the current coordinates.
3734    Called in response to the (Phase/Atoms tab) AtomCompute or
3735    Constraints/Edit Constr. "Show ISODISTORT modes" menu item, which
3736    should be enabled only when Phase['ISODISTORT'] is defined.
3737    '''
3738    def _onClose(event):
3739        dlg.EndModal(wx.ID_CANCEL)
3740    def fmtHelp(item,fullname):
3741        helptext = "A new variable"
3742        if item[-3]:
3743            helptext += " named "+str(item[-3])
3744        helptext += " is a linear combination of the following parameters:\n"
3745        first = True
3746        for term in item[:-3]:
3747            line = ''
3748            var = str(term[1])
3749            m = term[0]
3750            if first:
3751                first = False
3752                line += ' = '
3753            else:
3754                if m >= 0:
3755                    line += ' + '
3756                else:
3757                    line += ' - '
3758                m = abs(m)
3759            line += '%.3f*%s '%(m,var)
3760            varMean = G2obj.fmtVarDescr(var)
3761            helptext += "\n" + line + " ("+ varMean + ")"
3762        helptext += '\n\nISODISTORT full name: '+str(fullname)
3763        return helptext
3764
3765    Histograms,Phases = G2frame.GetUsedHistogramsAndPhasesfromTree() # init for constraint
3766    isophases = [p for p in Phases if 'ISODISTORT' in Phases[p]]
3767   
3768    if not isophases:
3769        G2frame.ErrorDialog('no ISODISTORT phases',
3770                            'Unexpected error: no ISODISTORT phases')
3771        return
3772    if phase and phase not in isophases:
3773        G2frame.ErrorDialog('not ISODISTORT phase',
3774                            'Unexpected error: selected phase is not an ISODISTORT phase')
3775        print('Unexpected error: selected phase is not an ISODISTORT phase',
3776                  phase,isophases)
3777    elif not phase and len(isophases) == 1:
3778        phase = isophases[0]
3779    elif not phase:
3780        dlg = wx.SingleChoiceDialog(G2frame,'Select phase from ISODISTORT phases',
3781                                        'Select Phase',isophases)
3782        if dlg.ShowModal() == wx.ID_OK:
3783            sel = dlg.GetSelection()
3784            phase = isophases[sel]
3785        else:
3786            return
3787    # if len(data.get('Histograms',[])) == 0:
3788    #     G2frame.ErrorDialog(
3789    #         'No data',
3790    #         'Sorry, this computation requires that a histogram first be added to the phase'
3791    #         )
3792    #     return
3793   
3794    covdata = G2frame.GPXtree.GetItemPyData(
3795        G2gd.GetGPXtreeItemId(G2frame,G2frame.root,'Covariance'))
3796    # make a lookup table for named NewVar Phase constraints
3797    sub = G2gd.GetGPXtreeItemId(G2frame,G2frame.root,'Constraints') 
3798    Constraints = G2frame.GPXtree.GetItemPyData(sub)
3799    constDict = {}
3800    for c in Constraints['Phase']:
3801        if c[-1] != 'f' or not c[-3]: continue
3802        constDict[str(c[-3])] = c
3803
3804    parmDict,varyList = G2frame.MakeLSParmDict()
3805
3806    dlg = wx.Dialog(G2frame,wx.ID_ANY,'ISODISTORT mode values',#size=(630,400),
3807                       style=wx.DEFAULT_DIALOG_STYLE|wx.RESIZE_BORDER)
3808    mainSizer = wx.BoxSizer(wx.VERTICAL)
3809    if 'ISODISTORT' not in Phases[phase]:
3810        G2frame.ErrorDialog('not ISODISTORT phase',
3811                            'Unexpected error: selected phase is not an ISODISTORT phase')
3812        return
3813    data = Phases[phase]
3814    ISO = data['ISODISTORT']
3815    mainSizer.Add(wx.StaticText(dlg,wx.ID_ANY,
3816                                'ISODISTORT mode computation for cordinates in phase '+
3817                                str(data['General'].get('Name'))))
3818    aSizer = wx.BoxSizer(wx.HORIZONTAL)
3819    panel1 = wxscroll.ScrolledPanel(
3820        dlg, wx.ID_ANY,#size=(100,200),
3821        style = wx.TAB_TRAVERSAL|wx.SUNKEN_BORDER)
3822    subSizer1 = wx.FlexGridSizer(cols=2,hgap=5,vgap=2)
3823    panel2 = wxscroll.ScrolledPanel(
3824        dlg, wx.ID_ANY,#size=(100,200),
3825        style = wx.TAB_TRAVERSAL|wx.SUNKEN_BORDER)
3826    subSizer2 = wx.FlexGridSizer(cols=3,hgap=5,vgap=2)
3827    subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,'Parameter name  '))
3828    subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,' value'),0,wx.ALIGN_RIGHT)
3829    subSizer2.Add((-1,-1))
3830    subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,'Mode name  '))
3831    subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,' value'),0,wx.ALIGN_RIGHT)
3832    # ISODISTORT displacive modes
3833    if 'G2VarList' in ISO:
3834        deltaList = []
3835        for gv,Ilbl in zip(ISO['G2VarList'],ISO['IsoVarList']):
3836            dvar = gv.varname()
3837            var = dvar.replace('::dA','::A')
3838            albl = Ilbl[:Ilbl.rfind('_')]
3839            v = Ilbl[Ilbl.rfind('_')+1:]
3840            pval = ISO['ParentStructure'][albl][['dx','dy','dz'].index(v)]
3841            if var in parmDict:
3842                cval = parmDict[var][0]
3843            else:
3844                dlg.EndModal(wx.ID_CANCEL)
3845                G2frame.ErrorDialog('Atom not found',"No value found for parameter "+str(var))
3846                return
3847            deltaList.append(cval-pval)
3848        modeVals = np.inner(ISO['Var2ModeMatrix'],deltaList)
3849        for lbl,xyz,var,val,norm,G2mode in zip(
3850                ISO['IsoVarList'],deltaList,
3851                ISO['IsoModeList'],modeVals,ISO['NormList'],ISO['G2ModeList'] ):
3852            #GSASIIpath.IPyBreak()
3853            if str(G2mode) in constDict:
3854                ch = G2G.HelpButton(panel2,fmtHelp(constDict[str(G2mode)],var))
3855                subSizer2.Add(ch,0,wx.LEFT|wx.RIGHT|WACV|wx.ALIGN_CENTER,1)
3856            else:
3857                subSizer2.Add((-1,-1))
3858            subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,str(lbl)))
3859            try:
3860                value = G2py3.FormatSigFigs(xyz)
3861            except TypeError:
3862                value = str(xyz)           
3863            subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,value),0,wx.ALIGN_RIGHT)
3864            subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,str(var)))
3865            try:
3866                value = G2py3.FormatSigFigs(val/norm)
3867                if 'varyList' in covdata:
3868                    if str(G2mode) in covdata['varyList']:
3869                        sig = covdata['sig'][covdata['varyList'].index(str(G2mode))]
3870                        value = G2mth.ValEsd(val/norm,sig/norm)
3871            except TypeError:
3872                value = '?'
3873            subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,value),0,wx.ALIGN_RIGHT)
3874            #GSASIIpath.IPyBreak()
3875    # ISODISTORT occupancy modes
3876    if 'G2OccVarList' in ISO:
3877        deltaList = []
3878        for gv,Ilbl in zip(ISO['G2OccVarList'],ISO['OccVarList']):
3879            var = gv.varname()
3880            albl = Ilbl[:Ilbl.rfind('_')]
3881            pval = ISO['BaseOcc'][albl]
3882            if var in parmDict:
3883                cval = parmDict[var][0]
3884            else:
3885                dlg.EndModal(wx.ID_CANCEL)
3886                G2frame.ErrorDialog('Atom not found',"No value found for parameter "+str(var))
3887                return
3888            deltaList.append(cval-pval)
3889        modeVals = np.inner(ISO['Var2OccMatrix'],deltaList)
3890        for lbl,delocc,var,val,norm,G2mode in zip(
3891                ISO['OccVarList'],deltaList,
3892                ISO['OccModeList'],modeVals,ISO['OccNormList'],ISO['G2OccModeList']):
3893            if str(G2mode) in constDict:
3894                ch = G2G.HelpButton(panel2,fmtHelp(constDict[str(G2mode)],var))
3895                subSizer2.Add(ch,0,wx.LEFT|wx.RIGHT|WACV|wx.ALIGN_CENTER,1)
3896            else:
3897                subSizer2.Add((-1,-1))
3898            subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,str(lbl)))
3899            try:
3900                value = G2py3.FormatSigFigs(delocc)
3901            except TypeError:
3902                value = str(delocc)
3903            subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,value),0,wx.ALIGN_RIGHT)
3904            #subSizer.Add((10,-1))
3905            subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,str(var)))
3906            try:
3907                value = G2py3.FormatSigFigs(val/norm)
3908                if 'varyList' in covdata:
3909                    if str(G2mode) in covdata['varyList']:
3910                        sig = covdata['sig'][covdata['varyList'].index(str(G2mode))]
3911                        value = G2mth.ValEsd(val/norm,sig/norm)
3912            except TypeError:
3913                value = '?'
3914            subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,value),0,wx.ALIGN_RIGHT)
3915
3916    # finish up ScrolledPanel
3917    panel1.SetSizer(subSizer1)
3918    panel2.SetSizer(subSizer2)
3919    panel1.SetAutoLayout(1)
3920    panel1.SetupScrolling()
3921    panel2.SetAutoLayout(1)
3922    panel2.SetupScrolling()
3923    # Allow window to be enlarged but not made smaller
3924    dlg.SetSizer(mainSizer)
3925    w1,l1 = subSizer1.GetSize()
3926    w2,l2 = subSizer2.GetSize()
3927    panel1.SetMinSize((w1+10,200))
3928    panel2.SetMinSize((w2+20,200))
3929    aSizer.Add(panel1,1, wx.ALL|wx.EXPAND,1)
3930    aSizer.Add(panel2,2, wx.ALL|wx.EXPAND,1)
3931    mainSizer.Add(aSizer,1, wx.ALL|wx.EXPAND,1)
3932
3933    # make OK button
3934    btnsizer = wx.BoxSizer(wx.HORIZONTAL)
3935    btn = wx.Button(dlg, wx.ID_CLOSE) 
3936    btn.Bind(wx.EVT_BUTTON,_onClose)
3937    btnsizer.Add(btn)
3938    mainSizer.Add(btnsizer, 0, wx.ALIGN_CENTER|wx.ALL, 5)
3939
3940    mainSizer.Fit(dlg)
3941    dlg.SetMinSize(dlg.GetSize())
3942    dlg.ShowModal()
3943    dlg.Destroy()
Note: See TracBrowser for help on using the repository browser.