[Pw_forum] how to calculate the polarization of a crystal by berry phase

[Jian ZHOU]

Dear all, 

In the example10 of pwscf, it run a berry phase calculation and get
the Born effective charge of Pb.

However, I am now trying to calculate the polarization of a Crystal,
for example for a tetragonal BaTiO3. I firstly do a scf calculation
with the following input file:

 &control
    calculation  = 'scf'
    restart_mode = 'from_scratch'
    pseudo_dir   = '/home/zj/espresso-2.1.4/pseudo/'
    outdir       = '/root/zj/compute14/pwscf/tet_batio3'
    tstress = .true.,
 /
 &system
    ibrav=6
    celldm(1)= 7.5435
    celldm(3)=1.0110
    nat=5
    ntyp=3
    ecutwfc=30.0
    occupations = 'fixed'
    degauss=0.00
 /
 &electrons
    conv_thr = 1e-9,
    mixing_beta = 0.3d0
 /
ATOMIC_SPECIES
  O     15.9994   O.pw91-van_ak.UPF
  Ti    47.867    Ti.pw91-nsp-van.UPF
  Ba    137.327   Ba.pw91-nsp-van.UPF
ATOMIC_POSITIONS
O 0 0.5 0.48500000
O 0.5 0 0.48500000
O 0.5 0.5 -0.02500000
Ti 0.5 0.5 0.50135000
Ba 0 0 0
K_POINTS {automatic}
  4 4 4 1 1 1


And then do a berry phase calculation with the input file:

 &control
    calculation   = 'nscf'
    pseudo_dir    = '/home/zj/espresso-2.1.4/pseudo/'
    outdir        = '/root/zj/compute14/pwscf/tet_batio3'
    lberry        = .true.
    gdir          = 3
    nppstr        = 7
 /
 &system
    ibrav=6
    celldm(1)= 7.5435
    celldm(3)=1.0110
    nat=5
    ntyp=3
    ecutwfc=30.0
    occupations = 'fixed'
    degauss=0.00
 /
 &electrons
    conv_thr      = 1e-9
    mixing_beta = 0.3d0
 /
ATOMIC_SPECIES
  O     15.9994   O.pw91-van_ak.UPF
  Ti    47.867    Ti.pw91-nsp-van.UPF
  Ba    137.327   Ba.pw91-nsp-van.UPF
ATOMIC_POSITIONS
O 0 0.5 0.48500000
O 0.5 0 0.48500000
O 0.5 0.5 -0.02500000
Ti 0.5 0.5 0.50135000
Ba 0 0 0
K_POINTS {automatic}
  4 4 7 1 1 1


The final output is :

                               SUMMARY OF PHASES
                               ~~~~~~~~~~~~~~~~~

                          Ionic Phase: -0.44095 (mod 2)
                     Electronic Phase:  0.63353 (mod 2)
                          TOTAL PHASE:  0.19258 (mod 2)


                             VALUES OF POLARIZATION
                             ~~~~~~~~~~~~~~~~~~~~~~

        The calculation of phases done along the direction of vector 3
        of the reciprocal lattice gives the following contribution to
        the polarization vector (in different units, and being Omega
        the volume of the unit cell):

           P =   1.4687319  (mod  15.2529570)  (e/Omega).bohr

           P =   0.0033843  (mod   0.0351467)  e/bohr^2

           P =   0.1934880  (mod   2.0093958)  C/m^2

        The polarization direction is:  ( 0.00000 , 0.00000 , 1.00000 )


It is said that the final polarization is  0.1934880 C/m^2. And the
experiment value is about 0.26 C/m^2.  A more high cutoff energy and
more k points is tested, and found the polarization does not change
much.

So, I am not sure If the procedure is correct. If it is correct, why
the value is much smaller than the experiment one.

Thank you.

Best wishes,

Jian