[gpaw-users] How to calculate dipole from density?

[Anders Hellman]

Dear GPAW-users,

I would like to calculate the dipole field at a certain point (point, line, sheet, or volume), where I could control the amount of charge density I want to include. 

For instance, I would like to calculate the dipole field 1Å away from the surface and perhaps only include the firs two surface layers. Now what I think it amount to is to read in the total charge density and thereafter chosen a reference point and than integrate the part of the charge density (wighted with for instance z if it is only the z-part I am interested in). As far as I know this is not doable in the GPAW calculator, by simply write a .get_???.  

However, although it sounds simple I have not been successful. Actually, I have not been able to even get out the dipole of CO (i.e. get the same value as .get_dipole_moment()) by integrating over the whole density. Below there is a crude script that is supposed to do that. Please, if there is anyone that have experience here that could help me I would be very thankful.

best,
Anders

The crude script looks like;

from ase import * 
from gpaw import *

gridrefinement=1
atoms,calc=restart('out.gpw') #Here I read in the gpw-file for CO in a box
density=calc.get_all_electron_density(gridrefinement=gridrefinement)
#density=calc.get_pseudo_density()
gridspace=calc.get_grid_spacings()
ngpts=calc.get_number_of_grid_points()

Z=np.zeros((ngpts[2]))
for i in range(ngpts[2]):
   Z[i]=gridspace[2]*i

print 'Summing up density', np.sum(np.sum(np.sum(density*gridspace[2])*gridspace[1])*gridspace[0])/gridrefinement**3

masscenter=atoms.get_center_of_mass()

for i in range(ngpts[0]):
   for j in range(ngpts[1]):
       for k in range(ngpts[2]):
           density[i][j][k]=density[i][j][k]*(Z[k]-masscenter[2])

print 'Summing up dipole', np.sum(np.sum(np.sum(density*gridspace[2])*gridspace[1])*gridspace[0])/gridrefinement**3