5.2. Running QE on PACE#
In this lecture, we will learn about running DFT calculations on PACE. We will first learn about SLURM submission scripts that are used to submit jobs on PACE. As a member of Medford group, you will have access to two supercomputer clusters here at Georgia Tech: PACE-Phoenix and PACE-Hive. If you are part of the VIP course, you will have access to PACE-ICE cluster. You shouldn’t use “#SBATCH -A GT-amedford6” if you are on PACE-ICE which is requesting resources from your PI’s account. If you are on PACE-Phoenix, you will need to request for resources from your PI’s account. This means that you are paying for the jobs you are running therefore, these resources should be used judiciously. An example run.sh script to submit a QE job is given below:
#!/bin/bash
#SBATCH -N 1 --ntasks-per-node=8
#SBATCH --mem-per-cpu=3G
#SBATCH -t12:00:00
#SBATCH -A gts-amedford6-joe
#SBATCH -J optimizer
#SBATCH -o Report-%j.out -e Report_err-%j.err
cd $SLURM_SUBMIT_DIR
source ~/.bashrc
module load anaconda3 #It is optional to load anaconda. You don't need to load it if it's already been loaded in the environment
# This path will change depending on the environment you wish to load
source /storage/coda1/p-amedford6/0/shared/rich_project_chbe-medford/medford-share/envs/espresso-6.7MaX-beef-ipi
python calc.py
To submit a job to the SLURM system, you have to use the qsub command. We usually call this file run.sh for convenience. The command for submitting the job will be qsub run.sh. When you do this, you are essentially handing the file to the SLURM queuing system which will scan the file for the blocks containing #SBATCH that specify the resources you are requesting for the job. Each line in this file starting with #SBATCH specifies something different. The line with tag -N are requesting certain resources (nodes and processors per node). You can request the resources based on requirement of your calculations. In this example script, we are requesting for 1 node and 8 processors per node. The next line (#SBATCH –mem-per-cpu=3G) specifies the memory you request per processor. Here it is 3 GB. Hence the total memory we request for is 24 GB. For Phoenix cluster, the default queue is inferno where the jobs are run. -J is the name you want to give the job, -o and -e are the filenames where the outputs and errors will be printed respectively. You can run pace-whoami check the maximum resources available for a particular queue.
If you are submitting a job to PACE-HIVE, you will need to specify the queue using #SBATCH -q. HIVE has different queues like hive, hive-himem, hive-interact which offer different resources and walltime. After the SBATCH block, you can source your environment and run the code. The next line with change directory (cd) command with an alias for current working directory (SLURM_SUBMIT_DIR) changes to the correct folder. More information about SLURM script and submission of jobs on PACE can be found in this link.
5.2.1. Running Quantum Espresso on PACE Phoenix#
For the purpose of DFT calculations, our group primarily uses Quantum Espresso (QE) which is a plane-wave method code. Various versions for QE are present in the shared folder of Medford group. The environment can be sourced in the bash file. An example command to source QE environment is:
source /storage/coda1/p-amedford6/0/shared/rich_project_chbe-medford/medford-share/envs/espresso-6.7Max-beef-ipi
Version 6.7Max is the environment for latest compiled QE code. We will first go through the details a simple code for running a DFT calculation using Quantum Espresso.
#example code for running Quantum Espresso
#import all the necessary modules
from ase import Atom, Atoms
from ase.optimize import BFGSLineSearch as BFGS
from ase.build import molecule
from espresso import iEspresso as Espresso
from ase.io import read, write
import numpy as np
import time
t0 = time.time()
images = []
#building the "molecule" object
image = molecule('H2O',vacuum=10.)
image.set_cell([10, 10, 10])
image.set_pbc([1,1,1])
image.center()
image.rattle(0.0001)
#Use QE as a single-point calculator (just SCF)
calc = Espresso(atoms=image,
pw=500.0,
xc='PBE',
kpts="gamma")
# retrieve Energy, force, stress
print('1st calculation with PBE')
print('Energy = {}'.format(image.get_potential_energy()))
print('Force = {}'.format(image.get_forces()))
print('Stress = {}'.format(image.get_stress()))
t1 = time.time()
print('Time = {} sec'.format(t1-t0))
print('-'*20)
# Run optimization (will change the atomic positions as it is relaxation calculation)
dyn = BFGS(image,logfile='opt_pbe.log', trajectory='h2o_pbe_optimization.traj')
dyn.run(fmax=0.005)
image.calc.close()
# Attach BEEF-vdw calculator
t2 = time.time()
calc = Espresso(atoms=image,
pw=350.0,
xc='BEEF-vdw')
# retrieve Energy, force, stress after optimization
print('2nd calculation with BEEF-vdw')
print('Energy = {}'.format(image.get_potential_energy()))
print('Force = {}'.format(image.get_forces()))
print('Stress = {}'.format(image.get_stress()))
E_ensemble = image.get_ensemble_energies()
print('Ensemble energies = {} +/- {}'.format(E_ensemble.mean(), E_ensemble.std()))
print('Time = {} sec'.format(time.time()-t2))
# Run optimization with BEEF
dyn = BFGS(image,logfile='opt_beef.log',trajectory='h2o_beef_optimization.traj')
dyn.run(fmax=0.005)
image.calc.close()
print('Totol time = {} sec'.format(time.time()-t0))
Now let’s take a look at this codeblock line by line to understand what is happening as we run this code. For setting up a DFT calculation, we first need to build the chemical system whose properties we want to calculate. In the example given, the object is a water molecule. By building the system, we are essentially specifying its atomic positions. Since it is a water molecule, it has to be put in a isolated box with 10 Angstrom vacuum in all three directions so that electron density of the molecule doesn’t overlap with that of its images. The molecule should be located at the center of the box. We can verify this by using ASE’s visualize method as we learnt in the previous lectures.
#building the "molecule" object
image = molecule('H2O',vacuum=10.)
image.set_cell([10, 10, 10])
image.set_pbc([1,1,1])
image.center()
image.rattle(0.0001)
After defining the system, we will set up a calculator where we define the parameters for the DFT calculation. These parameters will vary for different DFT codes. The important parameters that you will need to specify for your system in Quantum Espresso calculator are the following:
atoms: The Atoms object
pw: Plane-wave cutoff
kpts: K-points
xc: Exchange-correlation functional
spinpol: Spin-polarization
dipole: Dipole-correction
You will also need to decide if you want to run a single-point calculation or relaxation calculation. A single-point calculation might result in very high forces is most cases as the atoms in the initial guess might not be relaxed. The first step in the example script is to run a single-point calculation (also known as SCF calculation). SCF stands for self-consistent field. Note that the calculator also takes in the atoms object so that it knows where to get the data from.
calc = Espresso(atoms=image,
pw=500.0,
xc='PBE',
kpts="gamma")
print('1st calculation with PBE')
print('Energy = {}'.format(image.get_potential_energy()))
print('Force = {}'.format(image.get_forces()))
print('Stress = {}'.format(image.get_stress()))
t1 = time.time()
print('Time = {} sec'.format(t1-t0))
print('-'*20)
Some of the properties you want to obtain from a DFT calculation can be:
Energy and forces
Density of states
Electron density
Electrostatic potential
The next step in the example script is to run an optimization that minimizes the forces on atoms in the system. This is also known as geometric optimization because the geometry of system changes such that forces on atoms are below a certain threshold which is user defined.
# Run optimization (will change the atomic positions as it is relaxation calculation)
dyn = BFGS(image,logfile='opt_pbe.log', trajectory='h2o_pbe_optimization.traj')
dyn.run(fmax=0.005)
image.calc.close()
For optimization, we can use the optimizers available within ASE’s implementation. The force threshold in the example is 0.005 eV/Angstrom. The file “opt_pbe.log” will save the steps of optimization. You can find out the force and total ground state energy in this file. The optimization trajectory is saved in the “h2o_pbe_optimization.traj” file. You can visualize it to observe how the atoms move to obtain their relaxed atomic positions.
Next few lines in the example script show how to calculate energy, forces and stresses using the BEEF-vdw exchange correlation functional which is known to account for Van der waals interactions in chemical systems.
Strategy on selecting pseudopotentials for DFT calculations:
For running DFT calculations, we will need to select pseudopotentials. We have a large pool of pseudopotentials in our shared folder on PACE-HIVE and PACE-PHOENIX. We also have a module for selecting the pseudopotentials which can be imported using the following command:
from getpseudo import getpseudo
By default, we use the standard solid-state pseudopotential (SSSP) for the calculations. In future, we will default to a set of pseudopotentials which is compatible with both SPARC and Quantum Espresso. To check the available set of pseudopotentials, source the environment QE or SPARC environment (we will learn about SPARC in the future lessons). The commands for sourcing the environments are the following:
For QE:
source /storage/coda1/p-amedford6/0/shared/rich_project_chbe-medford/medford-share/envs/espresso-6.7MaX-beef-ipi
For SPARC:
source /storage/coda1/p-amedford6/0/shared/rich_project_chbe-medford/medford-share/envs/SPARC_latest_2021
Note that these environments are not actively maintained and it is advisable to check the documentation and use the existing files as a guide to build your own. After sourcing the environment, you can use the command “setpseudos” on terminal which will give a list of available pseudopotentials in the following format:
esp#pseudos
test_mostafa#210128#v1
test_mostafa#psp#GGA#ultrasoft#210217#v0
test_oncv#psp#pbe
test_mostafa#210128#v1#UPF
test_mostafa#210224#oncv#psps_soft
test_mostafa#210224#oncv#psps_standard
espresso_sg15_pp#brown#brenda
espresso_sg15_sparc#qimen_SG15#pseudopotentials#psp8#upf_LDA#upf#cut
espresso_sg15_sparc#qimen_SG15#pseudopotentials#psp8#upf_PBE#upf#cut
espresso_sg15_sparc#qimen_SG15#pseudopotentials#psp8#upf_PBE#psp#cut
espresso_sg15_sparc#qimen_SG15#pseudopotentials#psp8#upf_LDA#psp#cut
espresso_sg15_rel
For setting the path to pseudopotential directory in a DFT calculation, you can import the getpseudo module and set the path as an input argument to the calculator. For example, in SPARC, you can set the calculator as given in the following example:
calc = SPARC(atoms = sim, pseudo_dir=getpseudo("espresso_sg15_sparc#qimen_SG15#pseudopotentials#psp8#upf_PBE#psp#cut"), **parameters)
where, pseudo_dir is used to set the path to pseudopotential. “parameters” refers to all other input arguments to the SPARC calculator.
Similar to SPARC, for setting pseudopotential path in QE, you can refer to the following example:
calc = Espresso(atoms=sim,
pw = 500,
psppath = getpseudo("espresso_sg15_mcsh#related#psps_mcsh#related#psps_QE"),
convergence = {"energy": 0.0005, \
"mixing": 0.1, \
"maxsteps": 750, \
"diag": "david"},
smearing = "gauss",
sigma = 0.05,
nbands = -10,
xc = 'PBE',
kpts = (4,4,1),
spinpol = False,
dipole = {"status" : True},
nosym=True)
5.2.2. Exercises:#
Run a simple DFT calculation using QE on PACE-PHOENIX by following the steps given below:
Use ASE to generate a CO_\({2}\) molecule.
Run a single point calculation by sourcing the QE environment using a SLURM script and report the energy. Use PBE as the exchange correlation functional.
Run an optimization calculation for CO_\(2\) molecule by refering to the example given in this lecture for H_\({2}\)O molecule and report the energy of the molecule, forces and stresses on each atom.