Changeset 4366

Timestamp:

Mar 13, 2020 5:44:04 AM (4 years ago)

Author:

vondreele

Message:

Add plotting of bond distance & bond angle distributions found in potential energy restraint in RMCProfile

Location:

trunk

Files:

• GSASIIphsGUI.py (modified) (5 diffs)
• GSASIIplot.py (modified) (2 diffs)
• GSASIIpwd.py (modified) (2 diffs)

trunk/GSASIIphsGUI.py

@@ -5362 +5362 @@
                     i += 1
                 else:
-                    break                
+                    break
+            i = 1
+            while True:
+                if os.path.isfile(pName+'_anglehist_%d.csv'%i):
+                    os.remove(pName+'_anglehist_%d.csv'%i)
+                    i += 1
+                else:
+                    break
+                
             G2frame.OnFileSave(event)
             print (' GSAS-II project saved')
@@ -5479 +5487 @@
                     Partials = np.array(Partials).T
                     if 'Q' in label:
+                        continue            #skip these partials
                         XY = [[X.T,Y.T] for iy,Y in enumerate(Partials) if 'Va' not in Names[iy+1]]
                     else:
@@ -5491 +5500 @@
                             title = 'Neutron '+Labels[label][2]+pName
                         else:
+                            continue        #skip for now - x-ray partials are missing header record
                             title = 'X-ray '+Labels[label][2].replace('G','g')+pName
                             ylabel = 'g(R)'
@@ -5546 +5556 @@
                     labelY=r'$\mathsf{\chi^2}$',newPlot=True,Title='RMCP Chi^2 for '+pName,
                     lines=True,names=Names[3:])
+
+#get atoms from rmc6f file 
+            rmc6fName = pName+'.rmc6f'
+            rmc6f = open(rmc6fName,'r')
+            rmc6fAtoms = []
+            while True:
+                line = rmc6f.readline()
+                if 'Number of atoms:' in line:
+                    Natm = int(line.split(':')[1])
+                if 'Atoms:' in line:
+                    break
+            for iAtm in range(Natm):
+                line = rmc6f.readline().split()
+                rmc6fAtoms.append([line[1],float(line[3]),float(line[4]),float(line[5])])
+            rmc6f.close()
+#alt bond histograms - from rmc6 & bond files
+               
+            bondName = pName+'.bonds'
+            if os.path.exists(os.path.join(path,bondName)):                
+                nBond = len(RMCPdict['Potentials']['Stretch'])
+                bondList = []
+                bonds = open(bondName,'r')
+                while True:
+                    line = bonds.readline()
+                    if '............' in line:
+                        break       #positions at start of 1st bond list
+                for iBond in range(nBond):
+                    bondList.append([])
+                    Bonds = RMCPdict['Potentials']['Stretch'][iBond]
+                    for iAtm in range(Natm):     #same as in rmc6f above
+                        line = bonds.readline().split('::')
+                        if Bonds[0] in line[0]:
+                            items = line[1].replace(';','').split()[:-1:2]
+                            items = np.array([int(item) for item in items])
+                            bondList[iBond].append([iAtm,items])
+                bonds.close()
+                for iBond in range(nBond):
+                    Bonds = RMCPdict['Potentials']['Stretch'][iBond]
+                    title = '%s-%s'%(Bonds[0],Bonds[1])
+                    bondDist = G2pwd.GetRMCBonds(generalData,RMCPdict,rmc6fAtoms,bondList[iBond])
+                    print('%s mean %.3f(%d)'%(title,np.mean(bondDist),int(1000*np.std(bondDist))))
+                    G2plt.PlotBarGraph(G2frame,bondDist,Xname=r'$Bond, \AA$',Title=title+' from Potential Energy Restraint',
+                        PlotName='%s Bond for %s'%(title,pName))
+                
+#alt angle histograms - from rmc6 & triplets files
+            tripName = pName+'.triplets'
+            if os.path.exists(os.path.join(path,tripName)):
+                nAng = len(RMCPdict['Potentials']['Angles'])
+                tripList = []
+                triples = open(tripName,'r')
+                while True:
+                    line = triples.readline()
+                    if '............' in line:
+                        break       #positions at start of 1st triple list
+                for iAng in range(nAng):
+                    tripList.append([])
+                    Angles = RMCPdict['Potentials']['Angles'][iAng]
+                    for iAtm in range(Natm):     #same as in rmc6f above
+                        line = triples.readline().split('::')
+                        if Angles[1] in line[0]:
+                            items = line[1].replace(';','').split()[:-1:2]
+                            items = np.array([int(item) for item in items]).reshape((-1,2))
+                            tripList[iAng].append([iAtm,items])
+                        line = triples.readline()       #to skip a line
+                triples.close()
+                for iAng in range(nAng):
+                    Angles = RMCPdict['Potentials']['Angles'][iAng]
+                    angles = G2pwd.GetRMCAngles(generalData,RMCPdict,rmc6fAtoms,tripList[iAng])
+                    title = '%s-%s-%s'%(Angles[0],Angles[1],Angles[2])
+                    print('%s mean %.2f(%d)'%(title,np.mean(angles),int(100*np.std(angles))))
+                    G2plt.PlotBarGraph(G2frame,angles,Xname=r'$Angle, \AA$',Title=title+' from Potential Energy Restraint',
+                        PlotName='%s Angle for %s'%(title,pName))
+                                
 #bond odf plots                
             nPot = len(RMCPdict['Potentials']['Stretch'])
