Merge branch 'main' into kvec_isodistort_zyp

Author: Kvieta1990

GSASII/GSASIIlattice.py

@@ -225,8 +225,8 @@ def cellDijFill(pfx,phfx,SGData,parmDict):
     elif SGData['SGLaue'] in ['3R', '3mR']:
         A = [parmDict[pfx+'A0']+parmDict[phfx+'D11'],parmDict[pfx+'A0']+parmDict[phfx+'D11'],
             parmDict[pfx+'A0']+parmDict[phfx+'D11'],
-            parmDict[pfx+'A3']+parmDict[phfx+'D23'],parmDict[pfx+'A3']+parmDict[phfx+'D23'],
-            parmDict[pfx+'A3']+parmDict[phfx+'D23']]
+            parmDict[pfx+'A3']+parmDict[phfx+'D12'],parmDict[pfx+'A3']+parmDict[phfx+'D12'],
+            parmDict[pfx+'A3']+parmDict[phfx+'D12']]
     elif SGData['SGLaue'] in ['m3m','m3']:
         A = [parmDict[pfx+'A0']+parmDict[phfx+'D11'],parmDict[pfx+'A0']+parmDict[phfx+'D11'],
              parmDict[pfx+'A0']+parmDict[phfx+'D11'],0,0,0]
