Computermancy in Theoretical Condensed Matter
Explanation of the BC2 ESM external electric field and potential barrier by Yoshinari Takimoto

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

Water Adsorption geometry details calculator

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

Bugfix for SIESTA+ESM 3.1 in poison.f

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 correct use of SIESTA+ESM

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

MacroAve and finding the Work Function

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

The first line is the code you’re using, the second line can be either Potential or Charge, the third line is the name (before the file extension ‘.VH’) of your file to macroscopically average, the fourth line is how many step functions you want to use to smooth it (1 for surfaces, 2 for interfaces) and the next two lines are the width of those two functions. They ought to be equal to something like the planar lattice spacings of your material. The next line is how many electrons in the slab and the last line is the interpolation scheme, either splines or linear fit.
The reason you need to do all this preparation to get a work function in siesta is due to the fact that the energy reference in SIESTA is kind of arbitrary. The zero of the energy is not really the true vacuum level, it’s just set to be the average value of the deformation potential (vh-vna) on the mesh. Performing this calculation for a surface is easier than an interface. For an interface, see this paper: Phys. Rev. B 67, 155327 (2003) and this comment on SIESTA-L: http://www.mail-archive.com/siesta-l@listserv.uam.es/msg03754.htmlTo get a work function of a surface you need to take the difference between the average potential in the vacuum region and the Fermi energy of the slab. Remember that the zero-energy reference level is with respect to the average value of the deformation potential. Plot the deformation potential and then average it, you should then be able to line up everything with respect to the true vacuum level. MARK MY WORDS, unless you are using very thick slabs and ghost atoms above and below the surface the closest you are going to be able to get to an experimental work function is give or take 1eV. It is very sensitive to the surface relaxation and the choice of exchange-correlation functional.
When macroave finishes running you’ll get two files, the .PAV (planar average) and the .MAV (macroscopic average). I loaded the macroscopic averages up into my python console with pylab and plotted them like this (in DreamPie but any python 2.6 or higher console will do as long as you’ve installed pylab):
»> from pylab import *
»> vna = loadtxt(‘/home/abraham/work/macroave_6tl_pbe/vna.MAV’)
»> vh = loadtxt(‘/home/abraham/work/macroave_6tl_pbe/vh.MAV’)
»> plot(vh[:,0],vh[:,1])
»> plot(vh[:,0],vh[:,1]-vna[:,1])
Specifics of TiO2 convergence and Meshcutoff Question

I responded to this inquiry on the SIESTA mailing list:

