Electrostatic Potential (Quantum Espresso) - Courtesy of Andrew Doyle
####################################################### ## Tutorial by: Andrew ## Topic: Electrostatic potential in ase-espresso ## ## dependencies: ase, espresso, pickle ## ## Generates the data files necessary to analyze ## the charge density and electrostatic potential ## in and around a Pt(100) slab ####################################################### # Import necessary packages from ase import Atoms from ase.io import write from ase.lattice.surface import fcc100 from espresso import espresso import pickle # Generate the representation of the atoms in ASE slab = fcc100('Pt', size=(2,2,3), vacuum=10.0) # Create an espresso calculator object # (Tailor for desired settings and accuracy) calc = espresso(pw=400 , # Plane wave cutoff dw=4000, # Density wave cutoff kpts=(3,3,1), # (rather sparse) k-point (Brillouin) sampling xc='RPBE', #Exchange-correlation functional dipole={'status':True}, # Includes dipole correction (necessary for asymmetric slabs) outdir='pt100_analysis') # Espresso-generated files will be put here # Attach the calculator object to the atoms slab.set_calculator(calc) # Perform a static calculation to generate wavefunctions slab.get_potential_energy() # Create an image of the atoms write('slab.png',slab) # Get the electrostatic potential potential = calc.extract_ionic_and_hartree_potential() # Create and populate output file for electrostatic potential potential_file=open('potential.pickle','w') pickle.dump(potential,potential_file) potential_file.close()
Now we are done with the intensive calculations. The following scripts will interpret the results and can be run easily on any computer in a matter of seconds.
####################################################### ## Tutorial by: Andrew ## Topic: Plotting electrostatic potential data from ## ase-espresso ## ## dependencies: pickle, pylab, numpy ## ## Plots an XY-averaged electrostatic potential as ## a function of Z position through and outside ## a Pt(100) slab. ####################################################### # Import necessary packages import pickle import pylab as plt import numpy as np # Retrieves the data generated in the previous script potential_file = open('potential.pickle') potential = pickle.load(potential_file) origin=potential[0] # [X0,Y0,Z0] - The positions of the origin of the cell. Here they are [0,0,0] cell = potential[1] # [3x3 Matrix] Contains the [x,y,z] components of the 3 cell vectors cell_x = cell[0][0] # The x component of the first cell vector (Remember: This cell is rectangular) cell_y = cell[1][1] # The y component of the second cell vector (Remember: This cell is rectangular) cell_z = cell[2][2] # The z component of the third cell vector (Remember: This cell is rectangular) v = potential[2] # Contains all real-space potential field points. ####################################################### ## Understanding the architecture of the rest of the ## data can be tricky. Using the variables presented ## here, V(x,y,z)=v[x][y][z]. ## ## First we aim to extract all the information ## available to us, based on numpy's architecture. ####################################################### n_x = len(v) # Returns the number of entries in the first subset (number of X entries) n_y = len(v[0]) # Returns the number of entries in the second subset (number of Y entries - choice of '0' was arbitrary) n_z = len(v[0][0]) # Returns the number of entries in the third subset (number of Z entries - choice of '0' was arbitrary) potential = [] # Initializes a list for average potentials # Iterate over all space to average and process potential information for z in range(n_z): temp = 0 # Will serve as a partial sum for averaging purposes for x in range(n_x): for y in range(n_y): temp+=v[x][y][z] # Sums the potential at all points in an X-Y plane temp = temp / (n_x*n_y) # Divides by the total number of points potential.append( temp ) # Stores the value for average potential at this Z position # Translates the indices of points along the Z axis into real-space positions z_coordinate = np.linspace( 0 , cell_z , n_z ) # Plots the data, formats plot plt.plot(z_coordinate, potential) plt.xlabel(r'z ($\AA$)',size=18) plt.ylabel(r'$\langle V \ \rangle_{XY}$ (Volts)',size=18) plt.title('Electrostatics of Pt(100)',size=18) # Display plot plt.show() # If you'd prefer to save the plot, highlight the line above and uncomment the line below #plt.savefig('potential_analysis.png')
The end result should look something like this:
Plot self-consistent field convergence data
This script, when run in a quantum espresso job folder, will plot the estimated SCF accuracy and total DFT energy as a function of electronic step. This can be useful for diagnosing weird convergence behavior, or for impatiently waiting for the next ionic step to finish.
#!/home/vossj/suncat/bin/python import os import pylab as plt import sys import matplotlib.pyplot as plt2 from mpl_toolkits.axes_grid1 import host_subplot import mpl_toolkits.axisartist as AA rydberg = 13.6057 #rydberg to eV conversion # default to outdir/log if no log is given if len(sys.argv) > 1: log = sys.argv[1] else: log = 'outdir/log' filename = 'scf_temp.txt' filename2 = 'etmp.txt' os.system('grep scf %s > %s'%(log,filename)) os.system("grep 'total energy' %s > %s"%(log,filename2)) file = open(filename) lines = file.read().split('\n') lines = filter(None,lines) file2 = open(filename2) lines2 = file2.read().split('\n') lines2 = filter(None,lines2) scf = [] energy = [] iter = [] iter2 = [] n = 0 # append scf convergence data from log file to scf for line in lines: items = filter(None, line.split(' ')) try: scf.append( float(items[4])*rydberg ) n += 1 iter.append(n) except: continue # append total dft energy data from log to energy n = 0 for line in lines2: items = filter(None, line.split(' ')) try: n += 1 energy.append(float(items[3])*rydberg ) iter2.append(n) except: continue # find the plane-wave cutoff from the log file os.system("grep 'kinetic-energy cutoff' %s > %s"%(log,filename)) file = open(filename) items = filter(None, file.read().split(' ')) pw = float(items[3]) * 13.606 os.system('rm %s'%filename) os.system('rm %s'%filename2) # plot everything nicely host = host_subplot(111,axes_class = AA.Axes) plt2.subplots_adjust(right=0.75) par1 = host.twinx() host.set_xlabel('Number of Iterations') host.set_ylabel('Estimated SCF Accuracy (eV)',color='red') par1.set_ylabel('DFT Total-Energy (eV)',color='green') p1, = host.plot(iter,scf,label='SCF Acc') p2, = par1.plot(iter2,energy,label='DFT Energy') host.set_yscale('log') plt.semilogy(iter,scf) plt.xlabel('Number of Iterations') plt.ylabel('Estimated SCF Accuracy (eV)') plt.title(r'$PW_{Cut}$=%d eV'%pw) plt.tight_layout() plt.show()
The output should look something like this, for a nicely converging run with several ionic steps. Multiple ionic step directories will contain all ionic steps in the plot.
