source: trunk/GSASIIconstrGUI.py @ 4577

Last change on this file since 4577 was 4577, checked in by vondreele, 14 months ago

A much better fix to the C++ bug. Happens in Windows when gpx file opened directly or dropped om the G2 icon

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 160.6 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASIIconstrGUI - constraint GUI routines
3########### SVN repository information ###################
4# $Date: 2020-09-24 22:21:47 +0000 (Thu, 24 Sep 2020) $
5# $Author: vondreele $
6# $Revision: 4577 $
7# $URL: trunk/GSASIIconstrGUI.py $
8# $Id: GSASIIconstrGUI.py 4577 2020-09-24 22:21:47Z 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: 4577 $")
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.ALIGN_CENTER_VERTICAL|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|wx.ALIGN_CENTER_VERTICAL,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|wx.ALIGN_CENTER_VERTICAL,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|wx.ALIGN_CENTER_VERTICAL,5)
304            dataGridSizer.Add(scale,0,wx.RIGHT|wx.ALIGN_CENTER_VERTICAL,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|wx.ALIGN_CENTER_VERTICAL,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|wx.BOTTOM|wx.TOP, 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,wx.ALIGN_CENTER_VERTICAL|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|wx.ALIGN_CENTER_VERTICAL|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|wx.ALIGN_CENTER_VERTICAL|wx.ALIGN_CENTER,1)             # delete button
1128            if helptext:
1129                ch = G2G.HelpButton(pageDisplay,helptext)
1130                constSizer.Add(ch,0,wx.LEFT|wx.RIGHT|wx.ALIGN_CENTER_VERTICAL|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|wx.ALIGN_CENTER_VERTICAL|wx.ALIGN_CENTER,1)
1136            else:
1137                constSizer.Add((-1,-1))               
1138            constSizer.Add(wx.StaticText(pageDisplay,wx.ID_ANY,typeString),
1139                           0,wx.ALIGN_CENTER_VERTICAL|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,wx.ALIGN_CENTER_VERTICAL)
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,wx.ALIGN_CENTER_VERTICAL)
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    invTrans = nl.inv(Trans)
1585    opId = oldPhase['pId']
1586    npId = newPhase['pId']
1587    cx,ct,cs,cia = newPhase['General']['AtomPtrs']
1588    nAtoms = newPhase['Atoms']
1589    oSGData = oldPhase['General']['SGData']
1590    nSGData = newPhase['General']['SGData']
1591    #oAcof = G2lat.cell2A(oldPhase['General']['Cell'][1:7])
1592    #nAcof = G2lat.cell2A(newPhase['General']['Cell'][1:7])
1593    item = G2gd.GetGPXtreeItemId(G2frame,G2frame.root,'Constraints')
1594    if not item:
1595        print('Error: no constraints in Data Tree')
1596        return
1597    constraints = G2frame.GPXtree.GetItemPyData(item)
1598    xnames = ['dAx','dAy','dAz']
1599    # constraints on matching atom params between phases
1600    for ia,code in enumerate(atCodes):
1601        atom = nAtoms[ia]
1602        if not ia and atom[cia] == 'A':
1603            wx.MessageDialog(G2frame,
1604                'Anisotropic thermal motion constraints are not developed at the present time',
1605                'Anisotropic thermal constraint?',style=wx.ICON_INFORMATION).ShowModal()
1606        siteSym = G2spc.SytSym(atom[cx:cx+3],nSGData)[0]
1607        CSX = G2spc.GetCSxinel(siteSym)
1608#        CSU = G2spc.GetCSuinel(siteSym)
1609        item = code.split('+')[0]
1610        iat,opr = item.split(':')
1611        Nop = abs(int(opr))%100-1
1612        if '-' in opr:
1613            Nop *= -1
1614        Opr = oldPhase['General']['SGData']['SGOps'][abs(Nop)][0]
1615        if Nop < 0:         #inversion
1616            Opr *= -1
1617        XOpr = np.inner(Opr,Trans)
1618        for i,ix in enumerate(list(CSX[0])):
1619            if not ix:
1620                continue
1621            name = xnames[i]
1622            IndpCon = [1.0,G2obj.G2VarObj('%d::%s:%d'%(npId,name,ia))]
1623            DepCons = []
1624            for iop,opval in enumerate(XOpr[i]):
1625                if opval:
1626                    DepCons.append([opval,G2obj.G2VarObj('%d::%s:%s'%(opId,xnames[iop],iat))])
1627            if len(DepCons) == 1:
1628                constraints['Phase'].append([DepCons[0],IndpCon,None,None,'e'])
1629            elif len(DepCons) > 1:
1630                for Dep in DepCons:
1631                    Dep[0] *= -1
1632                constraints['Phase'].append([IndpCon]+DepCons+[0.0,None,'c'])
1633        for name in ['Afrac','AUiso']:
1634            IndpCon = [1.0,G2obj.G2VarObj('%d::%s:%d'%(npId,name,ia))]
1635            DepCons = [1.0,G2obj.G2VarObj('%d::%s:%s'%(opId,name,iat))]
1636            constraints['Phase'].append([DepCons,IndpCon,None,None,'e'])
1637           
1638        # unfinished Anisotropic constraint generation
1639#        Uids = [[0,0,'AU11'],[1,1,'AU22'],[2,2,'AU33'],[0,1,'AU12'],[0,2,'AU13'],[1,2,'AU23']]
1640#        DepConsDict = dict(zip(Us,[[],[],[],[],[],[]]))
1641#        for iu,Uid in enumerate(Uids):
1642#            UMT = np.zeros((3,3))
1643#            UMT[Uid[0],Uid[1]] = 1
1644#            nUMT = G2lat.prodMGMT(UMT,invTrans)
1645#            nUT = G2lat.UijtoU6(nUMT)
1646#            for iu,nU in enumerate(nUT):
1647#                if abs(nU) > 1.e-8:
1648#                    parm = '%d::%s;%s'%(opId,Us[iu],iat)
1649#                    DepConsDict[Uid[2]].append([abs(nU%1.),G2obj.G2VarObj(parm)])
1650#        nUcof = atom[iu:iu+6]
1651#        conStrings = []
1652#        for iU,Usi in enumerate(Us):
1653#            parm = '%d::%s;%d'%(npId,Usi,ia)
1654#            parmDict[parm] = nUcof[iU]
1655#            varyList.append(parm)
1656#            IndpCon = [1.0,G2obj.G2VarObj(parm)]
1657#            conStr = str([IndpCon,DepConsDict[Usi]])
1658#            if conStr in conStrings:
1659#                continue
1660#            conStrings.append(conStr)
1661#            if len(DepConsDict[Usi]) == 1:
1662#                if DepConsDict[Usi][0]:
1663#                    constraints['Phase'].append([IndpCon,DepConsDict[Usi][0],None,None,'e'])
1664#            elif len(DepConsDict[Usi]) > 1:       
1665#                for Dep in DepConsDict[Usi]:
1666#                    Dep[0] *= -1
1667#                constraints['Phase'].append([IndpCon]+DepConsDict[Usi]+[0.0,None,'c'])
1668           
1669        #how do I do Uij's for most Trans?
1670
1671    # constraints on lattice parameters between phases
1672#    T = nl.inv(Trans).T
1673    T = Trans.T
1674    conMat = [
1675        [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]],
1676        [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]],
1677        [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]],
1678        [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]],
1679        [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]],
1680        [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]]
1681        ]
1682    # Gnew = conMat * A:
1683#         T00**2*a0  T01**2*a1 T02**2*a2 T00*T01*a3    T00*T02*a4    T01*T02*a5
1684#         T10**2*a0  T11**2*a1 T12**2*a2 T10*T11*a3    T10*T12*a4    T11*T12*a5
1685#         T20**2*a0  T21**2*a1 T22**2*a2 T20*T21*a3    T20*T22*a4    T21*T22*a5
1686#         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
1687#         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
1688#         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
1689    # Generated as symbolic code using:
1690    # import sym
1691    # A0, A1, A2, A3, A4, A5 = sym.symbols('A0, A1, A2, A3, A4, A5')
1692    # G = sym.Matrix([ [A0,    A3/2.,  A4/2.], [A3/2.,  A1,    A5/2.], [A4/2., A5/2.,    A2]])
1693    # transformation matrix
1694    # T00, T10, T20, T01, T11, T21, T02, T12, T22 = sym.symbols('T00, T10, T20, T01, T11, T21, T02, T12, T22')
1695    # Tr = sym.Matrix([ [T00, T10, T20], [T01, T11, T21], [T02, T12, T22],])
1696    # Gnew = (Tr.T*G)*Tr
1697   
1698    #print('old A',G2lat.cell2A(oldPhase['General']['Cell'][1:7]))
1699    #print('new A',G2lat.cell2A(newPhase['General']['Cell'][1:7]))
1700   
1701#this is still incorrect for hex/trig/ortho/tetragonal --> monoclinic
1702   
1703#    for iAnew,Asi in enumerate(['A0','A1','A2','A3','A4','A5']): # loop through A[i] for new cell
1704#        Nparm = str(npId) + '::' + Asi
1705#        if Nparm != SetUniqAj(npId,iAnew,nSGData):
1706#            continue # skip: Ai constrained from Aj or must be zero
1707#        multDict = {}
1708#        for iAorg in range(6):
1709#            cA = conMat[iAnew][iAorg] # coeff for A[i] in constraint matrix
1710#            if abs(cA) < 1.e-8: continue
1711#            parm = SetUniqAj(opId,iAorg,oSGData) # translate to unique A[i] in original cell
1712#            if not parm: continue # must be zero
1713#            # sum coeff
1714#            if parm in multDict:
1715#                multDict[parm] += cA
1716#            else:
1717#                multDict[parm] = cA
1718#        # any non-zero multipliers?
1719#        maxMult = 0
1720#        for i in multDict:
1721#            maxMult = max(maxMult,abs(multDict[i]))
1722#        if maxMult <= 0:  # Nparm computes as zero; Fix this parameter
1723#            constraints['Phase'] += [[
1724#                [0.0,G2obj.G2VarObj(Nparm)],
1725#                None,None,'h']]
1726#        elif len(multDict) == 1:        # create equivalence
1727#            key = list(multDict.keys())[0]
1728#            constraints['Phase'] += [[
1729#                [1.0,G2obj.G2VarObj(key)],
1730#                [multDict[key],G2obj.G2VarObj(Nparm)],
1731#                None,None,'e']]
1732#        else:                           # create constraint
1733#            constr = [[-1.0,G2obj.G2VarObj(Nparm)]]
1734#            for key in multDict:
1735#                constr += [[multDict[key],G2obj.G2VarObj(key)]]
1736#            constr += [0.0,None,'c']
1737#            constraints['Phase'] += [constr]
1738   
1739    # constraints on HAP Scale, etc.
1740    for hId,hist in enumerate(UseList):    #HAP - seems OK
1741        ohapkey = '%d:%d:'%(opId,hId)
1742        nhapkey = '%d:%d:'%(npId,hId)
1743        IndpCon = [1.0,G2obj.G2VarObj(ohapkey+'Scale')]
1744        DepCons = [detTrans,G2obj.G2VarObj(nhapkey+'Scale')]
1745        constraints['HAP'].append([DepCons,IndpCon,None,None,'e'])
1746        for name in ['Size;i','Mustrain;i']:
1747            IndpCon = [1.0,G2obj.G2VarObj(ohapkey+name)]
1748            DepCons = [1.0,G2obj.G2VarObj(nhapkey+name)]
1749            constraints['HAP'].append([IndpCon,DepCons,None,None,'e'])
1750       
1751################################################################################
1752#### Rigid bodies
1753################################################################################
1754resRBsel = None
1755def UpdateRigidBodies(G2frame,data):
1756    '''Called when Rigid bodies tree item is selected.
1757    Displays the rigid bodies in the data window
1758    '''
1759    if not data.get('RBIds') or not data:
1760        data.update({'Vector':{'AtInfo':{}},'Residue':{'AtInfo':{}},
1761            'RBIds':{'Vector':[],'Residue':[]}})       #empty/bad dict - fill it
1762           
1763    global resList,resRBsel
1764    Indx = {}
1765    resList = []
1766    plotDefaults = {'oldxy':[0.,0.],'Quaternion':[0.,0.,0.,1.],'cameraPos':30.,'viewDir':[0,0,1],}
1767
1768    G2frame.rbBook = G2G.GSNoteBook(parent=G2frame.dataWindow)
1769    G2frame.dataWindow.GetSizer().Add(G2frame.rbBook,1,wx.ALL|wx.EXPAND)
1770    VectorRB = wx.ScrolledWindow(G2frame.rbBook)
1771    VectorRBDisplay = wx.Panel(VectorRB)
1772    G2frame.rbBook.AddPage(VectorRB,'Vector rigid bodies')
1773    ResidueRB = wx.ScrolledWindow(G2frame.rbBook)
1774    ResidueRBDisplay = wx.Panel(ResidueRB)
1775    G2frame.rbBook.AddPage(ResidueRB,'Residue rigid bodies')
1776   
1777    def OnPageChanged(event):
1778        global resList
1779        resList = []
1780        if event:       #page change event!
1781            page = event.GetSelection()
1782        else:
1783            try:
1784                page = G2frame.rbBook.GetSelection()
1785            except:
1786                if GSASIIpath.GetConfigValue('debug'): print('DBG_gpx open error:C++ Run time error - skipped')
1787                return
1788        G2frame.rbBook.ChangeSelection(page)
1789        text = G2frame.rbBook.GetPageText(page)
1790        if text == 'Vector rigid bodies':
1791            G2gd.SetDataMenuBar(G2frame,G2frame.dataWindow.VectorBodyMenu)
1792            G2frame.Bind(wx.EVT_MENU, AddVectorRB, id=G2G.wxID_VECTORBODYADD)
1793            G2frame.Bind(wx.EVT_MENU, ExtractPhaseRB, id=G2G.wxID_VECTORBODYIMP)
1794            G2frame.Bind(wx.EVT_MENU, AddVectTrans, id=G2G.wxID_VECTORBODYEXTD)
1795            G2frame.Bind(wx.EVT_MENU, SaveVectorRB, id=G2G.wxID_VECTORBODYSAV)
1796            G2frame.Bind(wx.EVT_MENU, ReadVectorRB, id=G2G.wxID_VECTORBODYRD)
1797            G2frame.Page = [page,'vrb']
1798            UpdateVectorRB()
1799        elif text == 'Residue rigid bodies':
1800            G2gd.SetDataMenuBar(G2frame,G2frame.dataWindow.RigidBodyMenu)
1801            G2frame.Bind(wx.EVT_MENU, AddResidueRB, id=G2G.wxID_RIGIDBODYADD)
1802            G2frame.Bind(wx.EVT_MENU, ExtractPhaseRB, id=G2G.wxID_RIGIDBODYIMP)
1803            G2frame.Bind(wx.EVT_MENU, OnImportRigidBody, id=G2G.wxID_RIGIDBODYIMPORT)
1804            G2frame.Bind(wx.EVT_MENU, OnSaveRigidBody, id=G2G.wxID_RIGIDBODYSAVE)
1805            G2frame.Bind(wx.EVT_MENU, OnDefineTorsSeq, id=G2G.wxID_RESIDUETORSSEQ) #enable only if residue RBs exist?
1806            G2frame.Bind(wx.EVT_MENU, DumpVectorRB, id=G2G.wxID_RESBODYSAV)
1807            G2frame.Bind(wx.EVT_MENU, LoadVectorRB, id=G2G.wxID_RESBODYRD)
1808            G2frame.Page = [page,'rrb']
1809            UpdateResidueRB()
1810        else:
1811            G2gd.SetDataMenuBar(G2frame)
1812            #G2frame.Page = [page,'rrb']
1813           
1814    def getMacroFile(macName):
1815        defDir = os.path.join(os.path.split(__file__)[0],'GSASIImacros')
1816        dlg = wx.FileDialog(G2frame,message='Choose '+macName+' rigid body macro file',
1817            defaultDir=defDir,defaultFile="",wildcard="GSAS-II macro file (*.mac)|*.mac",
1818            style=wx.FD_OPEN | wx.FD_CHANGE_DIR)
1819        try:
1820            if dlg.ShowModal() == wx.ID_OK:
1821                macfile = dlg.GetPath()
1822                macro = open(macfile,'Ur')
1823                head = macro.readline()
1824                if macName not in head:
1825                    print (head)
1826                    print ('**** ERROR - wrong restraint macro file selected, try again ****')
1827                    macro = []
1828            else: # cancel was pressed
1829                macro = []
1830        finally:
1831            dlg.Destroy()
1832        return macro        #advanced past 1st line
1833       
1834    def getTextFile():
1835        dlg = wx.FileDialog(G2frame,'Choose rigid body text file', G2frame.LastGPXdir, '',
1836            "GSAS-II text file (*.txt)|*.txt|XYZ file (*.xyz)|*.xyz|"
1837            "Sybyl mol2 file (*.mol2)|*.mol2|PDB file (*.pdb;*.ent)|*.pdb;*.ent",
1838            wx.FD_OPEN | wx.FD_CHANGE_DIR)
1839        try:
1840            if dlg.ShowModal() == wx.ID_OK:
1841                txtfile = dlg.GetPath()
1842                ext = os.path.splitext(txtfile)[1]
1843                text = open(txtfile,'Ur')
1844            else: # cancel was pressed
1845                ext = ''
1846                text = []
1847        finally:
1848            dlg.Destroy()
1849        if 'ent' in ext:
1850            ext = '.pdb'
1851        return text,ext.lower()
1852       
1853    def OnImportRigidBody(event):
1854        page = G2frame.rbBook.GetSelection()
1855        if 'Vector' in G2frame.rbBook.GetPageText(page):
1856            pass
1857        elif 'Residue' in G2frame.rbBook.GetPageText(page):
1858            ImportResidueRB()
1859           
1860    def OnSaveRigidBody(event):
1861        page = G2frame.rbBook.GetSelection()
1862        if 'Vector' in G2frame.rbBook.GetPageText(page):
1863            pass
1864        elif 'Residue' in G2frame.rbBook.GetPageText(page):
1865            SaveResidueRB()
1866           
1867    def DumpVectorRB(event):
1868        global resRBsel
1869        if resRBsel not in data['Residue']:
1870            return
1871        rbData = data['Residue'][resRBsel]
1872        pth = G2G.GetExportPath(G2frame)
1873        dlg = wx.FileDialog(G2frame, 'Choose file to save residue rigid body',
1874            pth, '', 'RRB files (*.resbody)|*.resbody',
1875            wx.FD_SAVE|wx.FD_OVERWRITE_PROMPT)
1876        try:
1877            if dlg.ShowModal() == wx.ID_OK:
1878                filename = dlg.GetPath()
1879                filename = os.path.splitext(filename)[0]+'.resbody'  # set extension
1880                fp = open(filename,'w')
1881                fp.write('Name: '+rbData['RBname']+'\n')
1882                fp.write('atNames: ')
1883                for i in rbData['atNames']:
1884                    fp.write(str(i)+" ") 
1885                fp.write('\n')
1886                for item in rbData['rbSeq']:
1887                    fp.write('rbSeq: ') 
1888                    fp.write('{:d} {:d} {:.1f}: '.format(*item[:3]))
1889                    for num in item[3]:
1890                        fp.write('{:d} '.format(num))
1891                    fp.write('\n')
1892                for i,sym in enumerate(rbData['rbTypes']):
1893                    fp.write("{:3s}".format(sym))
1894                    fp.write('{:8.5f}{:9.5f}{:9.5f}   '
1895                            .format(*rbData['rbXYZ'][i]))
1896                    fp.write('\n')
1897                fp.close()
1898                print ('Vector rigid body saved to: '+filename)
1899               
1900        finally:
1901            dlg.Destroy()
1902       
1903    def LoadVectorRB(event):
1904        AtInfo = data['Residue']['AtInfo']
1905        pth = G2G.GetExportPath(G2frame)
1906        dlg = wx.FileDialog(G2frame, 'Choose file to read vector rigid body',
1907            pth, '', 'RRB files (*.resbody)|*.resbody',
1908            wx.FD_OPEN)
1909        try:
1910            if dlg.ShowModal() == wx.ID_OK:
1911                filename = dlg.GetPath()
1912                filename = os.path.splitext(filename)[0]+'.resbody'  # set extension
1913                fp = open(filename,'r')
1914                l = fp.readline().strip()
1915                if 'Name' not in l:
1916                    fp.close()
1917                    G2frame.ErrorDialog('Read Error',
1918                        'File '+filename+' does not start with Name\nFirst line ='
1919                        +l+'\ninvalid file',parent=G2frame)
1920                    return
1921                name = l.split(':')[1].strip()
1922                line = fp.readline().strip().split(':')[1].split()
1923                atNames = [i for i in line]
1924                types = []
1925                coords = []
1926                l = fp.readline().strip()
1927                rbSeq = []
1928                while 'rbSeq' in l:
1929                    tag,vals,lst = l.split(':')
1930                    seq = []
1931                    for t,v in zip((int,int,float),vals.split()): 
1932                        seq.append(t(v))
1933                    seq.append([])
1934                    for num in lst.split():
1935                        seq[-1].append(int(num))
1936                    rbSeq.append(seq)
1937                    l = fp.readline().strip()                   
1938                while l:
1939                    nums = l.strip().split()
1940                    types.append(nums.pop(0))
1941                    t = types[-1]
1942                    if t not in AtInfo:
1943                        Info = G2elem.GetAtomInfo(t)
1944                        AtInfo[t] = [Info['Drad'],Info['Color']]
1945                    coords.append([float(nums.pop(0)) for j in range(3)])
1946                    l = fp.readline().strip()
1947                fp.close()
1948            else:
1949                return       
1950        finally:
1951            dlg.Destroy()
1952        coords = np.array(coords)
1953        rbid = ran.randint(0,sys.maxsize)
1954        data['Residue'][rbid] = {'RBname':name,
1955                'rbXYZ': coords,
1956                'rbRef':[0,1,2,False],
1957                'rbTypes':types, 'atNames':atNames,
1958                'useCount':0,
1959                'rbSeq':rbSeq, 'SelSeq':[0,0],}
1960        data['RBIds']['Residue'].append(rbid)
1961        UpdateResidueRB()
1962   
1963    def AddVectorRB(event):
1964        'Create a new vector rigid body'
1965        AtInfo = data['Vector']['AtInfo']
1966        dlg = G2G.MultiIntegerDialog(G2frame,'New Rigid Body',['No. atoms','No. translations'],[1,1])
1967        if dlg.ShowModal() == wx.ID_OK:
1968            nAtoms,nTrans = dlg.GetValues()
1969            rbid = ran.randint(0,sys.maxsize)
1970            vecMag = [1.0 for i in range(nTrans)]
1971            vecRef = [False for i in range(nTrans)]
1972            vecVal = [np.zeros((nAtoms,3)) for j in range(nTrans)]
1973            rbTypes = ['C' for i in range(nAtoms)]
1974            Info = G2elem.GetAtomInfo('C')
1975            AtInfo['C'] = [Info['Drad'],Info['Color']]
1976            data['Vector'][rbid] = {'RBname':'UNKRB','VectMag':vecMag,'rbXYZ':np.zeros((nAtoms,3)),
1977                'rbRef':[0,1,2,False],'VectRef':vecRef,'rbTypes':rbTypes,'rbVect':vecVal,'useCount':0}
1978            data['RBIds']['Vector'].append(rbid)
1979        dlg.Destroy()
1980        UpdateVectorRB()
1981
1982    def ExtractPhaseRB(event):
1983        'Extract a rigid body from a file with a phase'
1984        def SetupDrawing(atmData):
1985            '''Add the dicts needed for G2plt.PlotStructure to work to the
1986            reader .Phase object
1987            '''
1988            generalData = atmData['General']
1989            generalData['BondRadii'] = []
1990
1991            G2phG.SetDrawingDefaults(atmData['Drawing'])
1992            atmData['Drawing'].update(
1993                {'oldxy':[0.,0.],'Quaternion':[0.,0.,0.,1.],'cameraPos':150.,
1994                     'viewDir':[0,0,1],'atomPtrs': [2, 1, 6, 17],
1995                     })
1996            atmData['Drawing']['showRigidBodies'] = False
1997            generalData['Map'] = {'MapType':False, 'rho':[]}
1998            generalData['AtomTypes'] = []
1999            generalData['BondRadii'] = []
2000            generalData['AngleRadii'] = []
2001            generalData['vdWRadii'] = []
2002            generalData['Color'] = []
2003            generalData['Isotopes'] = {}
2004            generalData['Isotope'] = {}
2005            cx,ct,cs,cia = generalData['AtomPtrs']
2006            generalData['Mydir'] = G2frame.dirname
2007            for iat,atom in enumerate(atmData['Atoms']):
2008                atom[ct] = atom[ct].lower().capitalize()      #force elem symbol to standard form
2009                if atom[ct] not in generalData['AtomTypes'] and atom[ct] != 'UNK':
2010                    Info = G2elem.GetAtomInfo(atom[ct])
2011                    if not Info:
2012                        atom[ct] = 'UNK'
2013                        continue
2014                    atom[ct] = Info['Symbol'] # N.B. symbol might be changed by GetAtomInfo
2015                    generalData['AtomTypes'].append(atom[ct])
2016                    generalData['Z'] = Info['Z']
2017                    generalData['Isotopes'][atom[ct]] = Info['Isotopes']
2018                    generalData['BondRadii'].append(Info['Drad'])
2019                    generalData['AngleRadii'].append(Info['Arad'])
2020                    generalData['vdWRadii'].append(Info['Vdrad'])
2021                    if atom[ct] in generalData['Isotope']:
2022                        if generalData['Isotope'][atom[ct]] not in generalData['Isotopes'][atom[ct]]:
2023                            isotope = list(generalData['Isotopes'][atom[ct]].keys())[-1]
2024                            generalData['Isotope'][atom[ct]] = isotope
2025                    else:
2026                        generalData['Isotope'][atom[ct]] = 'Nat. Abund.'
2027                        if 'Nat. Abund.' not in generalData['Isotopes'][atom[ct]]:
2028                            isotope = list(generalData['Isotopes'][atom[ct]].keys())[-1]
2029                            generalData['Isotope'][atom[ct]] = isotope
2030                    generalData['Color'].append(Info['Color'])
2031                    if generalData['Type'] == 'magnetic':
2032                        if len(landeg) < len(generalData['AtomTypes']):
2033                            landeg.append(2.0)
2034            atmData['Drawing']['Atoms'] = []
2035            for atom in atmData['Atoms']:
2036                atmData['Drawing']['Atoms'].append(MakeDrawAtom(atmData,atom))
2037
2038        def onCancel(event,page=0):
2039            'complete or bail out from RB define, cleaning up'
2040            G2frame.rbBook.DeletePage(G2frame.rbBook.FindPage(pagename))
2041            G2frame.rbBook.SetSelection(page)
2042
2043        def Page1():
2044            '''Show the GUI for first stage of the rigid body with all atoms in
2045            phase in crystal coordinates. Select the atoms to go onto the
2046            next stage
2047            '''
2048            def ShowSelection(selections):
2049                'respond to change in atom selections'
2050                ct,cs = [1,8]
2051                generalData = rd.Phase['General']
2052                for i,atom in enumerate(rd.Phase['Drawing']['Atoms']):
2053                    if i in selections:
2054                        factor = 1
2055                    else:
2056                        factor = 2.5
2057                    atNum = generalData['AtomTypes'].index(atom[ct]) 
2058                    atom[cs] = list(np.array(generalData['Color'][atNum])//factor) 
2059                draw(*drawArgs)
2060            def onPage1OK(event):
2061                '1st section has been completed, move onto next'
2062                G2frame.G2plotNB.Delete(rd.Phase['General']['Name'])
2063                GetCoords(atmsel)
2064                Page2()
2065
2066            if 'macromolecular' == rd.Phase['General']['Type']:
2067                # for PDB imports, lets see if a quick reformat of atoms list will work
2068                rd.Phase['Atoms'] = [a[3:] for a in rd.Phase['Atoms']]
2069                rd.Phase['General']['AtomPtrs']  = [i-3 for i in rd.Phase['General']['AtomPtrs']]
2070                rd.Phase['General']['Type'] = 'nuclear'
2071            SetupDrawing(rd.Phase) # add information to reader object to allow plotting
2072            atomlist = [atom[0] for atom in rd.Phase['Atoms']]
2073            atmsel = list(range(len(rd.Phase['Atoms'])))
2074            # broken -- # why no bonds?
2075            #for atm in rd.Phase['Drawing']['Atoms']:
2076            #    atm[6] = 'balls & sticks'
2077
2078            draw,drawArgs = G2plt.PlotStructure(G2frame,rd.Phase,True)
2079            ShowSelection(atmsel)
2080
2081            if G2frame.rbBook.FindPage(pagename) is not None:
2082                G2frame.rbBook.DeletePage(G2frame.rbBook.FindPage(pagename))
2083
2084            RBImp = wx.ScrolledWindow(G2frame.rbBook)
2085            RBImpPnl = wx.Panel(RBImp)
2086            G2frame.rbBook.AddPage(RBImp,pagename)
2087            G2frame.rbBook.SetSelection(G2frame.rbBook.FindPage(pagename))
2088
2089            mainSizer = G2G.G2MultiChoiceWindow(RBImpPnl,
2090                            'Select atoms to import',
2091                            atomlist,atmsel,OnChange=ShowSelection)
2092
2093            # OK/Cancel buttons       
2094            btnsizer = wx.StdDialogButtonSizer()
2095            OKbtn = wx.Button(RBImpPnl, wx.ID_OK, 'Continue')
2096            OKbtn.SetDefault()
2097            btnsizer.AddButton(OKbtn)
2098            OKbtn.Bind(wx.EVT_BUTTON,onPage1OK)
2099            btn = wx.Button(RBImpPnl, wx.ID_CANCEL)
2100            btn.Bind(wx.EVT_BUTTON,onCancel)
2101            btnsizer.AddButton(btn)
2102            btnsizer.Realize()
2103            mainSizer.Add(btnsizer,0,wx.ALIGN_CENTER,50)
2104
2105            RBImpPnl.SetSizer(mainSizer,True)
2106
2107            mainSizer.Layout()   
2108            Size = mainSizer.GetMinSize()
2109            Size[0] += 40
2110            Size[1] = max(Size[1],G2frame.GetSize()[1]-200) + 20
2111            RBImpPnl.SetSize(Size)
2112            RBImp.SetScrollbars(10,10,Size[0]/10-4,Size[1]/10-1)
2113            RBImp.Scroll(0,0)
2114
2115        def Page2():
2116            '''Show the GUI for the second stage, where selected atoms are
2117            now in Cartesian space, manipulate the axes and export selected
2118            atoms to a vector or residue rigid body.
2119            '''
2120            def UpdateDraw(event=None):
2121                'Called when info changes in grid, replots'
2122                UpdateVectorBody(rbData)
2123                DrawCallback()
2124               
2125            def onSetAll(event):
2126                'Set all atoms as selected'
2127                grid.completeEdits()
2128                for i in range(len(rd.Phase['RBselection'])):
2129                    rd.Phase['RBselection'][i] = 1 # table needs 0/1 for T/F
2130                grid.ForceRefresh()
2131                UpdateDraw()
2132               
2133            def onToggle(event):
2134                'Toggles selection state for all atoms'
2135                grid.completeEdits()
2136                for i in range(len(rd.Phase['RBselection'])):
2137                    rd.Phase['RBselection'][i] = int(not rd.Phase['RBselection'][i])
2138                grid.ForceRefresh()
2139                UpdateDraw()
2140               
2141            def onSetOrigin(event):
2142                'Resets origin to midpoint between all selected atoms'
2143                grid.completeEdits()
2144                center = np.array([0.,0.,0.])
2145                count = 0
2146                for i in range(len(rd.Phase['RBselection'])):
2147                    if rd.Phase['RBselection'][i]:
2148                        count += 1
2149                        center += rd.Phase['RBcoords'][i]
2150                if count:
2151                    rd.Phase['RBcoords'] -= center/count
2152                grid.ForceRefresh()
2153                UpdateDraw()
2154               
2155            def onSetX(event):
2156                grid.completeEdits()
2157                center = np.array([0.,0.,0.])
2158                count = 0
2159                for i in range(len(rd.Phase['RBselection'])):
2160                    if rd.Phase['RBselection'][i]:
2161                        count += 1
2162                        center += rd.Phase['RBcoords'][i]
2163                if not count:
2164                    G2G.G2MessageBox(G2frame,'No atoms selected',
2165                                    'Selection required')
2166                    return
2167                XYZP = center/count
2168                if np.sqrt(sum(XYZP**2)) < 0.1:
2169                    G2G.G2MessageBox(G2frame,
2170                            'The selected atom(s) are too close to the origin',
2171                            'near origin')
2172                    return
2173                if bntOpts['direction'] == 'y':
2174                    YP = XYZP / np.sqrt(np.sum(XYZP**2))
2175                    ZP = np.cross((1,0,0),YP)
2176                    if sum(ZP*ZP) < .1: # pathological condition: Y' along X
2177                        ZP = np.cross((0,0,1),YP)
2178                    XP = np.cross(YP,ZP)
2179                elif bntOpts['direction'] == 'z':
2180                    ZP = XYZP / np.sqrt(np.sum(XYZP**2))
2181                    XP = np.cross((0,1,0),ZP)
2182                    if sum(XP*XP) < .1: # pathological condition: X' along Y
2183                        XP = np.cross((0,0,1),ZP)
2184                    YP = np.cross(ZP,XP)
2185                else:
2186                    XP = XYZP / np.sqrt(np.sum(XYZP**2))
2187                    YP = np.cross((0,0,1),XP)
2188                    if sum(YP*YP) < .1: # pathological condition: X' along Z
2189                        YP = np.cross((0,1,0),XP)
2190                    ZP = np.cross(XP,YP)
2191                trans = np.array((XP,YP,ZP))
2192                # update atoms in place
2193                rd.Phase['RBcoords'][:] = np.inner(trans,rd.Phase['RBcoords']).T
2194                grid.ForceRefresh()
2195                UpdateDraw()
2196
2197            def onSetPlane(event): 
2198                '''Compute least-squares plane for selected atoms;
2199                move atoms so that LS plane aligned with x-y plane,
2200                with minimum required change to x
2201                '''
2202                grid.completeEdits()
2203                selList = [i==1 for i in rd.Phase['RBselection']]
2204                XYZ = rd.Phase['RBcoords'][selList]
2205                if len(XYZ) < 3: 
2206                    G2G.G2MessageBox(G2frame,'A plane requires three or more atoms',
2207                                     'Need more atoms')
2208                    return
2209                # fit 3 ways (in case of singularity) and take result with lowest residual
2210                X,Y,Z = [XYZ[:,i] for i in (0,1,2)]
2211                XZ = copy.copy(XYZ)
2212                XZ[:,1] = 1
2213                (a,d,b), resd, rank, sing = nl.lstsq(XZ, -Y)
2214                resid_min = resd
2215                normal = a,1,b
2216                YZ = copy.copy(XYZ)
2217                YZ[:,0] = 1
2218                (d,a,b), resd, rank, sing = nl.lstsq(YZ, -X)
2219                if resid_min > resd:
2220                    resid_min = resd
2221                    normal = 1,a,b
2222                XY = copy.copy(XYZ)
2223                XY[:,2] = 1
2224                (a,b,d), resd, rank, sing = nl.lstsq(XY, -Z)
2225                if resid_min > resd:
2226                    resid_min = resd
2227                    normal = a,b,1
2228                # solve for  ax + bx + z + c = 0 or equivalently ax + bx + c = -z
2229                # try:
2230                # except:
2231                #     G2G.G2MessageBox(G2frame,
2232                #             'Error computing plane; are atoms in a line?',
2233                #             'Computation error')
2234                #     return
2235                if bntOpts['plane'] == 'xy':
2236                    # new coordinate system is
2237                    #   ZP, z' normal to plane
2238                    #   YP, y' = z' cross x (= [0,1,-b])
2239                    #   XP, x' = (z' cross x) cross z'
2240                    # this puts XP as close as possible to X with XP & YP in plane
2241                    ZP = np.array(normal)
2242                    ZP /= np.sqrt(np.sum(ZP**2))
2243                    YP = np.cross(ZP,[1,0,0])
2244                    if sum(YP*YP) < .1: # pathological condition: z' along x
2245                        YP = np.cross(ZP,[0,1,0])
2246                    YP /= np.sqrt(np.sum(YP**2))
2247                    XP = np.cross(YP,ZP)
2248                elif bntOpts['plane'] == 'yz':
2249                    # new coordinate system is
2250                    #   XP, x' normal to plane
2251                    #   ZP, z' = x' cross y
2252                    #   YP, y' = (x' cross y) cross x'
2253                    # this puts y' as close as possible to y with z' & y' in plane
2254                    XP = np.array(normal)
2255                    XP /= np.sqrt(np.sum(XP**2))
2256                    ZP = np.cross(XP,[0,1,0])
2257                    if sum(ZP*ZP) < .1: # pathological condition: x' along y
2258                        ZP = np.cross(XP,(0,0,1))
2259                    ZP /= np.sqrt(np.sum(ZP**2))
2260                    YP = np.cross(ZP,XP)
2261                elif bntOpts['plane'] == 'xz':
2262                    # new coordinate system is
2263                    #   YP, y' normal to plane
2264                    #   ZP, z' = x cross y'
2265                    #   XP, y' = (x cross y') cross z'
2266                    # this puts XP as close as possible to X with XP & YP in plane
2267                    YP = np.array(normal)
2268                    YP /= np.sqrt(np.sum(YP**2))
2269                    ZP = np.cross([1,0,0],YP)
2270                    if sum(ZP*ZP) < .1: # pathological condition: y' along x
2271                        ZP = np.cross([0,1,0],YP)
2272                    ZP /= np.sqrt(np.sum(ZP**2))
2273                    XP = np.cross(YP,ZP)
2274                else:
2275                    print('unexpected plane',bntOpts['plane'])
2276                    return
2277                trans = np.array((XP,YP,ZP))
2278                # update atoms in place
2279                rd.Phase['RBcoords'][:] = np.inner(trans,rd.Phase['RBcoords']).T
2280                grid.ForceRefresh()
2281                UpdateDraw()
2282
2283            def onWriteXYZ(event):
2284                '''Writes selected atoms in a .xyz file for use in Avogadro, etc.
2285                '''
2286                grid.completeEdits()
2287                center = np.array([0.,0.,0.])
2288                count = 0
2289                for i in range(len(rd.Phase['RBselection'])):
2290                    if rd.Phase['RBselection'][i]:
2291                        count += 1
2292                        center += rd.Phase['RBcoords'][i]
2293                if count:
2294                    center /= count
2295                else:
2296                    print('nothing selected')
2297                    return
2298                obj = G2IO.ExportBaseclass(G2frame,'XYZ','.xyz')
2299                #obj.InitExport(None)
2300                if obj.ExportSelect():    # set export parameters; ask for file name
2301                    return
2302                fp = obj.OpenFile()
2303                obj.Write(str(count))
2304                obj.Write('')
2305                for i in range(len(rd.Phase['RBselection'])):
2306                    if rd.Phase['RBselection'][i]:
2307                        line = ' ' + rd.Phase['RBtypes'][i]
2308                        for xyz in rd.Phase['RBcoords'][i]:
2309                            line += ' ' + str(xyz)
2310                        obj.Write(line)
2311                obj.CloseFile()
2312                #GSASIIpath.IPyBreak()
2313               
2314            def onAddVector(event):
2315                '''Adds selected atoms as a new vector rigid body.
2316                Closes out the importer tab when done.
2317                '''
2318                grid.completeEdits()
2319                rb = MakeVectorBody(os.path.split(filename)[1])
2320                UpdateVectorBody(rb,True)
2321                if len(rb['rbTypes']) < 3: return # must have at least 3 atoms
2322                rbid = ran.randint(0,sys.maxsize)
2323                data['Vector'][rbid] = rb
2324                data['RBIds']['Vector'].append(rbid)
2325                AtInfo = data['Vector']['AtInfo']
2326                for t in rb['rbTypes']:
2327                    if t in data['Vector']['AtInfo']: continue
2328                    Info = G2elem.GetAtomInfo(t)
2329                    data['Vector']['AtInfo'][t] = [Info['Drad'],Info['Color']]
2330                G2frame.G2plotNB.Delete('Rigid body')
2331                onCancel(event,0)
2332               
2333            def onAddResidue(event):
2334                '''Adds selected atoms as a new residue rigid body.
2335                Closes out the importer tab when done.
2336                '''
2337                grid.completeEdits()
2338                name = os.path.split(filename)[1]
2339                rbXYZ = []
2340                rbTypes = []
2341                atNames = []
2342                for i in rd.Phase['RBindex']:
2343                    if rd.Phase['RBselection'][i]:
2344                        rbXYZ.append(rd.Phase['RBcoords'][i])
2345                        rbTypes.append(rd.Phase['RBtypes'][i])
2346                        atNames.append(rd.Phase['RBlbls'][i])
2347                if len(rbTypes) < 3: return # must have at least 3 atoms
2348                rbXYZ = np.array(rbXYZ)
2349                rbid = ran.randint(0,sys.maxsize)
2350                data['Residue'][rbid] = {'RBname':name,'rbXYZ':rbXYZ,
2351                    'rbTypes':rbTypes,'atNames':atNames,'rbRef':[0,1,2,False],
2352                    'rbSeq':[],'SelSeq':[0,0],'useCount':0}
2353                data['RBIds']['Residue'].append(rbid)
2354                AtInfo = data['Residue']['AtInfo']
2355                for t in rbTypes:
2356                    if t in data['Residue']['AtInfo']: continue
2357                    Info = G2elem.GetAtomInfo(t)
2358                    data['Residue']['AtInfo'][t] = [Info['Drad'],Info['Color']]
2359
2360                print ('Rigid body added')
2361                G2frame.G2plotNB.Delete('Rigid body')
2362                onCancel(event,1)
2363
2364            if G2frame.rbBook.FindPage(pagename) is not None:
2365                G2frame.rbBook.DeletePage(G2frame.rbBook.FindPage(pagename))
2366            RBImp = wx.ScrolledWindow(G2frame.rbBook)
2367            RBImpPnl = wx.Panel(RBImp)
2368            G2frame.rbBook.AddPage(RBImp,pagename)
2369            G2frame.rbBook.SetSelection(G2frame.rbBook.FindPage(pagename))
2370            generalData = rd.Phase['General']
2371            AtInfo = {}
2372            ct = 1
2373            for t in rd.Phase['RBtypes']:
2374                if t in AtInfo: continue
2375                Info = G2elem.GetAtomInfo(t)
2376                AtInfo[t] = [Info['Drad'],Info['Color']]
2377            plotDefaults = {'oldxy':[0.,0.],'Quaternion':[0.,0.,0.,1.],'cameraPos':30.,'viewDir':[0,0,1],}
2378
2379            rd.Phase['RBindex'] = list(range(len(rd.Phase['RBtypes'])))
2380            rd.Phase['RBselection'] = len(rd.Phase['RBtypes']) * [1]
2381            rbData = MakeVectorBody()
2382            DrawCallback = G2plt.PlotRigidBody(G2frame,'Vector',
2383                                    AtInfo,rbData,plotDefaults)
2384
2385            mainSizer = wx.BoxSizer(wx.HORIZONTAL)
2386            btnSizer = wx.BoxSizer(wx.VERTICAL)
2387            btnSizer.Add(
2388                wx.StaticText(RBImpPnl,wx.ID_ANY,'Reorder atoms by dragging'),
2389                0,wx.ALL)
2390            btnSizer.Add((-1,15))
2391            btn = wx.Button(RBImpPnl, wx.ID_ANY, 'Set All')
2392            btn.Bind(wx.EVT_BUTTON,onSetAll)
2393            btnSizer.Add(btn,0,wx.ALIGN_CENTER)
2394            btn = wx.Button(RBImpPnl, wx.ID_ANY, 'Toggle')
2395            btn.Bind(wx.EVT_BUTTON,onToggle)
2396            btnSizer.Add(btn,0,wx.ALIGN_CENTER)
2397            btnSizer.Add((-1,15))
2398            btnSizer.Add(
2399                wx.StaticText(RBImpPnl,wx.ID_ANY,'Reorient using selected\natoms...'),
2400                0,wx.ALL)
2401            btnSizer.Add((-1,5))
2402            btn = wx.Button(RBImpPnl, wx.ID_ANY, 'Set origin')
2403            btn.Bind(wx.EVT_BUTTON,onSetOrigin)
2404            btnSizer.Add(btn,0,wx.ALIGN_CENTER)
2405
2406            bntOpts = {'plane':'xy','direction':'x'}
2407            inSizer = wx.BoxSizer(wx.HORIZONTAL)
2408            btn = wx.Button(RBImpPnl, wx.ID_ANY, 'Place in plane')
2409            btn.Bind(wx.EVT_BUTTON,onSetPlane)
2410            inSizer.Add(btn)
2411            inSizer.Add(G2G.G2ChoiceButton(RBImpPnl,('xy','yz','xz'),
2412                                           None,None,bntOpts,'plane'))
2413            btnSizer.Add(inSizer,0,wx.ALIGN_CENTER)
2414           
2415            inSizer = wx.BoxSizer(wx.HORIZONTAL)
2416            btn = wx.Button(RBImpPnl, wx.ID_ANY, 'Define as')
2417            btn.Bind(wx.EVT_BUTTON,onSetX)
2418            inSizer.Add(btn)
2419            inSizer.Add(G2G.G2ChoiceButton(RBImpPnl,('x','y','z'),
2420                                           None,None,bntOpts,'direction'))
2421            btnSizer.Add(inSizer,0,wx.ALIGN_CENTER)
2422           
2423            btnSizer.Add((-1,15))
2424            btnSizer.Add(
2425                wx.StaticText(RBImpPnl,wx.ID_ANY,'Use selected atoms to\ncreate...'),
2426                0,wx.ALL)
2427            btnSizer.Add((-1,5))
2428            btn = wx.Button(RBImpPnl, wx.ID_ANY, 'export as xyz')
2429            btn.Bind(wx.EVT_BUTTON,onWriteXYZ)
2430            btnSizer.Add(btn,0,wx.ALIGN_CENTER)
2431            btnSizer.Add((-1,10))
2432            btn = wx.Button(RBImpPnl, wx.ID_ANY, 'a Vector Body')
2433            btn.Bind(wx.EVT_BUTTON,onAddVector)
2434            btnSizer.Add(btn,0,wx.ALIGN_CENTER)
2435            btn = wx.Button(RBImpPnl, wx.ID_ANY, 'a Residue Body')
2436            btn.Bind(wx.EVT_BUTTON,onAddResidue)
2437            btnSizer.Add(btn,0,wx.ALIGN_CENTER)
2438            btn = wx.Button(RBImpPnl, wx.ID_CANCEL)
2439            btn.Bind(wx.EVT_BUTTON,onCancel)
2440            btnSizer.Add((-1,10))
2441            btnSizer.Add(btn,0,wx.ALIGN_CENTER)
2442
2443            mainSizer.Add(btnSizer)
2444            mainSizer.Add((5,5))
2445            grid = DragableRBGrid(RBImpPnl,rd.Phase,UpdateDraw)
2446            mainSizer.Add(grid)
2447            RBImpPnl.SetSizer(mainSizer,True)
2448            mainSizer.Layout()   
2449            Size = mainSizer.GetMinSize()
2450            Size[0] += 40
2451            Size[1] = max(Size[1],G2frame.GetSize()[1]-200) + 20
2452            RBImpPnl.SetSize(Size)
2453            RBImp.SetScrollbars(10,10,Size[0]/10-4,Size[1]/10-1)
2454            RBImp.Scroll(0,0)
2455
2456        def GetCoords(atmsel):
2457            '''Create orthogonal coordinates for selected atoms.
2458            Place the origin at the center of the body
2459            '''
2460            atms = rd.Phase['Atoms']
2461            cell = rd.Phase['General']['Cell'][1:7]
2462            Amat,Bmat = G2lat.cell2AB(cell)
2463            rd.Phase['RBcoords'] = np.array([np.inner(Amat,atms[i][3:6]) for i in atmsel])
2464            rd.Phase['RBcoords'] -= rd.Phase['RBcoords'].mean(axis=0)  # origin to middle
2465            rd.Phase['RBtypes'] = [atms[i][1] for i in atmsel]
2466            rd.Phase['RBlbls'] = [atms[i][0] for i in atmsel]
2467
2468        def UpdateVectorBody(rb,useSelection=False):
2469            '''Put the atoms in order to pass for plotting or for storage as
2470            a vector rigid body.
2471
2472            :param dict rb: rigid body contents created in :func:`MakeVectorBody`
2473            :param bool useSelection: True if the rd.Phase['RBselection']
2474              values will be used to select which atoms are included in the
2475              rigid body. If False (default) they are included in rb
2476              and are used for plotting.         
2477            '''
2478            coordlist = []
2479            typeslist = []
2480            sellist = []
2481            for i in rd.Phase['RBindex']:
2482                use = True
2483                if useSelection and not rd.Phase['RBselection'][i]: use = False
2484                if use:
2485                    coordlist.append(rd.Phase['RBcoords'][i])
2486                    typeslist.append(rd.Phase['RBtypes'][i])
2487                    sellist.append(rd.Phase['RBselection'][i])
2488            coordlist = np.array(coordlist)
2489            rb['rbXYZ'] = coordlist
2490            rb['rbVect'] = [coordlist]
2491            rb['rbTypes'] = typeslist
2492            if not useSelection:
2493                rb['Selection'] = sellist
2494            elif 'Selection' in rb:
2495                del rb['Selection']
2496
2497        def MakeVectorBody(name=''):
2498            '''Make the basic vector rigid body dict (w/o coordinates) used for
2499            export and for plotting
2500            '''
2501            nTrans = 1
2502            vecMag = [1.0]
2503            vecRef = [False]
2504            rb = {'RBname':name,'VectMag':vecMag,
2505                    'rbRef':[0,1,2,False],'VectRef':vecRef,
2506                    'useCount':0}
2507            UpdateVectorBody(rb)
2508            return rb
2509
2510        # too lazy to figure out why wx crashes
2511        if wx.__version__.split('.')[0] != '4':
2512            wx.MessageBox('Sorry, wxPython 4.x is required to run this command',
2513                                  caption='Update Python',
2514                                  style=wx.ICON_EXCLAMATION)
2515            return
2516        if platform.python_version()[:1] == '2':
2517            wx.MessageBox('Sorry, Python >=3.x is required to run this command',
2518                                  caption='Update Python',
2519                                  style=wx.ICON_EXCLAMATION)
2520            return
2521
2522        # get importer type and a phase file of that type
2523        G2sc.LoadG2fil()
2524        choices = [rd.formatName for  rd in G2sc.Readers['Phase']] 
2525        dlg = G2G.G2SingleChoiceDialog(G2frame,'Select the format of the file',
2526                                     'select format',choices)
2527        dlg.CenterOnParent()
2528        try:
2529            if dlg.ShowModal() == wx.ID_OK:
2530                col = dlg.GetSelection()
2531            else:
2532                col = None
2533                return
2534        finally:
2535            dlg.Destroy()
2536        reader = G2sc.Readers['Phase'][col]
2537
2538        choices = reader.formatName + " file ("
2539        w = ""
2540        for extn in reader.extensionlist:
2541            if w != "": w += ";"
2542            w += "*" + extn
2543        choices += w + ")|" + w
2544        #choices += "|zip archive (.zip)|*.zip"
2545        if not reader.strictExtension:
2546            choices += "|any file (*.*)|*.*"
2547        typ = '( type '+reader.formatName+')'
2548        filelist = G2G.GetImportFile(G2frame,
2549                        message="Choose phase input file"+typ,
2550                        defaultFile="",wildcard=choices,style=wx.FD_OPEN)
2551        if len(filelist) != 1: return
2552
2553        # read in the phase file
2554        filename = filelist[0]
2555        rd = reader
2556        with open(filename, 'Ur') as fp:
2557            rd.ReInitialize()
2558            rd.errors = ""
2559            if not rd.ContentsValidator(filename):   # Report error
2560                G2fl.G2Print("Warning: File {} has a validation error".format(filename))
2561                return
2562            if len(rd.selections) > 1:
2563                print("File {} has {} phases. This is unexpected."
2564                                    .format(filename,len(rd.selections)))
2565                return
2566
2567            rd.objname = os.path.basename(filename)
2568            try:
2569                flag = rd.Reader(filename)
2570            except Exception as msg:
2571                G2fl.G2Print("Warning: read of file {} failed\n{}".format(
2572                    filename,rd.errors))
2573                if GSASIIpath.GetConfigValue('debug'):
2574                    print(msg)
2575                    import traceback
2576                    print (traceback.format_exc())
2577                    GSASIIpath.IPyBreak()
2578                return
2579
2580        pagename = 'Rigid body importer'
2581        Page1()
2582        return
2583
2584    def AddVectTrans(event):
2585        'Add a translation to an existing vector rigid body'
2586        choices = []
2587        rbIdlist = []
2588        for rbid in data['RBIds']['Vector']:
2589            if rbid != 'AtInfo':
2590                rbIdlist.append(rbid)
2591                choices.append(data['Vector'][rbid]['RBname'])
2592        if len(choices) == 0:
2593            G2G.G2MessageBox(G2frame,'No Vector Rigid Bodies found',
2594                                 'No VR Bodies')
2595            return
2596        elif len(choices) == 1:
2597            rbid = rbIdlist[0]
2598        else:
2599            dlg = G2G.G2SingleChoiceDialog(G2frame,'Select the rigid body to save',
2600                                  'select format',choices)
2601            try:
2602                if dlg.ShowModal() == wx.ID_OK:
2603                    rbid = rbIdlist[dlg.GetSelection()]
2604                else:
2605                    return
2606            finally:
2607                dlg.Destroy()
2608        data['Vector'][rbid]['VectMag'] += [1.0]
2609        data['Vector'][rbid]['VectRef'] += [False]
2610        nAtoms = len(data['Vector'][rbid]['rbXYZ'])
2611        data['Vector'][rbid]['rbVect'] += [np.zeros((nAtoms,3))]
2612        UpdateVectorRB()
2613       
2614    def SaveVectorRB(event):
2615        choices = []
2616        rbIdlist = []
2617        for rbid in data['RBIds']['Vector']:
2618            if rbid != 'AtInfo':
2619                rbIdlist.append(rbid)
2620                choices.append(data['Vector'][rbid]['RBname'])
2621        if len(choices) == 0:
2622            G2G.G2MessageBox(G2frame,'No Vector Rigid Bodies found',
2623                                 'No VR Bodies')
2624            return
2625        elif len(choices) == 1:
2626            rbid = rbIdlist[0]
2627        else:
2628            dlg = G2G.G2SingleChoiceDialog(G2frame,'Select the rigid body to save',
2629                                  'select format',choices)
2630            try:
2631                if dlg.ShowModal() == wx.ID_OK:
2632                    rbid = rbIdlist[dlg.GetSelection()]
2633                else:
2634                    return
2635            finally:
2636                dlg.Destroy()
2637             
2638        pth = G2G.GetExportPath(G2frame)
2639        dlg = wx.FileDialog(G2frame, 'Choose file to save vector rigid body',
2640            pth, '', 'VRB files (*.vecbody)|*.vecbody',
2641            wx.FD_SAVE|wx.FD_OVERWRITE_PROMPT)
2642        try:
2643            if dlg.ShowModal() == wx.ID_OK:
2644                filename = dlg.GetPath()
2645                filename = os.path.splitext(filename)[0]+'.vecbody'  # set extension
2646                fp = open(filename,'w')
2647                fp.write('Name: '+data['Vector'][rbid]['RBname']+'\n')
2648                fp.write('Trans: ')
2649                for i in data['Vector'][rbid]['VectMag']:
2650                    fp.write(str(i)+" ") 
2651                fp.write('\n')
2652                ntrans = len(data['Vector'][rbid]['VectMag'])
2653                for i,sym in enumerate(data['Vector'][rbid]['rbTypes']):
2654                    fp.write("{:3s}".format(sym))
2655                    for j in range(ntrans):
2656                        fp.write('{:8.5f}{:9.5f}{:9.5f}   '
2657                            .format(*data['Vector'][rbid]['rbVect'][j][i]))
2658                    fp.write('\n')
2659                fp.close()
2660                print ('Vector rigid body saved to: '+filename)
2661        finally:
2662            dlg.Destroy()
2663           
2664    def ReadVectorRB(event):
2665        AtInfo = data['Vector']['AtInfo']
2666        pth = G2G.GetExportPath(G2frame)
2667        dlg = wx.FileDialog(G2frame, 'Choose file to read vector rigid body',
2668            pth, '', 'VRB files (*.vecbody)|*.vecbody',
2669            wx.FD_OPEN)
2670        try:
2671            if dlg.ShowModal() == wx.ID_OK:
2672                filename = dlg.GetPath()
2673                filename = os.path.splitext(filename)[0]+'.vecbody'  # set extension
2674                fp = open(filename,'r')
2675                l = fp.readline().strip()
2676                if 'Name' not in l:
2677                    fp.close()
2678                    G2frame.ErrorDialog('Read Error',
2679                        'File '+filename+' does not start with Name\nFirst line ='
2680                        +l+'\ninvalid file',parent=G2frame)
2681                    return
2682                name = l.split(':')[1].strip()
2683                trans = fp.readline().strip().split(':')[1].split()
2684                vecMag = [float(i) for i in trans]
2685                ntrans = len(trans)
2686                vecs = [[] for i in range(ntrans)]
2687                types = []
2688                l = fp.readline().strip()
2689                while l:
2690                    nums = l.strip().split()
2691                    types.append(nums.pop(0))
2692                    t = types[-1]
2693                    if t not in AtInfo:
2694                        Info = G2elem.GetAtomInfo(t)
2695                        AtInfo[t] = [Info['Drad'],Info['Color']]
2696                    for i in range(ntrans):
2697                        vecs[i].append([float(nums.pop(0)) for j in range(3)])
2698                    l = fp.readline().strip()
2699                fp.close()
2700            else:
2701                return       
2702        finally:
2703            dlg.Destroy()
2704        natoms = len(types)
2705        vecs = [np.array(vecs[i]) for i in range(ntrans)]
2706        rbid = ran.randint(0,sys.maxsize)
2707        data['Vector'][rbid] = {'RBname':name,'VectMag':vecMag,
2708                'rbXYZ':np.zeros((natoms,3)),
2709                'rbRef':[0,1,2,False],'VectRef':ntrans*[False],
2710                'rbTypes':types,
2711                'rbVect':vecs,'useCount':0}
2712        data['RBIds']['Vector'].append(rbid)
2713        UpdateVectorRB()
2714       
2715    def AddResidueRB(event):
2716        global resRBsel
2717        AtInfo = data['Residue']['AtInfo']
2718        macro = getMacroFile('rigid body')
2719        if not macro:
2720            return
2721        macStr = macro.readline()
2722        while macStr:
2723            items = macStr.split()
2724            if 'I' == items[0]:
2725                resRBsel = ran.randint(0,sys.maxsize)
2726                rbName = items[1]
2727                rbTypes = []
2728                rbXYZ = []
2729                rbSeq = []
2730                atNames = []
2731                nAtms,nSeq,nOrig,mRef,nRef = [int(items[i]) for i in [2,3,4,5,6]]
2732                for iAtm in range(nAtms):
2733                    macStr = macro.readline().split()
2734                    atName = macStr[0]
2735                    atType = macStr[1]
2736                    atNames.append(atName)
2737                    rbXYZ.append([float(macStr[i]) for i in [2,3,4]])
2738                    rbTypes.append(atType)
2739                    if atType not in AtInfo:
2740                        Info = G2elem.GetAtomInfo(atType)
2741                        AtInfo[atType] = [Info['Drad'],Info['Color']]
2742                rbXYZ = np.array(rbXYZ)-np.array(rbXYZ[nOrig-1])
2743                for iSeq in range(nSeq):
2744                    macStr = macro.readline().split()
2745                    mSeq = int(macStr[0])
2746                    for jSeq in range(mSeq):
2747                        macStr = macro.readline().split()
2748                        iBeg = int(macStr[0])-1
2749                        iFin = int(macStr[1])-1
2750                        angle = 0.0
2751                        nMove = int(macStr[2])
2752                        iMove = [int(macStr[i])-1 for i in range(3,nMove+3)]
2753                        rbSeq.append([iBeg,iFin,angle,iMove])
2754                data['Residue'][resRBsel] = {'RBname':rbName,'rbXYZ':rbXYZ,'rbTypes':rbTypes,
2755                    'atNames':atNames,'rbRef':[nOrig-1,mRef-1,nRef-1,True],'rbSeq':rbSeq,
2756                    'SelSeq':[0,0],'useCount':0}
2757                data['RBIds']['Residue'].append(resRBsel)
2758                print ('Rigid body '+rbName+' added')
2759            macStr = macro.readline()
2760        macro.close()
2761        UpdateResidueRB()
2762       
2763    def ImportResidueRB():
2764        global resRBsel
2765        AtInfo = data['Residue']['AtInfo']
2766        text,ext = getTextFile()
2767        if not text:
2768            return
2769        resRBsel = ran.randint(0,sys.maxsize)
2770        rbTypes = []
2771        rbXYZ = []
2772        atNames = []
2773        txtStr = text.readline()
2774        if 'xyz' in ext:
2775            txtStr = text.readline()
2776            txtStr = text.readline()
2777        elif 'mol2' in ext:
2778            while 'ATOM' not in txtStr:
2779                txtStr = text.readline()
2780            txtStr = text.readline()
2781        elif 'pdb' in ext:
2782            while 'ATOM' not in txtStr[:6] and 'HETATM' not in txtStr[:6]:
2783                txtStr = text.readline()
2784        items = txtStr.split()
2785        nat = 1
2786        while len(items):
2787            if 'txt' in ext:
2788                atName = items[0]
2789                atType = items[1]
2790                rbXYZ.append([float(items[i]) for i in [2,3,4]])
2791            elif 'xyz' in ext:
2792                atType = items[0]
2793                rbXYZ.append([float(items[i]) for i in [1,2,3]])
2794                atName = '%s%d'%(atType,nat)
2795            elif 'mol2' in ext:
2796                atType = items[1]
2797                atName = items[1]+items[0]
2798                rbXYZ.append([float(items[i]) for i in [2,3,4]])
2799            elif 'pdb' in ext:
2800                atType = items[-1]
2801                if not items[2][-1].isnumeric():
2802                    atName = '%s%d'%(items[2],nat)
2803                else:
2804                    atName = '5s'%items[2]
2805                xyz = txtStr[30:55].split()                   
2806                rbXYZ.append([float(x) for x in xyz])
2807            atNames.append(atName)
2808            rbTypes.append(atType)
2809            if atType not in AtInfo:
2810                Info = G2elem.GetAtomInfo(atType)
2811                AtInfo[atType] = [Info['Drad'],Info['Color']]
2812            txtStr = text.readline()
2813            if 'mol2' in ext and 'BOND' in txtStr:
2814                break
2815            if 'pdb' in ext and ('ATOM' not in txtStr[:6] and 'HETATM' not in txtStr[:6]):
2816                break
2817            items = txtStr.split()
2818            nat += 1
2819        if len(atNames) < 3:
2820            G2G.G2MessageBox(G2frame,'Not enough atoms in rigid body; must be 3 or more')
2821        else:
2822            rbXYZ = np.array(rbXYZ)-np.array(rbXYZ[0])
2823            Xxyz = rbXYZ[1]
2824            X = Xxyz/np.sqrt(np.sum(Xxyz**2))
2825            Yxyz = rbXYZ[2]
2826            Y = Yxyz/np.sqrt(np.sum(Yxyz**2))
2827            Mat = G2mth.getRBTransMat(X,Y)
2828            rbXYZ = np.inner(Mat,rbXYZ).T
2829            data['Residue'][resRBsel] = {'RBname':'UNKRB','rbXYZ':rbXYZ,'rbTypes':rbTypes,
2830                'atNames':atNames,'rbRef':[0,1,2,False],'rbSeq':[],'SelSeq':[0,0],'useCount':0}
2831            data['RBIds']['Residue'].append(resRBsel)
2832            print ('Rigid body UNKRB added')
2833        text.close()
2834        UpdateResidueRB()
2835       
2836    def SaveResidueRB():
2837        global resRBsel
2838        pth = G2G.GetExportPath(G2frame)
2839        dlg = wx.FileDialog(G2frame, 'Choose PDB file for Atom XYZ', pth, '', 
2840            'PDB files (*.pdb)|*.pdb',wx.FD_SAVE|wx.FD_OVERWRITE_PROMPT)
2841        try:
2842            if dlg.ShowModal() == wx.ID_OK:
2843                filename = dlg.GetPath()
2844                filename = os.path.splitext(filename)[0]+'.pdb'  # make extension .pdb
2845                File = open(filename,'w')       
2846                rbData =  data['Residue'][resRBsel]
2847                for iat,xyz in enumerate(rbData['rbXYZ']):
2848                    File.write('ATOM %6d  %-4s%3s     1    %8.3f%8.3f%8.3f  1.00  0.00          %2s\n'%(
2849                        iat,rbData['atNames'][iat],rbData['RBname'][:3],xyz[0],xyz[1],xyz[2],rbData['rbTypes'][iat]))
2850                File.close()
2851                print ('Atom XYZ saved to: '+filename)
2852        finally:
2853            dlg.Destroy()
2854       
2855       
2856    def FindNeighbors(Orig,XYZ,atTypes,atNames,AtInfo):
2857        Radii = []
2858        for Atype in atTypes:
2859            Radii.append(AtInfo[Atype][0])
2860        Radii = np.array(Radii)
2861        Neigh = []
2862        Dx = XYZ-XYZ[Orig]
2863        dist = np.sqrt(np.sum(Dx**2,axis=1))
2864        sumR = Radii[Orig]+Radii
2865        IndB = ma.nonzero(ma.masked_greater(dist-0.85*sumR,0.))
2866        for j in IndB[0]:
2867            if j != Orig and atTypes[j] != 'H':
2868                Neigh.append(atNames[j])
2869        return Neigh
2870       
2871    def FindAllNeighbors(XYZ,atTypes,atNames,AtInfo):
2872        NeighDict = {}
2873        for iat,xyz in enumerate(atNames):
2874            NeighDict[atNames[iat]] = FindNeighbors(iat,XYZ,atTypes,atNames,AtInfo)
2875        return NeighDict
2876       
2877    def FindRiding(Orig,Pivot,NeighDict):
2878        riding = [Orig,Pivot]
2879        iAdd = 1
2880        new = True
2881        while new:
2882            newAtms = NeighDict[riding[iAdd]]
2883            for At in newAtms:
2884                new = False
2885                if At not in riding:
2886                    riding.append(At)
2887                    new = True
2888            iAdd += 1
2889            if iAdd < len(riding):
2890                new = True
2891        return riding[2:]
2892                       
2893    def OnDefineTorsSeq(event):
2894        global resRBsel
2895        rbData = data['Residue'][resRBsel]
2896        if not len(rbData):
2897            return
2898        atNames = rbData['atNames']
2899        AtInfo = data['Residue']['AtInfo']
2900        atTypes = rbData['rbTypes']
2901        XYZ = rbData['rbXYZ']
2902        neighDict = FindAllNeighbors(XYZ,atTypes,atNames,AtInfo)
2903        TargList = []           
2904        dlg = wx.SingleChoiceDialog(G2frame,'Select origin atom for torsion sequence','Origin atom',rbData['atNames'])
2905        if dlg.ShowModal() == wx.ID_OK:
2906            Orig = dlg.GetSelection()
2907            TargList = neighDict[atNames[Orig]]
2908        dlg.Destroy()
2909        if not len(TargList):
2910            return
2911        dlg = wx.SingleChoiceDialog(G2frame,'Select pivot atom for torsion sequence','Pivot atom',TargList)
2912        if dlg.ShowModal() == wx.ID_OK:
2913            Piv = atNames.index(TargList[dlg.GetSelection()])
2914            riding = FindRiding(atNames[Orig],atNames[Piv],neighDict)
2915            Riding = []
2916            for atm in riding:
2917                Riding.append(atNames.index(atm))
2918            rbData['rbSeq'].append([Orig,Piv,0.0,Riding])           
2919        dlg.Destroy()
2920        UpdateResidueRB()
2921
2922    def UpdateVectorRB(Scroll=0):
2923        AtInfo = data['Vector']['AtInfo']
2924        refChoice = {}
2925        if 'DELETED' in str(G2frame.GetStatusBar()):   #seems to be no other way to do this (wx bug)
2926            if GSASIIpath.GetConfigValue('debug'):
2927                print ('DBG_wx error: Rigid Body/Status not cleanly deleted after Refine')
2928            return
2929        SetStatusLine(' You may use e.g. "c60" or "s60" for a vector entry')
2930        def rbNameSizer(rbid,rbData):
2931
2932            def OnRBName(event):
2933                event.Skip()
2934                Obj = event.GetEventObject()
2935                rbData['RBname'] = Obj.GetValue()
2936               
2937            def OnDelRB(event):
2938                Obj = event.GetEventObject()
2939                rbid = Indx[Obj.GetId()]
2940                if rbid in data['Vector']:
2941                    del data['Vector'][rbid]
2942                    data['RBIds']['Vector'].remove(rbid)
2943                    rbData['useCount'] -= 1
2944                wx.CallAfter(UpdateVectorRB)
2945            def OnPlotRB(event):
2946                Obj = event.GetEventObject()
2947                G2plt.PlotRigidBody(G2frame,'Vector',AtInfo,rbData,plotDefaults)
2948           
2949            nameSizer = wx.BoxSizer(wx.HORIZONTAL)
2950            nameSizer.Add(wx.StaticText(VectorRBDisplay,-1,'Rigid body name: '),
2951                0,wx.ALIGN_CENTER_VERTICAL)
2952            RBname = wx.TextCtrl(VectorRBDisplay,-1,rbData['RBname'])
2953            Indx[RBname.GetId()] = rbid
2954            RBname.Bind(wx.EVT_TEXT_ENTER,OnRBName)
2955            RBname.Bind(wx.EVT_KILL_FOCUS,OnRBName)
2956            nameSizer.Add(RBname,0,wx.ALIGN_CENTER_VERTICAL)
2957            nameSizer.Add((5,0),)
2958            plotRB =  wx.Button(VectorRBDisplay,wx.ID_ANY,'Plot',
2959                                style=wx.BU_EXACTFIT)
2960            plotRB.Bind(wx.EVT_BUTTON, OnPlotRB)
2961            Indx[plotRB.GetId()] = rbid
2962            nameSizer.Add(plotRB,0,wx.ALIGN_CENTER_VERTICAL)
2963            nameSizer.Add((5,0),)
2964            if not rbData['useCount']:
2965                delRB = wx.Button(VectorRBDisplay,wx.ID_ANY,"Delete",
2966                                style=wx.BU_EXACTFIT)
2967                delRB.Bind(wx.EVT_BUTTON, OnDelRB)
2968                Indx[delRB.GetId()] = rbid
2969                nameSizer.Add(delRB,0,wx.ALIGN_CENTER_VERTICAL)
2970            return nameSizer
2971           
2972        def rbRefAtmSizer(rbid,rbData):
2973           
2974            def OnRefSel(event):
2975                Obj = event.GetEventObject()
2976                iref = Indx[Obj.GetId()]
2977                sel = Obj.GetValue()
2978                rbData['rbRef'][iref] = atNames.index(sel)
2979                FillRefChoice(rbid,rbData)
2980           
2981            refAtmSizer = wx.BoxSizer(wx.HORIZONTAL)
2982            atNames = [name+str(i) for i,name in enumerate(rbData['rbTypes'])]
2983            rbRef = rbData.get('rbRef',[0,1,2,False])
2984            rbData['rbRef'] = rbRef
2985            if rbData['useCount']:
2986                refAtmSizer.Add(wx.StaticText(VectorRBDisplay,-1,
2987                    'Orientation reference atoms A-B-C: %s, %s, %s'%(atNames[rbRef[0]], \
2988                     atNames[rbRef[1]],atNames[rbRef[2]])),0)
2989            else:
2990                refAtmSizer.Add(wx.StaticText(VectorRBDisplay,-1,
2991                    'Orientation reference atoms A-B-C: '),0,wx.ALIGN_CENTER_VERTICAL)
2992                for i in range(3):
2993                    choices = [atNames[j] for j in refChoice[rbid][i]]
2994                    refSel = wx.ComboBox(VectorRBDisplay,-1,value='',
2995                        choices=choices,style=wx.CB_READONLY|wx.CB_DROPDOWN)
2996                    refSel.SetValue(atNames[rbRef[i]])
2997                    refSel.Bind(wx.EVT_COMBOBOX, OnRefSel)
2998                    Indx[refSel.GetId()] = i
2999                    refAtmSizer.Add(refSel,0,wx.ALIGN_CENTER_VERTICAL)
3000            return refAtmSizer
3001                       
3002        def rbVectMag(rbid,imag,rbData):
3003           
3004            def OnRBVectorMag(event):
3005                event.Skip()
3006                Obj = event.GetEventObject()
3007                rbid,imag = Indx[Obj.GetId()]
3008                try:
3009                    val = float(Obj.GetValue())
3010                    if val <= 0.:
3011                        raise ValueError
3012                    rbData['VectMag'][imag] = val
3013                except ValueError:
3014                    pass
3015                Obj.SetValue('%8.4f'%(val))
3016                wx.CallAfter(UpdateVectorRB,VectorRB.GetScrollPos(wx.VERTICAL))
3017                G2plt.PlotRigidBody(G2frame,'Vector',AtInfo,data['Vector'][rbid],plotDefaults)
3018               
3019            def OnRBVectorRef(event):
3020                Obj = event.GetEventObject()
3021                rbid,imag = Indx[Obj.GetId()]
3022                rbData['VectRef'][imag] = Obj.GetValue()
3023                       
3024            magSizer = wx.BoxSizer(wx.HORIZONTAL)
3025            magSizer.Add(wx.StaticText(VectorRBDisplay,-1,'Translation magnitude: '),
3026                0,wx.ALIGN_CENTER_VERTICAL)
3027            magValue = wx.TextCtrl(VectorRBDisplay,-1,'%8.4f'%(rbData['VectMag'][imag]))
3028            Indx[magValue.GetId()] = [rbid,imag]
3029            magValue.Bind(wx.EVT_TEXT_ENTER,OnRBVectorMag)
3030            magValue.Bind(wx.EVT_KILL_FOCUS,OnRBVectorMag)
3031            magSizer.Add(magValue,0,wx.ALIGN_CENTER_VERTICAL)
3032            magSizer.Add((5,0),)
3033            magref = wx.CheckBox(VectorRBDisplay,-1,label=' Refine?') 
3034            magref.SetValue(rbData['VectRef'][imag])
3035            magref.Bind(wx.EVT_CHECKBOX,OnRBVectorRef)
3036            Indx[magref.GetId()] = [rbid,imag]
3037            magSizer.Add(magref,0,wx.ALIGN_CENTER_VERTICAL)
3038            return magSizer
3039           
3040        def rbVectors(rbid,imag,mag,XYZ,rbData):
3041
3042            def TypeSelect(event):
3043                AtInfo = data['Vector']['AtInfo']
3044                r,c = event.GetRow(),event.GetCol()
3045                if vecGrid.GetColLabelValue(c) == 'Type':
3046                    PE = G2elemGUI.PickElement(G2frame,oneOnly=True)
3047                    if PE.ShowModal() == wx.ID_OK:
3048                        if PE.Elem != 'None':
3049                            El = PE.Elem.strip().lower().capitalize()
3050                            if El not in AtInfo:
3051                                Info = G2elem.GetAtomInfo(El)
3052                                AtInfo[El] = [Info['Drad'],Info['Color']]
3053                            rbData['rbTypes'][r] = El
3054                            vecGrid.SetCellValue(r,c,El)
3055                    PE.Destroy()
3056                wx.CallAfter(UpdateVectorRB,VectorRB.GetScrollPos(wx.VERTICAL))
3057
3058            def ChangeCell(event):
3059                r,c =  event.GetRow(),event.GetCol()
3060                if r >= 0 and (0 <= c < 3):
3061                    try:
3062                        val = float(vecGrid.GetCellValue(r,c))
3063                        rbData['rbVect'][imag][r][c] = val
3064                    except ValueError:
3065                        pass
3066                G2plt.PlotRigidBody(G2frame,'Vector',AtInfo,data['Vector'][rbid],plotDefaults)
3067                wx.CallAfter(UpdateVectorRB,VectorRB.GetScrollPos(wx.VERTICAL))
3068
3069            vecSizer = wx.BoxSizer()
3070            Types = 3*[wg.GRID_VALUE_FLOAT+':10,5',]+[wg.GRID_VALUE_STRING,]+3*[wg.GRID_VALUE_FLOAT+':10,5',]
3071            colLabels = ['Vector x','Vector y','Vector z','Type','Cart x','Cart y','Cart z']
3072            table = []
3073            rowLabels = []
3074            for ivec,xyz in enumerate(rbData['rbVect'][imag]):
3075                table.append(list(xyz)+[rbData['rbTypes'][ivec],]+list(XYZ[ivec]))
3076                rowLabels.append(str(ivec))
3077            vecTable = G2G.Table(table,rowLabels=rowLabels,colLabels=colLabels,types=Types)
3078            vecGrid = G2G.GSGrid(VectorRBDisplay)
3079            vecGrid.SetTable(vecTable, True)
3080            if 'phoenix' in wx.version():
3081                vecGrid.Bind(wg.EVT_GRID_CELL_CHANGED, ChangeCell)
3082            else:
3083                vecGrid.Bind(wg.EVT_GRID_CELL_CHANGE, ChangeCell)
3084            if not imag:
3085                vecGrid.Bind(wg.EVT_GRID_CELL_LEFT_DCLICK, TypeSelect)
3086            attr = wx.grid.GridCellAttr()
3087            attr.IncRef()
3088            attr.SetEditor(G2G.GridFractionEditor(vecGrid))
3089            for c in range(3):
3090                attr.IncRef()
3091                vecGrid.SetColAttr(c, attr)
3092            for row in range(vecTable.GetNumberRows()):
3093                if imag:
3094                    vecGrid.SetCellStyle(row,3,VERY_LIGHT_GREY,True)                   
3095                for col in [4,5,6]:
3096                    vecGrid.SetCellStyle(row,col,VERY_LIGHT_GREY,True)
3097#            vecGrid.SetScrollRate(0,0)
3098            vecGrid.AutoSizeColumns(False)
3099            vecSizer.Add(vecGrid)
3100            return vecSizer
3101       
3102        def FillRefChoice(rbid,rbData):
3103            choiceIds = [i for i in range(len(rbData['rbTypes']))]
3104           
3105            rbRef = rbData.get('rbRef',[-1,-1,-1,False])
3106            for i in range(3):
3107                choiceIds.remove(rbRef[i])
3108            refChoice[rbid] = [choiceIds[:],choiceIds[:],choiceIds[:]]
3109            for i in range(3):
3110                refChoice[rbid][i].append(rbRef[i])
3111                refChoice[rbid][i].sort()     
3112           
3113        if VectorRB.GetSizer(): VectorRB.GetSizer().Clear(True)
3114        VectorRBSizer = wx.BoxSizer(wx.VERTICAL)
3115        first = True
3116        for rbid in data['RBIds']['Vector']:
3117            if rbid != 'AtInfo':
3118                rbData = data['Vector'][rbid]
3119                FillRefChoice(rbid,rbData)
3120                if not first:
3121                    G2G.HorizontalLine(VectorRBSizer,VectorRBDisplay)
3122                VectorRBSizer.Add(rbNameSizer(rbid,rbData),0)
3123                VectorRBSizer.Add(rbRefAtmSizer(rbid,rbData),0)
3124                XYZ = np.array([[0.,0.,0.] for Ty in rbData['rbTypes']])
3125                for imag,mag in enumerate(rbData['VectMag']):
3126                    XYZ += mag*rbData['rbVect'][imag]
3127                    VectorRBSizer.Add(rbVectMag(rbid,imag,rbData),0)
3128                    VectorRBSizer.Add(rbVectors(rbid,imag,mag,XYZ,rbData),0)
3129                VectorRBSizer.Add((5,5),0)
3130                data['Vector'][rbid]['rbXYZ'] = XYZ       
3131                first = False
3132        VectorRBSizer.Layout()   
3133        VectorRBDisplay.SetSizer(VectorRBSizer,True)
3134        Size = VectorRBSizer.GetMinSize()
3135        Size[0] += 40
3136        Size[1] = max(Size[1],450) + 20
3137        VectorRBDisplay.SetSize(Size)
3138        VectorRB.SetScrollbars(10,10,Size[0]/10-4,Size[1]/10-1)
3139        VectorRB.Scroll(0,Scroll)
3140       
3141    def UpdateResidueRB():
3142        '''Draw the contents of the Residue Rigid Body tab for Rigid Bodies tree entry
3143        '''
3144        global resRBsel
3145        def rbNameSizer(rbid,rbData):
3146
3147            def OnDelRB(event):
3148                Obj = event.GetEventObject()
3149                rbid = Indx[Obj.GetId()]
3150                if rbid in data['Residue']: 
3151                    del data['Residue'][rbid]
3152                    data['RBIds']['Residue'].remove(rbid)
3153                wx.CallAfter(UpdateResidueRB)
3154               
3155            def OnStripH(event):
3156                Obj = event.GetEventObject()
3157                rbid = Indx[Obj.GetId()]
3158                if rbid in data['Residue']:
3159                    newNames = []
3160                    newTypes = []
3161                    newXYZ = []
3162                    for i,atype in enumerate(rbData['rbTypes']):
3163                        if atype != 'H':
3164                            newNames.append(rbData['atNames'][i])
3165                            newTypes.append(rbData['rbTypes'][i])
3166                            newXYZ.append(rbData['rbXYZ'][i])
3167                    rbData['atNames'] = newNames
3168                    rbData['rbTypes'] = newTypes
3169                    rbData['rbXYZ'] = newXYZ
3170                G2plt.PlotRigidBody(G2frame,'Residue',AtInfo,rbData,plotDefaults)
3171                wx.CallAfter(UpdateResidueRB)
3172                   
3173            def OnPlotRB(event):
3174                Obj = event.GetEventObject()
3175                G2plt.PlotRigidBody(G2frame,'Residue',AtInfo,rbData,plotDefaults)
3176               
3177            # start of rbNameSizer
3178            nameSizer = wx.BoxSizer(wx.HORIZONTAL)
3179            nameSizer.Add(wx.StaticText(ResidueRBDisplay,-1,'Residue name: '),
3180                0,wx.ALIGN_CENTER_VERTICAL)
3181            nameSizer.Add(G2G.ValidatedTxtCtrl(ResidueRBDisplay,rbData,'RBname'),0,WACV)
3182            nameSizer.Add((5,0),)
3183            plotRB =  wx.Button(ResidueRBDisplay,wx.ID_ANY,'Plot',style=wx.BU_EXACTFIT)
3184            plotRB.Bind(wx.EVT_BUTTON, OnPlotRB)
3185            Indx[plotRB.GetId()] = rbid
3186            nameSizer.Add(plotRB,0,wx.ALIGN_CENTER_VERTICAL)
3187            nameSizer.Add((5,0),)
3188            if not rbData['useCount']:
3189                delRB = wx.Button(ResidueRBDisplay,wx.ID_ANY,"Delete",
3190                                style=wx.BU_EXACTFIT)
3191                delRB.Bind(wx.EVT_BUTTON, OnDelRB)
3192                Indx[delRB.GetId()] = rbid
3193                nameSizer.Add(delRB,0,wx.ALIGN_CENTER_VERTICAL)
3194                if 'H'  in rbData['rbTypes']:
3195                    stripH = wx.Button(ResidueRBDisplay,wx.ID_ANY,
3196                                "Strip H-atoms",
3197                                style=wx.BU_EXACTFIT)
3198                    stripH.Bind(wx.EVT_BUTTON, OnStripH)
3199                    Indx[stripH.GetId()] = rbid
3200                    nameSizer.Add(stripH,0,wx.ALIGN_CENTER_VERTICAL)
3201            return nameSizer
3202           
3203        def rbResidues(rbid,rbData):
3204           
3205            def TypeSelect(event):
3206                AtInfo = data['Residue']['AtInfo']
3207                r,c = event.GetRow(),event.GetCol()
3208                if resGrid.GetColLabelValue(c) == 'Type':
3209                    PE = G2elemGUI.PickElement(G2frame,oneOnly=True)
3210                    if PE.ShowModal() == wx.ID_OK:
3211                        if PE.Elem != 'None':
3212                            El = PE.Elem.strip().lower().capitalize()
3213                            if El not in AtInfo:
3214                                Info = G2elem.GetAtomInfo(El)
3215                                AtInfo[El] = [Info['Drad'],Info['Color']]
3216                            rbData['rbTypes'][r] = El
3217                            resGrid.SetCellValue(r,c,El)
3218                    PE.Destroy()
3219
3220            def ChangeCell(event):
3221                r,c =  event.GetRow(),event.GetCol()
3222                if c == 0:
3223                    rbData['atNames'][r] = resGrid.GetCellValue(r,c)
3224                if r >= 0 and (2 <= c <= 4):
3225                    try:
3226                        val = float(resGrid.GetCellValue(r,c))
3227                        rbData['rbXYZ'][r][c-2] = val
3228                    except ValueError:
3229                        pass
3230                       
3231            def OnRefSel(event):
3232                Obj = event.GetEventObject()
3233                iref,res,jref = Indx[Obj.GetId()]
3234                sel = Obj.GetValue()
3235                ind = atNames.index(sel)
3236                if rbData['rbTypes'][ind] == 'H':
3237                    G2G.G2MessageBox(G2frame,'You should not select an H-atom for rigid body orientation')
3238                rbData['rbRef'][iref] = ind
3239                FillRefChoice(rbid,rbData)
3240                for i,ref in enumerate(RefObjs[jref]):
3241                    ref.SetItems([atNames[j] for j in refChoice[rbid][i]])
3242                    ref.SetValue(atNames[rbData['rbRef'][i]])                   
3243                rbXYZ = rbData['rbXYZ']
3244                if not iref:     #origin change
3245                    rbXYZ -= rbXYZ[ind]
3246                #TODO - transform all atom XYZ by axis choices
3247                Xxyz = rbXYZ[rbData['rbRef'][1]]
3248                X = Xxyz/np.sqrt(np.sum(Xxyz**2))
3249                Yxyz = rbXYZ[rbData['rbRef'][2]]
3250                Y = Yxyz/np.sqrt(np.sum(Yxyz**2))
3251                Mat = G2mth.getRBTransMat(X,Y)
3252                rbXYZ = np.inner(Mat,rbXYZ).T
3253                rbData['rbXYZ'] = rbXYZ
3254                res.ClearSelection()
3255                resTable = res.GetTable()
3256                for r in range(res.GetNumberRows()):
3257                    row = resTable.GetRowValues(r)
3258                    row[2:4] = rbXYZ[r]
3259                    resTable.SetRowValues(r,row)
3260                res.ForceRefresh()
3261                G2plt.PlotRigidBody(G2frame,'Residue',AtInfo,rbData,plotDefaults)
3262               
3263            Types = 2*[wg.GRID_VALUE_STRING,]+3*[wg.GRID_VALUE_FLOAT+':10,5',]
3264            colLabels = ['Name','Type','Cart x','Cart y','Cart z']
3265            table = []
3266            rowLabels = []
3267            for ivec,xyz in enumerate(rbData['rbXYZ']):
3268                table.append([rbData['atNames'][ivec],]+[rbData['rbTypes'][ivec],]+list(xyz))
3269                rowLabels.append(str(ivec))
3270            vecTable = G2G.Table(table,rowLabels=rowLabels,colLabels=colLabels,types=Types)
3271            resGrid = G2G.GSGrid(ResidueRBDisplay)
3272            Indx[resGrid.GetId()] = rbid
3273            resList.append(resGrid)
3274            resGrid.SetTable(vecTable, True)
3275            if 'phoenix' in wx.version():
3276                resGrid.Bind(wg.EVT_GRID_CELL_CHANGED, ChangeCell)
3277            else:
3278                resGrid.Bind(wg.EVT_GRID_CELL_CHANGE, ChangeCell)
3279            resGrid.Bind(wg.EVT_GRID_CELL_LEFT_DCLICK, TypeSelect)
3280            for c in range(2,5):
3281                attr = wx.grid.GridCellAttr()
3282                attr.IncRef()
3283                attr.SetEditor(G2G.GridFractionEditor(resGrid))
3284                resGrid.SetColAttr(c, attr)
3285            for row in range(vecTable.GetNumberRows()):
3286                resGrid.SetReadOnly(row,1,True)
3287            resGrid.AutoSizeColumns(False)
3288            vecSizer = wx.BoxSizer()
3289            vecSizer.Add(resGrid)
3290           
3291            refAtmSizer = wx.BoxSizer(wx.HORIZONTAL)
3292            atNames = rbData['atNames']
3293            rbRef = rbData['rbRef']
3294            if rbData['rbRef'][3] or rbData['useCount']:
3295                refAtmSizer.Add(wx.StaticText(ResidueRBDisplay,-1,
3296                    'Orientation reference non-H atoms A-B-C: %s, %s, %s'%(atNames[rbRef[0]], \
3297                     atNames[rbRef[1]],atNames[rbRef[2]])),0)
3298            else:
3299                refAtmSizer.Add(wx.StaticText(ResidueRBDisplay,-1,
3300                    'Orientation reference non-H atoms A-B-C: '),0,wx.ALIGN_CENTER_VERTICAL)
3301                refObj = [0,0,0]
3302                for i in range(3):
3303                    choices = [atNames[j] for j in refChoice[rbid][i]]
3304                    refSel = wx.ComboBox(ResidueRBDisplay,-1,value='',
3305                        choices=choices,style=wx.CB_READONLY|wx.CB_DROPDOWN)
3306                    refSel.SetValue(atNames[rbRef[i]])
3307                    refSel.Bind(wx.EVT_COMBOBOX, OnRefSel)
3308                    Indx[refSel.GetId()] = [i,resGrid,len(RefObjs)]
3309                    refObj[i] = refSel
3310                    refAtmSizer.Add(refSel,0,wx.ALIGN_CENTER_VERTICAL)
3311                RefObjs.append(refObj)
3312           
3313            mainSizer = wx.BoxSizer(wx.VERTICAL)
3314            mainSizer.Add(refAtmSizer)
3315            mainSizer.Add(vecSizer)
3316            return mainSizer
3317           
3318        def Add2SeqSizer(seqSizer,angSlide,rbid,iSeq,Seq,atNames):
3319           
3320            def ChangeAngle(event):
3321                event.Skip()
3322                Obj = event.GetEventObject()
3323                rbid,Seq = Indx[Obj.GetId()][:2]
3324                val = Seq[2]
3325                try:
3326                    val = float(Obj.GetValue())
3327                    Seq[2] = val
3328                except ValueError:
3329                    pass
3330                Obj.SetValue('%8.2f'%(val))
3331                G2plt.PlotRigidBody(G2frame,'Residue',AtInfo,data['Residue'][rbid],plotDefaults)
3332               
3333            def OnRadBtn(event):
3334                Obj = event.GetEventObject()
3335                Seq,iSeq,angId = Indx[Obj.GetId()]
3336                data['Residue'][rbid]['SelSeq'] = [iSeq,angId]
3337                angSlide.SetValue(int(100*Seq[2]))
3338               
3339            def OnDelBtn(event):
3340                Obj = event.GetEventObject()
3341                rbid,Seq = Indx[Obj.GetId()]
3342                data['Residue'][rbid]['rbSeq'].remove(Seq)       
3343                wx.CallAfter(UpdateResidueRB)
3344           
3345            iBeg,iFin,angle,iMove = Seq
3346            ang = wx.TextCtrl(ResidueRBDisplay,wx.ID_ANY,
3347                    '%8.2f'%(angle),size=(70,-1),style=wx.TE_PROCESS_ENTER)
3348            if not iSeq:
3349                radBt = wx.RadioButton(ResidueRBDisplay,wx.ID_ANY,
3350                                           '',style=wx.RB_GROUP)
3351                data['Residue'][rbid]['SelSeq'] = [iSeq,ang.GetId()]
3352                radBt.SetValue(True)
3353            else:
3354                radBt = wx.RadioButton(ResidueRBDisplay,wx.ID_ANY,'')
3355            radBt.Bind(wx.EVT_RADIOBUTTON,OnRadBtn)                   
3356            seqSizer.Add(radBt)
3357            delBt =  wx.Button(ResidueRBDisplay,wx.ID_ANY,'Del',
3358                                style=wx.BU_EXACTFIT)
3359            delBt.Bind(wx.EVT_BUTTON,OnDelBtn)
3360            seqSizer.Add(delBt)
3361            bond = wx.StaticText(ResidueRBDisplay,wx.ID_ANY,
3362                        '%s %s'%(atNames[iBeg],atNames[iFin]),size=(50,20))
3363            seqSizer.Add(bond,0,wx.ALIGN_CENTER_VERTICAL)
3364            Indx[radBt.GetId()] = [Seq,iSeq,ang.GetId()]
3365            Indx[delBt.GetId()] = [rbid,Seq]
3366            Indx[ang.GetId()] = [rbid,Seq,ang]
3367            ang.Bind(wx.EVT_TEXT_ENTER,ChangeAngle)
3368            ang.Bind(wx.EVT_KILL_FOCUS,ChangeAngle)
3369            seqSizer.Add(ang,0,wx.ALIGN_CENTER_VERTICAL)
3370            atms = ''
3371            for i in iMove:   
3372                atms += ' %s,'%(atNames[i])
3373            moves = wx.StaticText(ResidueRBDisplay,wx.ID_ANY,
3374                            atms[:-1],size=(200,20))
3375            seqSizer.Add(moves,1,wx.ALIGN_CENTER_VERTICAL|wx.EXPAND|wx.RIGHT)
3376            return seqSizer
3377           
3378        def SlideSizer():
3379           
3380            def OnSlider(event):
3381                Obj = event.GetEventObject()
3382                rbData = Indx[Obj.GetId()]
3383                iSeq,angId = rbData['SelSeq']
3384                val = float(Obj.GetValue())/100.
3385                rbData['rbSeq'][iSeq][2] = val
3386                Indx[angId][2].SetValue('%8.2f'%(val))
3387                G2plt.PlotRigidBody(G2frame,'Residue',AtInfo,rbData,plotDefaults)
3388           
3389            slideSizer = wx.BoxSizer(wx.HORIZONTAL)
3390            slideSizer.Add(wx.StaticText(ResidueRBDisplay,-1,'Selected torsion angle:'),0)
3391            iSeq,angId = rbData['SelSeq']
3392            angSlide = wx.Slider(ResidueRBDisplay,-1,
3393                int(100*rbData['rbSeq'][iSeq][2]),0,36000,size=(200,20),
3394                style=wx.SL_HORIZONTAL)
3395            angSlide.Bind(wx.EVT_SLIDER, OnSlider)
3396            Indx[angSlide.GetId()] = rbData
3397            slideSizer.Add(angSlide,0)           
3398            return slideSizer,angSlide
3399           
3400        def FillRefChoice(rbid,rbData):
3401            choiceIds = [i for i in range(len(rbData['atNames']))]
3402            for seq in rbData['rbSeq']:
3403                for i in seq[3]:
3404                    try:
3405                        choiceIds.remove(i)
3406                    except ValueError:
3407                        pass
3408            rbRef = rbData['rbRef']
3409            for i in range(3):
3410                try:
3411                    choiceIds.remove(rbRef[i])
3412                except ValueError:
3413                    pass
3414            refChoice[rbid] = [choiceIds[:],choiceIds[:],choiceIds[:]]
3415            for i in range(3):
3416                refChoice[rbid][i].append(rbRef[i])
3417                refChoice[rbid][i].sort()
3418               
3419        def OnSelect(event):
3420            global resRBsel
3421            sel = select.GetSelection()
3422            if sel == 0: return # 1st entry is blank
3423            rbname = rbchoice[sel-1]
3424            resRBsel = RBnames[rbname]
3425            wx.CallLater(100,UpdateResidueRB)
3426           
3427        #----- beginning of UpdateResidueRB -----------------------------------------------
3428        AtInfo = data['Residue']['AtInfo']
3429        refChoice = {}
3430        RefObjs = []
3431
3432        GS = ResidueRBDisplay.GetSizer()
3433        if GS: 
3434            try:        #get around a c++ error in wx 4.0; doing is again seems to be OK
3435                GS.Clear(True)
3436            except:
3437                GS.Clear(True)
3438       
3439        RBnames = {}
3440        for rbid in data['RBIds']['Residue']:
3441            RBnames.update({data['Residue'][rbid]['RBname']:rbid,})
3442        if not RBnames:
3443            return
3444        rbchoice = list(RBnames.keys())
3445        if GS: 
3446            ResidueRBSizer = GS
3447        else:
3448            ResidueRBSizer = wx.BoxSizer(wx.VERTICAL)
3449        if len(RBnames) > 1:
3450            selSizer = wx.BoxSizer(wx.HORIZONTAL)
3451            selSizer.Add(wx.StaticText(ResidueRBDisplay,
3452                                label=' Select residue to view:'),0)
3453            rbchoice.sort()
3454            select = wx.ComboBox(ResidueRBDisplay,choices=['']+rbchoice)
3455            select.Bind(wx.EVT_COMBOBOX,OnSelect)
3456            selSizer.Add(select,0)
3457            ResidueRBSizer.Add(selSizer,0)
3458        if resRBsel not in data['RBIds']['Residue']:
3459            resRBsel = RBnames[rbchoice[0]]
3460        rbData = data['Residue'][resRBsel]
3461        FillRefChoice(resRBsel,rbData)
3462        ResidueRBSizer.Add(rbNameSizer(resRBsel,rbData),0)
3463        ResidueRBSizer.Add(rbResidues(resRBsel,rbData),0)
3464        if len(rbData['rbSeq']):
3465            ResidueRBSizer.Add((-1,15),0)
3466            slideSizer,angSlide = SlideSizer()
3467            seqSizer = wx.FlexGridSizer(0,5,4,8)
3468            for lbl in 'Sel','','Bond','Angle','Riding atoms':
3469                seqSizer.Add(wx.StaticText(ResidueRBDisplay,wx.ID_ANY,lbl))
3470            ResidueRBSizer.Add(seqSizer)
3471#            for iSeq,Seq in enumerate(rbData['rbSeq']):
3472#                ResidueRBSizer.Add(SeqSizer(angSlide,resRBsel,iSeq,Seq,rbData['atNames']))
3473            for iSeq,Seq in enumerate(rbData['rbSeq']):
3474                Add2SeqSizer(seqSizer,angSlide,resRBsel,iSeq,Seq,rbData['atNames'])
3475            ResidueRBSizer.Add(slideSizer,)
3476
3477        ResidueRBSizer.Add((5,25),)
3478        ResidueRBSizer.Layout()   
3479        ResidueRBDisplay.SetSizer(ResidueRBSizer,True)
3480        ResidueRBDisplay.SetAutoLayout(True)
3481        Size = ResidueRBSizer.GetMinSize()
3482        ResidueRBDisplay.SetSize(Size)
3483       
3484        Size[0] += 40
3485        Size[1] = max(Size[1],450) + 20
3486        ResidueRB.SetSize(Size)
3487        ResidueRB.SetScrollbars(10,10,Size[0]/10-4,Size[1]/10-1)
3488        G2frame.dataWindow.SendSizeEvent()
3489       
3490        ResidueRBDisplay.Show()
3491       
3492    def SetStatusLine(text):
3493        G2frame.GetStatusBar().SetStatusText(text,1)                                     
3494
3495    G2gd.SetDataMenuBar(G2frame,G2frame.dataWindow.RigidBodyMenu)
3496    SetStatusLine('')
3497    UpdateVectorRB()
3498    G2frame.rbBook.Bind(wx.aui.EVT_AUINOTEBOOK_PAGE_CHANGED, OnPageChanged)
3499    wx.CallAfter(OnPageChanged,None)
3500
3501def ShowIsoDistortCalc(G2frame,phase=None):
3502    '''Compute the ISODISTORT mode values from the current coordinates.
3503    Called in response to the (Phase/Atoms tab) AtomCompute or
3504    Constraints/Edit Constr. "Show ISODISTORT modes" menu item, which
3505    should be enabled only when Phase['ISODISTORT'] is defined.
3506    '''
3507    def _onClose(event):
3508        dlg.EndModal(wx.ID_CANCEL)
3509    def fmtHelp(item,fullname):
3510        helptext = "A new variable"
3511        if item[-3]:
3512            helptext += " named "+str(item[-3])
3513        helptext += " is a linear combination of the following parameters:\n"
3514        first = True
3515        for term in item[:-3]:
3516            line = ''
3517            var = str(term[1])
3518            m = term[0]
3519            if first:
3520                first = False
3521                line += ' = '
3522            else:
3523                if m >= 0:
3524                    line += ' + '
3525                else:
3526                    line += ' - '
3527                m = abs(m)
3528            line += '%.3f*%s '%(m,var)
3529            varMean = G2obj.fmtVarDescr(var)
3530            helptext += "\n" + line + " ("+ varMean + ")"
3531        helptext += '\n\nISODISTORT full name: '+str(fullname)
3532        return helptext
3533
3534    Histograms,Phases = G2frame.GetUsedHistogramsAndPhasesfromTree() # init for constraint
3535    isophases = [p for p in Phases if 'ISODISTORT' in Phases[p]]
3536   
3537    if not isophases:
3538        G2frame.ErrorDialog('no ISODISTORT phases',
3539                            'Unexpected error: no ISODISTORT phases')
3540        return
3541    if phase and phase not in isophases:
3542        G2frame.ErrorDialog('not ISODISTORT phase',
3543                            'Unexpected error: selected phase is not an ISODISTORT phase')
3544        print('Unexpected error: selected phase is not an ISODISTORT phase',
3545                  phase,isophases)
3546    elif not phase and len(isophases) == 1:
3547        phase = isophases[0]
3548    elif not phase:
3549        dlg = wx.SingleChoiceDialog(G2frame,'Select phase from ISODISTORT phases',
3550                                        'Select Phase',isophases)
3551        if dlg.ShowModal() == wx.ID_OK:
3552            sel = dlg.GetSelection()
3553            phase = isophases[sel]
3554        else:
3555            return
3556    # if len(data.get('Histograms',[])) == 0:
3557    #     G2frame.ErrorDialog(
3558    #         'No data',
3559    #         'Sorry, this computation requires that a histogram first be added to the phase'
3560    #         )
3561    #     return
3562   
3563    covdata = G2frame.GPXtree.GetItemPyData(
3564        G2gd.GetGPXtreeItemId(G2frame,G2frame.root,'Covariance'))
3565    # make a lookup table for named NewVar Phase constraints
3566    sub = G2gd.GetGPXtreeItemId(G2frame,G2frame.root,'Constraints') 
3567    Constraints = G2frame.GPXtree.GetItemPyData(sub)
3568    constDict = {}
3569    for c in Constraints['Phase']:
3570        if c[-1] != 'f' or not c[-3]: continue
3571        constDict[str(c[-3])] = c
3572
3573    parmDict,varyList = G2frame.MakeLSParmDict()
3574
3575    dlg = wx.Dialog(G2frame,wx.ID_ANY,'ISODISTORT mode values',#size=(630,400),
3576                       style=wx.DEFAULT_DIALOG_STYLE|wx.RESIZE_BORDER)
3577    mainSizer = wx.BoxSizer(wx.VERTICAL)
3578    if 'ISODISTORT' not in Phases[phase]:
3579        G2frame.ErrorDialog('not ISODISTORT phase',
3580                            'Unexpected error: selected phase is not an ISODISTORT phase')
3581        return
3582    data = Phases[phase]
3583    ISO = data['ISODISTORT']
3584    mainSizer.Add(wx.StaticText(dlg,wx.ID_ANY,
3585                                'ISODISTORT mode computation for cordinates in phase '+
3586                                str(data['General'].get('Name'))))
3587    aSizer = wx.BoxSizer(wx.HORIZONTAL)
3588    panel1 = wxscroll.ScrolledPanel(
3589        dlg, wx.ID_ANY,#size=(100,200),
3590        style = wx.TAB_TRAVERSAL|wx.SUNKEN_BORDER)
3591    subSizer1 = wx.FlexGridSizer(cols=2,hgap=5,vgap=2)
3592    panel2 = wxscroll.ScrolledPanel(
3593        dlg, wx.ID_ANY,#size=(100,200),
3594        style = wx.TAB_TRAVERSAL|wx.SUNKEN_BORDER)
3595    subSizer2 = wx.FlexGridSizer(cols=3,hgap=5,vgap=2)
3596    subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,'Parameter name  '))
3597    subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,' value'),0,wx.ALIGN_RIGHT)
3598    subSizer2.Add((-1,-1))
3599    subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,'Mode name  '))
3600    subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,' value'),0,wx.ALIGN_RIGHT)
3601    # ISODISTORT displacive modes
3602    if 'G2VarList' in ISO:
3603        deltaList = []
3604        for gv,Ilbl in zip(ISO['G2VarList'],ISO['IsoVarList']):
3605            dvar = gv.varname()
3606            var = dvar.replace('::dA','::A')
3607            albl = Ilbl[:Ilbl.rfind('_')]
3608            v = Ilbl[Ilbl.rfind('_')+1:]
3609            pval = ISO['ParentStructure'][albl][['dx','dy','dz'].index(v)]
3610            if var in parmDict:
3611                cval = parmDict[var][0]
3612            else:
3613                dlg.EndModal(wx.ID_CANCEL)
3614                G2frame.ErrorDialog('Atom not found',"No value found for parameter "+str(var))
3615                return
3616            deltaList.append(cval-pval)
3617        modeVals = np.inner(ISO['Var2ModeMatrix'],deltaList)
3618        for lbl,xyz,var,val,norm,G2mode in zip(
3619                ISO['IsoVarList'],deltaList,
3620                ISO['IsoModeList'],modeVals,ISO['NormList'],ISO['G2ModeList'] ):
3621            #GSASIIpath.IPyBreak()
3622            if str(G2mode) in constDict:
3623                ch = G2G.HelpButton(panel2,fmtHelp(constDict[str(G2mode)],var))
3624                subSizer2.Add(ch,0,wx.LEFT|wx.RIGHT|wx.ALIGN_CENTER_VERTICAL|wx.ALIGN_CENTER,1)
3625            else:
3626                subSizer2.Add((-1,-1))
3627            subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,str(lbl)))
3628            try:
3629                value = G2py3.FormatSigFigs(xyz)
3630            except TypeError:
3631                value = str(xyz)           
3632            subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,value),0,wx.ALIGN_RIGHT)
3633            subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,str(var)))
3634            try:
3635                value = G2py3.FormatSigFigs(val/norm)
3636                if 'varyList' in covdata:
3637                    if str(G2mode) in covdata['varyList']:
3638                        sig = covdata['sig'][covdata['varyList'].index(str(G2mode))]
3639                        value = G2mth.ValEsd(val/norm,sig/norm)
3640            except TypeError:
3641                value = '?'
3642            subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,value),0,wx.ALIGN_RIGHT)
3643            #GSASIIpath.IPyBreak()
3644    # ISODISTORT occupancy modes
3645    if 'G2OccVarList' in ISO:
3646        deltaList = []
3647        for gv,Ilbl in zip(ISO['G2OccVarList'],ISO['OccVarList']):
3648            var = gv.varname()
3649            albl = Ilbl[:Ilbl.rfind('_')]
3650            pval = ISO['BaseOcc'][albl]
3651            if var in parmDict:
3652                cval = parmDict[var][0]
3653            else:
3654                dlg.EndModal(wx.ID_CANCEL)
3655                G2frame.ErrorDialog('Atom not found',"No value found for parameter "+str(var))
3656                return
3657            deltaList.append(cval-pval)
3658        modeVals = np.inner(ISO['Var2OccMatrix'],deltaList)
3659        for lbl,delocc,var,val,norm,G2mode in zip(
3660                ISO['OccVarList'],deltaList,
3661                ISO['OccModeList'],modeVals,ISO['OccNormList'],ISO['G2OccModeList']):
3662            if str(G2mode) in constDict:
3663                ch = G2G.HelpButton(panel2,fmtHelp(constDict[str(G2mode)],var))
3664                subSizer2.Add(ch,0,wx.LEFT|wx.RIGHT|wx.ALIGN_CENTER_VERTICAL|wx.ALIGN_CENTER,1)
3665            else:
3666                subSizer2.Add((-1,-1))
3667            subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,str(lbl)))
3668            try:
3669                value = G2py3.FormatSigFigs(delocc)
3670            except TypeError:
3671                value = str(delocc)
3672            subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,value),0,wx.ALIGN_RIGHT)
3673            #subSizer.Add((10,-1))
3674            subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,str(var)))
3675            try:
3676                value = G2py3.FormatSigFigs(val/norm)
3677                if 'varyList' in covdata:
3678                    if str(G2mode) in covdata['varyList']:
3679                        sig = covdata['sig'][covdata['varyList'].index(str(G2mode))]
3680                        value = G2mth.ValEsd(val/norm,sig/norm)
3681            except TypeError:
3682                value = '?'
3683            subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,value),0,wx.ALIGN_RIGHT)
3684
3685    # finish up ScrolledPanel
3686    panel1.SetSizer(subSizer1)
3687    panel2.SetSizer(subSizer2)
3688    panel1.SetAutoLayout(1)
3689    panel1.SetupScrolling()
3690    panel2.SetAutoLayout(1)
3691    panel2.SetupScrolling()
3692    # Allow window to be enlarged but not made smaller
3693    dlg.SetSizer(mainSizer)
3694    w1,l1 = subSizer1.GetSize()
3695    w2,l2 = subSizer2.GetSize()
3696    panel1.SetMinSize((w1+10,200))
3697    panel2.SetMinSize((w2+20,200))
3698    aSizer.Add(panel1,1, wx.ALL|wx.EXPAND,1)
3699    aSizer.Add(panel2,2, wx.ALL|wx.EXPAND,1)
3700    mainSizer.Add(aSizer,1, wx.ALL|wx.EXPAND,1)
3701
3702    # make OK button
3703    btnsizer = wx.BoxSizer(wx.HORIZONTAL)
3704    btn = wx.Button(dlg, wx.ID_CLOSE) 
3705    btn.Bind(wx.EVT_BUTTON,_onClose)
3706    btnsizer.Add(btn)
3707    mainSizer.Add(btnsizer, 0, wx.ALIGN_CENTER|wx.ALL, 5)
3708
3709    mainSizer.Fit(dlg)
3710    dlg.SetMinSize(dlg.GetSize())
3711    dlg.ShowModal()
3712    dlg.Destroy()
Note: See TracBrowser for help on using the repository browser.