source: trunk/GSASIIconstrGUI.py @ 4838

Last change on this file since 4838 was 4838, checked in by vondreele, 8 months ago

comment out MakeDrawAtom? in G2constrGUI & make the call to G2mth.MakeDrawAtom? which does same thing
move MakeDrawAtom? & DrawAtomsReplaeyID from G2phsGUI to G2mth & fix all references to them in G2phsGUI

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