reference position bug during dipole moment calculation?

[Phil G.]

Dear all,

in order to get a bulk polarization of a ferroelectric slab system (2x2x1 
slab with 3D periodic boundary condition) I let CP2K start to calculate the 
dipole moments of the slab system (in fact the unit cell has the hexagonal 
shape, but it does not matter if I cut off a piece of the large system of 
the duplicates of the unit cell in xyz coordinate system and this cut piece 
satisfy the 3D pbc in the xyz coordinate system?). In order to check the 
influence of the choice of a reference point on the dipole moment 
calculation, I let start several calculations with different reference 
points, e.g. COAC, ZERO, user_defined -> reference_point x y z.

If using the reference_point ZERO, then I notice an non-explainable 
difference between the quantity of the dipole moments for up-polarized slab 
and down-polarized slab (the slabs are the same, but are only 180 degree 
rotated to each other) and it seems to be the false choice for the dipole 
moment calculation.
If I use COAC as the reference point, then the both dipole moment values 
are almost of the same quantity (absolute value) in z direction, but of 
opposite sign (which is correct). If using COM as the reference point, then 
I am wondering why not really the center of the slab was used as expected 
(consider the slab system with lengths a_x, a_y, and c_z along the three 
x,y,z coordinate axes, so I expect the COM to be nearly at {a_x / 2, a_y / 
2, c_z / 2}, but in the dipole moment file the reference point lies outside 
of the slab system. Why? 
If using the reference_point USER_DEFINED, REFERENCE_POINT X Y Z , in which 
X, Y, and Z are of an arbitrary choice, for example at the center of the 
top slab surface (X=a_x / 2 , Y=a_y / 2, Z=0) and the z=0 is exactly the 
surface and below it there is the 2x2x1 slab system (i.e. with negative c_z 
lying in the bottom slab surface), then after calculation I get the info 
from dipole moment file that the y value has the opposite sign. Is there a 
bug in the cp2k during the dipole moment calculation?

After obtaining the calculated dipole moment values, I am able to calculate 
the bulk polarization if the z-dipole moment (which is much larger than in 
x- and y- directions) is divided by the slab volume (2x2x1 unit cells in 
the x,y,z coordinate system). Is that correct? If yes, why the value is 
much smaller compared to the experimental value at room temperature? (e.g. 
the experimental value of the polarization is 0.71 C/m^2  [APL 12 (2009), 
209]  and the calculated value is 0.042 C/m^2 , i.e. it is at least one 
order of magnitude smaller).
Is there any mistake I have calculated or should I take account anything 
for correct polarization calculation?

The input file, output and dipole moment file after 200 geometry 
optimization steps are attached here below.

Kind regards,

Phil

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20180712/cdd0b8b2/attachment.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: input.inp
Type: chemical/x-gamess-input
Size: 4484 bytes
Desc: not available
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20180712/cdd0b8b2/attachment.inp>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: out
Type: application/octet-stream
Size: 5661957 bytes
Desc: not available
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20180712/cdd0b8b2/attachment.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: mom-moments-1_200.dat
Type: application/octet-stream
Size: 2925 bytes
Desc: not available
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20180712/cdd0b8b2/attachment-0001.obj>