[gpaw-users] Dipole Correction in GPAW

[Ask Hjorth Larsen]

Dear Hongliang

On Wed, 14 Apr 2010, Hongliang Xin wrote:

> Dear All,
> 
> I got a simple question about GPAW. I do not know how GPAW makes dipole
> correction in z-direction.
> I found that even for a simple system like atomic oxygen adsorption on
> silver surface, very large vacuum space (25A) is required to converge the
> total energy. For system with core-hole on oxygen, it may require even
> larger space.
> 
> Thanks,
> 
> Hongliang

You can use dipole correction like this:

-------------------------
from ase.data.molecules import molecule
from gpaw import GPAW
from gpaw.poisson import PoissonSolver
from gpaw.dipole_correction import DipoleCorrectionPoissonSolver

system = molecule('H2O')
a = 8.0
system.set_pbc([1, 1, 0])
system.set_cell([a/2, a/2, a])
system.center()


p = PoissonSolver()
calc = GPAW(poissonsolver=DipoleCorrectionPoissonSolver(p, 2), # 2 means z-axis
             h=0.25,
             mode='lcao')

system.set_calculator(calc)
system.get_potential_energy()
-----------------------

I have attached a larger script which also plots some potentials and 
compares to the uncorrected function.

This requires the development version of GPAW and should be considered 
experimental.

Please let us know how the correction works, if anything behaves strangely 
and so on.  Atomic forces may be a bit off for dipole-corrected 
calculations.

Best regards
Ask
-------------- next part --------------
import numpy as np
import pylab as pl
from ase.data.molecules import molecule
from gpaw import GPAW
from gpaw.poisson import PoissonSolver
from gpaw.dipole_correction import DipoleCorrectionPoissonSolver

p1 = PoissonSolver()
p2 = DipoleCorrectionPoissonSolver(PoissonSolver(), 2)

calc1 = GPAW(mode='lcao', basis='sz', poissonsolver=p1, h=0.25)
calc2 = GPAW(mode='lcao', basis='sz', poissonsolver=p2, h=0.25)


system = molecule('HCl')
system.center(vacuum=3.0)
system.center(vacuum=6.0, axis=2)
system.set_pbc([1, 1, 0])

for i, c in enumerate([calc1, calc2]):
    system.set_calculator(c)
    system.get_potential_energy()
    vt_g = c.hamiltonian.vt_sg.sum(0)#c.get_effective_potential()
    #pl.plot(vt_g.sum(0).sum(0), label='%d' % i)
#pl.legend()

pl.figure()
rhot_g = calc1.density.rhot_g
phit1_g = np.zeros_like(rhot_g)
phit2_g = np.zeros_like(rhot_g)

p1.solve(phit1_g, rhot_g)
p2.solve(phit2_g, rhot_g)

dphit_g = phit2_g - phit1_g

pl.title('avg')
celldiag = calc1.density.gd.cell_cv.diagonal()
area = celldiag[0] * celldiag[1]
pl.plot(phit1_g.sum(0).sum(0) / area, label='phit1')
pl.plot(phit2_g.sum(0).sum(0) / area, label='phit2')
pl.plot(dphit_g.sum(0).sum(0) / area, label='corr')

x=10
pl.figure()
pl.title(str(x))
pl.plot(phit1_g[x, x, :], label='phit1')
pl.plot(phit2_g[x, x, :], label='phit2')
pl.plot(dphit_g[x, x, :], label='corr')

x = 3

pl.figure()
pl.title(str(x))
pl.plot(phit1_g[x, x, :], label='phit1')
pl.plot(phit2_g[x, x, :], label='phit2')
pl.plot(dphit_g[x, x, :], label='corr')


x = 0
pl.figure()
pl.title(str(x))
pl.plot(phit1_g[x, x, :], label='phit1')
pl.plot(phit2_g[x, x, :], label='phit2')
pl.plot(dphit_g[x, x, :], label='corr')



pl.legend()

pl.show()