******************
According to what I have found, if we run the CG or Broyden method to search for atom positions with least forces, say, 0.02 eV/Ang, the default MeshCutoff (e.g., 100.0 Ry) gives a configuration. However, if I simply displace the origin of my lattice from somewhere to some atom, keeping the primitive vectors unchanged, the forces exerted on each atom will change, which is unphysical.
What I can do is to increase the MeshCutoff to try to get converged results, i.e., displacing the origin of the lattice does not change the forces on each atom, and finally in agreement with experiments. But it turns out to be MeshCutoff = 500.0 Ry, so large that it is difficult to extend to larger systems.
My question is that how can I trust the atom configuration with the default MeshCutoff? Is it always necessary to set MeshCutoff to an exceedingly large value to get converged results?
*****************
To which I replied:
Personally I would only trust an 100 Ry cutoff calculation to provide a starting point for a geometric configuration destined for further relaxation. My experience tells me if you use oxygen you should always be using 350 Ry or more, and typically for some materials I’ve been simulating recently, as high as 550Ry. Incidentally, depending on your simulation cell, any cutoff will set the *actual* grid used to be higher, sometimes much higher. If you take a stable geometry and perform a single-point energy (0 CG steps) convergence test with respect to the meshcutoff, try and see what happens when you increment the mesh cutoff by steps of about 100 Ry or so all the way to about 2000 Ry- you should find that the energy won’t change by much after a certain point, and the cutoff you should choose will balance being respectably close to that converged energy. Always record the meshcutoff used by siesta, not the one in your input file.
For example, a while ago with bulk TiO2, I found that the difference in total energy between a ~1400 Ry cutoff and a ~550 Ry cutoff was about the same as the difference between ~450 Ry and ~550 Ry. I chose the 550 Ry grid as it was ‘sufficiently converged’ for my purposes, did an egg-box test, printed the proof in a couple quick plots and saved it somewhere accessible. Having a tight grid that is balanced by its memory usage is a trick you kind of have to get a knack for. It is true that some problems really do require tight expensive grids and a lot of memory, and it might seem to squeeze a lot of the time benefits out of using LCAO-based approach (vdW xc-functionals come to mind). In these cases your parallelization scheme and how you select initial reference geometries will make a lot of difference! The “egg-box” convergence test is another one you will find yourself doing often to check the reliability of your numerical grid, look elsewhere on the list for that.
The bottom line is, I would always explicitly include the MeshCutoff as an input in an fdf of mine for research purposes. Even when I need a ‘quick and dirty’ estimate of something or perhaps a stepping stone to a reference geometry to iron out some of the more expensive force-fluctuation math. 
********************
As an addendum, I included this in an email to the questioner off the SIESTA-L record, and wanted to share it here as well:
Yi Gao
I have, on a terminal here, a multi-step relaxation of TiO2 Anatase with the vdw-KBM xc-functional (a more computationally expensive nonlocal functional) with an 8 x 8 x 8 k-point mesh whose forces are less than 0.002 eV/Ang, and it was done in 3 steps of about 4 hours each on 8 processors. 0.02 eV/Ang might be acceptable in literature, but having very converged simulations of bulk and surface geometries will be useful if you are looking very deeply into the energetics of these materials.
If you’re looking to build a slab or nanostructure from a converged TiO2 run you had damn better make sure your bulk crystal structure has an extremely tight energy convergence. It’s actually a lot easier to do with bulk-phase materials than surface slabs or nanowires.
DO NOT try to relax the whole crystal to an arbitrarily high precision in one go. Do 20-30 cg steps at a time, each time increasing or decreasing several parameters in the fdf. Again, eventually you’ll develop a knack for this as you gain experience with different materials.
I have found the following parameters work pretty well for TiO2:
MaxSCFIterations        100
DM.NumberPulay          5
DM.MixingWeight         0.007   #you can reduce this down to 0.003, 0.001 or lower if you’re having problems converging surfaces
DM.NumberKick          17        # consider reducing this number to about 12 or so once you’ve proceeded very far in establishing convergence
DM.KickMixingWeight    0.0733 
DM.Tolerance            0.0001 #this is the highest DM.tolerance I would trust at a converged geometry. To improve speed when you’re reducing the initial forces from order(10) to order (0.1) this parameter can be set as high as 0.001
DM.EnergyTolerance      0.0001 eV  #default is okay here
MeshCutoff            550 Ry    #if you’re planning on doing surface slab simulations with adsorbates at any point, somewhere around here is where you want to be.
also, pay attention to these parameters too:
MD.TypeOfRun          cg
MD.NumCGsteps         200
MD.MaxCGDispl         0.0001  Ang    #when you are very near convergence, such a low number here is necessary. Far from a convergence basin, this can be as high as 0.2 Ang. In the intermediate regime, use your best judgment
MD.MaxForceTol        0.001 eV/Ang  #achieving this level of force convergence may be unrealistic for some systems., but increment this from 0.08 to 0.04 to 0.02, decreasing MD.MaxCGDispl each time. Sometimes the best you can do is 0.02 eV/Ang.
MD.VariableCell         .true.              #for bulk calculations to minimize the uniaxial stress.
Do not forget to perform the convergence test yourself to convince yourself that this is a good idea! It is more computationally expensive this way than using a 300 Ry cutoff, of course, but it is more reliable!
SLURM and ASE work together!

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)