Implement penalty fxn for negative SpinRB distributions & fix issues found in penaltyDeriv.

Author: vondreele

GSASII/GSASIIobj.py

@@ -33,7 +33,7 @@
 restraintNames = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
                   ['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
                   ['ChemComp','Sites'],['Texture','HKLs'],['Moments','Moments'],
-                  ['General','General']]
+                  ['General','General'],['SpinRB','SpinRBs']]
 '''Names of restraint keys for the restraint dict and the location of
 the restraints in each dict
 '''

GSASII/GSASIIphsGUI.py

@@ -10138,9 +10138,6 @@ def OnSliceShell(event):
                 RBObj['sliceSh'] = not RBObj['sliceSh']
                 G2plt.PlotStructure(G2frame,data)
 
-            def OnNegDen(event):
-                RBObj['noNeg'][0] = not RBObj['noNeg'][0]
-
             symchoice = ['53m','m3m','-43m','6/mmm','-6m2','-3m','3m','32','-3','3','4/mmm','-42m',
                 'mmm','2/m','2','m','-1','1']
             RBObj['hide'] = RBObj.get('hide',[False for i in range(len(RBObj['atType']))])
@@ -10197,18 +10194,11 @@ def OnNegDen(event):
             plotLine.Add(fadeShell,0,WACV)
             plotLine.Add(sliceShell,0,WACV)
             plotLine.Add(addShell,0,WACV)
-            negDen = wx.CheckBox(RigidBodies,label='Use negative density penalty for all shells? Penalty: ')
-            negDen.SetValue(RBObj['noNeg'][0])
-            negDen.Bind(wx.EVT_CHECKBOX,OnNegDen)            
-            secLine.Add(negDen,0,WACV)
-            secLine.Add(G2G.ValidatedTxtCtrl(RigidBodies,RBObj['noNeg'],1,nDig=(10,4),
-                xmin=0.,xmax=1.,typeHint=float),0,WACV)
             sprbSizer.Add(plotLine)
             G2G.HorizontalLine(sprbSizer,RigidBodies)
             sprbSizer.Add(topLine)
             sprbSizer.Add(secLine)
             sprbSizer.Add(SHsizer())
-            
             return sprbSizer
 
         def ResrbSizer(RBObj,resIndx):

GSASII/GSASIIrestrGUI.py

@@ -206,6 +206,8 @@ def OnAddRestraint(event):
             AddChemCompRestraint(restrData['ChemComp'])
         elif 'Texture' in G2frame.restrBook.GetPageText(page):
             AddTextureRestraint(restrData['Texture'])
+        elif 'Spin RB' in G2frame.restrBook.GetPageText(page):
+            AddSpinRBRestraint(restrData['SpinRB'])
         elif 'Moments' in G2frame.restrBook.GetPageText(page):
             AddMomentRestraint(restrData['Moments'])
         elif 'General' in G2frame.restrBook.GetPageText(page):
@@ -890,6 +892,19 @@ def AddMomentRestraint(momentRestData):
                 print ('**** ERROR - not enough atoms for a average moment restraint - try again ****')
         UpdateMomentRestr(momentRestData)                
 
+    def AddSpinRBRestraint(spinRestData):
+        SpinRB = RBModels.get('Spin',[])
+        indx = [item['Ids'][0] for item in SpinRB]
+        Snames = G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,ct-1)                
+        dlg = G2G.G2MultiChoiceDialog(G2frame,'Select atoms for negative spin RB restraint in '+General['Name'],
+                'Select atoms',Snames)
+        if dlg.ShowModal() == wx.ID_OK:
+            sels = dlg.GetSelections()
+            for sel in sels:
+                minspin = [indx[sel],0.01,False,1.0]    #spin atom id, esd, unit flag, unit esd
+                if minspin not in spinRestData['SpinRBs']:
+                    spinRestData['SpinRBs'].append(minspin)            
+        UpdateSpinRBRestr(spinRestData)                
 
     def AddTextureRestraint(textureRestData):
         dlg = wx.TextEntryDialog(G2frame,'Enter h k l for pole figure restraint','Enter HKL','')
@@ -2121,6 +2136,76 @@ def OnCellChange(event):
             mainSizer.Add(wx.StaticText(TextureRestr,-1,'No texture restraints for this phase'),0,)
         G2phsGUI.SetPhaseWindow(TextureRestr,mainSizer,Scroll=0)
 
