XANES with GPAW

Last modified by sjhuotar@helsinki_fi on 2024/02/07 06:22

XANES with GPAW

GPAW is a density-functional theory (DFT) Python code based on the projector-augmented wave (PAW) method. GPAW is available free (GNU Public License), can be compiled on your own machine, and is available at CSC supercomputers (see below). Using the projector-augmented wave (PAW) method allows one to get rid of the core electrons and work with soft pseudo valence wave functions, i.e. the PAW method is an all-electron method in the frozen core approximation, and there is a one to one transformation between the pseudo and all-electron quantities.

Pseudo wave functions, pseudo electron densities and potentials are represented on uniform real-space orthorhombic grids. In each of the three directions, the boundary conditions can be either periodic or open.

All the exchange-correlation functionals from the libxc library can be used.

Parallelization is done by domain decomposition of grids, as well as distributing k-points, spins, and bands over all processors. This leads to very good scalability in supercomputers.

TPA XANES

GPAW can calculate XANES intensities using the transition potential approximation (TPA) i.e. half-core hole (HCH) approach, and also full-core hole (FCH) approach. The calculation can be done either traditionally using empty orbitals (which must be present in the calculation), or using a recursion method which reaches very high in energy without needing to include many empty states in the calculation.

In order to get any meaningful energy scale, Delta-Kohn-Sham corrections need to be calculated for each TPA calculation.

Molecules

Normal vs. recursion.

Periodic boundary conditions

k-point sampling. Complications / code instabilities.

Delta Kohn Sham

Example code

import ase
import gpaw

# Choosing .xyz file (first it finds in current directory)
import os, fnmatch
xyzfiles = fnmatch.filter(os.listdir('.'),'*.xyz')
basename = xyzfiles[0]
basename = basename[0:-4] # this trims ".xyz"

ase.parallel.parprint('Running MPI calculations for ' + basename + '.xyz ...')

# Set calculation parameters.
h = 0.2     # grid spacing, using default value, good convergence (in Angstrom)
xc = 'PBE'  # exchange-correlation functional (default LDA)
vacuum = 5  # how much vacuum to add around cluster (in Angstrom)

# if necessary (CSC runs): create core-hole potentials in current directory
if 1:
  from gpaw.test import gen
   gpaw.setup_paths.insert(0, '.')
   gen('O', name='hch1s', xcname=xc, corehole=(1, 0, 0.5))
   gen('O', name='fch1s', xcname=xc, corehole=(1, 0, 1.0))

# Read coordinates. NB. 1st atom in list will be core-excited
atom = ase.io.read(basename+'.xyz')
atom.center(vacuum=vacuum)

# Calculate total charge (for LiCl solutions only...)
atomic_numbers = atom.get_atomic_numbers()
n_Li = sum(atomic_numbers == 3)
n_Cl = sum(atomic_numbers == 17)
total_charge = n_Li - n_Cl

calc = gpaw.GPAW(h=h, txt = basename + '_hch.log', charge = total_charge, xc=xc, setups = {0:'hch1s'})
atom.set_calculator(calc)

# Do the HCH calculation
ase.parallel.parprint(' HCH calculation...')
atom.get_potential_energy()
homo_lumo = calc.get_homo_lumo()
lumo = homo_lumo[1] # energy of LUMO state

# Save the checkpoint file?
#calc.write(basename+'.gpw')

# Recursion method XAS calculation
ase.parallel.parprint(' XAS (recursion) calculation...')

from gpaw.xas import *
r = RecursionMethod(calc)
r.run(600)
r.run(1400,inverse_overlap='approximate')
r.write(basename + '.rec',mode='all')

# Delta-Kohn-Sham calculation
ase.parallel.parprint(' GS calculation...')
calc = gpaw.GPAW(h=h,txt = basename + '_gs.log', charge = total_charge, xc=xc)
atom.set_calculator(calc)
e1 = atom.get_potential_energy() + calc.get_reference_energy()

ase.parallel.parprint(' FCH calculation...')
calc = gpaw.GPAW(h=h,txt=basename + '_fch.log',xc=xc,charge= total_charge-1, spinpol = True, occupations = gpaw.FermiDirac(0.0, fixmagmom=True), setups={0: 'fch1s'})
atom[0].magmom = 1
atom.set_calculator(calc)
e2 = atom.get_potential_energy() + calc.get_reference_energy()

dks_energy = e2 - e1

# Write Delta-Kohn-Sham value
f = ase.parallel.paropen(basename + '_deltaks.dat', 'w')
f.write(str(dks_energy))
f.close

