Reading through the threads, I had a couple of questions related to this topic. If someone, either developer or no has any thought's I'd love to hear them. My understanding is that VASP & Quantum Espresso include their dipole corrections according to the following two papers J. Neugebauer and M. Scheffler, Phys. Rev. B 46, 16 067 ~1992, and Bengtsson APL vol 59 12301 (1999)
http://journals.aps.org/prb/pdf/10.1103 ... B.59.12301 which corrects the energies by a factor of 1/2 due to the dipole field having an internal vs external origin.
In both of these papers, the authors point out the periodic boundary conditions with a slab that has a dipole can lead to erroneous results, due to the presence of an artificial dipole caused by forcing periodicity in the potential. IE, If the slab is oriented in the XY plane, then far from the slab there should be a potential offset between the vacuum level below the slab and above the slab of ±2*pi*m, where m is the net dipole of the slab. m is defined as the integral from -inf to inf over rho_av(z')z'dz', where rho_av(z) is the average charge density (sum of electron and ionic charge density). Without periodic boundary conditions, the electrostatic potential should be flat far from the surface, but due to the PBC's there is an artificial linear term forcing V(z=0) = V(z=L) where L is the length of the box. To correct for this, a linear dipole correction term can be subtracted from the potential that looks like V_dipole_correction (z) = 4*pi*m(z/zmax - 1/2). The paper subsequently defines corrections to the energy and forces that need to be added to correct for this linear term. The papers also differ on the magnitude of the dipole correction.
The authors then present some calculations for a water molecule in a 3 x 3 x 6 Angstrom box, which are easily reproduced within Abinit, and the expected non-flat potential from the PBCs can clearly be observed. I've plotted the xy averaged classical electrostatic potentials ( which as I understand should be the sum of the hartree potential and the psuedopotentials in real space). To calculate the dipole moment, m, I integrated the z-averaged electron density * zdz from 0 to L, and used the Z-averaged psuedocharge positions with the ionic charges of 6, 1, 1 for Oxygen, hydrogen and hydrogen respectively. This resulted in a dipole m of ~0.6 electron bohr/Area_unit_cell = 0.018 e/bohr. Using that as the surface dipole and plugging into the dipole correction, you get a dipole correction of V_Dipole_correction = -4*pi*.018*(Z/max(Z) - .5).
As a quick check, I changed outscf.F90 to print out VPSP and VHA to plot V_Electrostatic. Using cut 3d, I made an xy-averaged electrostatic potential, and fit the points in the vacuum to a line. The slope of that line was the same slope as the dipole correction term to within 1 part / 10^6.
I then altered rhotov.F90 to introduce this correction into the calculations self-consistently. Rather than calculate m each time, I made a dipole correction term by fitting the first ~10% of the electrostatic potential (after centering my atoms in middle of the slab). I still need to correct for the energies & forces, but does anyone see anything terrible with what I'm doing?
Thanks,
James
Code: Select all
!! Compute Dipole Correction Along the Z-axis Start by assuming that the slab is centered at zero and nspden=1
!! and no parallelizaton over ffts.
!! Start with the electrostatic potential (local VHartree + VPSP)
ngfft1dp=dble(ngfft(1))
ngfft2dp=dble(ngfft(2))
ngfft3dp=dble(ngfft(3))
do ifft=1,nfft
ves(ifft) = vhartr(ifft) + vpsp(ifft)
end do
!! Calculate the Z-averaged Electrostaic potential
do iz=1,ngfft(3)
vtmp=0
do ix=1,ngfft(1)
do iy=1,ngfft(2)
vtmp=vtmp + (ves(iy + (ix-1)*ngfft(1) +(iz-1)*ngfft(1)*ngfft(2)))/(ngfft1dp*ngfft2dp)
end do
end do
vesz(iz)=vtmp
!write (*,*) 'Electrostatic potential at Z=ZZ is ', iz ,vesz(iz)
end do
!! Potential should be flat far from the Slab. Assume Center of Density is halfway in the unit cell
!! Artificial periodic dipole should produce a linear potential far from the slab
!! Fit the z-averaged potential to a line, Use the dipole correction of
!! V_DIP = 4*pi*m*(z/z_box -1/2) after Bengtsson PRB Vol 59 p.12301 (1999)
!! TODO Don't assume center of charge density is halfway in unit cell
sumx=0.0_dp
sumx2=0.0_dp
sumxy=0.0_dp
sumy=0.0_dp
sumy2=0.0_dp
do iz=1,ngfft(3)/10
sumy=sumy +vesz(iz)
sumy2 = sumy2 + vesz(iz) * vesz(iz)
sumxy = sumxy + vesz(iz) * rprimd(3,3)*iz*(1.0/ngfft3dp)
sumx = sumx + 1.0*iz*(rprimd(3,3)/ngfft3dp)
sumx2 = sumx2 + 1.0*iz*iz*rprimd(3,3)*rprimd(3,3)/(ngfft3dp**2)
end do
n=ngfft(3)/10
dz=rprimd(3,3)/ngfft3dp
m=((1.0*n) * sumxy - sumx*sumy)/((1.0*n)*sumx2 - sumx**2)
write (*,*) 'Calculated slope of Vesz =', m
write (*,*) ' n = ', n
write (*,*) 'sums', sumx, sumx2, sumxy, sumy, sumy2, ngfft3dp
m= m/(4.0*3.14159265359)
do iz=1,ngfft(3)
Vdz(iz)=-4.0*3.14159265359 *m* ((1.0*iz)*dz - 0.5*rprimd(3,3))
write (*,*) 'Vdz at ZZ = ', dz*1.0*(iz-1),Vdz(iz)
end do
do iz=1, ngfft(3)
do ix = 1, ngfft(1)
do iy= 1, ngfft(2)
vdip(iy + (ix-1)*ngfft(1) + (iz-1)*ngfft(1)*ngfft(2)) = Vdz(iz)
end do
end do
end do
write (*,*) 'istep =', istep
!do iz=1, nfft
! write (*,*) 'Vdip at V(x,y,z)', vdip(iz)
!end do
if (optres==0) then
! Compute potential residual
if (.not. wvlbigdft) then
!$OMP PARALLEL DO COLLAPSE(2)
do ispden=1,min(dtset%nspden,2)
do ifft=1,nfft
vnew(ifft,ispden)=vhartr(ifft)+vpsp(ifft)+vdip(ifft)+vxc(ifft,ispden)+vzeeman(ispden)+Vmagconstr(ifft,ispden)
vresidnew(ifft,ispden)=vnew(ifft,ispden)-vtrial(ifft,ispden)
end do
end do