@@ -2460,12 +2460,6 @@ def RBsymChk(RBsym,cubic,coefNames,L=18):
     :param list coefNames: sp. harm coefficient names to be checked/converted
     :param int L:  maximum spherical harmonic order no. for cubic generation if needed
     '''
-    # cubicsigns = {'C(3,1)c':[-1,],'C(4,1)c':[1,1,],'C(6,1)c':[1,-1],'C(6,2)c':[1,-1],'C(7,1)c':[1,-1],'C(8,1)c':[1,1,1],
-    #     'C(9,1)c':[1,-1],'C(9,2)c':[1,-1],'C(10,1)c':[1,-1,-1],'C(10,2)c':[1,1,-1]}
-    # cubicnames = {'C(3,1)c':['C(3,2)',],'C(4,1)c':['C(4,0)','C(4,4)'],'C(6,1)c':['C(6,0)','C(6,4)'],
-    #     'C(6,2)c':['C(6,2)','C(6,6)'],'C(7,1)c':['C(7,2)','C(7,6)'],'C(8,1)c':['C(8,0)','C(8,4)','C(8,8)'],
-    #     'C(9,1)c':['C(9,2)','C((9,6)'],'C(9,2)c':['C(9,4)','C(9,8)'],
-    #     'C(10,1)c':['C(10,0)','C(10,4)','C(10,8)'],'C(10,2)c':['C(10,2)','C(10,6)','C(10,10)']}
     newNames = []
     newSgns = []
     if cubic:       #sytsym is a cubic site
@@ -2760,38 +2754,53 @@ def SamAng(Tth,Gangls,Sangl,IFCoup):
         (BA*dBBdF-BB*dBAdF)/BD],[np.zeros_like(BD),np.zeros_like(BD),np.zeros_like(BD)])
 
     return psi,gam,dPSdA,dGMdA
-
+'''BOH is even order terms from Table 15.2.1 of H.-J. Bunge, "Texture Analysis in Materials Science",
+Cuvillier Verlag, Gottingen, 1993, p 501. BOH = Blmu for cubic [001] parallel to Z. Used in GetKcl & GetKclKsl
+Dummy odd order terms not needed for texture analysis
+'''
 BOH = {
 'L=2':[[],[],[]],
 'L=4':[[0.30469720,0.36418281],[],[]],
 'L=6':[[-0.14104740,0.52775103],[],[]],
 'L=8':[[0.28646862,0.21545346,0.32826995],[],[]],
+'L=9':[[],[],[]],
 'L=10':[[-0.16413497,0.33078546,0.39371345],[],[]],
 'L=12':[[0.26141975,0.27266871,0.03277460,0.32589402],
     [0.09298802,-0.23773812,0.49446631,0.0],[]],
+'L=13':[[],[],[]],
 'L=14':[[-0.17557309,0.25821932,0.27709173,0.33645360],[],[]],
+'L=15':[[],[],[]],
 'L=16':[[0.24370673,0.29873515,0.06447688,0.00377,0.32574495],
     [0.12039646,-0.25330128,0.23950998,0.40962508,0.0],[]],
+'L=17':[[],[],[]],
 'L=18':[[-0.16914245,0.17017340,0.34598142,0.07433932,0.32696037],
     [-0.06901768,0.16006562,-0.24743528,0.47110273,0.0],[]],
+'L=19':[[],[],[]],
 'L=20':[[0.23067026,0.31151832,0.09287682,0.01089683,0.00037564,0.32573563],
     [0.13615420,-0.25048007,0.12882081,0.28642879,0.34620433,0.0],[]],
+'L=21':[[],[],[]],
 'L=22':[[-0.16109560,0.10244188,0.36285175,0.13377513,0.01314399,0.32585583],
     [-0.09620055,0.20244115,-0.22389483,0.17928946,0.42017231,0.0],[]],
+'L=23':[[],[],[]],
 'L=24':[[0.22050742,0.31770654,0.11661736,0.02049853,0.00150861,0.00003426,0.32573505],
     [0.13651722,-0.21386648,0.00522051,0.33939435,0.10837396,0.32914497,0.0],
     [0.05378596,-0.11945819,0.16272298,-0.26449730,0.44923956,0.0,0.0]],
+'L=25':[[],[],[]],
 'L=26':[[-0.15435003,0.05261630,0.35524646,0.18578869,0.03259103,0.00186197,0.32574594],
     [-0.11306511,0.22072681,-0.18706142,0.05439948,0.28122966,0.35634355,0.0],[]],
+'L=27':[[],[],[]],
 'L=28':[[0.21225019,0.32031716,0.13604702,0.03132468,0.00362703,0.00018294,0.00000294,0.32573501],
     [0.13219496,-0.17206256,-0.08742608,0.32671661,0.17973107,0.02567515,0.32619598,0.0],
     [0.07989184,-0.16735346,0.18839770,-0.20705337,0.12926808,0.42715602,0.0,0.0]],
+'L=29':[[],[],[]],
 'L=30':[[-0.14878368,0.01524973,0.33628434,0.22632587,0.05790047,0.00609812,0.00022898,0.32573594],
     [-0.11721726,0.20915005,-0.11723436,-0.07815329,0.31318947,0.13655742,0.33241385,0.0],
     [-0.04297703,0.09317876,-0.11831248,0.17355132,-0.28164031,0.42719361,0.0,0.0]],
+'L=31':[[],[],[]],
 'L=32':[[0.20533892,0.32087437,0.15187897,0.04249238,0.00670516,0.00054977,0.00002018,0.00000024,0.32573501],
     [0.12775091,-0.13523423,-0.14935701,0.28227378,0.23670434,0.05661270,0.00469819,0.32578978,0.0],
     [0.09703829,-0.19373733,0.18610682,-0.14407046,0.00220535,0.26897090,0.36633402,0.0,0.0]],
+'L=33':[[],[],[]],
 'L=34':[[-0.14409234,-0.01343681,0.31248977,0.25557722,0.08571889,0.01351208,0.00095792,0.00002550,0.32573508],
     [-0.11527834,0.18472133,-0.04403280,-0.16908618,0.27227021,0.21086614,0.04041752,0.32688152,0.0],
     [-0.06773139,0.14120811,-0.15835721,0.18357456,-0.19364673,0.08377174,0.43116318,0.0,0.0]]
@@ -3006,7 +3015,8 @@ def SphHarmAng(L,M,P,Th,Ph):
 
     :returns ylmp value/array: as reals
     '''