# Write the actual correction value to be applied to spectrum
f = ase.parallel.paropen(basename + '_shift.dat', 'w')
f.write(str(dks_energy - lumo))
f.close

ase.parallel.parprint(' ... MPI part done!') 

Recursion calculation spectrum output (single cpu)

from gpaw.xas import *
import numpy, sys

# This part of the calc. cannot be parallelized (but it only takes 1 min with 1 core)
import gpaw.mpi as mpi
if mpi.size ==1:

  # Broadening parameters (in eV)
   delta = 0.01 # Lorentzian broadening FWHM (mandatory)
   fwhm =  None # Gaussian broadening FWHM
   dx = 0.001   # energy scale step size

  # Filename (first .xyz file found)
  import os, fnmatch
   recfiles = fnmatch.filter(os.listdir('.'),'*.xyz')
   basename = recfiles[0]
   basename = basename[0:-4] # removes '.xyz'

   sys.setrecursionlimit(10000)

   x_start=-20
   x_end=100
   x_rec = x_start + numpy.arange(0, x_end - x_start ,dx)

  print " Producing spectrum..."

   r = RecursionMethod(filename=basename + '.rec')
   y = r.get_spectra(x_rec, delta=delta, fwhm=fwhm)
   y2 = sum(y)

  # Shift the spectrum according to Delta-Kohn-Sham and LUMO energy
   f = open(basename + '_shift.dat')
   shift = float(f.read())
   f.close()
   x_rec = x_rec + shift

   x2 = numpy.vstack((x_rec,y2)); x2 = x2.conj().transpose();
   numpy.savetxt(basename + '_xas.dat',x2,fmt='%.6f',delimiter=' ')

  print " ...done!"

else:
  print "This part can only be run in serial."

Running GPAW at CSC

N.B. Switch to Vuori June 2012 done below.

Programs can't be easily run interactively on the CSC clusters. The jobs are submitted to "queues" where they are run when there is capacity.

The runfiles below are for the codes shown above (gpaw.py and xas.py). These are submitted by running sbatch runfile. Job status can be queried with the command squeue.

Vuori resource monitoring
Vuori running programs
Running gpaw

Naturally, you have to wait for the first part to finish before running the second part.

Jobs cannot be run in the home directory (which has very little space). Instead, copy the python files and runfiles to $WRKDIR (or a subfolder of it), and submit the work from that directory. Copy the results to your own computer using scp.

#!/bin/csh
## 32 cores with 1 GB each, max 1 hour runtime
#SBATCH -n 32
#SBATCH -t 01:00:00
#SBATCH --mem-per-cpu=1000
#SBATCH -J run_mpi_%J
#SBATCH -e run_mpi_err_%J
#SBATCH -o run_mpi_out_%J
#SBATCH -p parallel

module load gpaw
srun gpaw-python runfile.py 

Vuori single-cpu runfile

#!/bin/csh
## 32 cores with 1 GB each, max 1 hour runtime
#SBATCH -n 1
#SBATCH -t 01:00:00
#SBATCH --mem-per-cpu=8096
#SBATCH -J run_mpi_%J
#SBATCH -e run_mpi_err_%J
#SBATCH -o run_mpi_out_%J
#SBATCH -p serial

module load gpaw
srun gpaw-python xas.py

Relevant output files of the MPI part will be (assuming input file filename.xyz):

  • run_mpi_out_XXXX: this gives the stdout of the run (with a numeric code XXXX) and includes a resource usage summary (CPU time and memory usage)
  • run_mpi_err_XXXX: this gives the stderr of the run, and is useful to check for input errors or insufficient runtime
  • filename_hch.log: GPAW output for TPA calculation
  • filename_fch.log: GPAW output for full-core hole calculation (necessary for Delta Kohn Sham)
  • filename_gs.log: GPAW output for ground state calculation (necessary for Delta Kohn Sham)
  • filename_deltaks.dat: the Delta Kohn Sham value calculated
  • filename_shift.dat: the actual shift applied to the XANES spectrum (Delta Kohn Sham value - LUMO energy)
  • filename.gpw: a restart file for the TPA part (not saved by default, takes loads of space)
  • filename.rec: the recursion coefficients, necessary for the single-cpu part (takes space)

Relevant output files of the single-cpu part will be:

  • filename_xas.dat: this gives the broadened XANES spectrum

CSC computers are compiled with the stable version, GPAW 0.7.2.7150. It is possible to use the newest version (which may not be the actual newest version, contact Jussi Elovaara if necessary) by changing module loading command in the runfile to module load gpaw/svn.