+    def UpdateSpinRBRestr(spinRestData):
+            
+        def OnDeleteRestraint(event):
+            rows = GetSelectedRows(SpinRBs,'delete',G2frame)
+            G2frame.GetStatusBar().SetStatusText('',1)
+            if not rows:
+                G2frame.GetStatusBar().SetStatusText('First select restraints to be deleted',1)
+                return
+            rows.sort()
+            rows.reverse()
+            for row in rows:
+                spinRestData['SpinRBs'].remove(spinRestData['SpinRBs'][row])
+            wx.CallAfter(UpdateSpinRBRestr,spinRestData)                
+            
+        def OnCellChange(event):
+            r,c = event.GetRow(),event.GetCol()
+            try:
+                if c in [1,3]:   #esds
+                    new = float(spinRBTable.GetValue(r,c))
+                    if new < -1. or new > 2.:
+                        raise ValueError
+                else:
+                    new = spinRBTable.GetValue(r,c)
+                spinRestData['SpinRBs'][r][c] = new
+            except ValueError:
+                pass            
+            wx.CallAfter(UpdateSpinRBRestr,spinRestData)                
+
+        try:
+            if SpinRBRestr.GetSizer(): SpinRBRestr.GetSizer().Clear(True)
+        except:
+            pass
+        mainSizer = wx.BoxSizer(wx.VERTICAL)
+        mainSizer.Add((5,5),0)
+        mainSizer.Add(WtBox(SpinRBRestr,spinRestData),0)
+        mainSizer.Add(wx.StaticText(SpinRBRestr,-1, 
+            'NB: The spin RB restraints suppress negative spin RB densities for the selected atoms\n'
+            '    "unit esd" gives a bias toward a flatter spin RB density'),0)
+        mainSizer.Add((5,5),0)
+
+        for i in (G2G.wxID_RESTDELETE,G2G.wxID_RESRCHANGEVAL,G2G.wxID_RESTCHANGEESD):
+            G2frame.dataWindow.RestraintEdit.Enable(id=i,enable=False)
+        if len(spinRestData['SpinRBs']):
+            table = []
+            rowLabels = []
+            Types = [wg.GRID_VALUE_STRING,wg.GRID_VALUE_FLOAT+':10,2',
+                wg.GRID_VALUE_BOOL,wg.GRID_VALUE_FLOAT+':10,2']
+            colLabels = ['atom','neg esd','use unit?','unit esd']
+            for i,[atomid,esd1,ifesd2,esd2] in enumerate(spinRestData['SpinRBs']):
+                table.append(['%s'%G2mth.GetAtomItemsById(Atoms,AtLookUp,atomid,ct-1)[0],esd1,ifesd2,esd2])
+                rowLabels.append(str(i))
+            spinRBTable = G2G.Table(table,rowLabels=rowLabels,colLabels=colLabels,types=Types)
+            SpinRBs = G2G.GSGrid(SpinRBRestr)
+            SpinRBs.SetTable(spinRBTable, True)
+            SpinRBs.AutoSizeColumns(False)
+            for r in range(len(spinRestData['SpinRBs'])):
+                SpinRBs.SetReadOnly(r,0,True)
+                SpinRBs.SetCellStyle(r,0,VERY_LIGHT_GREY,True)
+                if not spinRBTable.GetValue(r,2):
+                    SpinRBs.SetReadOnly(r,3,True)
+                    SpinRBs.SetCellStyle(r,3,VERY_LIGHT_GREY,True)
+                    SpinRBs.SetCellTextColour(r,3,VERY_LIGHT_GREY)
+            SpinRBs.Bind(wg.EVT_GRID_LABEL_LEFT_CLICK,OnRowSelect)
+            SpinRBs.Bind(wg.EVT_GRID_CELL_CHANGED, OnCellChange)
+            G2frame.dataWindow.RestraintEdit.Enable(id=G2G.wxID_RESTDELETE,enable=True)
+            G2frame.Bind(wx.EVT_MENU, OnDeleteRestraint, id=G2G.wxID_RESTDELETE)
+            mainSizer.Add(SpinRBs)
+        else:
+            mainSizer.Add(wx.StaticText(SpinRBRestr,-1,'No spin RB restraints for this phase'),0,)
+        G2phsGUI.SetPhaseWindow(SpinRBRestr,mainSizer,Scroll=0)        
 
     def UpdateGeneralRestr(generalRestData):
         '''Display any generalized restraint expressions'''
