How to calculate Polarization Properties

First: What's Polarization?

For more information check  the following link:    http://susi.theochem.tuwien.ac.at/reg_user/textbooks/WIEN2k_lecture-notes_2013/BerryPI.pdf

You can check this video:

Second: How to install BerryPI

The package is included in the version wien2k_2014 and recent versions.

To know how to install manually the package ckeck the following link:

https://github.com/rubel75/BerryPI/tree/master/Installation



Third: Spontaneous Polarization in BaTiO3

BaTiO3 is a ferroelectric material that exhibits a spontaneous polarization. For the calculation of spontaneous polarization of BaTiO3 two reference structures have been chosen. One is tetragonal non-centrosymmetric (lambda1), where the atoms are displaced from their equilibrium centrosymmetric positions in Z direction, and another structure is a centrosymmetric structure (lambda0).

1. We begin with the low symmetry structure lambda1 (non-centrosymmetric). The structure file is located in ../BerryPI/tutorials/tutorial1/lambda1. The head of the file lambda1.struct is shown below

BaTiO3min                                                                      
P                            4 99 P4mm                                         
           RELA                                                              
7.547566  7.547566  7.626934 90.000000 90.000000 90.000000                   
ATOM  -1: X=0.00000000 Y=0.00000000 Z=0.00000000
        MULT= 1          ISPLIT=-2
Ba1        NPT=  781  R0=0.00001000 RMT=    2.5000   Z:  56.00000              
LOCAL ROT MATRIX:    1.0000000 0.0000000 0.0000000
                   0.0000000 1.0000000 0.0000000
                   0.0000000 0.0000000 1.0000000
ATOM  -2: X=0.50000000 Y=0.50000000 Z=0.51517436
        MULT= 1          ISPLIT=-2
Ti1        NPT=  781  R0=0.00005000 RMT=    1.6400   Z:  22.00000              
LOCAL ROT MATRIX:    1.0000000 0.0000000 0.0000000
                   0.0000000 1.0000000 0.0000000
                   0.0000000 0.0000000 1.0000000
ATOM  -3: X=0.50000000 Y=0.50000000 Z=0.97356131
        MULT= 1          ISPLIT=-2
O 1        NPT=  781  R0=0.00010000 RMT=    1.5400   Z:   8.00000              
LOCAL ROT MATRIX:    1.0000000 0.0000000 0.0000000
                   0.0000000 1.0000000 0.0000000
                   0.0000000 0.0000000 1.0000000
ATOM  -4: X=0.50000000 Y=0.00000000 Z=0.48343742
        MULT= 2          ISPLIT= 8
    -4: X=0.00000000 Y=0.50000000 Z=0.48343742
O 2        NPT=  781  R0=0.00010000 RMT=    1.5400   Z:   8.00000              
LOCAL ROT MATRIX:    1.0000000 0.0000000 0.0000000
                   0.0000000 1.0000000 0.0000000
                   0.0000000 0.0000000 1.0000000
 8      NUMBER OF SYMMETRY OPERATIONS
...

Note Z-coordinates of atoms that deviate from "ideal" centrosymmetric positions.

2. Initialize calculation

init_lapw -b -vxc 13 -ecut -6 -numk 230

This should produce a k-mesh of 6x6x6. More details on sensitivity of the result to convergence parameters can be found at the end of this tutorial.

3. Run SCF sycle

run_lapw

You are welcome to specify additional convergence criteria using -ec and -cc flags. It is the user's responsibility to check the convergence with respect to the quantity of interest (polarization in this case). Laurence Marks's suggestion is to change TOT to FOR in case.in2c and then monitor :FCHECK in case.scfm, checking that the calculation is converged to some adequate level (perhaps 2.0 or better). More details on convergence studies can be found at the end of this tutorial.

4. Run BerryPI

berrypi -k 6:6:6

