Changeset 4366 for trunk/GSASIIphsGUI.py


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

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

File:
1 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       
Note: See TracChangeset for help on using the changeset viewer.