[gpaw-users] Problems with STM-simulation (copper fcc(111) surface)

[Tiia Tikkanen]

Hello,

I've been trying to do STM-simulation of copper surface (fcc111) with  
GPAW. I'd like to know who's responsible for the STM-code? I have  
managed to do calculations for aluminium surfaces, but similar  
calculations with copper surfaces fail.

Here's one problem:

I get an STM-image, that looks about the right, but the axis of the  
picture seem to be switched. The unit cell is approximately 5.1 Å *  
4.4 Å, but the STM-picture is 4.4 Å * 5.1 Å. Otherwise the picture  
looks about right, it only seems to be streched? See picture:  
https://docs.google.com/drawings/d/10NoQLOOdSe7wulxxXE6cB04Tv1vPbjO45DSCCing3GI/edit?pli=1&hl=en_US

Do you have an idea why the STM-picture turns out weird? I would  
really appreciate your help with this problem.

Here are the scripts I used (set.py to set the surface and tip,  
stmplot.py to do the STM-simulation and plot the image):

# set.py
#
# Sets the surface and the tip
# This is similar to the example code from the GPAW homepages:
#  
https://wiki.fysik.dtu.dk/gpaw/tutorials/negfstm/negfstm.html#setting-up-surface-and-tip

from ase import Atoms, Atom
from gpaw import GPAW, Mixer
from ase.lattice.surface import fcc111
from gpaw.transport.jstm import dump_hs, dump_lead_hs

calc = GPAW(h=0.2,
             mixer=Mixer(0.03, 5, weight=140.0),
             width=0.1,
             mode='lcao',
             basis='szp(dzp)',
             txt='dumphs.txt',
             usesymm=False)

# surface calculation
a = 4.0
srf = fcc111('Cu', size=(2, 2, 9), vacuum=8.0, orthogonal = True)
srf.pbc= (1, 1, 0)

srf.set_calculator(calc)
srf.get_potential_energy()
calc.write('srf', mode='all')
srf.write('cu111_surface.traj')

dump_hs(calc, 'srf', region='surface', cvl=4*3*9)

# tip calculation
a = 0.75 # lattice constant
tip = Atoms('H12', pbc=(1, 1, 0), cell=[5, 5, 12 * a + 7])
tip.positions[:,2] = [i * a for i in range(12)]
tip.positions[:] += (tip.cell / 2.0)[0, :] + (tip.cell / 2.0)[1, :]
tip.translate([0, 0, 6])

tip.set_calculator(calc)
tip.get_potential_energy()
calc.write('tip')
dump_hs(calc,Best regards,

Tiia Tikkanen 'tip', region='tip', cvl=4) # dump overlap and  
hamiltonian matrix

calc.set(kpts=(1, 1, 7)) # for the lead calculations we use
#kpoints in the z-direction

# surface principal layer calculation
srf_p = fcc111('Cu', size=(2, 2, 3), orthogonal = True)
srf_p.pbc = (1, 1, 1)

srf_p.set_calculator(calc)
srf_p.get_potential_energy()
dump_lead_hs(calc, 'srf_p') # dump overlap and hamiltonian matrix

# tip principal layer calculation
tip_p = Atoms('H4', pbc=(1,1,1), cell=[5, 5, 4*a])
tip_p.positions[:,2] = [i * a for i in range(4)]
tip_p.positions[:] += (tip_p.cell / 2.0)[0, :] + (tip_p.cell / 2.0)[1, :]

tip_p.set_calculator(calc)
tip_p.get_potential_energy()
dump_lead_hs(calc, 'tip_p') # dump overlap and hamiltonian matrix
#######

# stmplot.py
#
# Calculation of the STM-image at constant height.

import pickle

from gpaw import GPAW
from gpaw.transport.jstm import STM
from ase.io import write

import matplotlib
matplotlib.use('Agg')
import pylab

basename = 'Cu_dmin5'
savingformat = '.png'

#tip:
tip = GPAW('tip', txt=None)
h1, s1 = pickle.load(open('tip_hs.pckl'))
h10, s10 = pickle.load(open('tip_p_hs.pckl'))

#srf:
srf = GPAW('srf', txt=None)
h2, s2 = pickle.load(open('srf_hs.pckl'))
h20, s20 = pickle.load(open('srf_p_hs.pckl'))

stm = STM(tip, srf,
           hs1=(h1[0], s1[0]),
           hs10=(h10[0], s10[0]),
           hs2=(h2[0], s2[0]),
           hs20=(h20[0], s20[0]),
           bias = 0.05,
           de = 0.05/7.,
           logfile='scan.log')

stm.set(dmin=5)
stm.initialize()
stm.scan()

# plot
kuva = stm.plot(label='I[nA]')

pylab.savefig(basename+'surface'+savingformat)

#######

Best regards,
Tiia Tikkanen