Here we use the same k-mesh density as in the SCF calculation. It does not have to be that way. The density can be increased and tested for convergence. However, the experience shows that the polarization converges rapidly with increasing k-mesh.

The results of the calculation are

...
CALCULATION OF ELECTRONIC POLARIZATION
=======================================================================================
Value                           |  spin   |    dir(1)    |    dir(2)    |    dir(3)
---------------------------------------------------------------------------------------
Berry phase (rad) [-pi ... +pi]    sp(1)  [-3.149704e-11, -5.662171e-13,  1.526306e+00]
Berry phase (rad)                  up+dn  [-6.299408e-11, -1.132434e-12,  3.052611e+00]
Berry phase (rad) [-pi ... +pi]    up+dn  [-6.299408e-11, -1.132434e-12,  3.052611e+00]
Electronic polarization (C/m2)     sp(1)  [-9.964856e-12, -1.791366e-13,  4.879618e-01]
=======================================================================================

CALCULATION OF IONIC POLARIZATION
=======================================================================================
Elem.|  Fractional coord.  |  spin | Zion |    dir(1)    |    dir(2)    |    dir(3)
---------------------------------------------------------------------------------------
                                        +------------ Ionic phase (rad) ------------+
Ba  (0.0000, 0.0000, 0.0000)  sp(1) 10.00 [ 0.000000e+00,  0.000000e+00,  0.000000e+00]
Ti  (0.5000, 0.5000, 0.5152)  sp(1) 12.00 [ 3.769911e+01,  3.769911e+01,  3.884323e+01]
O  (0.5000, 0.5000, 0.9736)  sp(1)  6.00 [ 1.884956e+01,  1.884956e+01,  3.670240e+01]
O  (0.5000, 0.0000, 0.4834)  sp(1)  6.00 [ 1.884956e+01,  0.000000e+00,  1.822516e+01]
O  (0.0000, 0.5000, 0.4834)  sp(1)  6.00 [ 0.000000e+00,  1.884956e+01,  1.822516e+01]
---------------------------------------------------------------------------------------
Total ionic phase (rad)       sp(1)       [ 7.539822e+01,  7.539822e+01,  1.119960e+02]
Total ionic phase wrap. (rad) sp(1)       [ 0.000000e+00,  0.000000e+00, -1.101384e+00]
Ionic polarization (C/m2)     sp(1)       [ 0.000000e+00,  0.000000e+00, -1.760570e-01]
=======================================================================================

SUMMARY OF POLARIZATION CALCULATION
=======================================================================================
Value                           |  spin   |    dir(1)    |    dir(2)    |    dir(3)
---------------------------------------------------------------------------------------
Electronic polarization (C/m2)     sp(1)  [-9.964856e-12, -1.791366e-13,  4.879618e-01]
Ionic polarization (C/m2)          sp(1)  [ 0.000000e+00,  0.000000e+00, -1.760570e-01]
Tot. spin polariz.=Pion+Pel (C/m2) sp(1)  [-9.964856e-12, -1.791366e-13,  3.119048e-01]
---------------------------------------------------------------------------------------
TOTAL POLARIZATION (C/m2)          both   [-9.964856e-12, -1.791366e-13,  3.119048e-01]
=======================================================================================
...

It is tempting to interpret the results as the presence of a spontaneous polarization in Z-direction of the magnitude 0.31 C/m2. However this conclusion is premature at this point, since any observable polarization is a difference with respect to a reference structure (in our case, the centrosymmetric structure).

5. Next we proceed with the case lambda0 (centrosymmetric). The structure file is located in ../BerryPI/tutorials/tutorial1/lambda0. The head of the file lambda0.struct is shown below

BaTiO3min                                                                      
P                            4 99 P4mm                                         
           RELA                                                              
7.547566  7.547566  7.626934 90.000000 90.000000 90.000000                   
ATOM  -1: X=0.00000000 Y=0.00000000 Z=0.00000000
        MULT= 1          ISPLIT=-2
Ba1        NPT=  781  R0=0.00001000 RMT=    2.5000   Z:  56.00000              
LOCAL ROT MATRIX:    1.0000000 0.0000000 0.0000000
                   0.0000000 1.0000000 0.0000000
                   0.0000000 0.0000000 1.0000000
ATOM  -2: X=0.50000000 Y=0.50000000 Z=0.50000000
        MULT= 1          ISPLIT=-2
Ti1        NPT=  781  R0=0.00005000 RMT=    1.6400   Z:  22.00000              
LOCAL ROT MATRIX:    1.0000000 0.0000000 0.0000000
                   0.0000000 1.0000000 0.0000000
                   0.0000000 0.0000000 1.0000000
ATOM  -3: X=0.50000000 Y=0.50000000 Z=0.00000000
        MULT= 1          ISPLIT=-2
O 1        NPT=  781  R0=0.00010000 RMT=    1.5400   Z:   8.00000              
LOCAL ROT MATRIX:    1.0000000 0.0000000 0.0000000
                   0.0000000 1.0000000 0.0000000
                   0.0000000 0.0000000 1.0000000
ATOM  -4: X=0.50000000 Y=0.00000000 Z=0.50000000
        MULT= 2          ISPLIT= 8
    -4: X=0.00000000 Y=0.50000000 Z=0.50000000
O 2        NPT=  781  R0=0.00010000 RMT=    1.5400   Z:   8.00000              
LOCAL ROT MATRIX:    1.0000000 0.0000000 0.0000000
                   0.0000000 1.0000000 0.0000000
                   0.0000000 0.0000000 1.0000000
 8      NUMBER OF SYMMETRY OPERATIONS
...

Note Z-coordinates of atoms that correspond to the centrosymmetric positions. The symmetry operations are identical (!) to the lambda1 case, in spite of a higher symmetry of the lambda0 structure. It is done intentionally to keep both calculations as identically as possible. If we let Wien2k realize the inversion symmetry, it will switch to the "real" mode and the electronic phase will disappear.

6. Now we need to perform polarization calculation for this structure. In order to preserve settings, the following commands are executed, presumably you are still located in ../BerryPI/tutorials/tutorial1/lambda1 directory

cp * ../lambda0 # Copy all files from lambda1 to lambda0 directory

cd ../lambda0 # Change the current directory to lambda0

rm lambda1.struct # Remove the lambda1.struct file

rename_files lambda1 lambda0 # Rename all lambda1.* files to lambda0.* files

x kgen # Restore original k-mesh taking into account the symmetry with 230 k-points (shifted)

x dstart # Initialize the electron density for the new structure

run_lapw # Perform standard SCF calculation

At this point we intentionally avoid full initialization using init_lapw... to prevent Wien2k from realizing the full symmetry of the centrosymmetric structure.

7. Run BerryPI

berrypi -k 6:6:6

Once the calculation is completed the results will be outputed in the form

...
CALCULATION OF ELECTRONIC POLARIZATION
=======================================================================================
Value                           |  spin   |    dir(1)    |    dir(2)    |    dir(3)
---------------------------------------------------------------------------------------
Berry phase (rad) [-pi ... +pi]    sp(1)  [ 3.133816e-14,  4.003657e-14,  3.773516e-12]
Berry phase (rad)                  up+dn  [ 6.267632e-14,  8.007314e-14,  7.547033e-12]
Berry phase (rad) [-pi ... +pi]    up+dn  [ 6.267632e-14,  8.007314e-14,  7.547033e-12]
Electronic polarization (C/m2)     sp(1)  [ 9.914590e-15,  1.266654e-14,  1.206398e-12]
=======================================================================================

CALCULATION OF IONIC POLARIZATION
=======================================================================================
Elem.|  Fractional coord.  |  spin | Zion |    dir(1)    |    dir(2)    |    dir(3)
---------------------------------------------------------------------------------------
                                         +------------ Ionic phase (rad) ------------+
Ba  (0.0000, 0.0000, 0.0000)  sp(1) 10.00 [ 0.000000e+00,  0.000000e+00,  0.000000e+00]
Ti  (0.5000, 0.5000, 0.5000)  sp(1) 12.00 [ 3.769911e+01,  3.769911e+01,  3.769911e+01]
O  (0.5000, 0.5000, 0.0000)  sp(1)  6.00 [ 1.884956e+01,  1.884956e+01,  0.000000e+00]
O  (0.5000, 0.0000, 0.5000)  sp(1)  6.00 [ 1.884956e+01,  0.000000e+00,  1.884956e+01]
O  (0.0000, 0.5000, 0.5000)  sp(1)  6.00 [ 0.000000e+00,  1.884956e+01,  1.884956e+01]
---------------------------------------------------------------------------------------
Total ionic phase (rad)       sp(1)       [ 7.539822e+01,  7.539822e+01,  7.539822e+01]
Total ionic phase wrap. (rad) sp(1)       [ 0.000000e+00,  0.000000e+00,  0.000000e+00]
Ionic polarization (C/m2)     sp(1)       [ 0.000000e+00,  0.000000e+00,  0.000000e+00]
=======================================================================================

SUMMARY OF POLARIZATION CALCULATION
=======================================================================================
Value                           |  spin   |    dir(1)    |    dir(2)    |    dir(3)
---------------------------------------------------------------------------------------
Electronic polarization (C/m2)     sp(1)  [ 9.914590e-15,  1.266654e-14,  1.206398e-12]
Ionic polarization (C/m2)          sp(1)  [ 0.000000e+00,  0.000000e+00,  0.000000e+00]
Tot. spin polariz.=Pion+Pel (C/m2) sp(1)  [ 9.914590e-15,  1.266654e-14,  1.206398e-12]
---------------------------------------------------------------------------------------
TOTAL POLARIZATION (C/m2)          both   [ 9.914590e-15,  1.266654e-14,  1.206398e-12]
=======================================================================================
...

8. Spontaneous polarization is obtained by taking a difference in polarization between P(lambda1) and P(lambda0)

 Ps = [-9.964856e-12, -1.791366e-13,  3.119048e-01] -
      [ 9.914590e-15,  1.266654e-14,  1.206398e-12]
    = [0.0,  0.0,   0.312] C/m2

It should be noted that calculations are sometimes affected by "pi"-wrapping artefacts that produce spurious results. The reason is that it is impossible to distinguish between +pi and -pi on the 2pi circle. It is therefore advised to perform an additional calculation for an intermediate step (lambda0.5) that corresponds to interpolated coordinates between lambda0 and lambda1 structures. If no artefacts are present, the obtained value of Ps will be about one-half of the value determined above. Users are encouraged to perform similar check for this tutorial.

Here we deal with a structure where the lattice vectors are mutually orthogonal (all angles are 90 deg). In this case the electronic and ionic phases are additive, so do the polarizations. However, this is generally not a case, since the electronic polarization is computed in k-space and the ionic component is evaluated in the real space. Both sets of lattice vectors generally do not coincide and more manual work is required to transform the results into the common coordinate system (see Tutorial 3).

Results of convergence tests:

RKmax	k-mesh in SCF	k-mesh in BerryPI	-ec	-cc	:FCHECK	Ps (C/m2)
6	6x6x6	6x6x6	0.001	no	1.05	3.112792e-01
7	6x6x6	6x6x6	0.001	no	-8.66	3.121672e-01
7	6x6x6	6x6x6	0.0001	no	1.92	3.120008e-01
7	6x6x6	6x6x6	0.0001	0.001	0.17	3.120284e-01
7	6x6x6	4x4x4	0.0001	0.001	0.17	3.090667e-01
7	6x6x6	2x2x2	0.0001	0.001	0.17	3.081526e-01
7	6x6x6	1x1x1	0.0001	0.001	0.17	3.086679e-01
7	6x6x6	8x8x8	0.0001	0.001	0.17	3.134040e-01
7	6x6x6	10x10x10	0.0001	0.001	0.17	3.140543e-01
7	4x4x4	10x10x10	0.0001	0.001	0.35	3.137155e-01
7	10x10x9	10x10x10	0.0001	0.001	0.32	3.140688e-01
7	6x6x6	12x12x12	0.0001	0.001	0.17	3.144094e-01
7	6x6x6	14x14x14	0.0001	0.001	0.17	3.146237e-01
7	6x6x6	16x16x16	0.0001	0.001	0.17	3.147647e-01
8	6x6x6	16x16x16	0.0001	0.001	-0.45	3.147154e-01

Calculations in the table above were performed using identical structure file lambda1.struct. It is important to perform full relaxation of internal degrees of freedom, since inaccuracy in atomic positions will directly influence the value of polarization.

Note: completed using WIEN2k 14.1, BerryPI 1.3.3, Python 2.7.3, Numpy 1.9.2

Fourth: Calculation of Born effective charge of GaAs

For the calculation of Born effective charge of As in GaAs one of the 4 As atoms has been displaced along Z-axis form it's equilibrium position by +0.01(lambda1) and -0.01(lambda2) in fractional coordinates.


1 Case lambda1

Calculation of total phase (sum of electronic and ionic phase) for state where the As atom has been displaced by +0.01 (fractional coordinate) from its equilibrium position.

1.1 Change the current directory to ~/tutorial/tutorial2/lambda1 

1.2 Perform WIEN2k initialization 

$ init_lapw -b -vxc 13 -ecut -6 -numk 230

Here "-vxc 13" stands for PBE-GGA as exchange correlation function."-ecut -6" means the separation  energy of -6 Ry has been chosen to separate core electron from valance electron. "-numk 230" means that 230 k points has been chosen in Brillouion zone which generates 6*6*6 size k-mesh in the symmetric Brillouion zone

1.3 Execute WIEN2k scf calculation
 
$ run_lapw

in order to optimize the electron density.

1.4 Run BerryPI using python 
 
$ berrypi -k6:6:6

Here BerryPI program is running for the case (GaAs) located in the current directory.
6:6:6 means the calculation is being done using 6*6*6 k-mesh in the full Brillouin zone (non symmetric) with a total of 216 k points.

Note: k-mesh in BerryPI should not necessarily be identical to that used in the SCF cycle

1.5 Once the calculation is completed the result will be printed like this

                    ---PHASES/2*PI IN [0 to 2]RANGE---

TOTAL PHASE/(2*PI):  [0.99999999991459032, 1.0000000091132539, 0.97814709926563959]

                   ---PHASES/2*PI IN [-1 to +1]RANGE---
 
TOTAL PHASE/(2*PI):  [0.99999999991459032, -0.99999999088674607, 0.97814709926563959]

Here three total phase (sum of electronic and ionic phase)  values corresponds to X,Y and Z components of total phase, respectively. As the structure has only been perturbed in Z direction, only Z component of total phase has to be considered.  

Note: The total phase has been reported twice for different pi wrapping approaches. In this particular case both the phase values are the identical which is generally may not be the case. 


2 Case lambda2

Calculation of total phase (sum of electronic and ionic phase) for state where the As atom has been displaced by -0.01 (fractional coordinate) from its equilibrium position.

2.1 Copy all files from lambda1 to lambda2 directory
 
$ cp * ../lambda2

2.2 Change the current directory to lambda2 

$ cd ../lambda2

2.3 Remove the lambda1.struct file. 

$ rm lambda1.struct

2.4 Rename all lambda1.* files to lambda2.* files 

$ rename_files lambda1 lambda2


2.5 Restore original k-mesh taking into account the symmetry 

$ x kgen

with 230 k-points (Shifted) 

2.6 Initialize the electron density according to new structure
 
$ x dstart 

2.7 Run standard scf cycle 

$ run_lapw

2.8 Run BerryPI 
  
$ berrypi -k6:6:6

2.9 Once the calculation is completed the results will be printed like this
               ---PHASES/2*PI IN [0 to 2]RANGE---

 TOTAL PHASE/(2*PI):        [1.0000000009195591, 1.0000000089644974, 1.0218558597965477]

              ---PHASES/2*PI IN [-1 to +1]RANGE---

TOTAL PHASE/(2*PI):        [-0.9999999990804409, -0.99999999103550263, -0.97814414020345231]

Note: This time, the two reported phases (sum of electronic and ionic phase) are different from each other and only one of them needs to be considered as explained below.


3. Born effective charge

Calculation of Born effective charge using Z component of total phase value obtained for lambda1 and lambda2 case. The Born effective charge can be calculated using the following formula,

Z*_zz=(1/2*pi)*d_phi/d_rho 
Where,
Z*_zz is the born effective charge in Z direction for applied perturbation in Z direction.
d_phi is difference in total phase (sum of electronic and ionic phase) between the two structure(in the units of 2*pi)= 0.97814709926563959-1.0218558597965477= -0.04370876053090811(in the units of 2*pi)
d_rho is the relative displacement of the particular atom(in fractional coordinate)= 0.01-(-0.01)= 0.02

So, 
Z*_zz=(1/2*pi)* (-0.04370876053090811*2*pi)/0.02= -2.1854380265454054

Note: Here 2 values of total phase (sum of electronic and ionic phase) for both the cases are reported. When calculating the difference in total phase between two structures, one has to be careful of pi wrapping artifact. For example in this current study if the difference in phase is taken from the 2*pi modulo values this will lead to a phase difference of -0.04370876053090811. On the other hand if the difference is taken from -pi to +pi domain values this lead to a phase difference of 1.956291239469092. It has be realized that both this difference values represents the difference between both the cases but they just consider different paths when the phase difference is considered in a 2*pi(360 degree) circle[see "BerryPI: A software for studying polarization of crystalline solids with WIEN2k density functional all-electron package" for more details]. Taking the biggest difference (1.956291239469092) in phase will lead to inappropriate results. So while taking the difference in phase one has to take the smallest number which is -0.04370876053090811 in this case. This should also be taken care of while calculating the difference in polarization. In this particular case putting the smallest difference will lead to a born effective charge of -2.1854380265454054 for As.

Fifth: Calculation of Born effective charge of GaN

Check the following video:

Sixth: Polarization in GaN

Here we compute the polarization difference between wurtzite (W) and zinc-blende (ZB) structure of GaN. This polarization difference is responsible for the presence of a build-in electric field in heterostructures of group-III nitrides. The two structures (wurtzite and zinc-blende) are shown below.



Our basic steps will include obtaining SCF electron density for these structures followed by calculation of polarization and their subsequent difference
Ps = P(ZB) - P(W)


Note: Ps is a property of the interface rather than GaN alone.

Wurtzite GaN

Let's create a work directory GaN-W

Creating the structure file GaN-W.struct

GaN 4 atoms
H   LATTICE,NONEQUIV.ATOMS:  2186_P63mc
MODE OF CALC=RELA unit=bohr
  5.963131  5.963131  9.722374 90.000000 90.000000120.000000
ATOM  -1: X=0.33333333 Y=0.66666666 Z=0.50000000
          MULT= 2          ISPLIT= 4
      -1: X=0.66666667 Y=0.33333333 Z=0.00000000
Ga1        NPT=  781  R0=0.00005000 RMT=    1.9300   Z: 31.0
LOCAL ROT MATRIX:    1.0000000 0.0000000 0.0000000
                     0.0000000 1.0000000 0.0000000
                     0.0000000 0.0000000 1.0000000
ATOM  -2: X=0.33333333 Y=0.66666666 Z=0.87639300
          MULT= 2          ISPLIT= 4
      -2: X=0.66666667 Y=0.33333333 Z=0.37639300
N 1        NPT=  781  R0=0.00010000 RMT=    1.5000   Z:  7.0
LOCAL ROT MATRIX:    1.0000000 0.0000000 0.0000000
                     0.0000000 1.0000000 0.0000000
                     0.0000000 0.0000000 1.0000000
  12      NUMBER OF SYMMETRY OPERATIONS
 1 0 0 0.00000000
 0 1 0 0.00000000
 0 0 1 0.00000000
       1
...


It is the best to get the file from a repository BerryPI/tutorials/tutorial4/GaN-W/GaN-W.struct. Here the lattice parameters as well as atomic positions are fully optimized.

Initialize calculation


init_lapw -b -vxc 5 -rkmax 7 -numk 300

Here we opt for LDA exchange-correlation functional and the resultant 8x8x4 k-mesh, which is dense enough for the purpose of this calculations.

Run SCF calculation


run_lapw

It might be a good idea to speed up the process by getting more cores involved export OMP_NUM_THREADS=4

Run BerryPI


berrypi -k 8:8:8

The results of the calculation are
...
SUMMARY OF POLARIZATION CALCULATION
=======================================================================================
Value                           |  spin   |    dir(1)    |    dir(2)    |    dir(3)
---------------------------------------------------------------------------------------
Electronic polarization (C/m2)     sp(1)  [ 2.223396e-10,  1.205149e-07, -5.612307e-02]
Ionic polarization (C/m2)          sp(1)  [ 1.750356e-10, -2.049427e-07, -4.386009e-01]
Tot. spin polariz.=Pion+Pel (C/m2) sp(1)  [ 3.973752e-10, -8.442784e-08, -4.947239e-01]
---------------------------------------------------------------------------------------
TOTAL POLARIZATION (C/m2)          both   [ 3.973752e-10, -8.442784e-08, -4.947239e-01]
=======================================================================================
...


As we can see the polarization (-0.4947 C/m2) is non-zero only along Z-axis, which corresponds to [0001] direction.

Zinc-blende GaN

Now we can switch to another work directory GaN-ZB

Creating the structure file GaN-ZB.struct

GaN cubic
H   LATTICE,NONEQUIV.ATOMS:  6156_P3m1
MODE OF CALC=RELA unit=bohr
  5.963131  5.963131 14.606628 90.000000 90.000000120.000000
ATOM  -1: X=0.00000000 Y=0.00000000 Z=0.00000000
          MULT= 1          ISPLIT= 4
Ga1        NPT=  781  R0=0.00005000 RMT=    1.9300   Z: 31.0
LOCAL ROT MATRIX:    1.0000000 0.0000000 0.0000000
                     0.0000000 1.0000000 0.0000000
                     0.0000000 0.0000000 1.0000000
ATOM  -2: X=0.00000000 Y=0.00000000 Z=0.25000000
          MULT= 1          ISPLIT= 4
N 1        NPT=  781  R0=0.00010000 RMT=    1.5000   Z:  7.0
LOCAL ROT MATRIX:    1.0000000 0.0000000 0.0000000
                     0.0000000 1.0000000 0.0000000
                     0.0000000 0.0000000 1.0000000
ATOM  -3: X=0.33333333 Y=0.66666667 Z=0.33333333
          MULT= 1          ISPLIT= 4
Ga2        NPT=  781  R0=0.00005000 RMT=    1.9300   Z: 31.0
LOCAL ROT MATRIX:    1.0000000 0.0000000 0.0000000
                     0.0000000 1.0000000 0.0000000
                     0.0000000 0.0000000 1.0000000
ATOM  -4: X=0.33333333 Y=0.66666667 Z=0.58333333
          MULT= 1          ISPLIT= 4
N 2        NPT=  781  R0=0.00010000 RMT=    1.5000   Z:  7.0
LOCAL ROT MATRIX:    1.0000000 0.0000000 0.0000000
                     0.0000000 1.0000000 0.0000000
                     0.0000000 0.0000000 1.0000000
ATOM  -5: X=0.66666667 Y=0.33333333 Z=0.66666667
          MULT= 1          ISPLIT= 4
Ga3        NPT=  781  R0=0.00005000 RMT=    1.9300   Z: 31.0
LOCAL ROT MATRIX:    1.0000000 0.0000000 0.0000000
                     0.0000000 1.0000000 0.0000000
                     0.0000000 0.0000000 1.0000000
ATOM  -6: X=0.66666667 Y=0.33333333 Z=0.91666667
          MULT= 1          ISPLIT= 4
N 3        NPT=  781  R0=0.00010000 RMT=    1.5000   Z:  7.0
LOCAL ROT MATRIX:    1.0000000 0.0000000 0.0000000
                     0.0000000 1.0000000 0.0000000
                     0.0000000 0.0000000 1.0000000
   6      NUMBER OF SYMMETRY OPERATIONS
 1 0 0 0.00000000
 0 1 0 0.00000000
 0 0 1 0.00000000
       1
...


It is the best to get the file from a repository BerryPI/tutorials/tutorial4/GaN-ZB/GaN-ZB.struct. The structure has 6 atoms vs 4 atoms for the wurtzite structure. Here the lattice parameter a = b = 5.963131 Bohr is the same as for the wurtzite structure above (otherwise how do we imagine the heterointerface), however we cannot use the same c, even re-scaled. The ideal zinc-blende c = a*sqrt(6) = 14.606628 Bohr. Next steps are identical to the wirtzite case.

Initialize calculation


init_lapw -b -vxc 5 -rkmax 7 -numk 200

Here we opt for 200 k-points that results in 8x8x2 k-mesh. This is a consequence of smaller Brillouin zone for the ZB structure in comparison to W structure.

Run SCF calculation


run_lapw

Run BerryPI


berrypi -k 8:8:8

Here we keep k-mesh used for the Berry phase integration identical to the wurtzite case. The results of the calculation are
...
SUMMARY OF POLARIZATION CALCULATION
=======================================================================================
Value                           |  spin   |    dir(1)    |    dir(2)    |    dir(3)
---------------------------------------------------------------------------------------
Electronic polarization (C/m2)     sp(1)  [-1.219101e-03,  1.219101e-03, -1.727041e-04]
Ionic polarization (C/m2)          sp(1)  [ 1.165061e-10,  1.165088e-10, -4.644818e-01]
Tot. spin polariz.=Pion+Pel (C/m2) sp(1)  [-1.219101e-03,  1.219101e-03, -4.646545e-01]
---------------------------------------------------------------------------------------
TOTAL POLARIZATION (C/m2)          both   [-1.219101e-03,  1.219101e-03, -4.646545e-01]
=======================================================================================
...


This result can be quite puzzling, since we would expect the polarization to vanish in the ZB structure for symmetry reasons. The explanation is provided by David Vanderbilt and R. D. King-Smith Phys. Rev. B 48, 4442 (1993), Sec. III.E.

Polarization difference

Finally, we can compute the polarization difference
Ps = -4.646545e-01 - (-4.947239e-01) = 0.030 C/m2


The obtained result agrees with Ps = -0.029 C/m2 reported by Bernardini et al. Phys. Rev. B 56, R10024 (1997); it is also consistent with Ps = -0.022 C/m2 measured by Lähnemann et al. Phys. Rev. B 86, 081302(R) (2012) for the heterointerface W/ZB GaN. The sign of P and Ps depends on the orientation of our structure and will flip upon flipping the structure along c-axis.

Note: completed using WIEN2k 13.1, WIEN2WANNIER 1.0-beta2, BerryPI 1.2, Python 2.7.4, Numpy 1.6.2