Ionic contribution to the macroscopic polarization
Posted: Sat May 22, 2021 10:17 am
I might've spotted a bug in subroutine `POINT_CHARGE_DIPOL` in `dipol.F` which is only used in `pead.F` to calculate the ionic contribution to the macroscopic polarization.
In subroutine `POINT_CHARGE_DIPOL`, the following code:
essentially average out the position between ions with their periodic image (i.e. + a cell length) if they are sufficiently close (controlled by `TINY`) to the periodic boundary. (Note that the position of the periodic boundary is now defined relative to `POSCEN` which can be changed using `DIPOL` tag. This means one has to use the same `DIPOL` for both centro and non-centrosymmetric calculations.)
The logic here is that if the atom is at the periodic boundary, we should also count its equivalent position and add it to the ionic contribution. However, this is not necessarily true for the polarization calculations since all we need is the difference of polarization between a centrosymmetric phase and a non-centrosymmetric one. If the atom is considered to be "at the periodic boundary" for the centrosymmetric phase and "not at the periodic boundary" for the non-centrosymmetric phase, this procedure will cause problems.
I think these few lines of code have caused a lot of confusion in the community especially because the result of this routine doesn't match a simple manual calculation using $(P_{ion}-P_{DIPOL}) * ZVAL_{ion}$ where $P$ stand for "Position".
I have tried to comment these lines out and the result matches the expected values perfectly.
In subroutine `POINT_CHARGE_DIPOL`, the following code:
Code: Select all
!jF: the atoms at the cell boundary (position -0.5) cause a problem since
! one also had to add a contribution from the equivalent position at the
! opposite cell boundary (position +0.5) effectively adding up to zero
! (due to "position + -position") what should be fixed by following code.
! Such problems only occur for bulk cells. For isolated molecules etc.
! there are no atoms at the cell boundary, just vacuum :-). It is also no
! problem for the quadrupole moment (due to "position**2"). If you still
! fail to get the correct dipole moment play around with parameter TINY
! which defines the "thickness of the cell boundary region" ...
IF (ABS(ABS(X)-0.5_q)<TINY/LATT_CUR%ANORM(1)) X=0
IF (ABS(ABS(Y)-0.5_q)<TINY/LATT_CUR%ANORM(2)) Y=0
IF (ABS(ABS(Z)-0.5_q)<TINY/LATT_CUR%ANORM(3)) Z=0The logic here is that if the atom is at the periodic boundary, we should also count its equivalent position and add it to the ionic contribution. However, this is not necessarily true for the polarization calculations since all we need is the difference of polarization between a centrosymmetric phase and a non-centrosymmetric one. If the atom is considered to be "at the periodic boundary" for the centrosymmetric phase and "not at the periodic boundary" for the non-centrosymmetric phase, this procedure will cause problems.
I think these few lines of code have caused a lot of confusion in the community especially because the result of this routine doesn't match a simple manual calculation using $(P_{ion}-P_{DIPOL}) * ZVAL_{ion}$ where $P$ stand for "Position".
I have tried to comment these lines out and the result matches the expected values perfectly.