-    ylmp = spsp.sph_harm(M,L,rpd*Th,rpd*Ph)   #wants radians; order then degree
+    ylmp = spsp.sph_harm(M,L,rpd*Th,rpd*Ph)   #wants radians; order then degree; not normalized
+    #### TODO: this will be deprecated in future scipy; new one sph_harm_y in scipy 1.15.1
 
     if M > 0:
         return (-1)**M*P*np.real(ylmp)*SQ2
@@ -3029,7 +3039,7 @@ def CubicSHarm(L,M,Th,Ph):
 
     '''
     if L == 0:
-        return SphHarmAng(L,M,1,Th,Ph)
+        return SphHarmAng(0,0,1,Th,Ph)
     elif L == 3:
         return SphHarmAng(3,2,-1,Th,Ph)
     elif L == 4:

GSASII/GSASIIphsGUI.py

@@ -2919,6 +2919,16 @@ def OnTransform(event):
         Vvec = np.zeros(3)
         ifMag = False
         BNSlatt = ''
+        ifRB = False
+        RBModels = data['RBModels']
+        for item in RBModels:
+            if RBModels[item]:
+                ifRB = True
+        if ifRB:
+            msg = 'Rigid Bodies present'
+            text = 'Rigid bodies must be deleted before transform can be done'
+            wx.MessageBox(text,caption=msg,style=wx.ICON_EXCLAMATION)
+            return
         while True:
             dlg = TransformDialog(G2frame,data,Trans,Uvec,Vvec,ifMag,BNSlatt)
             try:
@@ -4490,7 +4500,8 @@ def chkUij(Uij,CSI): #needs to do something!!!
                 if 'Atoms' in data['Drawing'] and replot:
                     ci = colLabels.index('I/A')
                     G2mth.DrawAtomsReplaceByID(data,ci+8,atomData[r],ID)
-                    G2plt.PlotStructure(G2frame,data)
+                    G2frame.GetStatusBar().SetStatusText('Structure changed: Do "Edit Atoms/Update draw atoms" to refresh structure drawing',1)
+                    # G2plt.PlotStructure(G2frame,data)
                 if SGData['SpGrp'] != 'P 1':    #no need to update P 1 structures!
                     wx.CallAfter(Paint)
 
@@ -4515,7 +4526,8 @@ def AtomTypeSelect(event):
                 ID = atomData[r][ci+8]
                 if 'Atoms' in data['Drawing']:
                     G2mth.DrawAtomsReplaceByID(data,ci+8,atomData[r],ID)
-                    G2plt.PlotStructure(G2frame,data)
+                    G2frame.GetStatusBar().SetStatusText('Structure changed: Do "Edit Atoms/Update draw atoms" to refresh structure drawing',1)
+                    # G2plt.PlotStructure(G2frame,data)
                 SetupGeneral()
             else:
                 event.Skip()
@@ -4565,7 +4577,8 @@ def RowSelect(event):
                     G2frame.GetStatusBar().SetStatusText('Use right mouse click to brng up Atom editing options',1)
                     Atoms.ClearSelection()
                     Atoms.SelectRow(r,True)
-            G2plt.PlotStructure(G2frame,data)
+            G2frame.GetStatusBar().SetStatusText('Structure changed: Do "Edit Atoms/Update draw atoms" to refresh structure drawing',1)
+            # G2plt.PlotStructure(G2frame,data)
 
         def ChangeSelection(event):
             r,c =  event.GetRow(),event.GetCol()
@@ -4663,8 +4676,8 @@ def Paint(Scroll=0):
                                 Atoms.SetCellStyle(row,ci,WHITE,False)
                                 saveCSI = CSI[0][i]
                             Atoms.SetCellTextColour(row,ci,BLACK)
-                if 'F' in rbExcl:
-                    Atoms.SetCellStyle(row,colF,VERY_LIGHT_GREY,True)
+                # if 'F' in rbExcl:
+                #     Atoms.SetCellStyle(row,colF,VERY_LIGHT_GREY,True)
                 if 'X' in rbExcl:
                     for c in range(0,colX+3):
                         if c != colR:
@@ -5191,7 +5204,7 @@ def AtomModify(event):
             if dlg.ShowModal() == wx.ID_OK:
                 parm = dlg.GetValue()
                 for r in indx:
-                    if not Atoms.IsReadOnly(r,0):   #not if in RB!
+#                    if not Atoms.IsReadOnly(r,0):   #not if in RB!
                         atomData[r][cid] = parm
                 SetupGeneral()
                 FillAtomsGrid(Atoms)
@@ -5562,13 +5575,15 @@ def OnFracSplit(event):
         cx,ct,cs,ci = G2mth.getAtomPtrs(data)
         indx = getAtomSelections(Atoms,ct-1)
         if indx:
+            B2 = {}
             atomData = data['Atoms']
             PE = G2elemGUI.PickElement(G2frame,oneOnly=True)
             if PE.ShowModal() == wx.ID_OK:
                 Atype2 = PE.Elem.strip()
                 AtomInfo2 = G2elem.GetAtomInfo(Atype2)
                 FF2 = G2elem.GetFFtable([Atype2])[Atype2]
-                B2 = AtomInfo2['Isotopes']['Nat. Abund.']['SL'][0]
+                for iso in AtomInfo2['Isotopes']:
+                    B2[iso] = AtomInfo2['Isotopes'][iso]['SL'][0]
             PE.Destroy()
             SQ = 0.0
             PE2 = G2G.SingleFloatDialog(G2frame,'form factor','Enter sinth/lam for atom frac split',SQ)
@@ -5578,24 +5593,36 @@ def OnFracSplit(event):
             ff2 = G2elem.ScatFac(FF2,SQ**2)
         print(' X-ray site fractions for sin(th)/lam = %.3f'%SQ)
         for ind in indx:
-                Aname1 = atomData[ind][ct-1]
-                Atype1 = G2elem.FixValence(atomData[ind][ct])
-                Afrac = atomData[ind][cx+3]
-                Amult = float(atomData[ind][cx+5])
-                AtomInfo1 = G2elem.GetAtomInfo(Atype1)
-                if Atype1 == Atype2:
-                    print('ERROR - 2nd atom type must be different from selected atom')
-                    continue
-                FF1 = G2elem.GetFFtable([Atype1])[Atype1]
-                ff1 = G2elem.ScatFac(FF1,SQ**2)
-                B1 = AtomInfo1['Isotopes']['Nat. Abund.']['SL'][0]
+            Aname1 = atomData[ind][ct-1]
+            Atype1 = G2elem.FixValence(atomData[ind][ct])
+            if Atype1 == 'Q':
+                PE = G2elemGUI.PickElement(G2frame,oneOnly=False)
+                if PE.ShowModal() == wx.ID_OK:
+                    Atype1 = PE.Elem.strip()
+                PE.Destroy()
+            Afrac = atomData[ind][cx+3]
+            Amult = float(atomData[ind][cx+5])
+            AtomInfo1 = G2elem.GetAtomInfo(Atype1)
+            if Atype1 == Atype2:
+                print('ERROR - 2nd atom type must be different from selected atom')
+                continue
+            FF1 = G2elem.GetFFtable([Atype1])[Atype1]
+            ff1 = G2elem.ScatFac(FF1,SQ**2)
+            if ff1 != ff2:
                 frac1 = (ff1*Afrac-ff2)/(ff1-ff2)
-                bfrac1 = (B1*Afrac-B2)/(B1-B2)
                 print(' For %s: X-ray based site fractions %s = %.3f, %.3f/cell; %s = %.3f, %.3f/cell'     \
-                      %(Aname1,Atype1,frac1,frac1*Amult,Atype2,(1.-frac1),(1.-frac1)*Amult))
-                print('        neutron based site fractions %s = %.3f, %.3f/cell; %s = %.3f, %.3f/cell\n'      \
-                      %(Atype1,bfrac1,bfrac1*Amult,Atype2,(1.-bfrac1),(1.-bfrac1)*Amult))
-
+                  %(Aname1,Atype1,frac1,frac1*Amult,Atype2,(1.-frac1),(1.-frac1)*Amult))
+            
+            B1 = AtomInfo1['Isotopes']['Nat. Abund.']['SL'][0]
+            for iso in B2:
+                if B1 != B2[iso]:
+                    bfrac1 = (B1*Afrac-B2[iso])/(B1-B2[iso])
+                    atype2 = Atype2
+                    if 'nat' not in iso:
+                        atype2 += '_%s'%iso
+                    print('        neutron based site fractions %s = %.3f, %.3f/cell; %s = %.3f, %.3f/cell'      \
+                          %(Atype1,bfrac1,bfrac1*Amult,atype2,(1.-bfrac1),(1.-bfrac1)*Amult))
+    
     def OnValidProtein(event):
 
         def pickHandler(resName):
@@ -13867,14 +13894,16 @@ def RepaintRBInfo(rbType,rbIndx,Scroll=0):
                 return
             Indx.clear()
             rbObj = data['RBModels'][rbType][rbIndx]
-            data['Drawing']['viewPoint'][0] = data['Atoms'][AtLookUp[rbObj['Ids'][0]]][cx:cx+3]
             Quad = rbObj['Orient'][0]
             data['Drawing']['Quaternion'] = G2mth.invQ(Quad)
             if rbType == 'Residue':
+                data['Drawing']['viewPoint'][0] = rbObj['Orig'][0]
                 G2frame.bottomSizer =  ResrbSizer(rbObj,rbIndx)
             elif rbType == 'Spin':
+                data['Drawing']['viewPoint'][0] = data['Atoms'][AtLookUp[rbObj['Ids'][0]]][cx:cx+3]
                 G2frame.bottomSizer =  SpnrbSizer(rbObj,rbIndx)
             else: #Vector
+                data['Drawing']['viewPoint'][0] = rbObj['Orig'][0]
                 G2frame.bottomSizer =  VecrbSizer(rbObj,rbIndx)
             mainSizer.Add(G2frame.bottomSizer)
             mainSizer.Layout()
@@ -14252,7 +14281,7 @@ def OnTorAngle(invalid,value,tc):
                 recompute atom distances
                 '''
                 [tor,torSlide] = Indx[tc.GetId()]
