Hello Daojian,
no the dipole corrections have not been coded.
What you can do is use isolated boundary conditions (wf goes to 0) with the icoulomb http://www.abinit.org/documentation/hel ... l#icoulomb input flag, and HGH or GTH pseudopotentials. There is also a condition on the orthogonality of the rprim. Of course for a slab this is not what you want. The generalisation of icoulomb to 1d and 2d systems is underway.
Matthieu
On Tue, Mar 2, 2010 at 8:39 AM, Daojian Cheng <chengdaojian@gmail.com> wrote:
> Dear All,
>
>
> Can you tell me that the dipole correction for surface supercell calculations
> is available in Abinit or not? If it is available, how to use it?
>
> Thank you very much for your reply.
>
>
> Daojian
>
dipole correction
Moderator: bguster
dipole correction
Matthieu Verstraete
University of Liege, Belgium
University of Liege, Belgium
Re: dipole correction
no the dipole corrections have not been coded.
Question to the developers: Are there any plans to implement the dipole correction? If yes, is there a time line? Just curious.
Thanks, R.
Re: dipole correction
I know this is an old thread, but I wanted to bring it back up considering the number of people who post about surfaces and 2D materials. Are there any plans to incorporate dipole correction terms for super cells, or would they be particularly easy for someone to implement into the code?
Thanks,
James
Thanks,
James
Re: dipole correction
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
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