-
Notifications
You must be signed in to change notification settings - Fork 25
Tutorial 1: 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).
- 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.
- Initialize calculation
init_lapw -b -vxc 13 -ecut -6 -numk 230 -nohdlo
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. The -nohdlo switch is needed to avoid generation of high-derivative local orbitals by LSTART. Wien2wannier is incompatible with high-derivative local orbitals.
- 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.
- 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).
- 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.
- 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/lambda1directory
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.
- Run BerryPI
berrypi -k 6 6 6
Once the calculation is completed the results will be outputted 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]
=======================================================================================
...
- 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