@@ -2224,7 +2309,6 @@ def OnDelGenRestraint(event):              #does this work??
                 GridSiz.Add(wx.StaticText(GeneralRestr,wx.ID_ANY,txt))
             mainSizer.Add(GridSiz)
         G2phsGUI.SetPhaseWindow(GeneralRestr,mainSizer,Scroll=0)
-
 
     def OnPageChanged(event):
         page = event.GetSelection()
@@ -2287,6 +2371,12 @@ def OnPageChanged(event):
             G2frame.dataWindow.RestraintEdit.Enable(G2G.wxID_RESRCHANGEVAL,True)
             textureRestData = restrData['Texture']
             UpdateTextureRestr(textureRestData)
+        elif text == 'Spin RB':
+            G2gd.SetDataMenuBar(G2frame,G2frame.dataWindow.RestraintMenu)
+            G2frame.dataWindow.RestraintEdit.Enable(G2G.wxID_RESTRAINTADD,True)
+            G2frame.dataWindow.RestraintEdit.Enable(G2G.wxID_RESRCHANGEVAL,True)
+            spinRBRestData = restrData['SpinRB']
+            UpdateSpinRBRestr(spinRBRestData)
         elif text == 'Moments':
             G2gd.SetDataMenuBar(G2frame,G2frame.dataWindow.RestraintMenu)
             G2frame.dataWindow.RestraintEdit.Enable(G2G.wxID_RESTRAINTADD,True)
@@ -2374,11 +2464,14 @@ def OnSelectPage(event):
         restrData['Rama'] = {'wtFactor':1.0,'Coeff':{},'Ramas':[],'Use':True}
     if 'Texture' not in restrData:
         restrData['Texture'] = {'wtFactor':1.0,'HKLs':[],'Use':True}
+    if 'SpinRB' not in restrData:
+        restrData['SpinRB'] = {'wtFactor':1.0,'SpinRBs':[],'Use':True}
     if 'ChemComp' not in restrData:
         restrData['ChemComp'] = {'wtFactor':1.0,'Sites':[],'Use':True}
     if 'General' not in restrData:
         restrData['General'] = {'wtFactor':1.0,'General':[], 'Use':True}
     General = phasedata['General']
+    RBModels = phasedata['RBModels']
     if General['Type'] == 'magnetic' and 'Moments' not in restrData:
         restrData['Moments'] = {'wtFactor':1.0,'Moments':[],'Use':True}
     Cell = General['Cell'][1:7]          #skip flag & volume    
@@ -2476,6 +2569,12 @@ def OnSelectPage(event):
         Pages.append(txt)
 
     UpdateBondRestr(restrData['Bond'])
+    
+    if 'Spin' in phasedata['RBModels'] and len(phasedata['RBModels']['Spin']):
+        txt = 'Spin RB'
+        SpinRBRestr = wx.ScrolledWindow(G2frame.restrBook)
+        G2frame.restrBook.AddPage(SpinRBRestr,txt)
+        Pages.append(txt)    
 
     G2frame.restrBook.Bind(wx.aui.EVT_AUINOTEBOOK_PAGE_CHANGED, OnPageChanged)
 

GSASII/GSASIIstrIO.py

@@ -2147,11 +2147,13 @@ def PrintRestraints(cell,SGData,AtPtrs,Atoms,AtLookup,textureData,phaseRest,pFil
                         if num:
                             sum = np.sum(Z)
                         pFile.write ('   %d %d %d  %d %8.3f %8.3f %8d   %s    %8.3f\n'%(hkl[0],hkl[1],hkl[2],grid,esd1,sum,num,str(ifesd2),esd2))
+                elif name =='SpinRB':
+                    continue
                 elif name == 'General':
                     pFile.write('  target   sig     obs  expression        variables \n')
                     for expObj,target,esd in itemRest[rest]:
                         val = '?'
-                        calcobj = G2obj.ExpressionCalcObj(expObj)
+                        #calcobj = G2obj.ExpressionCalcObj(expObj)
                         # need to get parmDict to evaluate the value
                         #calcobj.SetupCalc(parmDict)
                         #val = ' {:8.3g} '.format(calcobj.EvalExpression())
@@ -4330,6 +4332,7 @@ def WriteVecRBModel(RBModel,sigDict={},irb=None):
 
 atmPattrn = re.compile("::A[xyz]:")
 fmtSplit =  re.compile('%([0-9]+)\\.([0-9]+)(.*)')
+
 def fmtESD(varname,SigDict,fmtcode,ndig=None,ndec=None):
     '''Format an uncertainty value as requested, but surround the
     number by () if the parameter is set by an equivalence