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=0
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.