Changeset 4366


Ignore:
Timestamp:
Mar 13, 2020 5:44:04 AM (19 months ago)
Author:
vondreele
Message:

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

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIphsGUI.py

    r4355 r4366  
    53625362                    i += 1
    53635363                else:
    5364                     break               
     5364                    break
     5365            i = 1
     5366            while True:
     5367                if os.path.isfile(pName+'_anglehist_%d.csv'%i):
     5368                    os.remove(pName+'_anglehist_%d.csv'%i)
     5369                    i += 1
     5370                else:
     5371                    break
     5372               
    53655373            G2frame.OnFileSave(event)
    53665374            print (' GSAS-II project saved')
     
    54795487                    Partials = np.array(Partials).T
    54805488                    if 'Q' in label:
     5489                        continue            #skip these partials
    54815490                        XY = [[X.T,Y.T] for iy,Y in enumerate(Partials) if 'Va' not in Names[iy+1]]
    54825491                    else:
     
    54915500                            title = 'Neutron '+Labels[label][2]+pName
    54925501                        else:
     5502                            continue        #skip for now - x-ray partials are missing header record
    54935503                            title = 'X-ray '+Labels[label][2].replace('G','g')+pName
    54945504                            ylabel = 'g(R)'
     
    55465556                    labelY=r'$\mathsf{\chi^2}$',newPlot=True,Title='RMCP Chi^2 for '+pName,
    55475557                    lines=True,names=Names[3:])
     5558
     5559#get atoms from rmc6f file
     5560            rmc6fName = pName+'.rmc6f'
     5561            rmc6f = open(rmc6fName,'r')
     5562            rmc6fAtoms = []
     5563            while True:
     5564                line = rmc6f.readline()
     5565                if 'Number of atoms:' in line:
     5566                    Natm = int(line.split(':')[1])
     5567                if 'Atoms:' in line:
     5568                    break
     5569            for iAtm in range(Natm):
     5570                line = rmc6f.readline().split()
     5571                rmc6fAtoms.append([line[1],float(line[3]),float(line[4]),float(line[5])])
     5572            rmc6f.close()
     5573#alt bond histograms - from rmc6 & bond files
     5574               
     5575            bondName = pName+'.bonds'
     5576            if os.path.exists(os.path.join(path,bondName)):               
     5577                nBond = len(RMCPdict['Potentials']['Stretch'])
     5578                bondList = []
     5579                bonds = open(bondName,'r')
     5580                while True:
     5581                    line = bonds.readline()
     5582                    if '............' in line:
     5583                        break       #positions at start of 1st bond list
     5584                for iBond in range(nBond):
     5585                    bondList.append([])
     5586                    Bonds = RMCPdict['Potentials']['Stretch'][iBond]
     5587                    for iAtm in range(Natm):     #same as in rmc6f above
     5588                        line = bonds.readline().split('::')
     5589                        if Bonds[0] in line[0]:
     5590                            items = line[1].replace(';','').split()[:-1:2]
     5591                            items = np.array([int(item) for item in items])
     5592                            bondList[iBond].append([iAtm,items])
     5593                bonds.close()
     5594                for iBond in range(nBond):
     5595                    Bonds = RMCPdict['Potentials']['Stretch'][iBond]
     5596                    title = '%s-%s'%(Bonds[0],Bonds[1])
     5597                    bondDist = G2pwd.GetRMCBonds(generalData,RMCPdict,rmc6fAtoms,bondList[iBond])
     5598                    print('%s mean %.3f(%d)'%(title,np.mean(bondDist),int(1000*np.std(bondDist))))
     5599                    G2plt.PlotBarGraph(G2frame,bondDist,Xname=r'$Bond, \AA$',Title=title+' from Potential Energy Restraint',
     5600                        PlotName='%s Bond for %s'%(title,pName))
     5601               
     5602#alt angle histograms - from rmc6 & triplets files
     5603            tripName = pName+'.triplets'
     5604            if os.path.exists(os.path.join(path,tripName)):
     5605                nAng = len(RMCPdict['Potentials']['Angles'])
     5606                tripList = []
     5607                triples = open(tripName,'r')
     5608                while True:
     5609                    line = triples.readline()
     5610                    if '............' in line:
     5611                        break       #positions at start of 1st triple list
     5612                for iAng in range(nAng):
     5613                    tripList.append([])
     5614                    Angles = RMCPdict['Potentials']['Angles'][iAng]
     5615                    for iAtm in range(Natm):     #same as in rmc6f above
     5616                        line = triples.readline().split('::')
     5617                        if Angles[1] in line[0]:
     5618                            items = line[1].replace(';','').split()[:-1:2]
     5619                            items = np.array([int(item) for item in items]).reshape((-1,2))
     5620                            tripList[iAng].append([iAtm,items])
     5621                        line = triples.readline()       #to skip a line
     5622                triples.close()
     5623                for iAng in range(nAng):
     5624                    Angles = RMCPdict['Potentials']['Angles'][iAng]
     5625                    angles = G2pwd.GetRMCAngles(generalData,RMCPdict,rmc6fAtoms,tripList[iAng])
     5626                    title = '%s-%s-%s'%(Angles[0],Angles[1],Angles[2])
     5627                    print('%s mean %.2f(%d)'%(title,np.mean(angles),int(100*np.std(angles))))
     5628                    G2plt.PlotBarGraph(G2frame,angles,Xname=r'$Angle, \AA$',Title=title+' from Potential Energy Restraint',
     5629                        PlotName='%s Angle for %s'%(title,pName))
     5630                               
    55485631#bond odf plots               
    55495632            nPot = len(RMCPdict['Potentials']['Stretch'])
     
    55595642                        numx,numy = odfData[:2]
    55605643                        G2plt.Plot3dXYZ(G2frame,int(numx),int(numy),odfData[2:],
    5561                             newPlot=False,Title='Bond %s-%s'%(bond[0],bond[1]),Centro=True) 
     5644                            newPlot=False,Title='Number of %s-%s Bonds'%(bond[0],bond[1]),Centro=True) 
    55625645                    OutFile.close()
    55635646       
  • trunk/GSASIIplot.py

    r4356 r4366  
    53805380################################################################################
    53815381           
    5382 def PlotBarGraph(G2frame,Xarray,Xname='',Title='',PlotName=None):
     5382def PlotBarGraph(G2frame,Xarray,Xname='',Yname='Number',Title='',PlotName=None,ifBinned=False):
    53835383    'Needs a description'
    53845384   
     
    54055405        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
    54065406    Page.Choice = None
    5407     nBins= max(10,len(Xarray)//10)
    5408     Bins,Dbins = np.histogram(Xarray,nBins)
    5409     wid = Dbins[1]-Dbins[0]
     5407    if ifBinned:
     5408        Dbins,Bins = Xarray
     5409        nBins = len(Dbins)
     5410        wid = Dbins[1]-Dbins[0]
     5411    else:
     5412        nBins= max(10,len(Xarray)//10)
     5413        Bins,Dbins = np.histogram(Xarray,nBins)
     5414        wid = Dbins[1]-Dbins[0]
     5415        Dbins = Dbins[:-1]
    54105416    Plot.set_title(Title)
    54115417    Plot.set_xlabel(Xname,fontsize=14)
    5412     Plot.set_ylabel(r'$Number$',fontsize=14)
    5413     Plot.bar(Dbins[:-1],Bins,width=wid,align='edge',facecolor='red',edgecolor='black')
     5418    Plot.set_ylabel(Yname,fontsize=14)
     5419    Plot.bar(Dbins,Bins,width=wid,align='edge',facecolor='red',edgecolor='black')
    54145420    Page.canvas.draw()
    54155421
  • trunk/GSASIIpwd.py

    r4347 r4366  
    23832383            fl.write('  > ANGLE_SEARCH :: %.1f%%\n'%RMCPdict['Potentials']['Angle search'])
    23842384            for angle in RMCPdict['Potentials']['Angles']:
    2385                 fl.write('  > ANGLE :: %s %s %s %.2f eV %.2f deg %.2f Ang %.2f Ang\n'%
     2385                fl.write('  > ANGLE :: %s %s %s %.2f eV %.2f deg %.2f %.2f Ang\n'%
    23862386                    (angle[1],angle[0],angle[2],angle[6],angle[3],angle[4],angle[5]))
    23872387    if RMCPdict['useBVS']:
     
    25032503    fl.close()
    25042504    return fname
     2505
     2506def GetRMCBonds(general,RMCPdict,Atoms,bondList):
     2507    bondDist = []
     2508    Cell = general['Cell'][1:7]
     2509    Supercell =  RMCPdict['SuperCell']
     2510    Trans = np.eye(3)*np.array(Supercell)
     2511    Cell = G2lat.TransformCell(Cell,Trans)[:6]
     2512    Amat,Bmat = G2lat.cell2AB(Cell)
     2513    indices = (-1,0,1)
     2514    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
     2515    for bonds in bondList:
     2516        Oxyz = np.array(Atoms[bonds[0]][1:])
     2517        Txyz = np.array([Atoms[tgt-1][1:] for tgt in bonds[1]])       
     2518        Dx = np.array([Txyz-Oxyz+unit for unit in Units])
     2519        Dx = np.sqrt(np.sum(np.inner(Dx,Amat)**2,axis=2))
     2520        for dx in Dx.T:
     2521            bondDist.append(np.min(dx))
     2522    return np.array(bondDist)
     2523   
     2524def GetRMCAngles(general,RMCPdict,Atoms,angleList):
     2525    bondAngles = []
     2526    Cell = general['Cell'][1:7]
     2527    Supercell =  RMCPdict['SuperCell']
     2528    Trans = np.eye(3)*np.array(Supercell)
     2529    Cell = G2lat.TransformCell(Cell,Trans)[:6]
     2530    Amat,Bmat = G2lat.cell2AB(Cell)
     2531    indices = (-1,0,1)
     2532    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
     2533    for angle in angleList:
     2534        Oxyz = np.array(Atoms[angle[0]][1:])
     2535        TAxyz = np.array([Atoms[tgt-1][1:] for tgt in angle[1].T[0]])       
     2536        TBxyz = np.array([Atoms[tgt-1][1:] for tgt in angle[1].T[1]])       
     2537        DAxV = np.inner(np.array([TAxyz-Oxyz+unit for unit in Units]),Amat)
     2538        DAx = np.sqrt(np.sum(DAxV**2,axis=2))
     2539        DBxV = np.inner(np.array([TBxyz-Oxyz+unit for unit in Units]),Amat)
     2540        DBx = np.sqrt(np.sum(DBxV**2,axis=2))
     2541        iDAx = np.argmin(DAx,axis=0)
     2542        iDBx = np.argmin(DBx,axis=0)
     2543        for i,[iA,iB] in enumerate(zip(iDAx,iDBx)):
     2544            DAv = DAxV[iA,i]/DAx[iA,i]
     2545            DBv = DBxV[iB,i]/DBx[iB,i]
     2546            bondAngles.append(npacosd(np.sum(DAv*DBv)))
     2547    return np.array(bondAngles)
    25052548   
    25062549################################################################################
Note: See TracChangeset for help on using the changeset viewer.