Tutorial 3: Non orthogonal lattice vectors

 

 

GaAs (2-atom basis)

Here we perform calculation of a Born effective charge in GaAs. When lattice angles are not 90 deg, a complication arises due to the real space conventional lattice vectors (used to express ionic polarization) being non-collinear with real space primitive lattice vectors (used to express electronic polarization). As a result, you cannot simply add ionic and electronic polarizations. BerryPI performs an additional transformation to Cartesian coordinates (linear algebra) as shown below.

1. Creating a template structure file GaAs1.struct

GaAs
F   LATTICE,NONEQUIV.ATOMS:  2
MODE OF CALC=RELA unit=bohr
 10.683093 10.683093 10.683093 90.000000 90.000000 90.000000
ATOM  -1: X=0.00000000 Y=0.00000000 Z=0.00000000
          MULT= 1          ISPLIT= 8
Ga         NPT=  781  R0=0.00005000 RMT=    2.0000   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.25100000 Y=0.25200000 Z=0.25300000
          MULT= 1          ISPLIT= 8
As         NPT=  781  R0=0.00005000 RMT=    2.0000   Z: 33.0
LOCAL ROT MATRIX:    1.0000000 0.0000000 0.0000000
                     0.0000000 1.0000000 0.0000000
                     0.0000000 0.0000000 1.0000000
   1      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/tutorial3/GaAs1/GaAs1.struct and modify ATOM -2: X=0.25100000 Y=0.25200000 Z=0.25300000, since some spaces may not be correctly displayed above. We intend to displace As and give it the least symmetric displacement from ideal 1/4 1/4 1/4. This will give us a maximum flexibility in the choice of a displacement vector later.

2. Initialize calculation

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

3. Edit the structure file GaAs1.struct

Introduce a displacement of +0.001 fractional units along Z-axis

...
ATOM  -2: X=0.25000000 Y=0.25000000 Z=0.25100000
...

Alternatively, you can get the file at BerryPI/tutorials/tutorial3/GaAs1/GaAs1.struct

4. Initialize the electron density only

x dstart

(!) Do not rerun init_lapw as it may realize a higher symmetry. The intention is to keep symmetry unchanged between subsequent runs.

5. Run SCF sycle

run_lapw

6. Run BerryPI

berrypi -k 6 6 6

The results of the calculation with [0 .. +2pi] phase wrapping are

SUMMARY OF POLARIZATION CALCULATION IN CARTESIAN COORDINATES
=======================================================================================
Value                           |  spin   |       X      |       Y      |       Z
---------------------------------------------------------------------------------------
Electronic polarization (C/m2)     sp(1)  [ 1.010064e+00,  1.010064e+00,  9.756000e-01]
Ionic polarization (C/m2)          sp(1)  [ 1.503956e+00,  1.503956e+00,  1.534036e+00]
Tot. spin polariz.=Pion+Pel (C/m2) sp(1)  [ 2.514021e+00,  2.514021e+00,  2.509636e+00]
---------------------------------------------------------------------------------------
TOTAL POLARIZATION (C/m2)          both   [ 2.514021e+00,  2.514021e+00,  2.509636e+00]
=======================================================================================

7. Perform a similar run for the displacement of As by -0.001 along the Z-axis

• Copy case: cp -r ../GaAs1 ../GaAs2; cd ../GaAs2

• Rename files: rename_files GaAs1 GaAs2

• Edit structure file GaAs2.struct:

  ...
  ATOM  -2: X=0.25000000 Y=0.25000000 Z=0.24900000
  ...
  

  Alternatively, you can get the file at BerryPI/tutorials/tutorial3/GaAs2/GaAs2.struct

• Restore k-mesh: x kgen with 800 points (shifted)

• Reinitialize the electron density: x dstart

• Run SCF cycle: rm *scf; rm *broy*; run_lapw

• Run BerryPI: berrypi -k 6 6 6

• Results of the calculation with [0 .. +2pi] phase wrapping are

SUMMARY OF POLARIZATION CALCULATION IN CARTESIAN COORDINATES
=======================================================================================
Value                           |  spin   |       X      |       Y      |       Z
---------------------------------------------------------------------------------------
Electronic polarization (C/m2)     sp(1)  [ 1.010061e+00,  1.010061e+00,  1.044524e+00]
Ionic polarization (C/m2)          sp(1)  [ 1.503956e+00,  1.503956e+00,  1.473877e+00]
Tot. spin polariz.=Pion+Pel (C/m2) sp(1)  [ 2.514018e+00,  2.514018e+00,  2.518402e+00]
---------------------------------------------------------------------------------------
TOTAL POLARIZATION (C/m2)          both   [ 2.514018e+00,  2.514018e+00,  2.518402e+00]
=======================================================================================

Please note that BerryPI will output results with [-pi .. +pi] and [0 .. +2pi] phase wrapping. Here we chose [0 .. +2pi] phase wrapping since [-pi .. +pi] wrapping leads to a situation where the electronic phase "jumps" by almost 2*pi (it is too much for such a small change in atomic positions) between two structures:

Berry phase (rad) [-pi ... +pi]    up+dn  [ 3.056875e+00,  3.056875e+00, -3.010334e+00]
Berry phase (rad) [-pi ... +pi]    up+dn  [-3.010348e+00, -3.010348e+00,  3.056869e+00]

while the ionic phase changes "smoothly"

Total ionic phase wrap. (rad) sp(1)       [-1.570796e+00, -1.570796e+00, -1.476549e+00]
Total ionic phase wrap. (rad) sp(1)       [-1.570796e+00, -1.570796e+00, -1.665044e+00]

It is an artefact related to 2*pi wrapping.

8. Total polarization difference

dPtot = Ptot(GaAs1) - Ptot(GaAs2)
    = [ 2.514021e+00,  2.514021e+00,  2.509636e+00]
     -[ 2.514018e+00,  2.514018e+00,  2.518402e+00]
    = [ 3.3e-06,       3.3e-06,      -8.766000e-03] C/m2

9. Change in the position of As atom in Cartesian coordinates is

dr(i) = sum_j [du(j)*R_conv(j,i)]

where i,j = 1,2,3 are components of conventional lattice vectors R_conv and u are fractional coordinates of the displaced atom. In WIEN2k all fractional coordinates are related to conventional (not primitive!) lattice vectors. For some (but not all!) lattices conventional and primitive lattices are the same. Here we find for the cubic lattice

dr = [0, 0, 0.251-0.249]*10.68309
 = [0, 0, 0.02137] bohr
 = [0, 0, 1.13065e-12] m

Here 10.68309 bohr is the length of conventional lattice vectors (see output CALCULATION OF IONIC POLARIZATION section):

...
The ionic positions and associated polarization vector is presented in this coord. system:
   dir(1) = [ 1.068309e+01,  0.000000e+00,  0.000000e+00] bohr
   dir(2) = [ 0.000000e+00,  1.068309e+01,  0.000000e+00] bohr
   dir(3) = [ 0.000000e+00,  0.000000e+00,  1.068309e+01] bohr

12. Finally, the effective charge is

Z*_As(3,3) = (Omega/e) * dPtot(3)/dr(3)
         = (304.81128*(5.29177e-11)^3/1.60217657e-19) * (-8.766e-03)/1.13065e-12
           = -2.186

Here the factor (5.29177e-11)^3 accounts for the Bohr3 -> m3 units conversion.

Obtained result is remarkably close the the value of Z* = -2.185 found in tutorial 2 using the supercell and also reproduces exactly the result of Wang and Vanderbilt Phys. Rev. B 75, 115116 (2007). The negative sign indicates a more electronegative element. Other components of the effective charge tensor Z*(i,j) (i,j=1,2,3) can be found in a similar way. For this particular structure, the symmetry implies identical diagonal components Z*(i,i) and almost zero off-diagonal components Z*(i,j). Eventually, the effective charge tensor takes the form

        | -2.19  0     0    |
Z*_As = |  0    -2.19  0    |
        |  0     0    -2.19 |

Note: completed using WIEN2k 19.2, BerryPI v. Dec 13 2020, Python 3.8.3, Numpy 1.19.1

 https://github.com/spichardo/BerryPI/tree/master/tutorials/tutorial3