source: trunk/GSASIIconstrGUI.py @ 4608

Last change on this file since 4608 was 4608, checked in by vondreele, 21 months ago

split VectorRB listings up into separate views for each RB in both Rigid Bodies & RB Models displays.
Various bug fixes for Vector RBs

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