source: trunk/GSASIIconstrGUI.py @ 3790

Last change on this file since 3790 was 3790, checked in by toby, 3 years ago

make sure off-diagonal terms that must compute as 0 are fixed

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 103.0 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASIIconstrGUI - constraint GUI routines
3########### SVN repository information ###################
4# $Date: 2019-01-20 03:42:46 +0000 (Sun, 20 Jan 2019) $
5# $Author: toby $
6# $Revision: 3790 $
7# $URL: trunk/GSASIIconstrGUI.py $
8# $Id: GSASIIconstrGUI.py 3790 2019-01-20 03:42:46Z toby $
9########### SVN repository information ###################
10'''
11*GSASIIconstrGUI: Constraint GUI routines*
12------------------------------------------
13
14Used to define constraints and rigid bodies.
15
16'''
17from __future__ import division, print_function
18import sys
19import wx
20import wx.grid as wg
21import random as ran
22import numpy as np
23import numpy.ma as ma
24import numpy.linalg as nl
25import os.path
26import GSASIIpath
27GSASIIpath.SetVersionNumber("$Revision: 3790 $")
28import GSASIIElem as G2elem
29import GSASIIElemGUI as G2elemGUI
30import GSASIIstrIO as G2stIO
31import GSASIImapvars as G2mv
32import GSASIImath as G2mth
33import GSASIIlattice as G2lat
34import GSASIIdataGUI as G2gd
35import GSASIIctrlGUI as G2G
36import GSASIIplot as G2plt
37import GSASIIobj as G2obj
38import GSASIIspc as G2spc
39VERY_LIGHT_GREY = wx.Colour(235,235,235)
40
41class ConstraintDialog(wx.Dialog):
42    '''Window to edit Constraint values
43    '''
44    def __init__(self,parent,title,text,data,separator='*',varname="",varyflag=False):
45        wx.Dialog.__init__(self,parent,-1,'Edit '+title, 
46            pos=wx.DefaultPosition,style=wx.DEFAULT_DIALOG_STYLE)
47        self.data = data[:]
48        self.newvar = [varname,varyflag]
49        panel = wx.Panel(self)
50        mainSizer = wx.BoxSizer(wx.VERTICAL)
51        topLabl = wx.StaticText(panel,-1,text)
52        mainSizer.Add((10,10),1)
53        mainSizer.Add(topLabl,0,wx.ALIGN_CENTER_VERTICAL|wx.LEFT,10)
54        mainSizer.Add((10,10),1)
55        dataGridSizer = wx.FlexGridSizer(cols=3,hgap=2,vgap=2)
56        self.OkBtn = wx.Button(panel,wx.ID_OK)
57        for id in range(len(self.data)):
58            lbl1 = lbl = str(self.data[id][1])
59            if lbl[-1] != '=': lbl1 = lbl + ' ' + separator + ' '
60            name = wx.StaticText(panel,wx.ID_ANY,lbl1,
61                                 style=wx.ALIGN_RIGHT)
62            scale = G2G.ValidatedTxtCtrl(panel,self.data[id],0,
63                                          typeHint=float,
64                                          OKcontrol=self.DisableOK)
65            dataGridSizer.Add(name,0,wx.LEFT|wx.RIGHT|wx.ALIGN_CENTER_VERTICAL,5)
66            dataGridSizer.Add(scale,0,wx.RIGHT,3)
67            if ':' in lbl:
68                dataGridSizer.Add(
69                    wx.StaticText(panel,-1,G2obj.fmtVarDescr(lbl)),
70                    0,wx.RIGHT|wx.ALIGN_CENTER_VERTICAL,3)
71            else:
72                dataGridSizer.Add((-1,-1))
73        if title == 'New Variable':
74            name = wx.StaticText(panel,wx.ID_ANY,"New variable's\nname (optional)",
75                                 style=wx.ALIGN_CENTER)
76            scale = G2G.ValidatedTxtCtrl(panel,self.newvar,0,
77                                          typeHint=str,notBlank=False)
78            dataGridSizer.Add(name,0,wx.LEFT|wx.RIGHT|wx.ALIGN_CENTER_VERTICAL,5)
79            dataGridSizer.Add(scale,0,wx.RIGHT|wx.ALIGN_CENTER_VERTICAL,3)
80            self.refine = wx.CheckBox(panel,label='Refine?')
81            self.refine.SetValue(self.newvar[1]==True)
82            self.refine.Bind(wx.EVT_CHECKBOX, self.OnCheckBox)
83            dataGridSizer.Add(self.refine,0,wx.RIGHT|wx.ALIGN_CENTER_VERTICAL,3)
84           
85        mainSizer.Add(dataGridSizer,0,wx.EXPAND)
86        self.OkBtn.Bind(wx.EVT_BUTTON, self.OnOk)
87        self.OkBtn.SetDefault()
88        cancelBtn = wx.Button(panel,wx.ID_CANCEL)
89        cancelBtn.Bind(wx.EVT_BUTTON, self.OnCancel)
90        btnSizer = wx.BoxSizer(wx.HORIZONTAL)
91        btnSizer.Add((20,20),1)
92        btnSizer.Add(self.OkBtn)
93        btnSizer.Add((20,20),1)
94        btnSizer.Add(cancelBtn)
95        btnSizer.Add((20,20),1)
96
97        mainSizer.Add(btnSizer,0,wx.EXPAND|wx.BOTTOM|wx.TOP, 10)
98        panel.SetSizer(mainSizer)
99        panel.Fit()
100        self.Fit()
101        self.CenterOnParent()
102       
103    def DisableOK(self,setting):
104        if setting:
105            self.OkBtn.Enable()
106        else:
107            self.OkBtn.Disable()
108
109    def OnCheckBox(self,event):
110        self.newvar[1] = self.refine.GetValue()
111
112    def OnOk(self,event):
113        parent = self.GetParent()
114        parent.Raise()
115        self.EndModal(wx.ID_OK)             
116
117    def OnCancel(self,event):
118        parent = self.GetParent()
119        parent.Raise()
120        self.EndModal(wx.ID_CANCEL)             
121
122    def GetData(self):
123        return self.data
124       
125################################################################################
126#####  Constraints
127################################################################################           
128       
129def UpdateConstraints(G2frame,data):
130    '''Called when Constraints tree item is selected.
131    Displays the constraints in the data window
132    '''
133       
134    def FindEquivVarb(name,nameList):
135        'Creates a list of variables appropriate to constrain with name'
136        outList = []
137        #phlist = []
138        items = name.split(':')
139        namelist = [items[2],]
140        if 'dA' in name:
141            namelist = ['dAx','dAy','dAz']
142        elif 'AU' in name:
143            namelist = ['AUiso','AU11','AU22','AU33','AU12','AU13','AU23']
144        elif 'AM' in name:
145            namelist = ['AMx','AMy','AMz']
146        elif items[-1] in ['A0','A1','A2','A3','A4','A5']:
147            namelist = ['A0','A1','A2','A3','A4','A5']
148        elif items[-1] in ['D11','D22','D33','D12','D13','D23']:
149            namelist = ['D11','D22','D33','D12','D13','D23']
150        elif 'Tm' in name:
151            namelist = ['Tmin','Tmax']
152        elif 'RB' in name:
153            rbfx = 'RB'+items[2][2]
154            if 'T' in name and 'Tr' not in name:
155                namelist = [rbfx+'T11',rbfx+'T22',rbfx+'T33',rbfx+'T12',rbfx+'T13',rbfx+'T23']
156            if 'L' in name:
157                namelist = [rbfx+'L11',rbfx+'L22',rbfx+'L33',rbfx+'L12',rbfx+'L13',rbfx+'L23']
158            if 'S' in name:
159                namelist = [rbfx+'S12',rbfx+'S13',rbfx+'S21',rbfx+'S23',rbfx+'S31',rbfx+'S32',rbfx+'SAA',rbfx+'SBB']
160            if 'U' in name:
161                namelist = [rbfx+'U',]
162
163        for item in nameList:
164            keys = item.split(':')
165            #if keys[0] not in phlist:
166            #    phlist.append(keys[0])
167            if items[1] == '*' and keys[2] in namelist: # wildcard -- select only sequential options
168                keys[1] = '*'
169                mitem = ':'.join(keys)
170                if mitem == name: continue
171                if mitem not in outList: outList.append(mitem)
172            elif keys[2] in namelist and item != name:
173                outList.append(item)
174        return outList
175       
176    def SelectVarbs(page,FrstVarb,varList,legend,constType):
177        '''Select variables used in constraints after one variable has
178        been selected. This routine determines the appropriate variables to be
179        used based on the one that has been selected and asks for more to be added.
180
181        It then creates the constraint and adds it to the constraints list.
182       
183        Called from OnAddEquivalence, OnAddFunction & OnAddConstraint (all but
184        OnAddHold)
185
186        :param list page: defines calling page -- type of variables to be used
187        :parm GSASIIobj.G2VarObj FrstVarb: reference to first selected variable
188        :param list varList: list of other appropriate variables to select from
189        :param str legend: header for selection dialog
190        :param str constType: type of constraint to be generated
191        :returns: a constraint, as defined in
192          :ref:`GSASIIobj <Constraint_definitions_table>`
193        '''
194        choices = [[i]+list(G2obj.VarDescr(i)) for i in varList]
195        meaning = G2obj.getDescr(FrstVarb.name)
196        if not meaning:
197            meaning = "(no definition found!)"
198        l = str(FrstVarb).split(':')
199        # make lists of phases & histograms to iterate over
200        phaselist = [l[0]]
201        if l[0]:
202            phaselbl = ['phase #'+l[0]]
203            if len(Phases) > 1:
204                phaselist += ['all'] 
205                phaselbl += ['all phases']
206        else:
207            phaselbl = ['']
208        histlist = [l[1]]
209        if l[1] == '*':
210            pass
211        elif l[1]:
212            histlbl = ['histogram #'+l[1]]
213            if len(Histograms) > 1:
214                histlist += ['all']
215                histlbl += ['all histograms']
216                typ = Histograms[G2obj.LookupHistName(l[1])[0]]['Instrument Parameters'][0]['Type'][1]
217                i = 0
218                for hist in Histograms:
219                    if Histograms[hist]['Instrument Parameters'][0]['Type'][1] == typ: i += 1
220                if i > 1:
221                    histlist += ['all='+typ]
222                    histlbl += ['all '+typ+' histograms']
223        else:
224            histlbl = ['']
225        # make a list of equivalent parameter names
226        nameList = [FrstVarb.name]
227        for var in varList:
228            nam = var.split(":")[2]
229            if nam not in nameList: nameList += [nam]
230        # add "wild-card" names to the list of variables
231        if l[1] == '*':
232            pass
233        elif page[1] == 'phs':
234            if 'RB' in FrstVarb.name:
235                pass
236            elif FrstVarb.atom is None:
237                for nam in nameList:
238                    for ph,plbl in zip(phaselist,phaselbl):
239                        if plbl: plbl = 'For ' + plbl
240                        var = ph+"::"+nam
241                        if var == str(FrstVarb) or var in varList: continue
242                        varList += [var]
243                        choices.append([var,plbl,meaning])
244            else:
245                for nam in nameList:
246                    for ph,plbl in zip(phaselist,phaselbl):
247                        if plbl: plbl = ' in ' + plbl
248                        for atype in ['']+TypeList:
249                            if atype:
250                                albl = "For "+atype+" atoms"
251                                akey = "all="+atype                       
252                            else:
253                                albl = "For all atoms"
254                                akey = "all"
255                            var = ph+"::"+nam+":"+akey
256                            if var == str(FrstVarb) or var in varList: continue
257                            varList += [var]
258                            choices.append([var,albl+plbl,meaning])
259        elif page[1] == 'hap':
260            if FrstVarb.name == "Scale":
261                meaning = "Phase fraction"
262            for nam in nameList:
263                for ph,plbl in zip(phaselist,phaselbl):
264                    if plbl: plbl = 'For ' + plbl
265                    for hst,hlbl in zip(histlist,histlbl):
266                        if hlbl:
267                            if plbl:
268                                hlbl = ' in ' + hlbl
269                            else:
270                                hlbl = 'For ' + hlbl                               
271                            var = ph+":"+hst+":"+nam
272                            if var == str(FrstVarb) or var in varList: continue
273                            varList += [var]
274                            choices.append([var,plbl+hlbl,meaning])
275        elif page[1] == 'hst':
276            if FrstVarb.name == "Scale":
277                meaning = "Scale factor"
278            for nam in nameList:
279                for hst,hlbl in zip(histlist,histlbl):
280                    if hlbl:
281                        hlbl = 'For ' + hlbl                               
282                        var = ":"+hst+":"+nam
283                        if var == str(FrstVarb) or var in varList: continue
284                        varList += [var]
285                        choices.append([var,hlbl,meaning])
286        elif page[1] == 'glb' or page[1] == 'sym':
287            pass
288        else:
289            raise Exception('Unknown constraint page '+ page[1])                   
290        if len(choices):
291            l1 = l2 = 1
292            for i1,i2,i3 in choices:
293                l1 = max(l1,len(i1))
294                l2 = max(l2,len(i2))
295            fmt = "{:"+str(l1)+"s} {:"+str(l2)+"s} {:s}"
296            atchoices = [fmt.format(*i1) for i1 in choices] # reformat list as str with columns
297            dlg = G2G.G2MultiChoiceDialog(
298                G2frame,legend,
299                'Constrain '+str(FrstVarb)+' with...',atchoices,
300                toggle=False,size=(625,400),monoFont=True)
301            dlg.CenterOnParent()
302            res = dlg.ShowModal()
303            Selections = dlg.GetSelections()[:]
304            dlg.Destroy()
305            if res != wx.ID_OK: return []
306            if len(Selections) == 0:
307                dlg = wx.MessageDialog(
308                    G2frame,
309                    'No variables were selected to include with '+str(FrstVarb),
310                    'No variables')
311                dlg.CenterOnParent()
312                dlg.ShowModal()
313                dlg.Destroy()
314                return []
315        else:
316            dlg = wx.MessageDialog(
317                G2frame,
318                'There are no appropriate variables to include with '+str(FrstVarb),
319                'No variables')
320            dlg.CenterOnParent()
321            dlg.ShowModal()
322            dlg.Destroy()
323            return []
324        # now process the variables provided by the user
325        varbs = [str(FrstVarb),] # list of selected variables
326        for sel in Selections:
327            var = varList[sel]
328            # phase(s) included
329            l = var.split(':')
330            if l[0] == "all":
331                phlist = [str(Phases[phase]['pId']) for phase in Phases]
332            else:
333                phlist = [l[0]]
334            # histogram(s) included
335            if l[1] == "all":
336                hstlist = [str(Histograms[hist]['hId']) for hist in Histograms]
337            elif '=' in l[1]:
338                htyp = l[1].split('=')[1]
339                hstlist = [str(Histograms[hist]['hId']) for hist in Histograms if 
340                           Histograms[hist]['Instrument Parameters'][0]['Type'][1] == htyp]
341            else:
342                hstlist = [l[1]]
343            if len(l) == 3:
344                for ph in phlist:
345                    for hst in hstlist:
346                        var = ph + ":" + hst + ":" + l[2]
347                        if var in varbs: continue
348                        varbs.append(var)
349            else: # constraints with atoms or rigid bodies
350                if len(l) == 5: # rigid body parameter
351                    var = ':'.join(l)
352                    if var in varbs: continue
353                    varbs.append(var)
354                elif l[3] == "all":
355                    for ph in phlist:
356                        key = G2obj.LookupPhaseName(ph)[0]
357                        for hst in hstlist: # should be blank
358                            for iatm,at in enumerate(Phases[key]['Atoms']):
359                                var = ph + ":" + hst + ":" + l[2] + ":" + str(iatm)
360                                if var in varbs: continue
361                                varbs.append(var)
362                elif '=' in l[3]:
363                    for ph in phlist:
364                        key = G2obj.LookupPhaseName(ph)[0]
365                        cx,ct,cs,cia = Phases[key]['General']['AtomPtrs']
366                        for hst in hstlist: # should be blank
367                            atyp = l[3].split('=')[1]
368                            for iatm,at in enumerate(Phases[key]['Atoms']):
369                                if at[ct] != atyp: continue
370                                var = ph + ":" + hst + ":" + l[2] + ":" + str(iatm)
371                                if var in varbs: continue
372                                varbs.append(var)
373                else:
374                    for ph in phlist:
375                        key = G2obj.LookupPhaseName(ph)[0]
376                        for hst in hstlist: # should be blank
377                            var = ph + ":" + hst + ":" + l[2] + ":" + l[3]
378                            if var in varbs: continue
379                            varbs.append(var)
380        if len(varbs) >= 1 or 'constraint' in constType:
381            constr = [[1.0,FrstVarb]]
382            for item in varbs[1:]:
383                constr += [[1.0,G2obj.G2VarObj(item)]]
384            if 'equivalence' in constType:
385                return [constr+[None,None,'e']]
386            elif 'function' in constType:
387                return [constr+[None,False,'f']]
388            elif 'constraint' in constType:
389                return [constr+[1.0,None,'c']]
390            else:
391                raise Exception('Unknown constraint type: '+str(constType))
392        else:
393            dlg = wx.MessageDialog(
394                G2frame,
395                'There are no selected variables to include with '+str(FrstVarb),
396                'No variables')
397            dlg.CenterOnParent()
398            dlg.ShowModal()
399            dlg.Destroy()
400        return []
401   
402    def ConstraintsLoad(data,newcons=[]):
403        '''Load all constraints. Constraints based on symmetry (etc.)
404        are generated by running :func:`GSASIIstrIO.GetPhaseData`.
405        '''
406        G2mv.InitVars()
407        #Find all constraints
408        constraintSet = []
409        for key in data:
410            if key.startswith('_'): continue
411            constraintSet += data[key]
412        if newcons:
413            constraintSet = constraintSet + newcons
414        constDictList,fixedList,ignored = G2stIO.ProcessConstraints(constraintSet)
415        # generate symmetry constraints to check for conflicts
416        rigidbodyDict = G2frame.GPXtree.GetItemPyData(   
417            G2gd.GetGPXtreeItemId(G2frame,G2frame.root,'Rigid bodies'))
418        rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
419        rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,Print=False)
420        Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,maxSSwave = G2stIO.GetPhaseData(
421            Phases,RestraintDict=None,rbIds=rbIds,Print=False) # generates atom symmetry constraints
422        return constDictList,phaseDict,fixedList
423           
424    def ConstraintsCheck(data,newcons=[]):
425        '''Load constraints & check them for errors. Since error checking
426        can cause changes in constraints in case of repairable conflicts
427        between equivalences, reload the constraints again after the check.
428        This could probably be done more effectively, but only reloading when
429        needed, but the reload should be quick.
430        '''
431        constDictList,phaseDict,fixedList = ConstraintsLoad(data,newcons)
432        msg = G2mv.EvaluateMultipliers(constDictList,phaseDict)
433        if msg:
434            return 'Unable to interpret multiplier(s): '+msg,''
435        res = G2mv.CheckConstraints('',constDictList,fixedList)
436        # reload constraints in case any were merged in MoveConfEquiv
437        ConstraintsLoad(data,newcons)
438        return res
439
440    def CheckAddedConstraint(newcons):
441        '''Check a new constraint that has just been input.
442        If there is an error display a message and discard the last entry
443
444        Since the varylist is not available, no warning messages
445        should be generated here
446
447        :returns: True if constraint should be added
448        '''
449       
450        errmsg,warnmsg = ConstraintsCheck(data,newcons)
451        if errmsg:
452            G2frame.ErrorDialog('Constraint Error',
453                'Error with newly added constraint:\n'+errmsg+
454                '\nIgnoring newly added constraint',parent=G2frame)
455            # reset error status
456            errmsg,warnmsg = ConstraintsCheck(data)
457            if errmsg:
458                print (errmsg)
459                print (G2mv.VarRemapShow([],True))
460            return False
461        elif warnmsg:
462            print ('Unexpected contraint warning:\n'+warnmsg)
463        return True
464
465    def CheckChangedConstraint():
466        '''Check all constraints after an edit has been made.
467        If there is an error display a message and reject the change.
468
469        Since the varylist is not available, no warning messages
470        should be generated.
471       
472        :returns: True if the edit should be retained
473        '''
474        errmsg,warnmsg = ConstraintsCheck(data)
475        if errmsg:
476            G2frame.ErrorDialog('Constraint Error',
477                'Error after editing constraint:\n'+errmsg+
478                '\nDiscarding last constraint edit',parent=G2frame)
479            # reset error status
480            errmsg,warnmsg = ConstraintsCheck(data)
481            if errmsg:
482                print (errmsg)
483                print (G2mv.VarRemapShow([],True))
484            return False
485        elif warnmsg:
486            print ('Unexpected contraint warning:\n'+warnmsg)
487        return True
488             
489    def PageSelection(page):
490        'Decode page reference'
491        if page[1] == "phs":
492            vartype = "phase"
493            varList = G2obj.removeNonRefined(phaseList)  # remove any non-refinable prms from list
494            constrDictEnt = 'Phase'
495        elif page[1] == "hap":
496            vartype = "Histogram*Phase"
497            varList = G2obj.removeNonRefined(hapList)  # remove any non-refinable prms from list
498            constrDictEnt = 'HAP'
499        elif page[1] == "hst":
500            vartype = "Histogram"
501            varList = G2obj.removeNonRefined(histList)  # remove any non-refinable prms from list
502            constrDictEnt = 'Hist'
503        elif page[1] == "glb":
504            vartype = "Global"
505            varList = G2obj.removeNonRefined(globalList)   # remove any non-refinable prms from list
506
507            constrDictEnt = 'Global'
508        elif page[1] == "sym":
509            return None,None,None
510        else:
511            raise Exception('Should not happen!')
512        return vartype,varList,constrDictEnt
513
514    def OnAddHold(event):
515        '''Create a new Hold constraint
516
517        Hold constraints allows the user to select one variable (the list of available
518        variables depends on which tab is currently active).
519        '''
520        page = G2frame.Page
521        vartype,varList,constrDictEnt = PageSelection(page)
522        if vartype is None: return
523        varList = G2obj.SortVariables(varList)
524        title1 = "Hold "+vartype+" variable"
525        if not varList:
526            G2frame.ErrorDialog('No variables','There are no variables of type '+vartype,
527                parent=G2frame)
528            return
529        l2 = l1 = 1
530        for i in varList:
531            l1 = max(l1,len(i))
532            loc,desc = G2obj.VarDescr(i)
533            l2 = max(l2,len(loc))
534        fmt = "{:"+str(l1)+"s} {:"+str(l2)+"s} {:s}"
535        varListlbl = [fmt.format(i,*G2obj.VarDescr(i)) for i in varList]
536        #varListlbl = ["("+i+") "+G2obj.fmtVarDescr(i) for i in varList]
537        legend = "Select variables to hold (Will not be varied, even if vary flag is set)"
538        dlg = G2G.G2MultiChoiceDialog(G2frame,
539            legend,title1,varListlbl,toggle=False,size=(625,400),monoFont=True)
540        dlg.CenterOnParent()
541        if dlg.ShowModal() == wx.ID_OK:
542            for sel in dlg.GetSelections():
543                Varb = varList[sel]
544                VarObj = G2obj.G2VarObj(Varb)
545                newcons = [[[0.0,VarObj],None,None,'h']]
546                if CheckAddedConstraint(newcons):
547                    data[constrDictEnt] += newcons
548        dlg.Destroy()
549        wx.CallAfter(OnPageChanged,None)
550       
551    def OnAddEquivalence(event):
552        '''add an Equivalence constraint'''
553        page = G2frame.Page
554        vartype,varList,constrDictEnt = PageSelection(page)
555        if vartype is None: return
556        title1 = "Create equivalence constraint between "+vartype+" variables"
557        title2 = "Select additional "+vartype+" variable(s) to be equivalent with "
558        if not varList:
559            G2frame.ErrorDialog('No variables','There are no variables of type '+vartype,
560                parent=G2frame)
561            return
562#        legend = "Select variables to make equivalent (only one of the variables will be varied when all are set to be varied)"
563        GetAddVars(page,title1,title2,varList,constrDictEnt,'equivalence')
564       
565    def OnAddAtomEquiv(event):
566        ''' Add equivalences between all parameters on atoms '''
567        page = G2frame.Page
568        vartype,varList,constrDictEnt = PageSelection(page)
569        if vartype is None: return
570        title1 = "Setup equivalent atom variables"
571        title2 = "Select additional atoms(s) to be equivalent with "
572        if not varList:
573            G2frame.ErrorDialog('No variables','There are no variables of type '+vartype,
574                parent=G2frame)
575            return
576#        legend = "Select atoms to make equivalent (only one of the atom variables will be varied when all are set to be varied)"
577        GetAddAtomVars(page,title1,title2,varList,constrDictEnt,'equivalence')
578       
579    def OnAddRiding(event):
580        ''' Add riding equivalences between all parameters on atoms '''
581        page = G2frame.Page
582        vartype,varList,constrDictEnt = PageSelection(page)
583        if vartype is None: return
584        title1 = "Setup riding atoms "
585        title2 = "Select additional atoms(s) to ride on "
586        if not varList:
587            G2frame.ErrorDialog('No variables','There are no variables of type '+vartype,
588                parent=G2frame)
589            return
590#        legend = "Select atoms to ride (only one of the atom variables will be varied when all are set to be varied)"
591        GetAddAtomVars(page,title1,title2,varList,constrDictEnt,'riding')
592   
593    def OnAddFunction(event):
594        '''add a Function (new variable) constraint'''
595        page = G2frame.Page
596        vartype,varList,constrDictEnt = PageSelection(page)
597        if vartype is None: return
598        title1 = "Setup new variable based on "+vartype+" variables"
599        title2 = "Include additional "+vartype+" variable(s) to be included with "
600        if not varList:
601            G2frame.ErrorDialog('No variables','There are no variables of type '+vartype,
602                parent=G2frame)
603            return
604#        legend = "Select variables to include in a new variable (the new variable will be varied when all included variables are varied)"
605        GetAddVars(page,title1,title2,varList,constrDictEnt,'function')
606                       
607    def OnAddConstraint(event):
608        '''add a constraint equation to the constraints list'''
609        page = G2frame.Page
610        vartype,varList,constrDictEnt = PageSelection(page)
611        if vartype is None: return
612        title1 = "Creating constraint on "+vartype+" variables"
613        title2 = "Select additional "+vartype+" variable(s) to include in constraint with "
614        if not varList:
615            G2frame.ErrorDialog('No variables','There are no variables of type '+vartype,
616                parent=G2frame)
617            return
618#        legend = "Select variables to include in a constraint equation (the values will be constrainted to equal a specified constant)"
619        GetAddVars(page,title1,title2,varList,constrDictEnt,'constraint')
620
621    def GetAddVars(page,title1,title2,varList,constrDictEnt,constType):
622        '''Get the variables to be added for OnAddEquivalence, OnAddFunction,
623        and OnAddConstraint. Then create and check the constraint.
624        '''
625        #varListlbl = ["("+i+") "+G2obj.fmtVarDescr(i) for i in varList]
626        if constType == 'equivalence':
627            omitVars = G2mv.GetDependentVars()
628        else:
629            omitVars = []
630        varList = G2obj.SortVariables([i for i in varList if i not in omitVars])
631        l2 = l1 = 1
632        for i in varList:
633            l1 = max(l1,len(i))
634            loc,desc = G2obj.VarDescr(i)
635            l2 = max(l2,len(loc))
636        fmt = "{:"+str(l1)+"s} {:"+str(l2)+"s} {:s}"
637        varListlbl = [fmt.format(i,*G2obj.VarDescr(i)) for i in varList]
638        dlg = G2G.G2SingleChoiceDialog(G2frame,'Select 1st variable:',
639            title1,varListlbl,monoFont=True,size=(625,400))
640        dlg.CenterOnParent()
641        if dlg.ShowModal() == wx.ID_OK:
642            if constType == 'equivalence':
643                omitVars = G2mv.GetDependentVars() + G2mv.GetIndependentVars()
644            sel = dlg.GetSelection()
645            FrstVarb = varList[sel]
646            VarObj = G2obj.G2VarObj(FrstVarb)
647            moreVarb = G2obj.SortVariables(FindEquivVarb(FrstVarb,[i for i in varList if i not in omitVars]))
648            newcons = SelectVarbs(page,VarObj,moreVarb,title2+FrstVarb,constType)
649            if len(newcons) > 0:
650                if CheckAddedConstraint(newcons):
651                    data[constrDictEnt] += newcons
652        dlg.Destroy()
653        wx.CallAfter(OnPageChanged,None)
654                       
655    def FindNeighbors(phase,FrstName,AtNames):
656        General = phase['General']
657        cx,ct,cs,cia = General['AtomPtrs']
658        Atoms = phase['Atoms']
659        atNames = [atom[ct-1] for atom in Atoms]
660        Cell = General['Cell'][1:7]
661        Amat,Bmat = G2lat.cell2AB(Cell)
662        atTypes = General['AtomTypes']
663        Radii = np.array(General['BondRadii'])
664        AtInfo = dict(zip(atTypes,Radii)) #or General['BondRadii']
665        Orig = atNames.index(FrstName.split()[1])
666        OType = Atoms[Orig][ct]
667        XYZ = G2mth.getAtomXYZ(Atoms,cx)       
668        Neigh = []
669        Dx = np.inner(Amat,XYZ-XYZ[Orig]).T
670        dist = np.sqrt(np.sum(Dx**2,axis=1))
671        sumR = AtInfo[OType]+0.5    #H-atoms only!
672        IndB = ma.nonzero(ma.masked_greater(dist-0.85*sumR,0.))
673        for j in IndB[0]:
674            if j != Orig:
675                Neigh.append(AtNames[j])
676        return Neigh
677       
678    def GetAddAtomVars(page,title1,title2,varList,constrDictEnt,constType):
679        '''Get the atom variables to be added for OnAddAtomEquiv. Then create and
680        check the constraints. Riding for H atoms only.
681        '''
682        Atoms = {G2obj.VarDescr(i)[0]:[] for i in varList if 'Atom' in G2obj.VarDescr(i)[0]}
683        for item in varList:
684            atName = G2obj.VarDescr(item)[0]
685            if atName in Atoms:
686                Atoms[atName].append(item)
687        AtNames = list(Atoms.keys())
688        AtNames.sort()
689        dlg = G2G.G2SingleChoiceDialog(G2frame,'Select 1st atom:',
690            title1,AtNames,monoFont=True,size=(625,400))
691        dlg.CenterOnParent()
692        FrstAtom = ''
693        if dlg.ShowModal() == wx.ID_OK:
694            sel = dlg.GetSelection()
695            FrstAtom = AtNames[sel]
696            if 'riding' in constType:
697                phaseName = (FrstAtom.split(' in ')[1]).strip()
698                phase = Phases[phaseName]
699                AtNames = FindNeighbors(phase,FrstAtom,AtNames)
700            else:
701                AtNames.remove(FrstAtom)
702        dlg.Destroy()
703        if FrstAtom == '':
704            print ('no atom selected')
705            return
706        dlg = G2G.G2MultiChoiceDialog(
707            G2frame,title2+FrstAtom,
708            'Constrain '+str(FrstAtom)+' with...',AtNames,
709            toggle=False,size=(625,400),monoFont=True)
710        if dlg.ShowModal() == wx.ID_OK:
711            Selections = dlg.GetSelections()[:]
712        else:
713            print ('no target atom selected')
714            dlg.Destroy()
715            return
716        dlg.Destroy()
717        for name in Atoms[FrstAtom]:
718            newcons = []
719            constr = []
720            if 'riding' in constType:
721                if 'AUiso' in name:
722                    constr = [[1.0,G2obj.G2VarObj(name)]]
723                elif 'AU11' in name:
724                    pass
725                elif 'AU' not in name:
726                    constr = [[1.0,G2obj.G2VarObj(name)]]
727            else:
728                constr = [[1.0,G2obj.G2VarObj(name)]]
729            pref = ':'+name.rsplit(':',1)[0].split(':',1)[1]    #get stuff between phase id & atom id
730            for sel in Selections:
731                name2 = Atoms[AtNames[sel]][0]
732                pid = name2.split(':',1)[0]                     #get phase id for 2nd atom
733                id = name2.rsplit(':',1)[-1]                    #get atom no. for 2nd atom
734                if 'riding' in constType:
735                    pref = pid+pref
736                    if 'AUiso' in pref:
737                        parts = pref.split('AUiso')
738                        constr += [[1.2,G2obj.G2VarObj('%s:%s'%(parts[0]+'AUiso',id))]]
739                    elif 'AU' not in pref:
740                        constr += [[1.0,G2obj.G2VarObj('%s:%s'%(pref,id))]]
741                else:
742                    constr += [[1.0,G2obj.G2VarObj('%s:%s'%(pid+pref,id))]]
743            if not constr:
744                continue
745            if 'frac' in pref and 'riding' not in constType:
746                newcons = [constr+[1.0,None,'c']]
747            else:
748                newcons = [constr+[None,None,'e']]
749            if len(newcons) > 0:
750                if CheckAddedConstraint(newcons):
751                    data[constrDictEnt] += newcons
752        wx.CallAfter(OnPageChanged,None)
753                       
754    def MakeConstraintsSizer(name,pageDisplay):
755        '''Creates a sizer displaying all of the constraints entered of
756        the specified type.
757
758        :param str name: the type of constraints to be displayed ('HAP',
759          'Hist', 'Phase', 'Global', 'Sym-Generated')
760        :param wx.Panel pageDisplay: parent panel for sizer
761        :returns: wx.Sizer created by method
762        '''
763        if name == 'Sym-Generated':         #show symmetry generated constraints
764            Sizer1 =  wx.BoxSizer(wx.VERTICAL)
765            Sizer1.Add(wx.StaticText(pageDisplay,wx.ID_ANY,
766                                    'Equivalences generated based on cell/space group input'))
767            Sizer1.Add((-1,5))
768            Sizer = wx.FlexGridSizer(0,2,0,0)
769            Sizer1.Add(Sizer)
770            for sym in G2mv.GetSymEquiv():
771                Sizer.Add(wx.StaticText(pageDisplay,wx.ID_ANY,'EQUIV'),
772                           0,wx.ALIGN_CENTER_VERTICAL|wx.ALIGN_CENTER|wx.RIGHT|wx.LEFT,3)
773                Sizer.Add(wx.StaticText(pageDisplay,wx.ID_ANY,sym))
774                Sizer.Add((-1,-1))
775                Sizer.Add((-1,2))
776            return Sizer1
777        constSizer = wx.FlexGridSizer(0,6,0,0)
778        maxlen = 50 # characters before wrapping a constraint
779        for Id,item in enumerate(data[name]):
780            refineflag = False
781            helptext = ""
782            eqString = ['',]
783            problemItem = False
784            for term in item[:-3]:
785                if str(term[1]) in G2mv.problemVars:
786                    problemItem = True
787            if item[-1] == 'h': # Hold on variable
788                constSizer.Add((-1,-1),0)              # blank space for edit button
789                typeString = 'FIXED'
790                var = str(item[0][1])
791                varMean = G2obj.fmtVarDescr(var)
792                eqString[-1] =  var +'   '
793                helptext = "Prevents variable:\n"+ var + " ("+ varMean + ")\nfrom being changed"
794            elif item[-1] == 'f' or item[-1] == 'e' or item[-1] == 'c': # not true on original-style (2011?) constraints
795                constEdit = wx.Button(pageDisplay,wx.ID_ANY,'Edit',style=wx.BU_EXACTFIT)
796                constEdit.Bind(wx.EVT_BUTTON,OnConstEdit)
797                constSizer.Add(constEdit,0,wx.LEFT|wx.RIGHT|wx.ALIGN_CENTER_VERTICAL|wx.ALIGN_CENTER,1)            # edit button
798                Indx[constEdit.GetId()] = [Id,name]
799                if item[-1] == 'f':
800                    helptext = "A new variable"
801                    if item[-3]:
802                        helptext += " named "+str(item[-3])
803                    helptext += " is created from a linear combination of the following variables:\n"
804                    for term in item[:-3]:
805                        var = str(term[1])
806                        if len(eqString[-1]) > maxlen:
807                            eqString.append(' ')
808                        m = term[0]
809                        if eqString[-1] != '':
810                            if m >= 0:
811                                eqString[-1] += ' + '
812                            else:
813                                eqString[-1] += ' - '
814                                m = abs(m)
815                        eqString[-1] += '{:g}*{:} '.format(m,var)
816                        varMean = G2obj.fmtVarDescr(var)
817                        helptext += "\n" + var + " ("+ varMean + ")"
818                    if '_Explain' in data:
819                        if data['_Explain'].get(item[-3]):
820                            helptext += '\n\n'
821                            helptext += data['_Explain'][item[-3]]
822                    # typeString = 'NEWVAR'
823                    # if item[-3]:
824                    #     eqString[-1] += ' = '+item[-3]
825                    # else:
826                    #     eqString[-1] += ' = New Variable'
827                    if item[-3]:
828                        typeString = item[-3] + ' = '
829                    else:
830                        typeString = 'New Variable = '
831                    #print 'refine',item[-2]
832                    refineflag = True
833                elif item[-1] == 'c':
834                    helptext = "The following variables constrained to equal a constant:"
835                    for term in item[:-3]:
836                        var = str(term[1])
837                        if len(eqString[-1]) > maxlen:
838                            eqString.append(' ')
839                        m = term[0]
840                        if eqString[-1] != '':
841                            if term[0] > 0:
842                                eqString[-1] += ' + '
843                            else:
844                                eqString[-1] += ' - '
845                                m = -term[0]
846                        eqString[-1] += '{:g}*{:} '.format(m,var)
847                        varMean = G2obj.fmtVarDescr(var)
848                        helptext += "\n" + var + " ("+ varMean + ")"
849                    typeString = 'CONST'
850                    eqString[-1] += ' = '+str(item[-3])
851                elif item[-1] == 'e':
852                    helptext = "The following variables are set to be equivalent, noting multipliers:"
853                    normval = item[:-3][1][0]
854                    for i,term in enumerate(item[:-3]):
855                        var = str(term[1])
856                        if term[0] == 0: term[0] = 1.0
857                        if len(eqString[-1]) > maxlen:
858                            eqString.append(' ')
859                        if i == 0: # move independent variable to end
860                            indepterm = term
861                            continue
862                        elif eqString[-1] != '':
863                            eqString[-1] += ' = '
864                        if normval/term[0] == 1:
865                            eqString[-1] += '{:}'.format(var)
866                        else:
867                            eqString[-1] += '{:g}*{:} '.format(normval/term[0],var)
868                        varMean = G2obj.fmtVarDescr(var)
869                        helptext += "\n" + var + " ("+ varMean + ")"
870                    if normval/indepterm[0] == 1:
871                        eqString[-1] += ' = {:} '.format(indepterm[1])
872                    else:
873                        eqString[-1] += ' = {:g}*{:} '.format(normval/indepterm[0],str(indepterm[1]))
874                    typeString = 'EQUIV'
875                else:
876                    print ('Unexpected constraint'+item)
877               
878            else:
879                print ('Removing old-style constraints')
880                data[name] = []
881                return constSizer
882            constDel = wx.Button(pageDisplay,wx.ID_ANY,'Delete',style=wx.BU_EXACTFIT)
883            constDel.Bind(wx.EVT_BUTTON,OnConstDel)
884            Indx[constDel.GetId()] = [Id,name]
885            constSizer.Add(constDel,0,wx.LEFT|wx.RIGHT|wx.ALIGN_CENTER_VERTICAL|wx.ALIGN_CENTER,1)             # delete button
886            if helptext:
887                ch = G2G.HelpButton(pageDisplay,helptext)
888                constSizer.Add(ch,0,wx.LEFT|wx.RIGHT|wx.ALIGN_CENTER_VERTICAL|wx.ALIGN_CENTER,1)
889            else:
890                constSizer.Add((-1,-1))
891            if refineflag:
892                ch = G2G.G2CheckBox(pageDisplay,'',item,-2)
893                constSizer.Add(ch,0,wx.LEFT|wx.RIGHT|wx.ALIGN_CENTER_VERTICAL|wx.ALIGN_CENTER,1)
894            else:
895                constSizer.Add((-1,-1))               
896            constSizer.Add(wx.StaticText(pageDisplay,wx.ID_ANY,typeString),
897                           0,wx.ALIGN_CENTER_VERTICAL|wx.ALIGN_CENTER|wx.RIGHT|wx.LEFT,3)
898            if problemItem: eqString[-1] += ' -- Conflict: see console'
899            if len(eqString) > 1:
900                Eq = wx.BoxSizer(wx.VERTICAL)
901                for s in eqString:
902                    line = wx.StaticText(pageDisplay,wx.ID_ANY,s)
903                    if problemItem: line.SetBackgroundColour(wx.YELLOW)
904                    Eq.Add(line,0,wx.ALIGN_CENTER_VERTICAL)
905                Eq.Add((-1,4))
906            else:
907                Eq = wx.StaticText(pageDisplay,wx.ID_ANY,eqString[0])
908                if problemItem: Eq.SetBackgroundColour(wx.YELLOW)
909            constSizer.Add(Eq,1,wx.ALIGN_CENTER_VERTICAL)
910        return constSizer
911               
912    def OnConstDel(event):
913        'Delete a constraint'
914        Obj = event.GetEventObject()
915        Id,name = Indx[Obj.GetId()]
916        del data[name][Id]
917        ConstraintsLoad(data)
918        wx.CallAfter(OnPageChanged,None)
919
920    def OnConstEdit(event):
921        '''Called to edit an individual contraint in response to a
922        click on its Edit button
923        '''
924        Obj = event.GetEventObject()
925        Id,name = Indx[Obj.GetId()]
926        if data[name][Id][-1] == 'f':
927            items = data[name][Id][:-3]
928            constType = 'New Variable'
929            if data[name][Id][-3]:
930                varname = data[name][Id][-3]
931            else:
932                varname = ""
933            lbl = 'Enter value for each term in constraint; sum = new variable'
934            dlg = ConstraintDialog(G2frame,constType,lbl,items,
935                                   varname=varname,varyflag=data[name][Id][-2])
936        elif data[name][Id][-1] == 'c':
937            items = data[name][Id][:-3]+[
938                [data[name][Id][-3],'fixed value =']]
939            constType = 'Constraint'
940            lbl = 'Edit value for each term in constant constraint sum'
941            dlg = ConstraintDialog(G2frame,constType,lbl,items)
942        elif data[name][Id][-1] == 'e':
943            items = data[name][Id][:-3]
944            constType = 'Equivalence'
945            lbl = 'The following terms are set to be equal:'
946            dlg = ConstraintDialog(G2frame,constType,lbl,items,'/')
947        else:
948            return
949        try:
950            prev = data[name][Id][:]
951            if dlg.ShowModal() == wx.ID_OK:
952                result = dlg.GetData()
953                for i in range(len(data[name][Id][:-3])):
954                    if type(data[name][Id][i]) is tuple: # fix non-mutable construct
955                        data[name][Id][i] = list(data[name][Id][i])
956                    data[name][Id][i][0] = result[i][0]
957                if data[name][Id][-1] == 'c':
958                    data[name][Id][-3] = str(result[-1][0])
959                elif data[name][Id][-1] == 'f':
960                    # process the variable name to put in global form (::var)
961                    varname = str(dlg.newvar[0]).strip().replace(' ','_')
962                    if varname.startswith('::'):
963                        varname = varname[2:]
964                    varname = varname.replace(':',';')
965                    if varname:
966                        data[name][Id][-3] = varname
967                    else:
968                        data[name][Id][-3] = ''
969                    data[name][Id][-2] = dlg.newvar[1]
970                if not CheckChangedConstraint():
971                    data[name][Id] = prev
972        except:
973            import traceback
974            print (traceback.format_exc())
975        finally:
976            dlg.Destroy()
977        wx.CallAfter(OnPageChanged,None)
978   
979    def UpdateConstraintPanel(panel,typ):
980        '''Update the contents of the selected Constraint
981        notebook tab. Called in :func:`OnPageChanged`
982        '''
983        if panel.GetSizer(): panel.GetSizer().Clear(True)
984        Siz = wx.BoxSizer(wx.VERTICAL)
985        Siz.Add((5,5),0)
986        Siz.Add(MakeConstraintsSizer(typ,panel),1,wx.EXPAND)
987        panel.SetSizer(Siz,True)
988        Size = Siz.GetMinSize()
989        Size[0] += 40
990        Size[1] = max(Size[1],450) + 20
991        panel.SetSize(Size)
992        panel.SetScrollbars(10,10,Size[0]/10-4,Size[1]/10-1)
993        panel.Show()
994
995    def OnPageChanged(event):
996        '''Called when a tab is pressed or when a "select tab" menu button is
997        used (see RaisePage), or to refresh the current tab contents (event=None)
998        '''
999        if event:       #page change event!
1000            page = event.GetSelection()
1001        else: # called directly, get current page
1002            page = G2frame.constr.GetSelection()
1003        G2frame.constr.ChangeSelection(page)
1004        text = G2frame.constr.GetPageText(page)
1005        G2frame.dataWindow.ConstraintEdit.Enable(G2G.wxID_EQUIVALANCEATOMS,False)
1006#        G2frame.dataWindow.ConstraintEdit.Enable(G2G.wxID_ADDRIDING,False)
1007        if text == 'Histogram/Phase':
1008            enableEditCons = [False]+4*[True]
1009            G2frame.Page = [page,'hap']
1010            UpdateConstraintPanel(HAPConstr,'HAP')
1011        elif text == 'Histogram':
1012            enableEditCons = [False]+4*[True]
1013            G2frame.Page = [page,'hst']
1014            UpdateConstraintPanel(HistConstr,'Hist')
1015        elif text == 'Phase':
1016            enableEditCons = 5*[True]
1017            G2frame.Page = [page,'phs']
1018            G2frame.dataWindow.ConstraintEdit.Enable(G2G.wxID_EQUIVALANCEATOMS,True)
1019#            G2frame.dataWindow.ConstraintEdit.Enable(G2G.wxID_ADDRIDING,True)
1020            if 'DELETED' in str(PhaseConstr):   #seems to be no other way to do this (wx bug)
1021                if GSASIIpath.GetConfigValue('debug'):
1022                    print ('DBG_wx error: PhaseConstr not cleanly deleted after Refine')
1023                return
1024            UpdateConstraintPanel(PhaseConstr,'Phase')
1025        elif text == 'Global':
1026            enableEditCons = [False]+4*[True]
1027            G2frame.Page = [page,'glb']
1028            UpdateConstraintPanel(GlobalConstr,'Global')
1029        else:
1030            enableEditCons = 5*[False]
1031            G2frame.Page = [page,'sym']
1032            UpdateConstraintPanel(SymConstr,'Sym-Generated')
1033        # remove menu items when not allowed
1034        for obj,flag in zip(G2frame.dataWindow.ConstraintEdit.GetMenuItems(),enableEditCons): 
1035            obj.Enable(flag)
1036        G2frame.dataWindow.SetDataSize()
1037
1038    def RaisePage(event):
1039        'Respond to a "select tab" menu button'
1040        try:
1041            i = (G2G.wxID_CONSPHASE,
1042                 G2G.wxID_CONSHAP,
1043                 G2G.wxID_CONSHIST,
1044                 G2G.wxID_CONSGLOBAL,
1045                 G2G.wxID_CONSSYM,
1046                ).index(event.GetId())
1047            G2frame.constr.SetSelection(i)
1048            wx.CallAfter(OnPageChanged,None)
1049        except ValueError:
1050            print('Unexpected event in RaisePage')
1051
1052    def SetStatusLine(text):
1053        G2frame.GetStatusBar().SetStatusText(text,1)
1054
1055    # UpdateConstraints execution starts here
1056    if not data:
1057        data.update({'Hist':[],'HAP':[],'Phase':[],'Global':[]})       #empty dict - fill it
1058    if 'Global' not in data:                                            #patch
1059        data['Global'] = []
1060    # DEBUG code ########################################
1061    #import GSASIIconstrGUI
1062    #reload(GSASIIconstrGUI)
1063    #reload(G2obj)
1064    #reload(G2stIO)
1065    #import GSASIIstrMain
1066    #reload(GSASIIstrMain)   
1067    #reload(G2mv)
1068    #reload(G2gd)
1069    ###################################################
1070    Histograms,Phases = G2frame.GetUsedHistogramsAndPhasesfromTree()
1071    if not len(Phases) or not len(Histograms):
1072        dlg = wx.MessageDialog(G2frame,'You need both phases and histograms to see Constraints',
1073            'No phases or histograms')
1074        dlg.CenterOnParent()
1075        dlg.ShowModal()
1076        dlg.Destroy()
1077        return
1078    G2obj.IndexAllIds(Histograms,Phases)
1079    ##################################################################################
1080    # patch: convert old-style (str) variables in constraints to G2VarObj objects
1081    for key,value in data.items():
1082        if key.startswith('_'): continue
1083        j = 0
1084        for cons in value:
1085            #print cons             # DEBUG
1086            for i in range(len(cons[:-3])):
1087                if type(cons[i][1]) is str:
1088                    cons[i][1] = G2obj.G2VarObj(cons[i][1])
1089                    j += 1
1090        if j:
1091            print (str(key) + ': '+str(j)+' variable(s) as strings converted to objects')
1092    ##################################################################################
1093    rigidbodyDict = G2frame.GPXtree.GetItemPyData(
1094        G2gd.GetGPXtreeItemId(G2frame,G2frame.root,'Rigid bodies'))
1095    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
1096    rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,Print=False)
1097    badPhaseParms = ['Ax','Ay','Az','Amul','AI/A','Atype','SHorder','mV0','mV1','mV2','AwaveType','FwaveType','PwaveType','MwaveType','Vol','isMag',]
1098    globalList = list(rbDict.keys())
1099    globalList.sort()
1100    try:
1101        AtomDict = dict([Phases[phase]['pId'],Phases[phase]['Atoms']] for phase in Phases)
1102    except KeyError:
1103        G2frame.ErrorDialog('Constraint Error','Constraints cannot be set until a cycle of least squares'+
1104                            ' has been run.\nWe suggest you refine a scale factor.')
1105        return
1106
1107    # create a list of the phase variables
1108    Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtable,BLtable,MFtable,maxSSwave = G2stIO.GetPhaseData(Phases,rbIds=rbIds,Print=False)
1109    phaseList = []
1110    for item in phaseDict:
1111        if item.split(':')[2] not in badPhaseParms:
1112            phaseList.append(item)
1113    phaseList.sort()
1114    phaseAtNames = {}
1115    phaseAtTypes = {}
1116    TypeList = []
1117    for item in phaseList:
1118        Split = item.split(':')
1119        if Split[2][:2] in ['AU','Af','dA','AM']:
1120            Id = int(Split[0])
1121            phaseAtNames[item] = AtomDict[Id][int(Split[3])][0]
1122            phaseAtTypes[item] = AtomDict[Id][int(Split[3])][1]
1123            if phaseAtTypes[item] not in TypeList:
1124                TypeList.append(phaseAtTypes[item])
1125        else:
1126            phaseAtNames[item] = ''
1127            phaseAtTypes[item] = ''
1128             
1129    # create a list of the hist*phase variables
1130    seqList = G2frame.testSeqRefineMode()
1131    if seqList: # for sequential refinement, only process 1st histgram in list
1132        histDict = {seqList[0]:Histograms[seqList[0]]}
1133    else:
1134        histDict = Histograms
1135    hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,histDict,Print=False,resetRefList=False)
1136    hapList = sorted([i for i in hapDict.keys() if i.split(':')[2] not in ('Type',)])
1137    if seqList: # convert histogram # to wildcard
1138        wildList = [] # list of variables with "*" for histogram number
1139        for i in hapList:
1140            s = i.split(':')
1141            if s[1] == "": continue
1142            s[1] = '*'
1143            sj = ':'.join(s)
1144            if sj not in wildList: wildList.append(sj)
1145        hapList = wildList
1146    histVary,histDict,controlDict = G2stIO.GetHistogramData(histDict,Print=False)
1147    histList = list(histDict.keys())
1148    histList.sort()
1149    if seqList: # convert histogram # to wildcard
1150        wildList = [] # list of variables with "*" for histogram number
1151        for i in histList:
1152            s = i.split(':')
1153            if s[1] == "": continue
1154            s[1] = '*'
1155            sj = ':'.join(s)
1156            if sj not in wildList: wildList.append(sj)
1157        histList = wildList       
1158    Indx = {}
1159    G2frame.Page = [0,'phs']
1160   
1161    G2gd.SetDataMenuBar(G2frame,G2frame.dataWindow.ConstraintMenu)
1162    SetStatusLine('')
1163   
1164    G2frame.Bind(wx.EVT_MENU, OnAddConstraint, id=G2G.wxID_CONSTRAINTADD)
1165    G2frame.Bind(wx.EVT_MENU, OnAddFunction, id=G2G.wxID_FUNCTADD)
1166    G2frame.Bind(wx.EVT_MENU, OnAddEquivalence, id=G2G.wxID_EQUIVADD)
1167    G2frame.Bind(wx.EVT_MENU, OnAddHold, id=G2G.wxID_HOLDADD)
1168    G2frame.Bind(wx.EVT_MENU, OnAddAtomEquiv, id=G2G.wxID_EQUIVALANCEATOMS)
1169#    G2frame.Bind(wx.EVT_MENU, OnAddRiding, id=G2G.wxID_ADDRIDING)
1170    # tab commands
1171    for id in (G2G.wxID_CONSPHASE,
1172               G2G.wxID_CONSHAP,
1173               G2G.wxID_CONSHIST,
1174               G2G.wxID_CONSGLOBAL,
1175               G2G.wxID_CONSSYM,
1176               ):
1177        G2frame.Bind(wx.EVT_MENU, RaisePage,id=id)
1178
1179    #G2frame.constr = G2G.GSNoteBook(parent=G2frame.dataWindow,size=G2frame.dataWindow.GetClientSize())
1180    G2frame.constr = G2G.GSNoteBook(parent=G2frame.dataWindow)
1181    G2frame.dataWindow.GetSizer().Add(G2frame.constr,1,wx.ALL|wx.EXPAND)
1182    # note that order of pages is hard-coded in RaisePage
1183    PhaseConstr = wx.ScrolledWindow(G2frame.constr)
1184    G2frame.constr.AddPage(PhaseConstr,'Phase')
1185    HAPConstr = wx.ScrolledWindow(G2frame.constr)
1186    G2frame.constr.AddPage(HAPConstr,'Histogram/Phase')
1187    HistConstr = wx.ScrolledWindow(G2frame.constr)
1188    G2frame.constr.AddPage(HistConstr,'Histogram')
1189    GlobalConstr = wx.ScrolledWindow(G2frame.constr)
1190    G2frame.constr.AddPage(GlobalConstr,'Global')
1191    SymConstr = wx.ScrolledWindow(G2frame.constr)
1192    G2frame.constr.AddPage(SymConstr,'Sym-Generated')
1193    wx.CallAfter(OnPageChanged,None)
1194    G2frame.constr.Bind(wx.aui.EVT_AUINOTEBOOK_PAGE_CHANGED, OnPageChanged)
1195    # validate all the constrants -- should not see any errors here normally
1196    errmsg,warnmsg = ConstraintsCheck(data)
1197    if errmsg:
1198        G2frame.ErrorDialog('Constraint Error',
1199                            'Error in constraints:\n'+errmsg+'\nCheck console output for more information',
1200                            parent=G2frame)
1201        print (errmsg)
1202        print (G2mv.VarRemapShow([],True))
1203    elif warnmsg:
1204        print ('Unexpected contraint warning:\n'+warnmsg)
1205       
1206################################################################################
1207# check scale & phase fractions, create constraint if needed
1208################################################################################
1209def CheckAllScalePhaseFractions(G2frame):
1210    '''Check if scale factor and all phase fractions are refined without a constraint
1211    for all used histograms, if so, offer the user a chance to create a constraint
1212    on the sum of phase fractions
1213    '''
1214    histograms, phases = G2frame.GetUsedHistogramsAndPhasesfromTree()
1215    for i,hist in enumerate(histograms):
1216        CheckScalePhaseFractions(G2frame,hist,histograms,phases)
1217       
1218def CheckScalePhaseFractions(G2frame,hist,histograms,phases):
1219    '''Check if scale factor and all phase fractions are refined without a constraint
1220    for histogram hist, if so, offer the user a chance to create a constraint
1221    on the sum of phase fractions
1222    '''
1223    if G2frame.testSeqRefineMode():
1224        histStr = '*'
1225    else:
1226        histStr = str(histograms[hist]['hId'])
1227    # Is this powder?
1228    if not hist.startswith('PWDR '): return
1229    # do this only if the scale factor is varied
1230    if not histograms[hist]['Sample Parameters']['Scale'][1]: return
1231    # are all phase fractions varied in all used histograms?
1232    phaseCount = 0
1233    for p in phases:
1234        if hist not in phases[p]['Histograms']: continue
1235        if phases[p]['Histograms'][hist]['Use'] and not phases[p]['Histograms'][hist]['Scale'][1]:
1236            return
1237        else:
1238            phaseCount += 1
1239   
1240    # all phase fractions and scale factor varied, now scan through constraints
1241    sub = G2gd.GetGPXtreeItemId(G2frame,G2frame.root,'Constraints') 
1242    Constraints = G2frame.GPXtree.GetItemPyData(sub)
1243    for c in Constraints.get('HAP',[]):
1244        if c[-1] != 'c': continue
1245        if not c[-3]: continue
1246        if len(c[:-3]) != phaseCount: continue
1247        # got a constraint equation with right number of terms, is it on phase fractions for
1248        # the correct histogram?
1249        if all([(i[1].name == 'Scale' and i[1].varname().split(':')[1] == histStr) for i in c[:-3]]):
1250            # got a constraint, this is OK
1251            return
1252    dlg = wx.MessageDialog(G2frame,'You are refining the scale factor and all phase fractions for histogram #'+
1253        histStr+'. This will produce an unstable refinement. '+
1254        'Do you want to constrain the sum of phase fractions?','Create constraint?',wx.OK|wx.CANCEL)
1255    if dlg.ShowModal() != wx.ID_OK:
1256        dlg.Destroy()
1257        return
1258    dlg.Destroy()
1259
1260    constr = []
1261    for p in phases:
1262        if hist not in phases[p]['Histograms']: continue
1263        if not phases[p]['Histograms'][hist]['Use']: continue
1264        constr += [[1.0,G2obj.G2VarObj(':'.join((str(phases[p]['pId']),histStr,'Scale')))]]
1265    constr += [1.0,None,'c']
1266    Constraints['HAP'] += [constr]
1267       
1268################################################################################
1269#### Make nuclear/magnetic phase transition constraints - called by OnTransform in G2phsGUI
1270################################################################################       
1271       
1272def TransConstraints(G2frame,oldPhase,newPhase,Trans,Vec,atCodes):
1273    '''Add constraints for new magnetic phase created via transformation of old
1274    nuclear one
1275    NB: A = [G11,G22,G33,2*G12,2*G13,2*G23]
1276    '''
1277   
1278    def SetUniqAj(pId,iA,SGData):
1279        SGLaue = SGData['SGLaue']
1280        SGUniq = SGData['SGUniq']
1281        if SGLaue in ['m3','m3m']:
1282            if iA in [0,1,2]:
1283                parm = '%d::%s'%(pId,'A0')
1284            else:
1285                parm = None
1286        elif SGLaue in ['4/m','4/mmm']:
1287            if iA in [0,1]:
1288                parm = '%d::%s'%(pId,'A0')
1289            elif iA == 2:
1290                parm = '%d::%s'%(pId,'A2')
1291            else:
1292                parm = None
1293        elif SGLaue in ['6/m','6/mmm','3m1', '31m', '3']:
1294            if iA in [0,1,3]:
1295                parm = '%d::%s'%(pId,'A0')
1296            elif iA == 2:
1297                parm = '%d::%s'%(pId,'A2')
1298            else:
1299                parm = None
1300        elif SGLaue in ['3R', '3mR']:
1301            if ia in [0,1,2]:
1302                parm = '%d::%s'%(pId,'A0')
1303            else:
1304                parm = '%d::%s'%(pId,'A3')
1305        elif SGLaue in ['mmm',]:
1306            if iA in [0,1,2]:
1307                parm = '%d::A%s'%(pId,iA)
1308            else:
1309                parm = None
1310        elif SGLaue == '2/m':
1311            if iA in [0,1,2]:
1312                parm = '%d::A%s'%(pId,iA)
1313            elif iA == 3 and SGUniq == 'c':
1314                parm = '%d::A%s'%(pId,iA)
1315            elif iA == 4 and SGUniq == 'b':
1316                parm = '%d::A%s'%(pId,iA)
1317            elif iA == 5 and SGUniq == 'a':
1318                parm = '%d::A%s'%(pId,iA)
1319            else:
1320                parm = None           
1321        else:
1322            parm = '%d::A%s'%(pId,iA)
1323        return parm
1324   
1325    Histograms,Phases = G2frame.GetUsedHistogramsAndPhasesfromTree()
1326    UseList = newPhase['Histograms']
1327    detTrans = np.abs(nl.det(Trans))
1328    opId = oldPhase['pId']
1329    npId = newPhase['pId']
1330    cx,ct,cs,cia = newPhase['General']['AtomPtrs']
1331    nAtoms = newPhase['Atoms']
1332    oSGData = oldPhase['General']['SGData']
1333    nSGData = newPhase['General']['SGData']
1334    #oAcof = G2lat.cell2A(oldPhase['General']['Cell'][1:7])
1335    #nAcof = G2lat.cell2A(newPhase['General']['Cell'][1:7])
1336    item = G2gd.GetGPXtreeItemId(G2frame,G2frame.root,'Constraints')
1337    if not item:
1338        print('Error: no constraints in Data Tree')
1339        return
1340    constraints = G2frame.GPXtree.GetItemPyData(item)
1341    xnames = ['dAx','dAy','dAz']
1342    # constraints on matching atom params between phases
1343    for ia,code in enumerate(atCodes):
1344        atom = nAtoms[ia]
1345        if not ia and atom[cia] == 'A':
1346            wx.MessageDialog(G2frame,
1347                'Anisotropic thermal motion constraints are not developed at the present time',
1348                'Anisotropic thermal constraint?',style=wx.ICON_INFORMATION).ShowModal()
1349        siteSym = G2spc.SytSym(atom[cx:cx+3],nSGData)[0]
1350        CSX = G2spc.GetCSxinel(siteSym)
1351#        CSU = G2spc.GetCSuinel(siteSym)
1352        item = code.split('+')[0]
1353        iat,opr = item.split(':')
1354        Nop = abs(int(opr))%100-1
1355        if '-' in opr:
1356            Nop *= -1
1357        Opr = oldPhase['General']['SGData']['SGOps'][abs(Nop)][0]
1358        if Nop < 0:         #inversion
1359            Opr *= -1
1360        XOpr = np.inner(Opr,Trans)
1361        for ix in list(set(CSX[0])):
1362            if not ix:
1363                continue
1364            name = xnames[ix-1]
1365            IndpCon = [1.0,G2obj.G2VarObj('%d::%s:%d'%(npId,name,ia))]
1366            DepCons = []
1367            for iop,opval in enumerate(XOpr[ix-1]):
1368                if opval:
1369                    DepCons.append([opval,G2obj.G2VarObj('%d::%s:%s'%(opId,xnames[iop],iat))])
1370            if len(DepCons) == 1:
1371                constraints['Phase'].append([DepCons[0],IndpCon,None,None,'e'])
1372            elif len(DepCons) > 1:
1373                for Dep in DepCons:
1374                    Dep[0] *= -1
1375                constraints['Phase'].append([IndpCon]+DepCons+[0.0,None,'c'])
1376        for name in ['Afrac','AUiso']:
1377            IndpCon = [1.0,G2obj.G2VarObj('%d::%s:%d'%(npId,name,ia))]
1378            DepCons = [1.0,G2obj.G2VarObj('%d::%s:%s'%(opId,name,iat))]
1379            constraints['Phase'].append([DepCons,IndpCon,None,None,'e'])
1380           
1381        # unfinished Anisotropic constraint generation
1382#        Uids = [[0,0,'AU11'],[1,1,'AU22'],[2,2,'AU33'],[0,1,'AU12'],[0,2,'AU13'],[1,2,'AU23']]
1383#        DepConsDict = dict(zip(Us,[[],[],[],[],[],[]]))
1384#        for iu,Uid in enumerate(Uids):
1385#            UMT = np.zeros((3,3))
1386#            UMT[Uid[0],Uid[1]] = 1
1387#            nUMT = G2lat.prodMGMT(UMT,invTrans)
1388#            nUT = G2lat.UijtoU6(nUMT)
1389#            for iu,nU in enumerate(nUT):
1390#                if abs(nU) > 1.e-8:
1391#                    parm = '%d::%s;%s'%(opId,Us[iu],iat)
1392#                    DepConsDict[Uid[2]].append([abs(nU%1.),G2obj.G2VarObj(parm)])
1393#        nUcof = atom[iu:iu+6]
1394#        conStrings = []
1395#        for iU,Usi in enumerate(Us):
1396#            parm = '%d::%s;%d'%(npId,Usi,ia)
1397#            parmDict[parm] = nUcof[iU]
1398#            varyList.append(parm)
1399#            IndpCon = [1.0,G2obj.G2VarObj(parm)]
1400#            conStr = str([IndpCon,DepConsDict[Usi]])
1401#            if conStr in conStrings:
1402#                continue
1403#            conStrings.append(conStr)
1404#            if len(DepConsDict[Usi]) == 1:
1405#                if DepConsDict[Usi][0]:
1406#                    constraints['Phase'].append([IndpCon,DepConsDict[Usi][0],None,None,'e'])
1407#            elif len(DepConsDict[Usi]) > 1:       
1408#                for Dep in DepConsDict[Usi]:
1409#                    Dep[0] *= -1
1410#                constraints['Phase'].append([IndpCon]+DepConsDict[Usi]+[0.0,None,'c'])
1411           
1412        #how do I do Uij's for most Trans?
1413
1414    # constraints on lattice parameters between phases
1415    T = nl.inv(Trans).T
1416    conMat = [
1417        [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]],
1418        [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]],
1419        [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]],
1420        [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]],
1421        [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]],
1422        [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]]
1423        ]
1424    # Gnew = conMat * A:
1425#         T00**2*a0  T01**2*a1 T02**2*a2 T00*T01*a3    T00*T02*a4    T01*T02*a5
1426#         T10**2*a0  T11**2*a1 T12**2*a2 T10*T11*a3    T10*T12*a4    T11*T12*a5
1427#         T20**2*a0  T21**2*a1 T22**2*a2 T20*T21*a3    T20*T22*a4    T21*T22*a5
1428#         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
1429#         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
1430#         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
1431    # Generated as symbolic code using:
1432    # import sym
1433    # A0, A1, A2, A3, A4, A5 = sym.symbols('A0, A1, A2, A3, A4, A5')
1434    # G = sym.Matrix([ [A0,    A3/2.,  A4/2.], [A3/2.,  A1,    A5/2.], [A4/2., A5/2.,    A2]])
1435    # transformation matrix
1436    # T00, T10, T20, T01, T11, T21, T02, T12, T22 = sym.symbols('T00, T10, T20, T01, T11, T21, T02, T12, T22')
1437    # Tr = sym.Matrix([ [T00, T10, T20], [T01, T11, T21], [T02, T12, T22],])
1438    # Gnew = (Tr.T*G)*Tr
1439   
1440    #print('old A',G2lat.cell2A(oldPhase['General']['Cell'][1:7]))
1441    #print('new A',G2lat.cell2A(newPhase['General']['Cell'][1:7]))
1442    for iAnew,Asi in enumerate(['A0','A1','A2','A3','A4','A5']): # loop through A[i] for new cell
1443        Nparm = str(npId) + '::' + Asi
1444        if Nparm != SetUniqAj(npId,iAnew,nSGData):
1445            continue # skip: Ai constrained from Aj or must be zero
1446        multDict = {}
1447        for iAorg in range(6):
1448            cA = conMat[iAnew][iAorg] # coeff for A[i] in constraint matrix
1449            if abs(cA) < 1.e-8: continue
1450            parm = SetUniqAj(opId,iAorg,oSGData) # translate to unique A[i] in original cell
1451            if not parm: continue # must be zero
1452            # sum coeff
1453            if parm in multDict:
1454                multDict[parm] += cA
1455            else:
1456                multDict[parm] = cA
1457        # any non-zero multipliers?
1458        maxMult = 0
1459        for i in multDict:
1460            maxMult = max(maxMult,abs(multDict[i]))
1461        if maxMult <= 0:  # Nparm computes as zero; Fix this parameter
1462            constraints['Phase'] += [[
1463                [0.0,G2obj.G2VarObj(Nparm)],
1464                None,None,'h']]
1465        elif len(multDict) == 1:        # create equivalence
1466            key = list(multDict.keys())[0]
1467            constraints['Phase'] += [[
1468                [1.0,G2obj.G2VarObj(key)],
1469                [multDict[key],G2obj.G2VarObj(Nparm)],
1470                None,None,'e']]
1471        else:                           # create constraint
1472            constr = [[-1.0,G2obj.G2VarObj(Nparm)]]
1473            for key in multDict:
1474                constr += [[multDict[key],G2obj.G2VarObj(key)]]
1475            constr += [0.0,None,'c']
1476            constraints['Phase'] += [constr]
1477   
1478    # constraints on HAP Scale, etc.
1479    for hId,hist in enumerate(UseList):    #HAP - seems OK
1480        ohapkey = '%d:%d:'%(opId,hId)
1481        nhapkey = '%d:%d:'%(npId,hId)
1482        IndpCon = [1.0,G2obj.G2VarObj(ohapkey+'Scale')]
1483        DepCons = [detTrans,G2obj.G2VarObj(nhapkey+'Scale')]
1484        constraints['HAP'].append([DepCons,IndpCon,None,None,'e'])
1485        for name in ['Size;i','Mustrain;i']:
1486            IndpCon = [1.0,G2obj.G2VarObj(ohapkey+name)]
1487            DepCons = [1.0,G2obj.G2VarObj(nhapkey+name)]
1488            constraints['HAP'].append([IndpCon,DepCons,None,None,'e'])
1489       
1490################################################################################
1491#### Rigid bodies
1492################################################################################
1493
1494def UpdateRigidBodies(G2frame,data):
1495    '''Called when Rigid bodies tree item is selected.
1496    Displays the rigid bodies in the data window
1497    '''
1498    if not data.get('RBIds') or not data:
1499        data.update({'Vector':{'AtInfo':{}},'Residue':{'AtInfo':{}},
1500            'RBIds':{'Vector':[],'Residue':[]}})       #empty/bad dict - fill it
1501           
1502    global resList,rbId
1503    Indx = {}
1504    resList = []
1505    plotDefaults = {'oldxy':[0.,0.],'Quaternion':[0.,0.,0.,1.],'cameraPos':30.,'viewDir':[0,0,1],}
1506
1507    G2frame.rbBook = G2G.GSNoteBook(parent=G2frame.dataWindow)
1508    G2frame.dataWindow.GetSizer().Add(G2frame.rbBook,1,wx.ALL|wx.EXPAND)
1509    VectorRB = wx.ScrolledWindow(G2frame.rbBook)
1510    VectorRBDisplay = wx.Panel(VectorRB)
1511    G2frame.rbBook.AddPage(VectorRB,'Vector rigid bodies')
1512    ResidueRB = wx.ScrolledWindow(G2frame.rbBook)
1513    ResidueRBDisplay = wx.Panel(ResidueRB)
1514    G2frame.rbBook.AddPage(ResidueRB,'Residue rigid bodies')
1515   
1516    def OnPageChanged(event):
1517        global resList
1518        resList = []
1519        if event:       #page change event!
1520            page = event.GetSelection()
1521        else:
1522            page = G2frame.rbBook.GetSelection()
1523        G2frame.rbBook.ChangeSelection(page)
1524        text = G2frame.rbBook.GetPageText(page)
1525        if text == 'Vector rigid bodies':
1526            G2gd.SetDataMenuBar(G2frame,G2frame.dataWindow.VectorBodyMenu)
1527            G2frame.Bind(wx.EVT_MENU, AddVectorRB, id=G2G.wxID_VECTORBODYADD)
1528            G2frame.Page = [page,'vrb']
1529            UpdateVectorRB()
1530        elif text == 'Residue rigid bodies':
1531            G2gd.SetDataMenuBar(G2frame,G2frame.dataWindow.RigidBodyMenu)
1532            G2frame.Bind(wx.EVT_MENU, AddResidueRB, id=G2G.wxID_RIGIDBODYADD)
1533            G2frame.Bind(wx.EVT_MENU, OnImportRigidBody, id=G2G.wxID_RIGIDBODYIMPORT)
1534            G2frame.Bind(wx.EVT_MENU, OnDefineTorsSeq, id=G2G.wxID_RESIDUETORSSEQ) #enable only if residue RBs exist?
1535            G2frame.Page = [page,'rrb']
1536            UpdateResidueRB()
1537           
1538    def getMacroFile(macName):
1539        defDir = os.path.join(os.path.split(__file__)[0],'GSASIImacros')
1540        dlg = wx.FileDialog(G2frame,message='Choose '+macName+' rigid body macro file',
1541            defaultDir=defDir,defaultFile="",wildcard="GSAS-II macro file (*.mac)|*.mac",
1542            style=wx.FD_OPEN | wx.FD_CHANGE_DIR)
1543        try:
1544            if dlg.ShowModal() == wx.ID_OK:
1545                macfile = dlg.GetPath()
1546                macro = open(macfile,'Ur')
1547                head = macro.readline()
1548                if macName not in head:
1549                    print (head)
1550                    print ('**** ERROR - wrong restraint macro file selected, try again ****')
1551                    macro = []
1552            else: # cancel was pressed
1553                macro = []
1554        finally:
1555            dlg.Destroy()
1556        return macro        #advanced past 1st line
1557       
1558    def getTextFile():
1559        dlg = wx.FileDialog(G2frame,'Choose rigid body text file', '.', '',
1560            "GSAS-II text file (*.txt)|*.txt|XYZ file (*.xyz)|*.xyz|"
1561            "Sybyl mol2 file (*.mol2)|*.mol2|PDB file (*.pdb;*.ent)|*.pdb;*.ent",
1562            wx.FD_OPEN | wx.FD_CHANGE_DIR)
1563        try:
1564            if dlg.ShowModal() == wx.ID_OK:
1565                txtfile = dlg.GetPath()
1566                ext = os.path.splitext(txtfile)[1]
1567                text = open(txtfile,'Ur')
1568            else: # cancel was pressed
1569                ext = ''
1570                text = []
1571        finally:
1572            dlg.Destroy()
1573        if 'ent' in ext:
1574            ext = '.pdb'
1575        return text,ext.lower()
1576       
1577    def OnImportRigidBody(event):
1578        page = G2frame.rbBook.GetSelection()
1579        if 'Vector' in G2frame.rbBook.GetPageText(page):
1580            pass
1581        elif 'Residue' in G2frame.rbBook.GetPageText(page):
1582            ImportResidueRB()
1583           
1584    def AddVectorRB(event):
1585        AtInfo = data['Vector']['AtInfo']
1586        dlg = G2G.MultiIntegerDialog(G2frame,'New Rigid Body',['No. atoms','No. translations'],[1,1])
1587        if dlg.ShowModal() == wx.ID_OK:
1588            nAtoms,nTrans = dlg.GetValues()
1589            rbId = ran.randint(0,sys.maxsize)
1590            vecMag = [1.0 for i in range(nTrans)]
1591            vecRef = [False for i in range(nTrans)]
1592            vecVal = [np.zeros((nAtoms,3)) for j in range(nTrans)]
1593            rbTypes = ['C' for i in range(nAtoms)]
1594            Info = G2elem.GetAtomInfo('C')
1595            AtInfo['C'] = [Info['Drad'],Info['Color']]
1596            data['Vector'][rbId] = {'RBname':'UNKRB','VectMag':vecMag,'rbXYZ':np.zeros((nAtoms,3)),
1597                'rbRef':[0,1,2,False],'VectRef':vecRef,'rbTypes':rbTypes,'rbVect':vecVal,'useCount':0}
1598            data['RBIds']['Vector'].append(rbId)
1599        dlg.Destroy()
1600        UpdateVectorRB()
1601       
1602    def AddResidueRB(event):
1603        AtInfo = data['Residue']['AtInfo']
1604        macro = getMacroFile('rigid body')
1605        if not macro:
1606            return
1607        macStr = macro.readline()
1608        while macStr:
1609            items = macStr.split()
1610            if 'I' == items[0]:
1611                rbId = ran.randint(0,sys.maxsize)
1612                rbName = items[1]
1613                rbTypes = []
1614                rbXYZ = []
1615                rbSeq = []
1616                atNames = []
1617                nAtms,nSeq,nOrig,mRef,nRef = [int(items[i]) for i in [2,3,4,5,6]]
1618                for iAtm in range(nAtms):
1619                    macStr = macro.readline().split()
1620                    atName = macStr[0]
1621                    atType = macStr[1]
1622                    atNames.append(atName)
1623                    rbXYZ.append([float(macStr[i]) for i in [2,3,4]])
1624                    rbTypes.append(atType)
1625                    if atType not in AtInfo:
1626                        Info = G2elem.GetAtomInfo(atType)
1627                        AtInfo[atType] = [Info['Drad'],Info['Color']]
1628                rbXYZ = np.array(rbXYZ)-np.array(rbXYZ[nOrig-1])
1629                for iSeq in range(nSeq):
1630                    macStr = macro.readline().split()
1631                    mSeq = int(macStr[0])
1632                    for jSeq in range(mSeq):
1633                        macStr = macro.readline().split()
1634                        iBeg = int(macStr[0])-1
1635                        iFin = int(macStr[1])-1
1636                        angle = 0.0
1637                        nMove = int(macStr[2])
1638                        iMove = [int(macStr[i])-1 for i in range(3,nMove+3)]
1639                        rbSeq.append([iBeg,iFin,angle,iMove])
1640                data['Residue'][rbId] = {'RBname':rbName,'rbXYZ':rbXYZ,'rbTypes':rbTypes,
1641                    'atNames':atNames,'rbRef':[nOrig-1,mRef-1,nRef-1,True],'rbSeq':rbSeq,
1642                    'SelSeq':[0,0],'useCount':0}
1643                data['RBIds']['Residue'].append(rbId)
1644                print ('Rigid body '+rbName+' added')
1645            macStr = macro.readline()
1646        macro.close()
1647        UpdateResidueRB()
1648       
1649    def ImportResidueRB():
1650        AtInfo = data['Residue']['AtInfo']
1651        text,ext = getTextFile()
1652        if not text:
1653            return
1654        rbId = ran.randint(0,sys.maxsize)
1655        rbTypes = []
1656        rbXYZ = []
1657        atNames = []
1658        txtStr = text.readline()
1659        if 'xyz' in ext:
1660            txtStr = text.readline()
1661            txtStr = text.readline()
1662        elif 'mol2' in ext:
1663            while 'ATOM' not in txtStr:
1664                txtStr = text.readline()
1665            txtStr = text.readline()
1666        elif 'pdb' in ext:
1667            while 'ATOM' not in txtStr[:6] and 'HETATM' not in txtStr[:6]:
1668                txtStr = text.readline()
1669                #print txtStr
1670        items = txtStr.split()
1671        while len(items):
1672            if 'txt' in ext:
1673                atName = items[0]
1674                atType = items[1]
1675                rbXYZ.append([float(items[i]) for i in [2,3,4]])
1676            elif 'xyz' in ext:
1677                atType = items[0]
1678                rbXYZ.append([float(items[i]) for i in [1,2,3]])
1679                atName = atType+str(len(rbXYZ))
1680            elif 'mol2' in ext:
1681                atType = items[1]
1682                atName = items[1]+items[0]
1683                rbXYZ.append([float(items[i]) for i in [2,3,4]])
1684            elif 'pdb' in ext:
1685                atType = items[-1]
1686                atName = items[2]
1687                xyz = txtStr[30:55].split()                   
1688                rbXYZ.append([float(x) for x in xyz])
1689            atNames.append(atName)
1690            rbTypes.append(atType)
1691            if atType not in AtInfo:
1692                Info = G2elem.GetAtomInfo(atType)
1693                AtInfo[atType] = [Info['Drad'],Info['Color']]
1694            txtStr = text.readline()
1695            if 'mol2' in ext and 'BOND' in txtStr:
1696                break
1697            if 'pdb' in ext and ('ATOM' not in txtStr[:6] and 'HETATM' not in txtStr[:6]):
1698                break
1699            items = txtStr.split()
1700        if len(atNames) < 3:
1701            G2G.G2MessageBox(G2frame,'Not enough atoms in rigid body; must be 3 or more')
1702        else:
1703            rbXYZ = np.array(rbXYZ)-np.array(rbXYZ[0])
1704            Xxyz = rbXYZ[1]
1705            X = Xxyz/np.sqrt(np.sum(Xxyz**2))
1706            Yxyz = rbXYZ[2]
1707            Y = Yxyz/np.sqrt(np.sum(Yxyz**2))
1708            Mat = G2mth.getRBTransMat(X,Y)
1709            rbXYZ = np.inner(Mat,rbXYZ).T
1710            data['Residue'][rbId] = {'RBname':'UNKRB','rbXYZ':rbXYZ,'rbTypes':rbTypes,
1711                'atNames':atNames,'rbRef':[0,1,2,False],'rbSeq':[],'SelSeq':[0,0],'useCount':0}
1712            data['RBIds']['Residue'].append(rbId)
1713            print ('Rigid body UNKRB added')
1714        text.close()
1715        UpdateResidueRB(rbId)
1716       
1717    def FindNeighbors(Orig,XYZ,atTypes,atNames,AtInfo):
1718        Radii = []
1719        for Atype in atTypes:
1720            Radii.append(AtInfo[Atype][0])
1721        Radii = np.array(Radii)
1722        Neigh = []
1723        Dx = XYZ-XYZ[Orig]
1724        dist = np.sqrt(np.sum(Dx**2,axis=1))
1725        sumR = Radii[Orig]+Radii
1726        IndB = ma.nonzero(ma.masked_greater(dist-0.85*sumR,0.))
1727        for j in IndB[0]:
1728            if j != Orig and atTypes[j] != 'H':
1729                Neigh.append(atNames[j])
1730        return Neigh
1731       
1732    def FindAllNeighbors(XYZ,atTypes,atNames,AtInfo):
1733        NeighDict = {}
1734        for iat,xyz in enumerate(atNames):
1735            NeighDict[atNames[iat]] = FindNeighbors(iat,XYZ,atTypes,atNames,AtInfo)
1736        return NeighDict
1737       
1738    def FindRiding(Orig,Pivot,NeighDict):
1739        riding = [Orig,Pivot]
1740        iAdd = 1
1741        new = True
1742        while new:
1743            newAtms = NeighDict[riding[iAdd]]
1744            for At in newAtms:
1745                new = False
1746                if At not in riding:
1747                    riding.append(At)
1748                    new = True
1749            iAdd += 1
1750            if iAdd < len(riding):
1751                new = True
1752        return riding[2:]
1753                       
1754    def OnDefineTorsSeq(event):
1755        global rbId
1756        rbData = data['Residue'][rbId]
1757        if not len(rbData):
1758            return
1759        atNames = rbData['atNames']
1760        AtInfo = data['Residue']['AtInfo']
1761        atTypes = rbData['rbTypes']
1762        XYZ = rbData['rbXYZ']
1763        neighDict = FindAllNeighbors(XYZ,atTypes,atNames,AtInfo)
1764        TargList = []           
1765        dlg = wx.SingleChoiceDialog(G2frame,'Select origin atom for torsion sequence','Origin atom',rbData['atNames'])
1766        if dlg.ShowModal() == wx.ID_OK:
1767            Orig = dlg.GetSelection()
1768            TargList = neighDict[atNames[Orig]]
1769        dlg.Destroy()
1770        if not len(TargList):
1771            return
1772        dlg = wx.SingleChoiceDialog(G2frame,'Select pivot atom for torsion sequence','Pivot atom',TargList)
1773        if dlg.ShowModal() == wx.ID_OK:
1774            Piv = atNames.index(TargList[dlg.GetSelection()])
1775            riding = FindRiding(atNames[Orig],atNames[Piv],neighDict)
1776            Riding = []
1777            for atm in riding:
1778                Riding.append(atNames.index(atm))
1779            rbData['rbSeq'].append([Orig,Piv,0.0,Riding])           
1780        dlg.Destroy()
1781        UpdateResidueRB(rbId)
1782
1783    def UpdateVectorRB(Scroll=0):
1784        AtInfo = data['Vector']['AtInfo']
1785        refChoice = {}
1786        if 'DELETED' in str(G2frame.GetStatusBar()):   #seems to be no other way to do this (wx bug)
1787            if GSASIIpath.GetConfigValue('debug'):
1788                print ('DBG_wx error: Rigid Body/Status not cleanly deleted after Refine')
1789            return
1790        SetStatusLine(' You may use e.g. "c60" or "s60" for a vector entry')
1791        def rbNameSizer(rbId,rbData):
1792
1793            def OnRBName(event):
1794                event.Skip()
1795                Obj = event.GetEventObject()
1796                rbData['RBname'] = Obj.GetValue()
1797               
1798            def OnDelRB(event):
1799                Obj = event.GetEventObject()
1800                rbId = Indx[Obj.GetId()]
1801                if rbId in data['Vector']:
1802                    del data['Vector'][rbId]
1803                    data['RBIds']['Vector'].remove(rbId)
1804                    rbData['useCount'] -= 1
1805                wx.CallAfter(UpdateVectorRB)
1806               
1807            def OnPlotRB(event):
1808                Obj = event.GetEventObject()
1809                Obj.SetValue(False)
1810                G2plt.PlotRigidBody(G2frame,'Vector',AtInfo,rbData,plotDefaults)
1811           
1812            nameSizer = wx.BoxSizer(wx.HORIZONTAL)
1813            nameSizer.Add(wx.StaticText(VectorRBDisplay,-1,'Rigid body name: '),
1814                0,wx.ALIGN_CENTER_VERTICAL)
1815            RBname = wx.TextCtrl(VectorRBDisplay,-1,rbData['RBname'])
1816            Indx[RBname.GetId()] = rbId
1817            RBname.Bind(wx.EVT_TEXT_ENTER,OnRBName)
1818            RBname.Bind(wx.EVT_KILL_FOCUS,OnRBName)
1819            nameSizer.Add(RBname,0,wx.ALIGN_CENTER_VERTICAL)
1820            nameSizer.Add((5,0),)
1821            plotRB = wx.CheckBox(VectorRBDisplay,-1,'Plot?')
1822            Indx[plotRB.GetId()] = rbId
1823            plotRB.Bind(wx.EVT_CHECKBOX,OnPlotRB)
1824            nameSizer.Add(plotRB,0,wx.ALIGN_CENTER_VERTICAL)
1825            nameSizer.Add((5,0),)
1826            if not rbData['useCount']:
1827                delRB = wx.CheckBox(VectorRBDisplay,-1,'Delete?')
1828                Indx[delRB.GetId()] = rbId
1829                delRB.Bind(wx.EVT_CHECKBOX,OnDelRB)
1830                nameSizer.Add(delRB,0,wx.ALIGN_CENTER_VERTICAL)
1831            return nameSizer
1832           
1833        def rbRefAtmSizer(rbId,rbData):
1834           
1835            def OnRefSel(event):
1836                Obj = event.GetEventObject()
1837                iref = Indx[Obj.GetId()]
1838                sel = Obj.GetValue()
1839                rbData['rbRef'][iref] = atNames.index(sel)
1840                FillRefChoice(rbId,rbData)
1841           
1842            refAtmSizer = wx.BoxSizer(wx.HORIZONTAL)
1843            atNames = [name+str(i) for i,name in enumerate(rbData['rbTypes'])]
1844            rbRef = rbData.get('rbRef',[0,1,2,False])
1845            rbData['rbRef'] = rbRef
1846            if rbData['useCount']:
1847                refAtmSizer.Add(wx.StaticText(VectorRBDisplay,-1,
1848                    'Orientation reference atoms A-B-C: %s, %s, %s'%(atNames[rbRef[0]], \
1849                     atNames[rbRef[1]],atNames[rbRef[2]])),0)
1850            else:
1851                refAtmSizer.Add(wx.StaticText(VectorRBDisplay,-1,
1852                    'Orientation reference atoms A-B-C: '),0,wx.ALIGN_CENTER_VERTICAL)
1853                for i in range(3):
1854                    choices = [atNames[j] for j in refChoice[rbId][i]]
1855                    refSel = wx.ComboBox(VectorRBDisplay,-1,value='',
1856                        choices=choices,style=wx.CB_READONLY|wx.CB_DROPDOWN)
1857                    refSel.SetValue(atNames[rbRef[i]])
1858                    refSel.Bind(wx.EVT_COMBOBOX, OnRefSel)
1859                    Indx[refSel.GetId()] = i
1860                    refAtmSizer.Add(refSel,0,wx.ALIGN_CENTER_VERTICAL)
1861            return refAtmSizer
1862                       
1863        def rbVectMag(rbId,imag,rbData):
1864           
1865            def OnRBVectorMag(event):
1866                event.Skip()
1867                Obj = event.GetEventObject()
1868                rbId,imag = Indx[Obj.GetId()]
1869                try:
1870                    val = float(Obj.GetValue())
1871                    if val <= 0.:
1872                        raise ValueError
1873                    rbData['VectMag'][imag] = val
1874                except ValueError:
1875                    pass
1876                Obj.SetValue('%8.4f'%(val))
1877                wx.CallAfter(UpdateVectorRB,VectorRB.GetScrollPos(wx.VERTICAL))
1878                G2plt.PlotRigidBody(G2frame,'Vector',AtInfo,data['Vector'][rbId],plotDefaults)
1879               
1880            def OnRBVectorRef(event):
1881                Obj = event.GetEventObject()
1882                rbId,imag = Indx[Obj.GetId()]
1883                rbData['VectRef'][imag] = Obj.GetValue()
1884                       
1885            magSizer = wx.BoxSizer(wx.HORIZONTAL)
1886            magSizer.Add(wx.StaticText(VectorRBDisplay,-1,'Translation magnitude: '),
1887                0,wx.ALIGN_CENTER_VERTICAL)
1888            magValue = wx.TextCtrl(VectorRBDisplay,-1,'%8.4f'%(rbData['VectMag'][imag]))
1889            Indx[magValue.GetId()] = [rbId,imag]
1890            magValue.Bind(wx.EVT_TEXT_ENTER,OnRBVectorMag)
1891            magValue.Bind(wx.EVT_KILL_FOCUS,OnRBVectorMag)
1892            magSizer.Add(magValue,0,wx.ALIGN_CENTER_VERTICAL)
1893            magSizer.Add((5,0),)
1894            magref = wx.CheckBox(VectorRBDisplay,-1,label=' Refine?') 
1895            magref.SetValue(rbData['VectRef'][imag])
1896            magref.Bind(wx.EVT_CHECKBOX,OnRBVectorRef)
1897            Indx[magref.GetId()] = [rbId,imag]
1898            magSizer.Add(magref,0,wx.ALIGN_CENTER_VERTICAL)
1899            return magSizer
1900           
1901        def rbVectors(rbId,imag,mag,XYZ,rbData):
1902
1903            def TypeSelect(event):
1904                AtInfo = data['Vector']['AtInfo']
1905                r,c = event.GetRow(),event.GetCol()
1906                if vecGrid.GetColLabelValue(c) == 'Type':
1907                    PE = G2elemGUI.PickElement(G2frame,oneOnly=True)
1908                    if PE.ShowModal() == wx.ID_OK:
1909                        if PE.Elem != 'None':
1910                            El = PE.Elem.strip().lower().capitalize()
1911                            if El not in AtInfo:
1912                                Info = G2elem.GetAtomInfo(El)
1913                                AtInfo[El] = [Info['Drad'],Info['Color']]
1914                            rbData['rbTypes'][r] = El
1915                            vecGrid.SetCellValue(r,c,El)
1916                    PE.Destroy()
1917                wx.CallAfter(UpdateVectorRB,VectorRB.GetScrollPos(wx.VERTICAL))
1918
1919            def ChangeCell(event):
1920                r,c =  event.GetRow(),event.GetCol()
1921                if r >= 0 and (0 <= c < 3):
1922                    try:
1923                        val = float(vecGrid.GetCellValue(r,c))
1924                        rbData['rbVect'][imag][r][c] = val
1925                    except ValueError:
1926                        pass
1927                G2plt.PlotRigidBody(G2frame,'Vector',AtInfo,data['Vector'][rbId],plotDefaults)
1928                wx.CallAfter(UpdateVectorRB,VectorRB.GetScrollPos(wx.VERTICAL))
1929
1930            vecSizer = wx.BoxSizer()
1931            Types = 3*[wg.GRID_VALUE_FLOAT+':10,5',]+[wg.GRID_VALUE_STRING,]+3*[wg.GRID_VALUE_FLOAT+':10,5',]
1932            colLabels = ['Vector x','Vector y','Vector z','Type','Cart x','Cart y','Cart z']
1933            table = []
1934            rowLabels = []
1935            for ivec,xyz in enumerate(rbData['rbVect'][imag]):
1936                table.append(list(xyz)+[rbData['rbTypes'][ivec],]+list(XYZ[ivec]))
1937                rowLabels.append(str(ivec))
1938            vecTable = G2G.Table(table,rowLabels=rowLabels,colLabels=colLabels,types=Types)
1939            vecGrid = G2G.GSGrid(VectorRBDisplay)
1940            vecGrid.SetTable(vecTable, True)
1941            if 'phoenix' in wx.version():
1942                vecGrid.Bind(wg.EVT_GRID_CELL_CHANGED, ChangeCell)
1943            else:
1944                vecGrid.Bind(wg.EVT_GRID_CELL_CHANGE, ChangeCell)
1945            if not imag:
1946                vecGrid.Bind(wg.EVT_GRID_CELL_LEFT_DCLICK, TypeSelect)
1947            attr = wx.grid.GridCellAttr()
1948            attr.SetEditor(G2G.GridFractionEditor(vecGrid))
1949            for c in range(3):
1950                vecGrid.SetColAttr(c, attr)
1951            for row in range(vecTable.GetNumberRows()):
1952                if imag:
1953                    vecGrid.SetCellStyle(row,3,VERY_LIGHT_GREY,True)                   
1954                for col in [4,5,6]:
1955                    vecGrid.SetCellStyle(row,col,VERY_LIGHT_GREY,True)
1956#            vecGrid.SetScrollRate(0,0)
1957            vecGrid.AutoSizeColumns(False)
1958            vecSizer.Add(vecGrid)
1959            return vecSizer
1960       
1961        def FillRefChoice(rbId,rbData):
1962            choiceIds = [i for i in range(len(rbData['rbTypes']))]
1963           
1964            rbRef = rbData.get('rbRef',[-1,-1,-1,False])
1965            for i in range(3):
1966                choiceIds.remove(rbRef[i])
1967            refChoice[rbId] = [choiceIds[:],choiceIds[:],choiceIds[:]]
1968            for i in range(3):
1969                refChoice[rbId][i].append(rbRef[i])
1970                refChoice[rbId][i].sort()     
1971           
1972        if VectorRB.GetSizer(): VectorRB.GetSizer().Clear(True)
1973        VectorRBSizer = wx.BoxSizer(wx.VERTICAL)
1974        for rbId in data['RBIds']['Vector']:
1975            if rbId != 'AtInfo':
1976                rbData = data['Vector'][rbId]
1977                FillRefChoice(rbId,rbData)
1978                VectorRBSizer.Add(rbNameSizer(rbId,rbData),0)
1979                VectorRBSizer.Add(rbRefAtmSizer(rbId,rbData),0)
1980                XYZ = np.array([[0.,0.,0.] for Ty in rbData['rbTypes']])
1981                for imag,mag in enumerate(rbData['VectMag']):
1982                    XYZ += mag*rbData['rbVect'][imag]
1983                    VectorRBSizer.Add(rbVectMag(rbId,imag,rbData),0)
1984                    VectorRBSizer.Add(rbVectors(rbId,imag,mag,XYZ,rbData),0)
1985                VectorRBSizer.Add((5,5),0)
1986                data['Vector'][rbId]['rbXYZ'] = XYZ       
1987        VectorRBSizer.Layout()   
1988        VectorRBDisplay.SetSizer(VectorRBSizer,True)
1989        Size = VectorRBSizer.GetMinSize()
1990        Size[0] += 40
1991        Size[1] = max(Size[1],450) + 20
1992        VectorRBDisplay.SetSize(Size)
1993        VectorRB.SetScrollbars(10,10,Size[0]/10-4,Size[1]/10-1)
1994        VectorRB.Scroll(0,Scroll)
1995       
1996    def UpdateResidueRB(rbId=0):
1997        AtInfo = data['Residue']['AtInfo']
1998        refChoice = {}
1999        RefObjs = []
2000
2001        def rbNameSizer(rbId,rbData):
2002
2003            def OnRBName(event):
2004                Obj = event.GetEventObject()
2005                rbData['RBname'] = Obj.GetValue()
2006                wx.CallAfter(UpdateResidueRB,rbId)
2007               
2008            def OnDelRB(event):
2009                Obj = event.GetEventObject()
2010                rbId = Indx[Obj.GetId()]
2011                if rbId in data['Residue']: 
2012                    del data['Residue'][rbId]
2013                    data['RBIds']['Residue'].remove(rbId)
2014                wx.CallAfter(UpdateResidueRB)
2015               
2016            def OnStripH(event):
2017                Obj = event.GetEventObject()
2018                rbId = Indx[Obj.GetId()]
2019                if rbId in data['Residue']:
2020                    newNames = []
2021                    newTypes = []
2022                    newXYZ = []
2023                    for i,atype in enumerate(rbData['rbTypes']):
2024                        if atype != 'H':
2025                            newNames.append(rbData['atNames'][i])
2026                            newTypes.append(rbData['rbTypes'][i])
2027                            newXYZ.append(rbData['rbXYZ'][i])
2028                    rbData['atNames'] = newNames
2029                    rbData['rbTypes'] = newTypes
2030                    rbData['rbXYZ'] = newXYZ
2031                G2plt.PlotRigidBody(G2frame,'Residue',AtInfo,rbData,plotDefaults)
2032                wx.CallAfter(UpdateResidueRB,rbId)
2033                   
2034            def OnPlotRB(event):
2035                Obj = event.GetEventObject()
2036                Obj.SetValue(False)
2037                G2plt.PlotRigidBody(G2frame,'Residue',AtInfo,rbData,plotDefaults)
2038           
2039            nameSizer = wx.BoxSizer(wx.HORIZONTAL)
2040            nameSizer.Add(wx.StaticText(ResidueRBDisplay,-1,'Residue name: '),
2041                0,wx.ALIGN_CENTER_VERTICAL)
2042            RBname = wx.TextCtrl(ResidueRBDisplay,-1,rbData['RBname'])
2043            Indx[RBname.GetId()] = rbId
2044            RBname.Bind(wx.EVT_TEXT_ENTER,OnRBName)
2045            RBname.Bind(wx.EVT_KILL_FOCUS,OnRBName)
2046            nameSizer.Add(RBname,0,wx.ALIGN_CENTER_VERTICAL)
2047            nameSizer.Add((5,0),)
2048            plotRB = wx.CheckBox(ResidueRBDisplay,-1,'Plot?')
2049            Indx[plotRB.GetId()] = rbId
2050            plotRB.Bind(wx.EVT_CHECKBOX,OnPlotRB)
2051            nameSizer.Add(plotRB,0,wx.ALIGN_CENTER_VERTICAL)
2052            nameSizer.Add((5,0),)
2053            if not rbData['useCount']:
2054                delRB = wx.CheckBox(ResidueRBDisplay,-1,'Delete?')
2055                Indx[delRB.GetId()] = rbId
2056                delRB.Bind(wx.EVT_CHECKBOX,OnDelRB)
2057                nameSizer.Add(delRB,0,wx.ALIGN_CENTER_VERTICAL)
2058                if 'H'  in rbData['rbTypes']:
2059                    stripH = wx.CheckBox(ResidueRBDisplay,-1,'Strip H-atoms?')
2060                    Indx[stripH.GetId()] = rbId
2061                    stripH.Bind(wx.EVT_CHECKBOX,OnStripH)
2062                    nameSizer.Add(stripH,0,wx.ALIGN_CENTER_VERTICAL)
2063            return nameSizer
2064           
2065        def rbResidues(rbId,rbData):
2066           
2067            def TypeSelect(event):
2068                AtInfo = data['Residue']['AtInfo']
2069                r,c = event.GetRow(),event.GetCol()
2070                if vecGrid.GetColLabelValue(c) == 'Type':
2071                    PE = G2elemGUI.PickElement(G2frame,oneOnly=True)
2072                    if PE.ShowModal() == wx.ID_OK:
2073                        if PE.Elem != 'None':
2074                            El = PE.Elem.strip().lower().capitalize()
2075                            if El not in AtInfo:
2076                                Info = G2elem.GetAtomInfo(El)
2077                                AtInfo[El] = [Info['Drad']['Color']]
2078                            rbData['rbTypes'][r] = El
2079                            vecGrid.SetCellValue(r,c,El)
2080                    PE.Destroy()
2081
2082            def ChangeCell(event):
2083                r,c =  event.GetRow(),event.GetCol()
2084                if r >= 0 and (0 <= c < 3):
2085                    try:
2086                        val = float(vecGrid.GetCellValue(r,c))
2087                        rbData['rbXYZ'][r][c] = val
2088                    except ValueError:
2089                        pass
2090                       
2091            def RowSelect(event):
2092                r,c =  event.GetRow(),event.GetCol()
2093                if c < 0:                   #only row clicks
2094                    for vecgrid in resList:
2095                        vecgrid.ClearSelection()
2096                    vecGrid.SelectRow(r,True)
2097
2098            def OnRefSel(event):
2099               
2100                Obj = event.GetEventObject()
2101                iref,res,jref = Indx[Obj.GetId()]
2102                sel = Obj.GetValue()
2103                ind = atNames.index(sel)
2104                if rbData['rbTypes'][ind] == 'H':
2105                    G2G.G2MessageBox(G2frame,'You should not select an H-atom for rigid body orientation')
2106                rbData['rbRef'][iref] = ind
2107                FillRefChoice(rbId,rbData)
2108                for i,ref in enumerate(RefObjs[jref]):
2109                    ref.SetItems([atNames[j] for j in refChoice[rbId][i]])
2110                    ref.SetValue(atNames[rbData['rbRef'][i]])                   
2111                rbXYZ = rbData['rbXYZ']
2112                if not iref:     #origin change
2113                    rbXYZ -= rbXYZ[ind]
2114                #TODO - transform all atom XYZ by axis choices
2115                Xxyz = rbXYZ[rbData['rbRef'][1]]
2116                X = Xxyz/np.sqrt(np.sum(Xxyz**2))
2117                Yxyz = rbXYZ[rbData['rbRef'][2]]
2118                Y = Yxyz/np.sqrt(np.sum(Yxyz**2))
2119                Mat = G2mth.getRBTransMat(X,Y)
2120                rbXYZ = np.inner(Mat,rbXYZ).T
2121                rbData['rbXYZ'] = rbXYZ
2122                res.ClearSelection()
2123                resTable = res.GetTable()
2124                for r in range(res.GetNumberRows()):
2125                    row = resTable.GetRowValues(r)
2126                    row[2:4] = rbXYZ[r]
2127                    resTable.SetRowValues(r,row)
2128                res.ForceRefresh()
2129                G2plt.PlotRigidBody(G2frame,'Residue',AtInfo,rbData,plotDefaults)
2130               
2131            Types = 2*[wg.GRID_VALUE_STRING,]+3*[wg.GRID_VALUE_FLOAT+':10,5',]
2132            colLabels = ['Name','Type','Cart x','Cart y','Cart z']
2133            table = []
2134            rowLabels = []
2135            for ivec,xyz in enumerate(rbData['rbXYZ']):
2136                table.append([rbData['atNames'][ivec],]+[rbData['rbTypes'][ivec],]+list(xyz))
2137                rowLabels.append(str(ivec))
2138            vecTable = G2G.Table(table,rowLabels=rowLabels,colLabels=colLabels,types=Types)
2139            vecGrid = G2G.GSGrid(ResidueRBDisplay)
2140            Indx[vecGrid.GetId()] = rbId
2141            resList.append(vecGrid)
2142            vecGrid.SetTable(vecTable, True)
2143            if 'phoenix' in wx.version():
2144                vecGrid.Bind(wg.EVT_GRID_CELL_CHANGED, ChangeCell)
2145            else:
2146                vecGrid.Bind(wg.EVT_GRID_CELL_CHANGE, ChangeCell)
2147            vecGrid.Bind(wg.EVT_GRID_CELL_LEFT_DCLICK, TypeSelect)
2148            vecGrid.Bind(wg.EVT_GRID_LABEL_LEFT_CLICK, RowSelect)
2149            attr = wx.grid.GridCellAttr()
2150            attr.SetEditor(G2G.GridFractionEditor(vecGrid))
2151            for c in range(3):
2152                vecGrid.SetColAttr(c, attr)
2153            for row in range(vecTable.GetNumberRows()):
2154                for col in range(5):
2155                    vecGrid.SetCellStyle(row,col,VERY_LIGHT_GREY,True)
2156            vecGrid.AutoSizeColumns(False)
2157            vecSizer = wx.BoxSizer()
2158            vecSizer.Add(vecGrid)
2159           
2160            refAtmSizer = wx.BoxSizer(wx.HORIZONTAL)
2161            atNames = rbData['atNames']
2162            rbRef = rbData['rbRef']
2163            if rbData['rbRef'][3] or rbData['useCount']:
2164                refAtmSizer.Add(wx.StaticText(ResidueRBDisplay,-1,
2165                    'Orientation reference non-H atoms A-B-C: %s, %s, %s'%(atNames[rbRef[0]], \
2166                     atNames[rbRef[1]],atNames[rbRef[2]])),0)
2167            else:
2168                refAtmSizer.Add(wx.StaticText(ResidueRBDisplay,-1,
2169                    'Orientation reference non-H atoms A-B-C: '),0,wx.ALIGN_CENTER_VERTICAL)
2170                refObj = [0,0,0]
2171                for i in range(3):
2172                    choices = [atNames[j] for j in refChoice[rbId][i]]
2173                    refSel = wx.ComboBox(ResidueRBDisplay,-1,value='',
2174                        choices=choices,style=wx.CB_READONLY|wx.CB_DROPDOWN)
2175                    refSel.SetValue(atNames[rbRef[i]])
2176                    refSel.Bind(wx.EVT_COMBOBOX, OnRefSel)
2177                    Indx[refSel.GetId()] = [i,vecGrid,len(RefObjs)]
2178                    refObj[i] = refSel
2179                    refAtmSizer.Add(refSel,0,wx.ALIGN_CENTER_VERTICAL)
2180                RefObjs.append(refObj)
2181           
2182            mainSizer = wx.BoxSizer(wx.VERTICAL)
2183            mainSizer.Add(refAtmSizer)
2184            mainSizer.Add(vecSizer)
2185            return mainSizer
2186           
2187        def SeqSizer(angSlide,rbId,iSeq,Seq,atNames):
2188           
2189            def ChangeAngle(event):
2190                event.Skip()
2191                Obj = event.GetEventObject()
2192                rbId,Seq = Indx[Obj.GetId()][:2]
2193                val = Seq[2]
2194                try:
2195                    val = float(Obj.GetValue())
2196                    Seq[2] = val
2197                except ValueError:
2198                    pass
2199                Obj.SetValue('%8.2f'%(val))
2200                G2plt.PlotRigidBody(G2frame,'Residue',AtInfo,data['Residue'][rbId],plotDefaults)
2201               
2202            def OnRadBtn(event):
2203                Obj = event.GetEventObject()
2204                Seq,iSeq,angId = Indx[Obj.GetId()]
2205                data['Residue'][rbId]['SelSeq'] = [iSeq,angId]
2206                angSlide.SetValue(int(100*Seq[2]))
2207               
2208            def OnDelBtn(event):
2209                Obj = event.GetEventObject()
2210                rbId,Seq = Indx[Obj.GetId()]
2211                data['Residue'][rbId]['rbSeq'].remove(Seq)       
2212                wx.CallAfter(UpdateResidueRB,rbId)
2213           
2214            seqSizer = wx.FlexGridSizer(0,5,2,2)
2215            seqSizer.AddGrowableCol(3,0)
2216            iBeg,iFin,angle,iMove = Seq
2217            ang = wx.TextCtrl(ResidueRBDisplay,-1,'%8.2f'%(angle),size=(50,20))
2218            if not iSeq:
2219                radBt = wx.RadioButton(ResidueRBDisplay,-1,'',style=wx.RB_GROUP)
2220                data['Residue'][rbId]['SelSeq'] = [iSeq,ang.GetId()]
2221            else:
2222                radBt = wx.RadioButton(ResidueRBDisplay,-1,'')
2223            radBt.Bind(wx.EVT_RADIOBUTTON,OnRadBtn)                   
2224            seqSizer.Add(radBt)
2225            delBt = wx.RadioButton(ResidueRBDisplay,-1,'')
2226            delBt.Bind(wx.EVT_RADIOBUTTON,OnDelBtn)
2227            seqSizer.Add(delBt)
2228            bond = wx.TextCtrl(ResidueRBDisplay,-1,'%s %s'%(atNames[iBeg],atNames[iFin]),size=(50,20))
2229            seqSizer.Add(bond,0,wx.ALIGN_CENTER_VERTICAL)
2230            Indx[radBt.GetId()] = [Seq,iSeq,ang.GetId()]
2231            Indx[delBt.GetId()] = [rbId,Seq]
2232            Indx[ang.GetId()] = [rbId,Seq,ang]
2233            ang.Bind(wx.EVT_TEXT_ENTER,ChangeAngle)
2234            ang.Bind(wx.EVT_KILL_FOCUS,ChangeAngle)
2235            seqSizer.Add(ang,0,wx.ALIGN_CENTER_VERTICAL)
2236            atms = ''
2237            for i in iMove:   
2238                atms += ' %s,'%(atNames[i])
2239            moves = wx.TextCtrl(ResidueRBDisplay,-1,atms[:-1],size=(200,20))
2240            seqSizer.Add(moves,1,wx.ALIGN_CENTER_VERTICAL|wx.EXPAND|wx.RIGHT)
2241            return seqSizer
2242           
2243        def SlideSizer():
2244           
2245            def OnSlider(event):
2246                Obj = event.GetEventObject()
2247                rbData = Indx[Obj.GetId()]
2248                iSeq,angId = rbData['SelSeq']
2249                val = float(Obj.GetValue())/100.
2250                rbData['rbSeq'][iSeq][2] = val
2251                Indx[angId][2].SetValue('%8.2f'%(val))
2252                G2plt.PlotRigidBody(G2frame,'Residue',AtInfo,rbData,plotDefaults)
2253           
2254            slideSizer = wx.BoxSizer(wx.HORIZONTAL)
2255            slideSizer.Add(wx.StaticText(ResidueRBDisplay,-1,'Selected torsion angle:'),0)
2256            iSeq,angId = rbData['SelSeq']
2257            angSlide = wx.Slider(ResidueRBDisplay,-1,
2258                int(100*rbData['rbSeq'][iSeq][2]),0,36000,size=(200,20),
2259                style=wx.SL_HORIZONTAL)
2260            angSlide.Bind(wx.EVT_SLIDER, OnSlider)
2261            Indx[angSlide.GetId()] = rbData
2262            slideSizer.Add(angSlide,0)           
2263            return slideSizer,angSlide
2264           
2265        def FillRefChoice(rbId,rbData):
2266            choiceIds = [i for i in range(len(rbData['atNames']))]
2267            for seq in rbData['rbSeq']:
2268                for i in seq[3]:
2269                    try:
2270                        choiceIds.remove(i)
2271                    except ValueError:
2272                        pass
2273            rbRef = rbData['rbRef']
2274            for i in range(3):
2275                try:
2276                    choiceIds.remove(rbRef[i])
2277                except ValueError:
2278                    pass
2279            refChoice[rbId] = [choiceIds[:],choiceIds[:],choiceIds[:]]
2280            for i in range(3):
2281                refChoice[rbId][i].append(rbRef[i])
2282                refChoice[rbId][i].sort()
2283               
2284        def OnSelect(event):
2285            rbname = rbchoice[select.GetSelection()]
2286            rbId = RBnames[rbname]
2287            wx.CallLater(100,UpdateResidueRB,rbId)
2288           
2289        GS = ResidueRBDisplay.GetSizer()
2290        if GS: 
2291            try:        #get around a c++ error in wx 4.0; doing is again seems to be OK
2292                GS.Clear(True)
2293            except:
2294                GS.Clear(True)
2295       
2296        RBnames = {}
2297        for rbid in data['RBIds']['Residue']:
2298            RBnames.update({data['Residue'][rbid]['RBname']:rbid,})
2299        if not RBnames:
2300            return
2301        rbchoice = list(RBnames.keys())
2302        ResidueRBSizer = wx.BoxSizer(wx.VERTICAL)
2303        if len(RBnames) > 1:
2304            selSizer = wx.BoxSizer(wx.HORIZONTAL)
2305            selSizer.Add(wx.StaticText(ResidueRBDisplay,label=' Select residue to view:'),0)
2306            rbchoice.sort()
2307            select = wx.ComboBox(ResidueRBDisplay,choices=rbchoice)
2308            select.Bind(wx.EVT_COMBOBOX,OnSelect)
2309            selSizer.Add(select,0)
2310            ResidueRBSizer.Add(selSizer,0)
2311        if not rbId:
2312            rbId = RBnames[rbchoice[0]]
2313        rbData = data['Residue'][rbId]
2314        FillRefChoice(rbId,rbData)
2315        ResidueRBSizer.Add(rbNameSizer(rbId,rbData),0)
2316        ResidueRBSizer.Add(rbResidues(rbId,rbData),0)
2317        ResidueRBSizer.Add((5,5),0)
2318        if rbData['rbSeq']:
2319            slideSizer,angSlide = SlideSizer()
2320        if len(rbData['rbSeq']):
2321            ResidueRBSizer.Add(wx.StaticText(ResidueRBDisplay,-1,
2322                'Sel  Del  Bond             Angle      Riding atoms'),
2323                0,wx.ALIGN_CENTER_VERTICAL)                       
2324        for iSeq,Seq in enumerate(rbData['rbSeq']):
2325            ResidueRBSizer.Add(SeqSizer(angSlide,rbId,iSeq,Seq,rbData['atNames']))
2326        if rbData['rbSeq']:
2327            ResidueRBSizer.Add(slideSizer,)
2328
2329        ResidueRBSizer.Add((5,25),)
2330        ResidueRBSizer.Layout()   
2331        ResidueRBDisplay.SetSizer(ResidueRBSizer,True)
2332        ResidueRBDisplay.SetAutoLayout(True)
2333        Size = ResidueRBSizer.GetMinSize()
2334        ResidueRBDisplay.SetSize(Size)
2335        ResidueRBDisplay.Show()
2336       
2337    def SetStatusLine(text):
2338        G2frame.GetStatusBar().SetStatusText(text,1)                                     
2339
2340    G2gd.SetDataMenuBar(G2frame,G2frame.dataWindow.RigidBodyMenu)
2341    SetStatusLine('')
2342    UpdateVectorRB()
2343    G2frame.rbBook.Bind(wx.aui.EVT_AUINOTEBOOK_PAGE_CHANGED, OnPageChanged)
2344    wx.CallAfter(OnPageChanged,None)
Note: See TracBrowser for help on using the repository browser.