-                torSlide.SetValue(int(value*10))
+                torSlide.ChangeValue(int(value*10))
                 UpdateTablePlot()
 
             def OnTorSlide(event):
@@ -14264,7 +14293,7 @@ def OnTorSlide(event):
                 Tors = data['testRBObj']['rbObj']['Torsions'][tor]
                 val = float(Obj.GetValue())/10.
                 Tors[0] = val
-                ang.SetValue(val)
+                ang.ChangeValue(val)
                 UpdateTablePlot()
 
             def UpdateTable(event=None):
@@ -14300,7 +14329,7 @@ def OnAzSlide(event):
                 rbObj['OrientVec'][0] = float(Obj.GetValue())/10.
                 for i in range(4):
                     val = rbObj['OrientVec'][i]
-                    G2frame.testRBObjSizers['OrientVecSiz'][i].SetValue(val)
+                    G2frame.testRBObjSizers['OrientVecSiz'][i].ChangeValue(val)
                 Q = G2mth.AVdeg2Q(rbObj['OrientVec'][0],
                                 np.inner(Amat,rbObj['OrientVec'][1:]))
                 rbObj['Orient'][0] = Q
@@ -14317,7 +14346,7 @@ def UpdateOrientation(*args,**kwargs):
                     np.inner(Amat,rbObj['OrientVec'][1:]))
                 rbObj['Orient'][0] = Q
                 try:
-                    G2frame.testRBObjSizers['OrientVecSiz'][4].SetValue(
+                    G2frame.testRBObjSizers['OrientVecSiz'][4].ChangeValue(
                         int(10*rbObj['OrientVec'][0]))
                 except:
                     pass

GSASII/GSASIIplot.py

@@ -2733,7 +2733,7 @@ def OnMotion(event):
     Plot.set_title(Title)
     try:
         Page.figure.colorbar(mcolors,shrink=0.75,label='Rotation',boundaries=range(91))
-    except TypeError:
+    except:
         print('mpl error - no colorbar shown')
     Page.canvas.draw()
 
@@ -7418,8 +7418,10 @@ def Draw(caller='',Fade=[],NPkey=False):
                 Backbone = Backbones[chain]
                 RenderBackbone(Backbone,BackboneColor,bondR)
         if drawingData['showVoids']:
-            for x,y,z in drawingData['Voids']:
-                RenderSphere(x,y,z,.05,(0.,0.,1.),True)
+            RC = len(drawingData['Voids'])*[[0.05,2*Bl]]
+            RenderDots(drawingData['Voids'],RC)
+            # for x,y,z in drawingData['Voids']:
+            #     RenderSphere(x,y,z,.05,(0.,0.,1.),True)
         if drawingData['unitCellBox']:
             RenderBox()
             if drawingData['Plane'][1]:

GSASII/GSASIIpwdGUI.py

@@ -1647,7 +1647,13 @@ def OnBackSave(event):
                 filename = os.path.splitext(filename)[0]+'.pwdrbck'
                 File = open(filename,'w')
                 File.write("#GSAS-II background parameter file; do not add/delete items!\n")
-                File.write(str(data[0])+'\n')
+                out_tmp = list()
+                for item in data[0]:
+                    if isinstance(item, np.float64) or isinstance(item, float):
+                        out_tmp.append(float(item))
+                    else:
+                        out_tmp.append(item)
+                File.write(str(out_tmp)+'\n')
                 for item in data[1]:
                     if item in ['nPeaks','background PWDR','nDebye'] or not len(data[1][item]):
                         File.write(item+':'+str(data[1][item])+'\n')

GSASII/GSASIIstrMath.py

@@ -453,11 +453,9 @@ def MakePolar(Orient,QB):
 
     dFFdS = {}
     atFlg = []
-#    R,Th,Ph = G2lat.H2ThPh2(np.reshape(HKL,(-1,3)),Bmat)
     SQR = np.repeat(SQ,HKL.shape[1])
     for iAt,Atype in enumerate(Tdata):
-        if 'Q' in Atype:
-#            R,Th,Ph = G2lat.H2ThPh2(np.reshape(HKL,(-1,3)),Bmat)
+        if 'Q' in Atype:        #spinning RB
             atFlg.append(1.0)
             SHdat = SHCdict[iAt]
             symAxis = np.array(SHdat['symAxis'])
@@ -527,7 +525,7 @@ def MakePolar(Orient,QB):
                         SHMk = G2lat.KslCalc(item,ThMk,PhMk)
                         BS = 1.0
                         if 'Q' in Atm:
-                            BS = sp.spherical_jn(l,1.0)    #Slater term here?
+                            BS = sp.spherical_jn(l,1.0)/(4.*np.pi)    #Slater term here?
                         else:
                             BS = sp.spherical_jn(l,QR*R)/(4.*np.pi)    #Bessel function
                             BSP = sp.spherical_jn(l,QR*(R+0.01))/(4.*np.pi)
@@ -547,7 +545,7 @@ def MakePolar(Orient,QB):
             dFFdS[Oiname] = dSHdOi
             dFFdS[Ojname] = dSHdOj
             dFFdS[Okname] = dSHdOk
-        elif iAt in SHCdict and 'X' in hType:
+        elif iAt in SHCdict and 'X' in hType:   #X-ray deformation
             radial = SHCdict[-iAt]['Radial']
             orKeys = [item for item in ORBtables[Atype] if item not in ['Slater','ZSlater','NSlater','SZE','popCore','popVal']]
             if 'B' in radial: