Explanation of the BC2 ESM external electric field and potential barrier by Yoshinari Takimoto

I made this script to calculate the configurational details of water molecules on a surface. The surface normal is in the +x direction but can be anything (I use TiO2 rutile (110) in my work). I hope you find the script helpful for analyzing water geometries in a facile way. Feel free to remix the work for any arbitrary surface unit cell you’d like, as long as you give me credit in the header of your modified script. Thanks!

#!/usr/bin/env python

#Water adsorption geometry details calculator

#Version 1.1 by Abraham Hmiel July 11 2012

#call like this:

#» python water_configuration.py slab.xyz slab+h2o.xyz

#use the 1st argument as the reference geometry

#use the second argument containing water molecules

#works best when your list of atoms follows: slab, then h2o at the end

#with O coming before the H atoms in the list for each water molecule like so:

#<SLAB>

#X1 Y1 Z1 O

#X2 Y2 Z2 H

#X3 Y3 Z3 H

#X1 Y1 Z1 O

#X2 Y2 Z2 H

#X3 Y3 Z3 H

#… and so on

#Licensed Under Creative Commons attribute remix sharealike license 3.0 Global

from ase import *

from ase.io import read

import sys

import numpy as np

import pylab

ind = 0

#import .xyz files from the command line arguments

for arg in sys.argv:

    ind = ind + 1

    if ind == 2:

        atoms_1 = read(arg)

        print “importing reference geometry .xyz file:”

        print arg

    elif ind == 3:

        atoms_2 = read(arg)

        print “importing water configuration .xyz file:”

        print arg

#prompt the user for the type of adsorbate and the surface unit cell

prompt = ‘> ‘

print “”

print “Pick the type of water adsorbate:”

print “1 - h2o associative adsorption”

print “2 - h2o dissociative adsorption”

print “3 - mixed h2o assoc/dissoc. adsorption”

adsorbate_type = raw_input(prompt)

if not (adsorbate_type == ‘1’ or adsorbate_type == ‘2’ or adsorbate_type == ‘3’):

    print “Please re-run the program with an appropriate choice of task.”

    exit(1)

print “”

print “Pick the surface unit cell:”

print “1 - 2x1”

print “2 - 4x2”

surface_unit_cell = raw_input(prompt)

if not (surface_unit_cell == ‘1’ or surface_unit_cell == ‘2’):

    print “Please re-run the program with an appropriate choice of unit cell.”

    exit(1)

num_Ti_inplane = int(4.0**float(surface_unit_cell))

#create matrices with only position data

ref_geom = atoms_1.get_positions()

h2o_geom = atoms_2.get_positions()

#switch the two arrays if they were given in the wrong order

if pylab.size(ref_geom) > pylab.size(h2o_geom):

    ref_geom = atoms_2.get_positions()

    h2o_geom = atoms_1.get_positions()

#construct a sub-array of only the surface of the h2o_geom array

#the assumption is that the h2o are appended to the end of the list,

#if this is not true then this code will not work!

numatoms_ref = ref_geom.shape[0]

numatoms_h2o = h2o_geom.shape[0]

surf_h2o = h2o_geom[0:numatoms_ref,:]

#calculate the difference in positions between the reference and

#the reconstructed h2o-adsorbed surface

#after - before = difference

diff_surf = surf_h2o - ref_geom

#then, ignore the frozen rows of atoms for the analysis

diff_surf = diff_surf[diff_surf.all(1)]

#this is the number of deleted rows

numfrozenrows = surf_h2o.shape[0] - diff_surf.shape[0]

#find the average position of the Ti planes:

#make a list of all the chemical symbols

#find the x-coordinate of all the Ti atoms

chemsym = atoms_1.get_chemical_symbols()

Ti_ref = np.zeros([numatoms_ref/3, 1])

Ti_h2o = np.zeros([numatoms_ref/3, 1])

numplanes = numatoms_ref/(3*num_Ti_inplane)

index_count=0

index_newarr_ref=0

index_newarr_h2o=0

for j in chemsym:

    if j==’Ti’:

        Ti_ref[index_newarr_ref] = ref_geom[index_count,0]

        Ti_h2o[index_newarr_h2o] = surf_h2o[index_count,0]

        index_newarr_ref = index_newarr_ref+1

        index_newarr_h2o = index_newarr_h2o+1

    index_count=index_count+1

index_count = 0

Ti_av_ref = np.zeros(numplanes)

Ti_av_h2o = np.zeros(numplanes)

Ti_planed_ref = np.zeros(numplanes-1)

Ti_planed_h2o = np.zeros(numplanes-1)

#average the x-coordinate of the atoms in the Ti plane,

#store in Ti_av_ref & Ti_av_h2o

for j in np.arange(numplanes):

    Ti_ref_x = 0.0

    Ti_h2o_x = 0.0

    for k in np.arange(num_Ti_inplane):

        Ti_ref_x = Ti_ref_x + Ti_ref[j*num_Ti_inplane+k]

        Ti_h2o_x = Ti_h2o_x + Ti_h2o[j*num_Ti_inplane+k]

    Ti_av_ref[j] = Ti_ref_x/num_Ti_inplane

    Ti_av_h2o[j] = Ti_h2o_x/num_Ti_inplane

#take differences between successive planes, store in Ti_planed_xxx

for j in np.arange(numplanes-1):

    Ti_planed_ref[j] = abs(Ti_av_ref[j]-Ti_av_ref[j+1])

    Ti_planed_h2o[j] = abs(Ti_av_h2o[j]-Ti_av_h2o[j+1])

#use the distance from the final Ti plane to find the distance

#to the oxygen atom, for each oxygen atom

numh2o = (numatoms_h2o-numatoms_ref)/3

o_tiplane_d = np.zeros(numh2o)

for j in np.arange(numh2o):

    o_tiplane_d[j] = abs(h2o_geom[numatoms_ref+(j*3),0]-Ti_av_h2o[-1])

#Do the following in the purely associative case

if adsorbate_type == ‘1’:

    oh_bond_d = np.zeros(numh2o*2)

    h2o_angle = np.zeros(numh2o)

    h2o_bond_angle = np.zeros(numh2o)

    oh_bond_sym = list()

    oh_bond_index = list()

    for j in np.arange(numh2o):

        #calculate length of the O-H bonds in the water molecules

        #oxygen must occur before hydrogen to work!

        for k in np.arange(2):

            oh_bond_d[j*2+k] = np.sqrt((h2o_geom[numatoms_ref+(j*3),0]- \

                                   h2o_geom[numatoms_ref+(j*3)+k+1,0])**2 + \

                                   (h2o_geom[numatoms_ref+(j*3),1]- \

                                   h2o_geom[numatoms_ref+(j*3)+k+1,1])**2 + \

                                   (h2o_geom[numatoms_ref+(j*3),2]- \

                                   h2o_geom[numatoms_ref+(j*3)+k+1,2])**2)

            oh_bond_index.append(numatoms_ref+(j*3)+1)

            oh_bond_index.append(numatoms_ref+(j*3)+k+2)

            oh_bond_sym.append(atoms_2[numatoms_ref+(j*3)].symbol)

            oh_bond_sym.append(atoms_2[numatoms_ref+(j*3)+k+1].symbol)

        #calculate the angle the h2o molecule makes with the yz plane

        #calculate the unit normal vector perpendicular to the plane

        #then calculate the dot product with [1, 0, 0], divide by norm,

        #take the arccos

        h2o_normal = np.cross(h2o_geom[numatoms_ref+(j*3)+1]- \

                              h2o_geom[numatoms_ref+(j*3)], \

                              h2o_geom[numatoms_ref+(j*3)+2]- \

                              h2o_geom[numatoms_ref+(j*3)])

        h2o_angle[j] = np.arccos(np.dot(h2o_normal,np.array([1.0, 0.0, 0.0]))/ \

                                 np.linalg.norm(h2o_normal))*180./3.1415927

        if h2o_angle[j] > 90.0:

            h2o_angle[j] = 180.0-h2o_angle[j]

        #calculate the H-O-H bond angle by vector subtraction

        #and taking the ‘inverse’ dot product

        vec1 = h2o_geom[numatoms_ref+(j*3)+1] - \

                        h2o_geom[numatoms_ref+(j*3)]

        vec2 = h2o_geom[numatoms_ref+(j*3)+2]- \

                        h2o_geom[numatoms_ref+(j*3)]

        #make sure to adjust the branch of arccos that will give a

        #result of around 105

        h2o_bond_angle[j] = 180. - (np.arccos(abs(np.dot(vec1, vec2))/ \

                            (np.linalg.norm(vec1)*np.linalg.norm(vec2))) * \

                            180./3.1415927)

elif adsorbate_type == ‘2’:

    #calculate each individual OH-bond distance

    oh_bond_d = np.zeros(numh2o*2)

    oh_bond_sym = list()

    oh_bond_index = list()

    #find the nearest-neighbor of each hydrogen atom

    index_count = 0

    for k in np.arange(numatoms_ref,numatoms_h2o):

        if atoms_2[k].symbol == ‘H’:

            #set up fiducial ‘maximum radius’

            h_x_dist = 10000

            n = atoms_2[k].get_position()

            for l in np.arange(numatoms_h2o-(num_Ti_inplane*4),numatoms_h2o):

                m = h2o_geom[l]

                h_x_dist_try = np.sqrt((m[0]-n[0])**2 + (m[1]-n[1])**2 + \

                                   (m[2]-n[2])**2)

                if h_x_dist_try < h_x_dist and h_x_dist_try > 0.00001:

                    h_x_dist = h_x_dist_try

                    l_ind = l

            #store indices and symbols as well as distance

            oh_bond_d[index_count] = h_x_dist

            index_count = index_count + 1

            oh_bond_sym.append(atoms_2[k].symbol)

            oh_bond_index.append(k+1)

            oh_bond_sym.append(atoms_2[l_ind].symbol)

            oh_bond_index.append(l_ind+1)

#mixed adsorption case:

#will do later

elif adsorbate_type == 3:

    #identify the molecule that is dissociatively adsorbed

    #look at h2o oxygens and find nearest neighbor hydrogens

    #if one’s distance is more than 1.2 Angstrom, that group

    #is dissociatively adsorbed.

    oh_bond_cutoff = 1.2 #Angstrom

    oh_bond_d = np.zeros(numh2o*2)

    h2o_angle = np.zeros(numh2o)

    h2o_bond_angle = np.zeros(numh2o)

    oh_bond_sym = list()

    oh_bond_index = list()

    for j in np.arange(numh2o):

        dist_check = np.zeros(2)

        for k in np.arange(2):

            dist_check[k] = np.sqrt((h2o_geom[numatoms_ref+(j*3),0]- \

                                   h2o_geom[numatoms_ref+(j*3)+k+1,0])**2 + \

                                   (h2o_geom[numatoms_ref+(j*3),1]- \

                                   h2o_geom[numatoms_ref+(j*3)+k+1,1])**2 + \

                                   (h2o_geom[numatoms_ref+(j*3),2]- \

                                   h2o_geom[numatoms_ref+(j*3)+k+1,2])**2)

        if dist_check[0] < oh_bond_cutoff and \

           dist_check[1] < oh_bond_cutoff:

            #associative adsorption

        #calculate the angle the h2o molecule makes with the yz plane

        #calculate the unit normal vector perpendicular to the plane

        #then calculate the dot product with [1, 0, 0], divide by norm,

        #take the arccos

                h2o_normal = np.cross(h2o_geom[numatoms_ref+(j*3)+1]- \

                              h2o_geom[numatoms_ref+(j*3)], \

                              h2o_geom[numatoms_ref+(j*3)+2]- \

                              h2o_geom[numatoms_ref+(j*3)])

                h2o_angle[j] = np.arccos(np.dot(h2o_normal,np.array( \

                               [1.0, 0.0, 0.0]))/  \

                                 np.linalg.norm(h2o_normal))*180./3.1415927

                if h2o_angle[j] > 90.0:

                    h2o_angle[j] = 180.0-h2o_angle[j]

        #calculate the H-O-H bond angle by vector subtraction

        #and taking the ‘inverse’ dot product

                vec1 = h2o_geom[numatoms_ref+(j*3)+1] - \

                        h2o_geom[numatoms_ref+(J*3)]

                vec2 = h2o_geom[numatoms_ref+(j*3)+2]- \

                        h2o_geom[numatoms_ref+(j*3)]

        #make sure to adjust the branch of arccos that will give a

        #result of around 105

                h2o_bond_angle[j] = 180. - (np.arccos(abs(np.dot(vec1, vec2))/ \

                            (np.linalg.norm(vec1)*np.linalg.norm(vec2))) * \

                            180./3.1415927)

        else:

                #find the nearest-neighbor of each hydrogen atom

                index_count = 0

                for k in np.arange(numatoms_ref,numatoms_h2o):

                    if atoms_2[k].symbol == ‘H’:

                        #set up fiducial ‘maximum radius’

                        h_x_dist = 10000

                        n = atoms_2[k].get_position()

                        for l in np.arange(numatoms_h2o- \

                                           (num_Ti_inplane*4),numatoms_h2o):

                            m = h2o_geom[l]

                            h_x_dist_try = np.sqrt((m[0]-n[0])**2 + \

                                                   (m[1]-n[1])**2 + \

                                                   (m[2]-n[2])**2)

                            if h_x_dist_try < h_x_dist and \

                                   h_x_dist_try > 0.00001:

                                h_x_dist = h_x_dist_try

                                l_ind = l

                        #store indices and symbols as well as distance

                        oh_bond_d[index_count] = h_x_dist

                        index_count = index_count + 1

                        oh_bond_sym.append(atoms_2[k].symbol)

                        oh_bond_index.append(k+1)

                        oh_bond_sym.append(atoms_2[l_ind].symbol)

                        oh_bond_index.append(l_ind+1)

#calculate the distance between the oxygen atom

#in the h2o and the nearest Ti5d atom on the surface

index_count = 0

oti_bond_d = np.zeros(numh2o)

oti_bond_sym = list()

oti_bond_index = list()

for k in np.arange(numatoms_ref,numatoms_h2o):

    if atoms_2[k].symbol == ‘O’:

        #set up fiducial ‘maximum radius’

        o_ti_dist = 10000

        n = atoms_2[k].get_position()

        for l in np.arange(numatoms_h2o-(num_Ti_inplane*4),numatoms_h2o):

            #print atoms_2[l]

            if atoms_2[l].symbol == ‘Ti’:

                m = surf_h2o[l]

                o_ti_dist_try = np.sqrt((m[0]-n[0])**2 + (m[1]-n[1])**2 + \

                                   (m[2]-n[2])**2)

                if o_ti_dist_try < o_ti_dist:

                    o_ti_dist = o_ti_dist_try

                    l_ind = l

        #store indices and symbols as well as distance

        oti_bond_d[index_count] = o_ti_dist

        index_count = index_count + 1

        oti_bond_sym.append(atoms_2[k].symbol)

        oti_bond_index.append(k+1)

        oti_bond_sym.append(atoms_2[l_ind].symbol)

        oti_bond_index.append(l_ind+1)

#print out the data

print “RESULTS:”

print “Relaxation distances between adsorbate case and “

print “clean surface, not counting frozen atoms (Angstrom):”

print ‘\n’.join(’ ‘.join(str(cell) for cell in row) for row in diff_surf)

print “”

print “Distance between averaged Ti planes in clean surface (Angstrom):”

print Ti_planed_ref

print “”

print “Distance between averaged Ti planes in hydrated surface (Angstrom):”

print Ti_planed_h2o

print “”

print “Distance between topmost Ti plane and Oxygen atom(s) in H2O (Angstrom):”

print o_tiplane_d

if adsorbate_type == ‘1’:

    print “”

    print “O-H bond lengths in H2O (Angstrom):”

    print oh_bond_d

    print “”

    print “O-H bond symbols:”

    print oh_bond_sym

    print “”

    print “O-H bond indices:”

    print oh_bond_index

    print “”

    print “Angle between H2O plane and surface normal (Degrees):”

    print h2o_angle

    print “H—O—H bond angle (Degrees):”

    print h2o_bond_angle

elif adsorbate_type == ‘2’:

    print “”

    print “O-H nearest-neighbor bond lengths in dissociated H2O (Angstrom):”

    print oh_bond_d

    print “”

    print “O-H bond symbols:”

    print oh_bond_sym

    print “”

    print “O-H bond indices:”

    print oh_bond_index

    print “”

elif adsorbate_type == 3:

    print “”

    print “O-H nearest-neighbor bond lengths (Angstrom):”

    print oh_bond_d

    print “”

    print “O-H bond symbols:”

    print oh_bond_sym

    print “”

    print “O-H bond indices:”

    print oh_bond_index

    print “”

    print “Angle between H2O plane and surface normal (Degrees):”

    print h2o_angle

    print “H—O—H bond angle (Degrees):”

    print h2o_bond_angle

print “”

print “O-Ti nearest-neighbor bond lengths from H2O to surface (Angstrom):”

print oti_bond_d

print “”

print “O-Ti bond symbols:”

print oti_bond_sym

print “”

print “O-Ti bond indices:”

print oti_bond_index

I recently had a correspondence with Dr. Yoshi Takimoto of the Sugino Group about issues with getting bad adsorption energies with their code, which ought to be a relatively simple matter. As it turns out, there is an integer casing error in the F90 code which modifies the solution of Poisson’s equation for the metal-metal boundary condition (bc2). To fix this, the lines in the poison.f (after applying the ESM patch) must be changed from:

              else if (ESM_bc.eq.2) then

                rg =

     .            -4._dp*PI*

     .            (exp(+sgp*(xx-ESM_x1))-exp(-sgp*(xx+3._dp*ESM_x1)))

     .            *tmp(1)/(1-exp(-4._dp*sgp*ESM_x1))

     .            +4._dp*PI*

     .            (exp(+sgp*(xx-3._dp*ESM_x1))-exp(-sgp*(xx+ESM_x1)))

     .            *tmp(2)/(1-exp(-4._dp*sgp*ESM_x1))

to:

              else if (ESM_bc.eq.2) then

                rg =

     .            -4._dp*PI*

     .            (exp(+sgp*(xx-ESM_x1))-exp(-sgp*(xx+3._dp*ESM_x1)))

     .            *tmp(1)/(1._dp-exp(-4._dp*sgp*ESM_x1))

     .            +4._dp*PI*

     .            (exp(+sgp*(xx-3._dp*ESM_x1))-exp(-sgp*(xx+ESM_x1)))

     .            *tmp(2)/(1._dp-exp(-4._dp*sgp*ESM_x1))

Note the absence of the “._dp” after the “1” in the denominator, which occurs twice. I am currently testing this bugfix. However, I feel that it might not solve the issue entirely, as the absurd adsorption energies are prevalent in the bc1 case (vacuum-vacuum) as well. Furthermore, what’s still puzzling is that the only energy that seems “off” is the slab calculation, not the water molecule or the slab+adsorbate. I feel there is still much work to be done here, as the periodic boundary condition case yields results that I expect.

