source: trunk/GSASIIconstrGUI.py @ 3904

Last change on this file since 3904 was 3904, checked in by vondreele, 3 years ago

fix lattice parm constraint bug

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