source: trunk/GSASIIconstrGUI.py

Last change on this file was 5002, checked in by vondreele, 8 weeks ago

allow constraints between like Debye & Bragg peak parameters.

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