Changeset 934


Ignore:
Timestamp:
May 28, 2013 4:29:16 PM (9 years ago)
Author:
vondreele
Message:

changes for MC/SA

Location:
trunk
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIconstrGUI.py

    r929 r934  
    969969            rbData['rbSeq'].append([Orig,Piv,0.0,Riding])           
    970970        dlg.Destroy()
    971         UpdateResidueRB()       
     971        UpdateResidueRB()
    972972
    973973    def UpdateVectorRB(Scroll=0):
     
    13541354                data['Residue'][rbId]['SelSeq'] = [iSeq,angId]
    13551355                angSlide.SetValue(int(100*Seq[2]))
    1356            
    1357             seqSizer = wx.FlexGridSizer(0,4,2,2)
     1356               
     1357            def OnDelBtn(event):
     1358                Obj = event.GetEventObject()
     1359                rbId,Seq = Indx[Obj.GetId()]
     1360                data['Residue'][rbId]['rbSeq'].remove(Seq)       
     1361                wx.CallAfter(UpdateResidueRB)
     1362           
     1363            seqSizer = wx.FlexGridSizer(0,5,2,2)
    13581364            seqSizer.AddGrowableCol(3,0)
    13591365            iBeg,iFin,angle,iMove = Seq
     
    13661372            radBt.Bind(wx.EVT_RADIOBUTTON,OnRadBtn)                   
    13671373            seqSizer.Add(radBt)
     1374            delBt = wx.RadioButton(ResidueRBDisplay,-1,'')
     1375            delBt.Bind(wx.EVT_RADIOBUTTON,OnDelBtn)
     1376            seqSizer.Add(delBt)
    13681377            bond = wx.TextCtrl(ResidueRBDisplay,-1,'%s %s'%(atNames[iBeg],atNames[iFin]),size=(50,20))
    13691378            seqSizer.Add(bond,0,wx.ALIGN_CENTER_VERTICAL)
    13701379            Indx[radBt.GetId()] = [Seq,iSeq,ang.GetId()]
     1380            Indx[delBt.GetId()] = [rbId,Seq]
    13711381            Indx[ang.GetId()] = [rbId,Seq,ang]
    13721382            ang.Bind(wx.EVT_TEXT_ENTER,ChangeAngle)
     
    14341444            if len(rbData['rbSeq']):
    14351445                ResidueRBSizer.Add(wx.StaticText(ResidueRBDisplay,-1,
    1436                     'Sel  Bond             Angle      Riding atoms'),
     1446                    'Sel  Del  Bond             Angle      Riding atoms'),
    14371447                    0,wx.ALIGN_CENTER_VERTICAL)                       
    14381448            for iSeq,Seq in enumerate(rbData['rbSeq']):
  • trunk/GSASIIgrid.py

    r933 r934  
    5959[ wxID_FOURCALC, wxID_FOURSEARCH, wxID_FOURCLEAR, wxID_PEAKSMOVE, wxID_PEAKSCLEAR,
    6060    wxID_CHARGEFLIP, wxID_PEAKSUNIQUE, wxID_PEAKSDELETE, wxID_PEAKSDA,
    61     wxID_PEAKSDISTVP, wxID_PEAKSVIEWPT, wxID_FINDEQVPEAKS,wxID_SHOWBONDS,
    62 ] = [wx.NewId() for item in range(13)]
     61    wxID_PEAKSDISTVP, wxID_PEAKSVIEWPT, wxID_FINDEQVPEAKS,wxID_SHOWBONDS,wxID_RUNMCSA,
     62] = [wx.NewId() for item in range(14)]
    6363
    6464[ wxID_PWDRADD, wxID_HKLFADD,wxID_PWDANALYSIS,wxID_DATADELETE,
     
    7979[ wxID_DRAWRESTRBOND, wxID_DRAWRESTRANGLE, wxID_DRAWRESTRPLANE, wxID_DRAWRESTRCHIRAL,
    8080] = [wx.NewId() for item in range(4)]
     81
     82[ wxID_ADDMCSAATOM,wxID_ADDMCSARB,wxID_CLEARMCSARB,
     83] = [wx.NewId() for item in range(3)]
    8184
    8285[ wxID_CLEARTEXTURE,wxID_REFINETEXTURE,
     
    16501653        self.GeneralCalc.Append(help='Clear map',id=wxID_FOURCLEAR, kind=wx.ITEM_NORMAL,
    16511654            text='Clear map')
     1655        self.GeneralCalc.Append(help='Run Monte Carlo - Simulated Annealing',id=wxID_RUNMCSA, kind=wx.ITEM_NORMAL,
     1656            text='Run MC/SA')
    16521657        self.PostfillDataMenu()
    16531658       
     
    17671772        self.DrawAtomRigidBody.Append(id=wxID_DRAWDEFINERB, kind=wx.ITEM_NORMAL,text='Define rigid body',
    17681773            help='Define rigid body with selected atoms')
     1774        self.PostfillDataMenu()
     1775
     1776        # Phase / MCSA tab
     1777        self.MCSAMenu = wx.MenuBar()
     1778        self.PrefillDataMenu(self.MCSAMenu,helpType='MC/SA')
     1779        self.MCSAEdit = wx.Menu(title='')
     1780        self.MCSAMenu.Append(menu=self.MCSAEdit, title='MC/SA')
     1781        self.MCSAEdit.Append(id=wxID_ADDMCSAATOM, kind=wx.ITEM_NORMAL,text='Add atom',
     1782            help='Add single atom to MC/SA model')
     1783        self.MCSAEdit.Append(id=wxID_ADDMCSARB, kind=wx.ITEM_NORMAL,text='Add rigid body',
     1784            help='Add rigid body to MC/SA model' )
     1785        self.MCSAEdit.Append(id=wxID_CLEARMCSARB, kind=wx.ITEM_NORMAL,text='Clear rigid bodies',
     1786            help='Clear all rigid bodies from MC/SA model' )
    17691787        self.PostfillDataMenu()
    17701788           
  • trunk/GSASIImath.py

    r915 r934  
    282282        XYZ[i] = np.inner(Bmat,X)+RBObj['Orig'][0]
    283283    return XYZ,Cart
     284
     285def UpdateMCSAxyz(Bmat,mcsaModels,RBData):
     286    xyz = []
     287    atTypes = []
     288    iatm = 0
     289    for model in mcsaModels[1:]:        #skip the MD model
     290        if model['Type'] == 'Atom':
     291            xyz.append(model['Pos'][0])
     292            atTypes.append(model['atType'])
     293            iatm += 1
     294        else:
     295            rideList = []
     296            RBRes = RBData[model['Type']][model['RBId']]
     297            Pos = np.array(model['Pos'][0])
     298            Qori = np.array(model['Ori'][0])
     299            if model['Type'] == 'Vector':
     300                vecs = RBRes['rbVect']
     301                mags = RBRes['VectMag']
     302                Cart = np.zeros_like(vecs[0])
     303                for vec,mag in zip(vecs,mags):
     304                    Cart += vec*mag
     305            elif model['Type'] == 'Residue':
     306                Cart = np.array(RBRes['rbXYZ'])
     307                for itor,seq in enumerate(RBRes['rbSeq']):
     308                    QuatA = AVdeg2Q(model['Tor'][0][itor],Cart[seq[0]]-Cart[seq[1]])
     309                    for ride in seq[3]:
     310                        Cart[ride] = prodQVQ(QuatA,Cart[ride]-Cart[seq[1]])+Cart[seq[1]]
     311                    rideList += seq[3]
     312            centList = set(range(len(Cart)))-set(rideList)
     313            if model['MolCent']:
     314                cent = np.zeros(3)
     315                for i in centList:
     316                    cent += Cart[i]
     317                Cart -= cent/len(centList)
     318            for i,x in enumerate(Cart):
     319                xyz.append(np.inner(Bmat,prodQVQ(Qori,x))+Pos)
     320                atType = RBRes['rbTypes'][i]
     321                atTypes.append(atType)
     322                iatm += 1
     323    return np.array(xyz),atTypes
    284324   
    285325def UpdateRBUIJ(Bmat,Cart,RBObj):
     
    13601400        XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
    13611401    return XY
     1402
     1403def mcsaSearch(data,reflType,reflData,covData,pgbar):
     1404    generalData = data['General']
     1405    mcsaControls = generalData['MCSA controls']
     1406    reflName = mcsaData['Data source']
     1407    phaseName = generalData['Name']
     1408    MCSAObjs = data['MCSAModels']
     1409           
     1410#             = {'':'','Annealing':[50.0,0.001,0.90,1000],
     1411#            'dmin':2.0,'Algolrithm':'Normal','Jump coeff':[0.95,0.5]} #'Normal','Random jump','Tremayne jump'
     1412
     1413    return {}
     1414
    13621415       
    13631416def getWave(Parms):
     
    14401493        V /= d
    14411494    else:
    1442         return [0.,0.,0.,1.]    #identity
     1495        V = np.array([0.,0.,1.])
    14431496    p = A/2.
    14441497    Q[0] = np.cos(p)
     
    14551508        V /= d
    14561509    else:
    1457         return [0.,0.,0.,1.]    #identity
     1510        V = np.array([0.,0.,1.])
    14581511    p = A/2.
    14591512    Q[0] = cosd(p)
     
    14671520    A = 2.*acosd(Q[0])
    14681521    V = np.array(Q[1:])
    1469     if nl.norm(Q[1:]):
    1470         V = Q[1:]/nl.norm(Q[1:])
     1522    d = np.sqrt(np.sum(V**2))
     1523    if d:
     1524        V /= d
    14711525    else:
    14721526        A = 0.
    1473         V = np.array([0.,0.,1.])
     1527        V = np.array([0.,0.,0.])
    14741528    return A,V
    14751529   
     
    14801534    A = 2.*np.arccos(Q[0])
    14811535    V = np.array(Q[1:])
    1482     if nl.norm(Q[1:]):
    1483         V = Q[1:]/nl.norm(Q[1:])
     1536    d = np.sqrt(np.sum(V**2))
     1537    if d:
     1538        V /= d
    14841539    else:
    14851540        A = 0.
    1486         V = np.array([0.,0.,1.])
     1541        V = np.array([0.,0.,0.])
    14871542    return A,V
    14881543   
  • trunk/GSASIIphsGUI.py

    r927 r934  
    8383    if 'RBModels' not in data:
    8484        data['RBModels'] = {}
     85    if 'MCSA' not in data:
     86        data['MCSA'] = {'Models':[{'Type':'MD','Coef':[1.0,False,[0.,3.],],'axis':[0,0,1]}],'Results':[],'AtInfo':{}}
    8587#end patch   
    8688
     
    120122        if 'MCSA controls' not in generalData:
    121123            generalData['MCSA controls'] = {'Data source':'','Annealing':[50.0,0.001,0.90,1000],
    122             'dmin':2.0,'Algolrithm':'Normal','Jump coeff':[0.95,0.5]} #'Normal','Random jump','Tremayne jump'
     124            'dmin':2.0,'Algolrithm':'Random jump','Jump coeff':[0.95,0.5],'nRuns':1} #'Random jump','Tremayne jump'
    123125# end of patches
    124126        generalData['NoAtoms'] = {}
     
    649651                    pass
    650652                dmin.SetValue("%.3f"%(MCSA['dmin']))          #reset in case of error
     653               
     654            def OnNoRuns(event):
     655                MCSA['nRuns'] = int(noRuns.GetValue())
    651656           
    652657            def OnAlist(event):
     
    690695                refList = ['Pawley reflections']
    691696            for item in data['Histograms'].keys():
    692                 if 'HKLF' in item:
     697                if 'HKLF' in item or 'PWDR' in item:
    693698                    refList.append(item)
    694699            mcsaSizer = wx.BoxSizer(wx.VERTICAL)
     
    707712            mcsaSizer.Add((5,5),)
    708713            line2Sizer = wx.BoxSizer(wx.HORIZONTAL)
    709             Achoice = ['Normal','Random jump','Tremayne jump']
     714            Rchoice = ['1','2','3','5','10','15','20']
     715            line2Sizer.Add(wx.StaticText(General,label=' No. MC/SA runs: '),0,wx.ALIGN_CENTER_VERTICAL)
     716            noRuns = wx.ComboBox(General,-1,value=str(MCSA.get('nRuns',1)),choices=Rchoice,
     717                style=wx.CB_READONLY|wx.CB_DROPDOWN)
     718            noRuns.Bind(wx.EVT_COMBOBOX,OnNoRuns)
     719            line2Sizer.Add(noRuns,0,wx.ALIGN_CENTER_VERTICAL)
     720            Achoice = ['Random jump','Tremayne jump']
    710721            line2Sizer.Add(wx.StaticText(General,label=' MC/SA algorithm: '),0,wx.ALIGN_CENTER_VERTICAL)
    711722            Alist = wx.ComboBox(General,-1,value=MCSA['Algolrithm'],choices=Achoice,
     
    12541265            print '**** ERROR - no rigid bodies defined ****'
    12551266            return
     1267        dlg = wx.SingleChoiceDialog(G2frame.dataFrame,'Select','Rigid body',rbNames.keys())
     1268        if dlg.ShowModal() == wx.ID_OK:
     1269            sel = dlg.GetSelection()
     1270            rbname = rbNames.keys()[sel]
     1271            rbType,rbId = rbNames[rbname]
     1272            RB = rbData[rbType][rbId]
    12561273       
    12571274    def OnAtomMove(event):
     
    27832800        if True: # Bob wants this tab to always resize -- let's try that.
    27842801            Size = mainSizer.Fit(G2frame.dataFrame)
    2785             Size[0] = max(Size[0]+35,500)           # leave some extra room and don't get too small
     2802  #          Size[0] = max(Size[0]+35,400)           # leave some extra room and don't get too small
     2803            Size[0] += 35                           #compensate for scroll bar
    27862804            Size[1] = max(Size[1]+35,350)                           #compensate for status bar
    27872805            drawOptions.SetScrollbars(10,10,Size[0]/10-4,Size[1]/10-1)
     
    39173935                wx.EndBusyCursor()
    39183936            FillRigidBodyGrid()
     3937           
     3938################################################################################
     3939##### MC/SA routines
     3940################################################################################
     3941
     3942    def UpdateMCSA(refresh=True):
     3943        Indx = {}
     3944       
     3945        def OnPosRef(event):
     3946            Obj = event.GetEventObject()
     3947            model,item,ix = Indx[Obj.GetId()]
     3948            model[item][1][ix] = Obj.GetValue()
     3949           
     3950        def OnPosVal(event):
     3951            Obj = event.GetEventObject()
     3952            model,item,ix = Indx[Obj.GetId()]
     3953            try:
     3954                model[item][0][ix] = float(Obj.GetValue())
     3955            except ValueError:
     3956                pass
     3957            Obj.SetValue("%.4f"%(model[item][0][ix]))
     3958            G2plt.PlotStructure(G2frame,data)
     3959           
     3960        def OnPosRange(event):
     3961            Obj = event.GetEventObject()
     3962            model,item,ix = Indx[Obj.GetId()]
     3963            Range = Obj.GetValue().split()
     3964            try:
     3965                rmin,rmax = [float(Range[i]) for i in range(2)]
     3966                if rmin >= rmax:
     3967                    raise ValueError
     3968            except (ValueError,IndexError):
     3969                rmin,rmax = model[item][2][ix]
     3970            model[item][2][ix] = [rmin,rmax]
     3971            Obj.SetValue('%.3f %.3f'%(rmin,rmax))                 
     3972               
     3973        def atomSizer(model):
     3974           
     3975            atomsizer = wx.FlexGridSizer(1,7,5,5)
     3976            atomsizer.Add(wx.StaticText(MCSA,-1,' Atom: '+model['name']+': '),0,wx.ALIGN_CENTER_VERTICAL)
     3977            for ix,item in enumerate(['x','y','z']):
     3978                posRef = wx.CheckBox(MCSA,-1,label=item+': ')
     3979                posRef.SetValue(model['Pos'][1][ix])
     3980                posRef.Bind(wx.EVT_CHECKBOX,OnPosRef)
     3981                Indx[posRef.GetId()] = [model,'Pos',ix]
     3982                atomsizer.Add(posRef,0,wx.ALIGN_CENTER_VERTICAL)
     3983                posVal = wx.TextCtrl(MCSA,-1,'%.4f'%(model['Pos'][0][ix]),style=wx.TE_PROCESS_ENTER)
     3984                posVal.Bind(wx.EVT_TEXT_ENTER,OnPosVal)
     3985                posVal.Bind(wx.EVT_KILL_FOCUS,OnPosVal)
     3986                Indx[posVal.GetId()] = [model,'Pos',ix]
     3987                atomsizer.Add(posVal,0,wx.ALIGN_CENTER_VERTICAL)
     3988            atomsizer.Add((5,5),0)
     3989            for ix,item in enumerate(['x','y','z']):
     3990                atomsizer.Add(wx.StaticText(MCSA,-1,' Range: '),0,wx.ALIGN_CENTER_VERTICAL)
     3991                rmin,rmax = model['Pos'][2][ix]
     3992                posRange = wx.TextCtrl(MCSA,-1,'%.3f %.3f'%(rmin,rmax),style=wx.TE_PROCESS_ENTER)
     3993                Indx[posRange.GetId()] = [model,'Pos',ix]
     3994                posRange.Bind(wx.EVT_TEXT_ENTER,OnPosRange)
     3995                posRange.Bind(wx.EVT_KILL_FOCUS,OnPosRange)
     3996                atomsizer.Add(posRange,0,wx.ALIGN_CENTER_VERTICAL)
     3997            return atomsizer
     3998           
     3999        def rbSizer(model):
     4000           
     4001            def OnOrVar(event):
     4002                Obj = event.GetEventObject()
     4003                model = Indx[Obj.GetId()]
     4004                model['Ovar'] = Obj.GetValue()
     4005           
     4006            def OnOriVal(event):
     4007                Obj = event.GetEventObject()
     4008                model,ix = Indx[Obj.GetId()]
     4009                A,V = G2mth.Q2AVdeg(model['Ori'][0])
     4010#                V = np.inner(Bmat,V)
     4011                if ix:
     4012                    Anew = A
     4013                    Vec = Obj.GetValue().split()
     4014                    try:
     4015                        Vnew = [float(Vec[i]) for i in range(3)]
     4016                    except ValueError:
     4017                        Vnew = V
     4018                else:
     4019                    Vnew = V
     4020                    try:
     4021                        Anew = float(Obj.GetValue())
     4022                    except ValueError:
     4023                        Anew = A
     4024#                V = np.inner(Amat,V)
     4025                Q = G2mth.AVdeg2Q(Anew,Vnew)
     4026                model['Ori'][0] = Q
     4027                G2plt.PlotStructure(G2frame,data)
     4028                UpdateMCSA()
     4029
     4030            def OnMolCent(event):
     4031                Obj = event.GetEventObject()
     4032                model = Indx[Obj.GetId()]
     4033                model['MolCent'] = Obj.GetValue()
     4034                G2plt.PlotStructure(G2frame,data)
     4035           
     4036            rbsizer = wx.BoxSizer(wx.VERTICAL)
     4037            rbsizer1 = wx.FlexGridSizer(1,7,5,5)
     4038            rbsizer1.Add(wx.StaticText(MCSA,-1,model['Type']+': '+model['name']+': '),0,wx.ALIGN_CENTER_VERTICAL)
     4039            for ix,item in enumerate(['x','y','z']):
     4040                posRef = wx.CheckBox(MCSA,-1,label=item+': ')
     4041                posRef.SetValue(model['Pos'][1][ix])
     4042                posRef.Bind(wx.EVT_CHECKBOX,OnPosRef)
     4043                Indx[posRef.GetId()] = [model,'Pos',ix]
     4044                rbsizer1.Add(posRef,0,wx.ALIGN_CENTER_VERTICAL)
     4045                posVal = wx.TextCtrl(MCSA,-1,'%.4f'%(model['Pos'][0][ix]),style=wx.TE_PROCESS_ENTER)
     4046                posVal.Bind(wx.EVT_TEXT_ENTER,OnPosVal)
     4047                posVal.Bind(wx.EVT_KILL_FOCUS,OnPosVal)
     4048                Indx[posVal.GetId()] = [model,'Pos',ix]
     4049                rbsizer1.Add(posVal,0,wx.ALIGN_CENTER_VERTICAL)
     4050            molcent = wx.CheckBox(MCSA,-1,label=' Use mol. center? ')
     4051            molcent.SetValue(model['MolCent'])
     4052            molcent.Bind(wx.EVT_CHECKBOX,OnMolCent)
     4053            Indx[molcent.GetId()] = model
     4054            rbsizer1.Add(molcent,0,wx.ALIGN_CENTER_VERTICAL)
     4055            for ix,item in enumerate(['x','y','z']):
     4056                rbsizer1.Add(wx.StaticText(MCSA,-1,' Range: '),0,wx.ALIGN_CENTER_VERTICAL)
     4057                rmin,rmax = model['Pos'][2][ix]
     4058                posRange = wx.TextCtrl(MCSA,-1,'%.3f %.3f'%(rmin,rmax),style=wx.TE_PROCESS_ENTER)
     4059                Indx[posRange.GetId()] = [model,'Pos',ix]
     4060                posRange.Bind(wx.EVT_TEXT_ENTER,OnPosRange)
     4061                posRange.Bind(wx.EVT_KILL_FOCUS,OnPosRange)
     4062                rbsizer1.Add(posRange,0,wx.ALIGN_CENTER_VERTICAL)
     4063               
     4064            rbsizer2 = wx.FlexGridSizer(1,6,5,5)
     4065            Orien,OrienV = G2mth.Q2AVdeg(model['Ori'][0])
     4066            Ori = [Orien,]+list(OrienV)
     4067            rbsizer2.Add(wx.StaticText(MCSA,-1,'Oa: '),0,wx.ALIGN_CENTER_VERTICAL)
     4068            angVal = wx.TextCtrl(MCSA,-1,'%.5f'%(Ori[0]),style=wx.TE_PROCESS_ENTER)
     4069            angVal.Bind(wx.EVT_TEXT_ENTER,OnOriVal)
     4070            angVal.Bind(wx.EVT_KILL_FOCUS,OnOriVal)
     4071            Indx[angVal.GetId()] = [model,0]
     4072            rbsizer2.Add(angVal,0,wx.ALIGN_CENTER_VERTICAL)
     4073            rbsizer2.Add(wx.StaticText(MCSA,-1,'Oi,Oj,Ok: '),0,wx.ALIGN_CENTER_VERTICAL)
     4074            vecVal = wx.TextCtrl(MCSA,-1,'%.3f %.3f %.3f'%(Ori[1],Ori[2],Ori[3]),style=wx.TE_PROCESS_ENTER)
     4075            vecVal.Bind(wx.EVT_TEXT_ENTER,OnOriVal)
     4076            vecVal.Bind(wx.EVT_KILL_FOCUS,OnOriVal)
     4077            Indx[vecVal.GetId()] = [model,1]
     4078            rbsizer2.Add(vecVal,0,wx.ALIGN_CENTER_VERTICAL)
     4079            rbsizer2.Add(wx.StaticText(MCSA,-1,' Vary? '),0,wx.ALIGN_CENTER_VERTICAL)
     4080            choice = [' ','A','AV']
     4081            orvar = wx.ComboBox(MCSA,-1,value=model['Ovar'],choices=choice,
     4082                style=wx.CB_READONLY|wx.CB_DROPDOWN)
     4083            orvar.Bind(wx.EVT_COMBOBOX, OnOrVar)
     4084            Indx[orvar.GetId()] = model
     4085            rbsizer2.Add(orvar,0,wx.ALIGN_CENTER_VERTICAL)
     4086            rbsizer2.Add(wx.StaticText(MCSA,-1,' Range: Oa: '),0,wx.ALIGN_CENTER_VERTICAL)
     4087            Rge = model['Ori'][2]
     4088            angRange = wx.TextCtrl(MCSA,-1,'%.3f %.3f'%(Rge[0][0],Rge[0][1]),style=wx.TE_PROCESS_ENTER)
     4089            Indx[angRange.GetId()] = [model,'Ori',0]
     4090            angRange.Bind(wx.EVT_TEXT_ENTER,OnPosRange)
     4091            angRange.Bind(wx.EVT_KILL_FOCUS,OnPosRange)
     4092            rbsizer2.Add(angRange,0,wx.ALIGN_CENTER_VERTICAL)
     4093            rbsizer2.Add(wx.StaticText(MCSA,-1,'Oi,Oj,Ok: '),0,wx.ALIGN_CENTER_VERTICAL)
     4094            for io,item in enumerate(['Oi','Oj','Ok']):
     4095                rmin,rmax = Rge[io+1]
     4096                vecRange = wx.TextCtrl(MCSA,-1,'%.3f %.3f '%(rmin,rmax),style=wx.TE_PROCESS_ENTER)
     4097                Indx[vecRange.GetId()] = [model,'Ori',io+1]
     4098                vecRange.Bind(wx.EVT_TEXT_ENTER,OnPosRange)
     4099                vecRange.Bind(wx.EVT_KILL_FOCUS,OnPosRange)
     4100                rbsizer2.Add(vecRange,0,wx.ALIGN_CENTER_VERTICAL)
     4101            rbsizer.Add(rbsizer1)   
     4102            rbsizer.Add(rbsizer2)   
     4103            if model['Type'] == 'Residue':
     4104                rbsizer3 = wx.FlexGridSizer(1,8,5,5)
     4105                for it,tor in enumerate(model['Tor'][0]):
     4106                    name = 'Tor('+str(it)+')'
     4107                    torRef = wx.CheckBox(MCSA,-1,label=' %s: '%(name))
     4108                    torRef.SetValue(model['Tor'][1][it])
     4109                    torRef.Bind(wx.EVT_CHECKBOX,OnPosRef)
     4110                    Indx[torRef.GetId()] = [model,'Tor',it]
     4111                    rbsizer3.Add(torRef,0,wx.ALIGN_CENTER_VERTICAL)
     4112                    torVal = wx.TextCtrl(MCSA,-1,'%.4f'%(tor),style=wx.TE_PROCESS_ENTER)
     4113                    torVal.Bind(wx.EVT_TEXT_ENTER,OnPosVal)
     4114                    torVal.Bind(wx.EVT_KILL_FOCUS,OnPosVal)
     4115                    Indx[torVal.GetId()] = [model,'Tor',it]
     4116                    rbsizer3.Add(torVal,0,wx.ALIGN_CENTER_VERTICAL)
     4117                    rbsizer3.Add(wx.StaticText(MCSA,-1,' Range: '),0,wx.ALIGN_CENTER_VERTICAL)
     4118                    rmin,rmax = model['Tor'][2][it]
     4119                    torRange = wx.TextCtrl(MCSA,-1,'%.3f %.3f'%(rmin,rmax),style=wx.TE_PROCESS_ENTER)
     4120                    Indx[torRange.GetId()] = [model,'Tor',it]
     4121                    torRange.Bind(wx.EVT_TEXT_ENTER,OnPosRange)
     4122                    torRange.Bind(wx.EVT_KILL_FOCUS,OnPosRange)
     4123                    rbsizer3.Add(torRange,0,wx.ALIGN_CENTER_VERTICAL)
     4124                rbsizer.Add(rbsizer3)
     4125               
     4126            return rbsizer
     4127           
     4128        def MDSizer(POData):
     4129           
     4130            def OnPORef(event):
     4131                POData['Coef'][1] = poRef.GetValue()
     4132               
     4133            def OnPOVal(event):
     4134                try:
     4135                    mdVal = float(poVal.GetValue())
     4136                    if mdVal > 0:
     4137                        POData['Coef'][0] = mdVal
     4138                except ValueError:
     4139                    pass
     4140                poVal.SetValue("%.3f"%(POData['Coef'][0]))
     4141               
     4142            def OnPORange(event):
     4143                Range = poRange.GetValue().split()
     4144                try:
     4145                    rmin,rmax = [float(Range[i]) for i in range(2)]
     4146                    if rmin >= rmax:
     4147                        raise ValueError
     4148                except (ValueError,IndexError):
     4149                    rmin,rmax = POData['Coef'][2]
     4150                POData['Coef'][2] = [rmin,rmax]
     4151                poRange.SetValue('%.3f %.3f'%(rmin,rmax))                 
     4152               
     4153            def OnPOAxis(event):
     4154                Saxis = poAxis.GetValue().split()
     4155                try:
     4156                    hkl = [int(Saxis[i]) for i in range(3)]
     4157                except (ValueError,IndexError):
     4158                    hkl = POData['axis']
     4159                if not np.any(np.array(hkl)):
     4160                    hkl = POData['axis']
     4161                POData['axis'] = hkl
     4162                h,k,l = hkl
     4163                poAxis.SetValue('%3d %3d %3d'%(h,k,l))                 
     4164               
     4165            poSizer = wx.BoxSizer(wx.HORIZONTAL)
     4166            poRef = wx.CheckBox(MCSA,-1,label=' March-Dollase ratio: ')
     4167            poRef.SetValue(POData['Coef'][1])
     4168            poRef.Bind(wx.EVT_CHECKBOX,OnPORef)
     4169            poSizer.Add(poRef,0,wx.ALIGN_CENTER_VERTICAL)
     4170            poVal = wx.TextCtrl(MCSA,-1,'%.3f'%(POData['Coef'][0]),style=wx.TE_PROCESS_ENTER)
     4171            poVal.Bind(wx.EVT_TEXT_ENTER,OnPOVal)
     4172            poVal.Bind(wx.EVT_KILL_FOCUS,OnPOVal)
     4173            poSizer.Add(poVal,0,wx.ALIGN_CENTER_VERTICAL)
     4174            poSizer.Add(wx.StaticText(MCSA,-1,' Range: '),0,wx.ALIGN_CENTER_VERTICAL)
     4175            rmin,rmax = POData['Coef'][2]
     4176            poRange = wx.TextCtrl(MCSA,-1,'%.3f %.3f'%(rmin,rmax),style=wx.TE_PROCESS_ENTER)
     4177            poRange.Bind(wx.EVT_TEXT_ENTER,OnPORange)
     4178            poRange.Bind(wx.EVT_KILL_FOCUS,OnPORange)
     4179            poSizer.Add(poRange,0,wx.ALIGN_CENTER_VERTICAL)                       
     4180            poSizer.Add(wx.StaticText(MCSA,-1,' Unique axis, H K L: '),0,wx.ALIGN_CENTER_VERTICAL)
     4181            h,k,l = POData['axis']
     4182            poAxis = wx.TextCtrl(MCSA,-1,'%3d %3d %3d'%(h,k,l),style=wx.TE_PROCESS_ENTER)
     4183            poAxis.Bind(wx.EVT_TEXT_ENTER,OnPOAxis)
     4184            poAxis.Bind(wx.EVT_KILL_FOCUS,OnPOAxis)
     4185            poSizer.Add(poAxis,0,wx.ALIGN_CENTER_VERTICAL)
     4186            return poSizer
     4187       
     4188        # UpdateMCSA executable code starts here
     4189        if refresh:
     4190            MCSA.DestroyChildren()
     4191        general = data['General']
     4192        Amat,Bmat = G2lat.cell2AB(general['Cell'][1:7])
     4193        RBData = G2frame.PatternTree.GetItemPyData(   
     4194            G2gd.GetPatternTreeItemId(G2frame,G2frame.root,'Rigid bodies'))
     4195        Indx = {}
     4196        atomStyle = 'balls & sticks'
     4197        if 'macro' in general['Type']:
     4198            atomStyle = 'sticks'
     4199        G2frame.dataFrame.SetStatusText('')
     4200        mainSizer = wx.BoxSizer(wx.VERTICAL)
     4201        if not data['MCSA']['Models']:
     4202            mainSizer.Add((5,5),0)
     4203            mainSizer.Add(wx.StaticText(MCSA,-1,'No MC/SA models:'),0,wx.ALIGN_CENTER_VERTICAL)
     4204            mainSizer.Add((5,5),0)
     4205        else:
     4206            mainSizer.Add((5,5),0)
     4207            mainSizer.Add(wx.StaticText(MCSA,-1,'MC/SA models:'),0,wx.ALIGN_CENTER_VERTICAL)
     4208            mainSizer.Add((5,5),0)
     4209            for model in data['MCSA']['Models']:
     4210                if model['Type'] == 'MD':
     4211                    mainSizer.Add(MDSizer(model))
     4212                elif model['Type'] == 'Atom':
     4213                    mainSizer.Add(atomSizer(model))
     4214                else:
     4215                    mainSizer.Add(rbSizer(model))
     4216                G2gd.HorizontalLine(mainSizer,MCSA)
     4217               
     4218        if not data['MCSA']['Results']:
     4219            mainSizer.Add((5,5),0)
     4220            mainSizer.Add(wx.StaticText(MCSA,-1,'No MC/SA results:'),0,wx.ALIGN_CENTER_VERTICAL)
     4221            mainSizer.Add((5,5),0)
     4222        else:
     4223            for result in data['MCSA']['Results']:
     4224                print result
     4225
     4226        MCSA.SetSizer(mainSizer)
     4227        if G2frame.dataFrame.PhaseUserSize is None:
     4228            mainSizer.FitInside(G2frame.dataFrame)
     4229            Size = mainSizer.GetMinSize()
     4230            Size[0] += 40
     4231            Size[1] = max(Size[1],290) + 35
     4232            MCSA.SetSize(Size)
     4233            MCSA.SetScrollbars(10,10,Size[0]/10-4,Size[1]/10-1)
     4234            Size[1] = min(Size[1],450)
     4235            G2frame.dataFrame.setSizePosLeft(Size)
     4236        else:
     4237            Size = G2frame.dataFrame.PhaseUserSize
     4238            MCSA.SetSize(Size)
     4239            MCSA.SetScrollbars(10,10,Size[0]/10-4,Size[1]/10-1)
     4240            G2frame.dataFrame.Update()
     4241
     4242    def OnRunMCSA(event):
     4243        generalData = data['General']
     4244        mcsaControls = generalData['MCSA controls']
     4245        reflName = mcsaControls['Data source']
     4246        phaseName = generalData['Name']
     4247        MCSAdata = data['MCSA']
     4248        covData = {}
     4249        if 'PWDR' in reflName:
     4250            PatternId = G2gd.GetPatternTreeItemId(G2frame,G2frame.root, reflName)
     4251            reflSets = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId,'Reflection Lists'))
     4252            reflData = reflSets[phaseName]
     4253            reflType = 'PWDR'
     4254        elif 'HKLF' in reflName:
     4255            PatternId = G2gd.GetPatternTreeItemId(G2frame,G2frame.root, reflName)
     4256            reflData = G2frame.PatternTree.GetItemPyData(PatternId)[1]
     4257            reflType = 'HKLF'
     4258        elif reflName == 'Pawley reflections':
     4259            reflData = data['Pawley ref']
     4260            covData = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,G2frame.root, 'Covariance'))
     4261            reflType = 'Pawley'
     4262        else:
     4263            print '**** ERROR - No data defined for MC/SA run'
     4264            return
     4265        MCSAmodels = MCSAdata['Models']
     4266        if not len(MCSAmodels):
     4267            print '**** ERROR - no models defined for MC/SA run****'
     4268            return
     4269        pgbar = wx.ProgressDialog('MC/SA','Residual Rcf =',101.0,
     4270            style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_CAN_ABORT)
     4271        screenSize = wx.ClientDisplayRect()
     4272        Size = pgbar.GetSize()
     4273        Size = (int(Size[0]*1.2),Size[1]) # increase size a bit along x
     4274        pgbar.SetPosition(wx.Point(screenSize[2]-Size[0]-305,screenSize[1]+5))
     4275        pgbar.SetSize(Size)
     4276        try:
     4277            MCSAdata['Results'] = G2mth.mcsaSearch(data,reflType,reflData,covData,pgbar)
     4278        finally:
     4279            pgbar.Destroy()
     4280        if not data['Drawing']:                 #if new drawing - no drawing data!
     4281            SetupDrawingData()
     4282        print ' MC/SA run finished: best model residual: %f.3'%(0.01)
     4283        UpdateMCSA()
     4284
     4285    def OnMCSAaddAtom(event):
     4286        dlg = G2elemGUI.PickElement(G2frame)
     4287        if dlg.ShowModal() == wx.ID_OK:
     4288            El = dlg.Elem.strip()
     4289            Info = G2elem.GetAtomInfo(El)
     4290        dlg.Destroy()
     4291       
     4292        atom = {'Type':'Atom','atType':El,'Pos':[[0.,0.,0.],
     4293            [False,False,False],[[0.,1.],[0.,1.],[0.,1.]]],
     4294            'name':El+'('+str(len(data['MCSA']['Models']))+')'}     
     4295        data['MCSA']['Models'].append(atom)
     4296        data['MCSA']['AtInfo'][El] = [Info['Drad'],Info['Color']]
     4297        G2plt.PlotStructure(G2frame,data)
     4298        UpdateMCSA()
     4299       
     4300    def OnMCSAaddRB(event):
     4301        rbData = G2frame.PatternTree.GetItemPyData(   
     4302            G2gd.GetPatternTreeItemId(G2frame,G2frame.root,'Rigid bodies'))
     4303        rbNames = {}
     4304        for rbVec in rbData['Vector']:
     4305            if rbVec != 'AtInfo':
     4306                rbNames[rbData['Vector'][rbVec]['RBname']] = ['Vector',rbVec]
     4307        for rbRes in rbData['Residue']:
     4308            if rbRes != 'AtInfo':
     4309                rbNames[rbData['Residue'][rbRes]['RBname']] = ['Residue',rbRes]
     4310        if not rbNames:
     4311            print '**** ERROR - no rigid bodies defined ****'
     4312            return
     4313        dlg = wx.SingleChoiceDialog(G2frame.dataFrame,'Select','Rigid body',rbNames.keys())
     4314        if dlg.ShowModal() == wx.ID_OK:
     4315            sel = dlg.GetSelection()
     4316            rbname = rbNames.keys()[sel]
     4317            rbType,rbId = rbNames[rbname]
     4318            RB = rbData[rbType][rbId]
     4319        body = {'name':RB['RBname']+'('+str(len(data['MCSA']['Models']))+')','RBId':rbId,'Type':rbType,
     4320            'Pos':[[0.,0.,0.],[False,False,False],[[0.,1.],[0.,1.],[0.,1.]]],'Ovar':'','MolCent':False,
     4321            'Ori':[[1.,0.,0.,0.],[False,False,False,False],[[-180.,180.],[-1.,1.],[-1.,1.],[-1.,1.]]]}
     4322        if rbType == 'Residue':
     4323            body['Tor'] = [[],[],[]]
     4324            for i,tor in enumerate(RB['rbSeq']):
     4325                body['Tor'][0].append(0.0)
     4326                body['Tor'][1].append(False)
     4327                body['Tor'][2].append([0.,360.])
     4328        data['MCSA']['Models'].append(body)
     4329        data['MCSA']['rbData'] = rbData
     4330        data['MCSA']['AtInfo'].update(rbData[rbType]['AtInfo'])
     4331        G2plt.PlotStructure(G2frame,data)
     4332        UpdateMCSA()
     4333       
     4334    def OnMCSAclear(event):
     4335        data['MCSA'] = {'Models':[{'Type':'MD','Coef':[1.0,False,[0.,3.],],'axis':[0,0,1]}],'Results':[],'AtInfo':{}}
     4336        G2plt.PlotStructure(G2frame,data)
     4337        UpdateMCSA()           
    39194338                       
    3920 
    39214339################################################################################
    39224340##### Pawley routines
     
    43834801        else:
    43844802            print 'Bad charge flip map - no peak search done'
    4385                
     4803                           
    43864804    def OnTextureRefine(event):
    43874805        print 'refine texture?'
     
    44924910            #G2plt.PlotStructure(G2frame,data)
    44934911            wx.CallAfter(G2plt.PlotStructure,G2frame,data)
     4912        elif text == 'MC/SA':
     4913            G2gd.SetDataMenuBar(G2frame,G2frame.dataFrame.MCSAMenu)
     4914            G2frame.dataFrame.Bind(wx.EVT_MENU, OnMCSAaddAtom, id=G2gd.wxID_ADDMCSAATOM)
     4915            G2frame.dataFrame.Bind(wx.EVT_MENU, OnMCSAaddRB, id=G2gd.wxID_ADDMCSARB)
     4916            G2frame.dataFrame.Bind(wx.EVT_MENU, OnMCSAclear, id=G2gd.wxID_CLEARMCSARB)
     4917            UpdateMCSA()                       
     4918            wx.CallAfter(G2plt.PlotStructure,G2frame,data)
    44944919        elif text == 'Texture':
    44954920            G2gd.SetDataMenuBar(G2frame,G2frame.dataFrame.TextureMenu)
     
    44974922            G2frame.dataFrame.Bind(wx.EVT_MENU, OnTextureClear, id=G2gd.wxID_CLEARTEXTURE)
    44984923            UpdateTexture()                       
    4499             G2plt.PlotTexture(G2frame,data,Start=True)
    4500            
     4924            G2plt.PlotTexture(G2frame,data,Start=True)           
    45014925        else:
    45024926            G2gd.SetDataMenuBar(G2frame)
     
    45194943    MapPeaks = G2gd.GSGrid(G2frame.dataDisplay)
    45204944    G2frame.dataDisplay.AddPage(MapPeaks,'Map peaks')
     4945    MCSA = wx.ScrolledWindow(G2frame.dataDisplay)
     4946    G2frame.dataDisplay.AddPage(MCSA,'MC/SA')
    45214947    Texture = wx.ScrolledWindow(G2frame.dataDisplay)
    45224948    G2frame.dataDisplay.AddPage(Texture,'Texture')
     
    45284954    G2frame.dataFrame.Bind(wx.EVT_MENU, OnChargeFlip, id=G2gd.wxID_CHARGEFLIP)
    45294955    G2frame.dataFrame.Bind(wx.EVT_MENU, OnFourClear, id=G2gd.wxID_FOURCLEAR)
     4956    G2frame.dataFrame.Bind(wx.EVT_MENU, OnRunMCSA, id=G2gd.wxID_RUNMCSA)   
    45304957    G2frame.dataDisplay.Bind(wx.aui.EVT_AUINOTEBOOK_PAGE_CHANGED, OnPageChanged)
    45314958
    45324959    G2gd.SetDataMenuBar(G2frame,G2frame.dataFrame.DataGeneral)
    4533     SetupGeneral() # this is done all over the place (including in
    4534     # UpdateGeneral). Why here too?
     4960    SetupGeneral()
    45354961   
    45364962    GeneralData = data['General']
  • trunk/GSASIIplot.py

    r924 r934  
    26182618    testRBObj = data.get('testRBObj',{})
    26192619    rbObj = testRBObj.get('rbObj',{})
     2620    MCSA = data.get('MCSA',{})
     2621    mcsaModels = []
     2622    if len(MCSA):
     2623        mcsaModels = MCSA['Models']
    26202624    drawAtoms = drawingData.get('Atoms',[])
    26212625    mapData = {}
     
    34643468                RenderBonds(x,y,z,rbBonds[ind],0.03,Gr)
    34653469                RenderLabel(x,y,z,name,0.2,Or)
     3470        if len(mcsaModels) > 1 and pageName == 'MC/SA':             #skip the default MD entry
     3471            XYZ,atTypes = G2mth.UpdateMCSAxyz(Bmat,mcsaModels,MCSA.get('rbData',{}))
     3472            rbBonds = FindPeaksBonds(XYZ)
     3473            for ind,[x,y,z] in enumerate(XYZ):
     3474                aType = atTypes[ind]
     3475                name = aType+str(ind)
     3476                color = np.array(MCSA['AtInfo'][aType][1])
     3477                RenderSphere(x,y,z,0.2,color/255.)
     3478                RenderBonds(x,y,z,rbBonds[ind],0.03,Gr)
     3479                RenderLabel(x,y,z,name,0.2,Or)
    34663480        if Backbones:
    34673481            for chain in Backbones:
     
    35543568        Bonds = FindBonds(XYZ)
    35553569    elif rbType == 'Residue':
    3556         atNames = [str(i)+':'+Ty for i,Ty in enumerate(rbData['atNames'])]
     3570#        atNames = [str(i)+':'+Ty for i,Ty in enumerate(rbData['atNames'])]
     3571        atNames = rbData['atNames']
    35573572        XYZ = np.copy(rbData['rbXYZ'])      #don't mess with original!
    35583573        Seq = rbData['rbSeq']
     
    36103625    def SetLights():
    36113626        glEnable(GL_DEPTH_TEST)
    3612         glShadeModel(GL_SMOOTH)
     3627        glShadeModel(GL_FLAT)
    36133628        glEnable(GL_LIGHTING)
    36143629        glEnable(GL_LIGHT0)
     
    37523767        glMultMatrixf(matRot.T)
    37533768        RenderUnitVectors(0.,0.,0.)
    3754         radius = 0.4
     3769        radius = 0.2
    37553770        for iat,atom in enumerate(XYZ):
    37563771            x,y,z = atom
     
    37583773            color = np.array(CL)/255.
    37593774            RenderSphere(x,y,z,radius,color)
    3760             RenderBonds(x,y,z,Bonds[iat],0.1,color)
     3775            RenderBonds(x,y,z,Bonds[iat],0.05,color)
    37613776            RenderLabel(x,y,z,atNames[iat],radius)
    37623777        Page.canvas.SwapBuffers()
  • trunk/GSASIIstrIO.py

    r926 r934  
    550550       
    551551################################################################################
    552 ##### Rigid Body Models
     552##### Rigid Body Models  and not General.get('doPawley')
    553553################################################################################
    554554       
  • trunk/GSASIIstrMath.py

    r927 r934  
    15421542                print 'TOF Undefined at present'
    15431543                raise Exception    #no TOF yet
    1544             #do atom derivatives -  for RB,F,X & U so far             
    1545             corr = dervDict['int']/refl[9]
    1546             if Ka2:
    1547                 corr2 = dervDict2['int']/refl[9]
    1548             for name in varylist+dependentVars:
    1549                 if '::RBV;' in name:
    1550                     pass
    1551                 else:
    1552                     try:
    1553                         aname = name.split(pfx)[1][:2]
    1554                         if aname not in ['Af','dA','AU','RB']: continue # skip anything not an atom or rigid body param
    1555                     except IndexError:
    1556                         continue
    1557                 if name in varylist:
    1558                     dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
    1559                     if Ka2:
    1560                         dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
    1561                 elif name in dependentVars:
    1562                     depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
    1563                     if Ka2:
    1564                         depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
     1544            if not Phase['General'].get('doPawley'):
     1545                #do atom derivatives -  for RB,F,X & U so far             
     1546                corr = dervDict['int']/refl[9]
     1547                if Ka2:
     1548                    corr2 = dervDict2['int']/refl[9]
     1549                for name in varylist+dependentVars:
     1550                    if '::RBV;' in name:
     1551                        pass
     1552                    else:
     1553                        try:
     1554                            aname = name.split(pfx)[1][:2]
     1555                            if aname not in ['Af','dA','AU','RB']: continue # skip anything not an atom or rigid body param
     1556                        except IndexError:
     1557                            continue
     1558                    if name in varylist:
     1559                        dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
     1560                        if Ka2:
     1561                            dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
     1562                    elif name in dependentVars:
     1563                        depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
     1564                        if Ka2:
     1565                            depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
    15651566#        print 'profile derv time: %.3fs'%(time.time()-time0)
    15661567    # now process derivatives in constraints
Note: See TracChangeset for help on using the changeset viewer.