Changeset 4366
 Timestamp:
 Mar 13, 2020 5:44:04 AM (4 years ago)
 Location:
 trunk
 Files:

 3 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIIphsGUI.py
r4355 r4366 5362 5362 i += 1 5363 5363 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 5365 5373 G2frame.OnFileSave(event) 5366 5374 print (' GSASII project saved') … … 5479 5487 Partials = np.array(Partials).T 5480 5488 if 'Q' in label: 5489 continue #skip these partials 5481 5490 XY = [[X.T,Y.T] for iy,Y in enumerate(Partials) if 'Va' not in Names[iy+1]] 5482 5491 else: … … 5491 5500 title = 'Neutron '+Labels[label][2]+pName 5492 5501 else: 5502 continue #skip for now  xray partials are missing header record 5493 5503 title = 'Xray '+Labels[label][2].replace('G','g')+pName 5494 5504 ylabel = 'g(R)' … … 5546 5556 labelY=r'$\mathsf{\chi^2}$',newPlot=True,Title='RMCP Chi^2 for '+pName, 5547 5557 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 5548 5631 #bond odf plots 5549 5632 nPot = len(RMCPdict['Potentials']['Stretch']) … … 5559 5642 numx,numy = odfData[:2] 5560 5643 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) 5562 5645 OutFile.close() 5563 5646 
trunk/GSASIIplot.py
r4356 r4366 5380 5380 ################################################################################ 5381 5381 5382 def PlotBarGraph(G2frame,Xarray,Xname='', Title='',PlotName=None):5382 def PlotBarGraph(G2frame,Xarray,Xname='',Yname='Number',Title='',PlotName=None,ifBinned=False): 5383 5383 'Needs a description' 5384 5384 … … 5405 5405 Page.canvas.mpl_connect('motion_notify_event', OnMotion) 5406 5406 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] 5410 5416 Plot.set_title(Title) 5411 5417 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') 5414 5420 Page.canvas.draw() 5415 5421 
trunk/GSASIIpwd.py
r4347 r4366 2383 2383 fl.write(' > ANGLE_SEARCH :: %.1f%%\n'%RMCPdict['Potentials']['Angle search']) 2384 2384 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'% 2386 2386 (angle[1],angle[0],angle[2],angle[6],angle[3],angle[4],angle[5])) 2387 2387 if RMCPdict['useBVS']: … … 2503 2503 fl.close() 2504 2504 return fname 2505 2506 def 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[tgt1][1:] for tgt in bonds[1]]) 2518 Dx = np.array([TxyzOxyz+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 2524 def 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[tgt1][1:] for tgt in angle[1].T[0]]) 2536 TBxyz = np.array([Atoms[tgt1][1:] for tgt in angle[1].T[1]]) 2537 DAxV = np.inner(np.array([TAxyzOxyz+unit for unit in Units]),Amat) 2538 DAx = np.sqrt(np.sum(DAxV**2,axis=2)) 2539 DBxV = np.inner(np.array([TBxyzOxyz+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) 2505 2548 2506 2549 ################################################################################
Note: See TracChangeset
for help on using the changeset viewer.