Computing the piezoelectric tensor for AlN

Version: 2016.0

In this tutorial you will learn how to obtain the piezoelectric tensor and its coefficients with ATK. Specifically, you will:

  1. compute the piezoelectric tensor;
  2. calculate the \({e}_{33}\) coefficient by following these steps:
  • calculate the changes in the relevant structural parameters due to strain;
  • compute the piezoelectric tensor in the clamped-ion approximation and extract the coefficient \({e}_{33}(0)\);
  • compute the Born effective charge.

introbar

Introduction

Piezoelectric materials exhibit an induced electric polarization upon the application of external macroscopic strain. The polarization can be reversed by applying an external electric field. These materials find application in a variety of Microelectromechanical Systems (MEMS).

In this tutorial we are going to study the piezoelectric constant of AlN. This synthetic ceramic material is used in a variety of MEMS such as Surface Acoustic Wave sensors (SAWs), and Film Bulk Acoustic Resonators (FBAR).

The piezoelectric tensor

In the absence of external fields, the total macroscopic polarization \(P\) of a solid is the sum of the spontaneous polarization \({P}_{eq}\) (strain independent) of the equilibrium structure, and of the piezoelectric polarization induced by strain \({P}_{p}\) (strain dependent).

\[{P} = {P}_{eq} + {P}_{p}\]

The piezoelectric tensor can be expressed as:

\[{\gamma}_{\delta \alpha} = \frac{\Delta {P}_{\delta}}{\Delta {\epsilon}_{\alpha}}\]

In ATK \({\gamma}_{\delta \alpha}\) is calculated using a finite-difference approach and \({P}\) is calculated using a Berry-phase approach. You can refer to the Polarization tutorial for more information on how \({P}\) is evaluated using modern theory of polarization.

Important

The current implementation requires the unit cell to be orthogonal (simple cubic, tetragonal, or orthorhombic) for the calculation of \(P\). There is no explicit check for this in the code, but the results cannot be expected to be correct if a non-orthogonal cell is used.

Computing the piezoelectric tensor

The first thing you should do is to import the AlN (Hexagonal) structure in the Stash, optimize its bulk structure and change the cell geometry to Orthorhombic. In this tutorial, you will use a AlN (Orthorhombic) bulk optimized at PBE using FHI pseudopotentials and DZP basis set.

Create a new VNL project, and download the NetCDF file containing the optimized structure: (AlN_orthorhombic.nc). Save it in the Project Folder and its content will become visible in the LabFloor.

  1. Go to Scripter and double click the analysis_from_file_icon Analysis from File and analysis_icon PiezoelectricTensor (click Analysis ‣ PiezoelectricTensor) blocks to add them to the script.

    snap1

  2. Double-click the analysis_from_file_icon Analysis from file block and select the object containing the optimized bulk structure (object id: gID001), as shown below.

    snap2

  3. In the analysis_icon PiezoelectricTensor block, set the k-point sampling as in the next figure and tick “Optimize strained geometries”.

    snap3

  4. Save the NetCDF file and send the script to the Job Manager to run it.

Analysis of the \({e}_{33}\) coefficient

In the LabFloor, select the PiezoelectricTensor object and click the “Text Representation” plugin in the left hand panel to inspect the piezoelectric tensor.

snap4

snap5

The predicted value for \({e}_{33}\) is 1.4849 C/m2. This value is in good agreement with the value found in the literature (1.46 C/m2). [+BFV97]

In the next section we are going to see how you can calculate \({e}_{33}\) in an alternative way.

Alternative way of calculating the piezoelectric coefficient \({e}_{33}\)

The piezoelectric coefficient \({e}_{33}\) can be also calculated as the sum of two terms [+BFV97] [+YINU12]:

  1. the clamped-ion term \({e}_{33}(0)\), that expresses the electronic response to strain.
  2. a second term describing the effect of the internal strain on the piezoelectric polarization.

The complete expression for \({e}_{33}\) thus reads:

\[{e}_{33} = {e}_{33}(0) + \frac{4 e Z^{*}}{\sqrt{3} a^{2}} \frac{du}{d{\epsilon}_{3}}\]

where \(e\) is the electronic charge, \({\epsilon}_{3}\) is the macroscopic applied strain and \(Z^{*}\) is the Born effective charge, which depends on the change in polarization upon the dispacement of an ion (or rather, a periodic sublattice of equivalent ions):

\[Z_{\nu,i,j}^{*} = \frac{\Omega}{|e|} \frac{\partial {P}_{t}^{i}}{\partial{r}_{j}^{\nu}}\]

where \(\Omega\) is the unit cell volume, \({P}_t^i\) is the total polarization along the Cartesian direction \(i\), and \({r}^\nu_j\) is the coordinate of ion \(\nu\) in direction \(j\). The Born effective charge is also referred to as effective charge or dynamical charge.

Note

The Born effective charge is a tensor. Ìn fact, when an ion is displaced in direction \(j\), this will clearly affect the polarization in same direction \(j\). However, it may also lead to a change in polarization in another direction \(i\) perpendicular to \(j\).