@@ -5559 +5642 @@
                         numx,numy = odfData[:2]
                         G2plt.Plot3dXYZ(G2frame,int(numx),int(numy),odfData[2:],
-                            newPlot=False,Title='Bond %s-%s'%(bond[0],bond[1]),Centro=True)  
+                            newPlot=False,Title='Number of %s-%s Bonds'%(bond[0],bond[1]),Centro=True)  
                     OutFile.close()
         

trunk/GSASIIplot.py

@@ -5380 +5380 @@
 ################################################################################
             
-def PlotBarGraph(G2frame,Xarray,Xname='',Title='',PlotName=None):
+def PlotBarGraph(G2frame,Xarray,Xname='',Yname='Number',Title='',PlotName=None,ifBinned=False):
     'Needs a description'
     
@@ -5405 +5405 @@
         Page.canvas.mpl_connect('motion_notify_event', OnMotion)
     Page.Choice = None
-    nBins= max(10,len(Xarray)//10)
-    Bins,Dbins = np.histogram(Xarray,nBins)
-    wid = Dbins[1]-Dbins[0]
+    if ifBinned:
+        Dbins,Bins = Xarray
+        nBins = len(Dbins)
+        wid = Dbins[1]-Dbins[0]
+    else:
+        nBins= max(10,len(Xarray)//10)
+        Bins,Dbins = np.histogram(Xarray,nBins)
+        wid = Dbins[1]-Dbins[0]
+        Dbins = Dbins[:-1]
     Plot.set_title(Title)
     Plot.set_xlabel(Xname,fontsize=14)
-    Plot.set_ylabel(r'$Number$',fontsize=14)
-    Plot.bar(Dbins[:-1],Bins,width=wid,align='edge',facecolor='red',edgecolor='black')
+    Plot.set_ylabel(Yname,fontsize=14)
+    Plot.bar(Dbins,Bins,width=wid,align='edge',facecolor='red',edgecolor='black')
     Page.canvas.draw()
 

trunk/GSASIIpwd.py

@@ -2383 +2383 @@
             fl.write('  > ANGLE_SEARCH :: %.1f%%\n'%RMCPdict['Potentials']['Angle search'])
             for angle in RMCPdict['Potentials']['Angles']:
-                fl.write('  > ANGLE :: %s %s %s %.2f eV %.2f deg %.2f Ang %.2f Ang\n'%
+                fl.write('  > ANGLE :: %s %s %s %.2f eV %.2f deg %.2f %.2f Ang\n'%
                     (angle[1],angle[0],angle[2],angle[6],angle[3],angle[4],angle[5]))
     if RMCPdict['useBVS']:
@@ -2503 +2503 @@
     fl.close()
     return fname
+
+def GetRMCBonds(general,RMCPdict,Atoms,bondList):
+    bondDist = []
+    Cell = general['Cell'][1:7]
+    Supercell =  RMCPdict['SuperCell']
+    Trans = np.eye(3)*np.array(Supercell)
+    Cell = G2lat.TransformCell(Cell,Trans)[:6]
+    Amat,Bmat = G2lat.cell2AB(Cell)
+    indices = (-1,0,1)
+    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
+    for bonds in bondList:
+        Oxyz = np.array(Atoms[bonds[0]][1:])
+        Txyz = np.array([Atoms[tgt-1][1:] for tgt in bonds[1]])        
+        Dx = np.array([Txyz-Oxyz+unit for unit in Units])
+        Dx = np.sqrt(np.sum(np.inner(Dx,Amat)**2,axis=2))
+        for dx in Dx.T:
+            bondDist.append(np.min(dx))
+    return np.array(bondDist)
+    
+def GetRMCAngles(general,RMCPdict,Atoms,angleList):
+    bondAngles = []
+    Cell = general['Cell'][1:7]
+    Supercell =  RMCPdict['SuperCell']
+    Trans = np.eye(3)*np.array(Supercell)
+    Cell = G2lat.TransformCell(Cell,Trans)[:6]
+    Amat,Bmat = G2lat.cell2AB(Cell)
+    indices = (-1,0,1)
+    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
+    for angle in angleList:
+        Oxyz = np.array(Atoms[angle[0]][1:])
+        TAxyz = np.array([Atoms[tgt-1][1:] for tgt in angle[1].T[0]])        
+        TBxyz = np.array([Atoms[tgt-1][1:] for tgt in angle[1].T[1]])        
+        DAxV = np.inner(np.array([TAxyz-Oxyz+unit for unit in Units]),Amat)
+        DAx = np.sqrt(np.sum(DAxV**2,axis=2))
+        DBxV = np.inner(np.array([TBxyz-Oxyz+unit for unit in Units]),Amat)
+        DBx = np.sqrt(np.sum(DBxV**2,axis=2))
+        iDAx = np.argmin(DAx,axis=0)
+        iDBx = np.argmin(DBx,axis=0)
+        for i,[iA,iB] in enumerate(zip(iDAx,iDBx)):
+            DAv = DAxV[iA,i]/DAx[iA,i]
+            DBv = DBxV[iB,i]/DBx[iB,i]
+            bondAngles.append(npacosd(np.sum(DAv*DBv)))
+    return np.array(bondAngles)
     
 ################################################################################