The ESM package (available here: http://sugino.issp.u-tokyo.ac.jp/esm/index.php?Programs#SIESTA )  for calculation of real interfaces under charging and electrical bias, has existed for a little bit of time now. This method enables one to describe the surface and interface, e.g. a slab system that is periodic in the direction parallel to the surface but having arbitrary shape in the perpendicular direction, sandwiched by semi-infinite media, such as vacuum, an electrode, or an electrolyte. Here, we call the medium “effective screening medium” (ESM).

The problem is (sorry, developers) that the ESM documentation can be misleading at times. If you don’t know what you are doing, you will obtain very bad adsorption energies if you follow through with the typical procedure. In fact, some of the differences in (slab+adsorbate)/(slab)(adsorbate) energies for the neutral atom energy, DUscf, DEna, Exc, and other energies will be horribly off from the periodic boundary condition case. I even got energies from by BSSE correction single point energies that were lower than my relaxed slabs! Talk about unphysical… my initial simulations of a water molecule on a Rutile TiO2 (110) 4x2 surface had adsorption energies of -2.1 eV, which is ridiculous. To get around this and recover output that actually makes sense, here’s what I did- you should follow this exactly, even if it doesn’t make sense to begin:

1. Make sure the geometry of your surface slab is centered about the line x=0 and make sure that the x direction is the surface normal vector. It’s important to leave enough vacuum between the slab, preferably you should be using a vacuum spacing of at least 13 Angstrom, take care that any adsorbate molecules are not close to the ESM boundary at the center of the unit cell.

2. Make sure the “ESM.ZeroAt” line is set to ‘center’ - for some reason I have gotten spurious energies when setting the ESM zero to the left side of the simulation cell. 

3. Do not use a spacer or barrier for a CG relaxation. 

4. It’s helpful to write the potential along X. Use the flag “ESM.WriteRhoPotX .true.” - this may be handy in that if all you’re interested in is the planar average of the electron density or electrostatic potential, you won’t have to make up a macroave input file and then run macroave.

5. When adding and subtracting charge or applying an electric field, use the boundary condition ‘bc2’ and try to stay in the linear response regime or else you will have convergence problems. Start with a small perturbation and slowly increase it, taking note of any changes in the convergence of the SCF cycle rather than starting from a specified field and charge and hoping to get lucky. Re-use the density matrix from a previous run of smaller field or charge to stay on track.

6. Do not use the SlabDipoleCorrection .true. flag at all if you are using bc1/bc2/bc3. This is only for periodic boundary conditions to induce a field in the vacuum region that in turn creates a discontinuity in the electrostatic potential so that the two sides of the slab can stay electronically distinct. It is not necessary to use this in ESM because of the broken symmetry in the x direction.

That’s all with the helpful hints for ESM + SIESTA-3.1 this time. I’ll probably have more, eventually

Finding the work function of a specific surface or interface is a relatively straightforward task in ab-initio quantum chemistry. The program MACROAVE, bundled with SIESTA, can help you do this and it works for abinit as well. 

Your first step is to prepare a converged (Forces < 0.02 eV/Ang and lower, if possible) slab calculation. Now, with the LCAO method you are never going to get a “right answer” because of the limitations of the basis set. This is because the electrostatic potential “fall-off” region between the slab and the vacuum is not going to be adequately represented with an atomic basis set unless you use a few tricks. The first thing to try is diffuse functions on the surface atoms. This overall is not going to help much, but you can stretch its efficacy much farther if you also include a row of frozen ghost atoms above the surface on both sides to complete the basis in the interstitial region.

Macroave is compiled in the /Util directory of your siesta build. You should » make macroave and then » make permute - permute can shift a surface slab cell’s potential or charge density file on the grid from the slab normal direction of z to x. This is necessary if you’re using SIESTA+ESM since you don’t get a choice of slab normal direction with ESM— it is always in the +x direction. Macroave will by default assume that the z direction is the slab normal direction. 

The next step is to turn on the parameters to write the potentials. You will need to set #saveelectrostaticpotential and #savevna to true in your input file, then if need be, use permute to shape the densities into the right normal direction. Then, shorten the names of the potential files to less than 10 characters or so. Macroave will always be looking for something with a .VH file extension so if you have a neutral atom potential file, name it vna.VH or something. Your input file for macroave should always be called macroave.in and it looks like this.

Siesta

Potential

vt

2

5.165

6.235

2304.0

spline

I responded to this inquiry on the SIESTA mailing list:

Hello all, for my first post back from my hiatus, I’ll show you how to use ASE in conjunction with a queuing script:

Let’s say we have an ASE python file like so, where the object is to calculate the total energy of a water dimer with the oxygens facing each other, having the O-O distance changing from far to close:

Read More

undefinedauguments:

Ultrathin Electronic Escarpment (Part III)

undefinedauguments:

Ultrathin Electronic Escarpment (Part II)