In the calculations below the derivative will be approximated using finite differences:

\[\displaystyle \frac{\partial P_t^z}{\partial r^\nu_z} \approx \frac{P_t^z(+\delta \hat{z},\nu)-P_t^z(-\delta \hat{z},\nu)}{2\delta},\]

where \({P}_t^z(\pm\delta \hat{z},\nu)\) is the polarization along the z-direction when atom \(\nu\) has been displaced by the amount \(\delta\) in the positive/negative z-direction.

In the following, we will obtain \({e}_{33}\) by calculating each term explicitly. Specifically, we will:

  1. calculate the variation of the interatomic distance due to strain \((\frac{du}{d{\epsilon}_{\alpha}})\);
  2. compute the piezoelectric constant in the clamped-ion approximation \({e}_{33}(0)\);
  3. compute the Born effective charge \({Z}^{*}\).

The variation of the interatomic distance due to strain

This variation is expressed as \(\frac{du}{d{\epsilon}_{\alpha}}\) where \({du}\) represents the difference in the interatomic distance \({u}\) (\({du}= {u'} - u\)) and \({d{\epsilon}_{3}}\) represents the applied strain (\({d{\epsilon}_{3}}= {c'} - c\)) in the direction parallel to the bond.

In order to calculate \({du}\) we will:

  1. relax the bulk structure and calculate \({u}\);
  2. apply a 1% strain to the unit cell in c direction, relax the internal parameters and get the interatomic distance under strain, \({u}^{'}\).

Calculating the interatomic distance \({u}\) for the relaxed system

  1. Go to Builder and click Add ‣ Form Database, locate “AlN (Hexagonal)” in the database, and add it to the Stash.

    snap6

  2. Send the bulk configuration to the Script Generator.

  3. In the Script Generator, double-click the calculator_icon New Calculator and optimization_icon OptimizeGeometry (click Optimization ‣ Optimize Geometry) blocks to add them to the script.

    snap6b

  4. In the calculator_icon New Calculator block set the following parameters:

    • exchange and correlation: GGA.PBE
    • Basis set: DZP
    • k-point sampling: (9,9,9)

    snap7

  5. In the optimization_icon OptimizeGeometry block untick “Constrain Lattice Vectors”.

    snap8

  6. Save the NetCDF file as AlN_hexagonal.nc and send the script to the Job Manager to run it.

Once the job is finshed, the optimized bulk structure should appear in the LabFloor. Select it and click the Viewer plugin to visualize it.

snap8a

snap8b

snap8c

The table below shows the relaxed structural parameters.

Table 5 Table: Relaxed cell parameters in
  ATK Bernardini et al. [+BFV97]
a 3.116 3.0766
c 5.007 4.9810
c/a 1.607 1.6190
u/c 0.381 0.380

The value of the relaxed \({c}\) lattice constant is 5.007 Å. In the next step, we will apply a 1% compressive strain in the c direction by diminishing the lattice constant in c direction in 0.05 Å.

Internal relaxation under strain in z direction

  1. Send the relaxed structure to the Builder and rename it as AlN_hexagonal_strain.nc.

  2. Go to Bulk Tools ‣ Lattice Parameters to decrease the lattice parameter in 1% in c direction.

    snap9

    snap10

  3. Send the structure to the Script Generator to create a script.

  4. In the Script Generator, double-click the calculator_icon New Calculator and optimization_icon OptimizeGeometry blocks.

  5. In the calculator_icon New Calculator, set the following parameters:

    • exchange and correlation: GGA.PBE
    • k-point sampling: (9,9,9)
  6. In the optimization_icon OptimizeGeometry block, constrain the lattice vectors as shown in the next figure.

    snap11

  7. Send the script to the Job Manager and run the script.

Once the calculation is finished, select the optimized bulk structure in the LabFloor and click the Viewer plugin to visualize the parameters.

snap12

snap13

The predicted \({u'}\) value is 1.89985 Å, and therefore the predicted \(\frac{du}{d{\epsilon}_{3}}\) value is -0.1866. This value is in agreement with the value found in the literature (-0.18) [+BFV97].

Computing the piezoelectric tensor in the clamped-ion approximation

  1. Send the relaxed structure to the Builder.

  2. In the Builder, go to Bulk Tools ‣ Supercell to transform the cell according to the figure below.

    snap14

  3. Send the bulk configuration to the Script Generator and add the calculator_icon New Calculator and the analysis_icon Analysis ‣ PiezoelectricTensor blocks to the script.

    snap15

  4. In the calculator_icon New Calculator block set the following parameters:

    • exchange and correlation: GGA.PBE
    • k-point sampling: (9,9,9)
  5. In the analysis_icon PiezoelectricTensor block set a [(19,9,9), (9,19,9), (9,9,19)] k-point sampling as shown below.

    snap16

  6. Save the NetCDF file as piezelectric_tensor_ci.nc and send the script to the Job Manager to run it.

Once the calculation is finished, the PiezoelectricTensor object will appear in the LabFloor. Use the mouse to select the object and click the “Text Representation” plugin to visualize the tensor.

snap17

snap18

The predicted \({e}_{33}(0)\) value is -0.4424 C/m2 , which is in good agreement with the value found in the literature [+BFV97].

Computing the Born effective charge

In order to compute the Born effective charge in ATK, you need to use the script provided here (born_charges.py).

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
#read saved configuration
configuration0 = nlread('piezoelectric_tensor_ci.nc', object_id='gID000')[0]

# Get the fractional coordinates of read configuration
R0 = configuration0.fractionalCoordinates()

# Get the elements
elements = configuration0.elements()

# Get the lattice
lattice = configuration0.bravaisLattice()

# From the lattice extract unit cell volume and length of lattice vector in z-direction
volume = lattice.unitCellVolume()

vectors = lattice.primitiveVectors()
c = vectors[2][2]

# Get the calculator
calculator = configuration0.calculator()

# Number of atoms
numberOfAtoms = len(R0)

# Fractional displacement in the +/- z direction
delta_z = 0.01

# Array for storing the calculated Born Charges
bornCharges = numpy.zeros(numberOfAtoms)


# Loop over atoms in unit cell
for nAtom in range(numberOfAtoms):
    # Loop over displacement in positive/negative z-direction
    
    # List with polarization values
    Pz = []
    for s in [1,-1]:
        # Make a copy of the initial coordinates
        R = R0.copy()
        
        # Displace z-coordinate if atom 'nAtom'
        R[nAtom,2] += s*delta_z
        
        # Make a new configuration with the displaced atom
        configuration = BulkConfiguration(bravais_lattice=lattice,
                                                elements=elements,
                                                fractional_coordinates=R)
        
        # Set the calculator using the saved configuration as initial state
        configuration.setCalculator(calculator,initial_state=configuration0)
        
        # Update the configuration (DFT calculation)
        configuration.update()
        
        # Calculate polarization. Only use fine k-sampling in the z-direction
        polarization = Polarization(configuration=configuration,
                                    kpoints_a=MonkhorstPackGrid(2,2,2),
                                    kpoints_b=MonkhorstPackGrid(2,2,2),
                                    kpoints_c=MonkhorstPackGrid(5,5,20),
                                    )
        
        # Print polarization
        nlprint(polarization)
        
        # Get the total cartesian polarization
        Pt = polarization.totalCartesianPolarization()
        
        # Append the z-component to the Pz list
        Pz.append(Pt[2])
    
    # Make a finite difference approximation for the derivative
    dP = (Pz[0]-Pz[1])/(2*delta_z*c)
    
    # Calculate Born charge
    born_charge = volume*dP
    
    # Add the value (in units of electron charge) to the list 'bornCharges'
    bornCharges[nAtom] = born_charge.inUnitsOf(elementary_charge) 



# Finally print out the results
print ''
print '+------------------------------+'
print '| Born effective charges (e)   |'
print '+------------------------------+'

for nAtom in range(numberOfAtoms):
    print '  %2s' %elements[nAtom].symbol() + '  :      %4.4f       ' %bornCharges[nAtom]    

print '+------------------------------+'
print '  Sum :      %4.4f       ' %numpy.sum(bornCharges)
print '+------------------------------+'
  • The script starts reading the results from the previous calculation. This will serve as a good starting guess for the DFT calculations to be performed later in the script (where the atoms are displaced).
  • Extract the fractional coordinates, the list of elements, the lattice, and the calculator from the configuration, and define a few convenient variables such as the unit cell volume and the length of the C lattice vector.
  • The parameter delta_z is the fractional displacement when calculating the derivatives.
  • There are two for loops in the script. The outer loop runs over the atoms in the unit cell, and the inner one (over the variable s) performs the displacements in the positive and negative z-direction.
  • At the end of the script, the results are printed out.

Running the script

  • Download the above born_charges.py script and save it in the Project Folder.
  • Drag the script to the Job Manager to run it.
  • When the job is finished, the calculated Born effective charges are written in the born_charges.log file.
  • As indicated by the calculated sum of the Born effective charges, the calculated values fulfill the acoustic sum rule \(\sum_s Z_{s,ij}^*=0\) with only a small error.

You have already got all the components you needed to compute to calculate \({e}_{33}\). By substituting the computed values in the formula for \({e}_{33}\), you will obtain \({e}_{33} =\) 1.464 C/m2, which is in good agreement with the value 1.4849 C/m2 calculated above and with the value found in the literature (1.464 C/m2) [+BFV97].

References

[+BFV97](1, 2, 3, 4, 5, 6) Fabio Bernardini, Vincenzo Fiorentini, and David Vanderbilt. Spontaneous polarization and piezoelectric constants of iii-v nitrides. Phys. Rev. B, 56:R10024–R10027, Oct 1997. doi:10.1103/PhysRevB.56.R10024.
[+YINU12]T. Yokoyama, Y. Iwazaki, T. Nishihara, and M. Ueda. Analysis on electromechanical coupling coefficients in aln-based bulk acoustic wave resonators based on first-principle calculations. In Ultrasonics Symposium (IUS), 2012 IEEE International, 551–554. Oct 2012. doi:10.1109/ULTSYM.2012.0137.