CASTEP tutorials > Calculating elastic constants for BN

Calculating elastic constants for BN

Purpose: Illustrates the use of CASTEP to calculate elastic constants.

Modules: Materials Visualizer, CASTEP


Prerequisites: Predicting the lattice parameters of AlAs from first principles


Recent developments in density functional theory (DFT) methods applicable to studies of large periodic systems have become essential in addressing problems in materials design and processing. The DFT tools can be used to guide and lead the design of new materials, allowing researchers to understand the underlying chemistry and physics of processes.


In this tutorial, you will learn how to use CASTEP to calculate elastic constants and other mechanical properties. In the first part you will optimize the structure of cubic BN and then you will calculate its elastic constants.

This tutorial covers:

In order to ensure that you can follow this tutorial exactly as intended, you should use the Settings Organizer dialog to ensure that all your project settings are set to their Accelrys default values. See the Creating a project tutorial for instructions on how to restore default project settings.

1. Getting started

Begin by starting Materials Studio and creating a new project.

Open the New Project dialog and enter BN_elastic as the project name, click the OK button.

The new project is created with BN_elastic listed in the Project Explorer. The next step is to import the BN structure.

Select File | Import... from the menu bar to open the Import Document dialog. Navigate to the folder Structures/semiconductors/ and select BN.msi. Click the Open button.

The crystal structure of BN is displayed.

Structure of cubic BN
Structure of BN cubic

This is the conventional representation of the BN structure. In order to reduce the computation time, you should convert to the primitive representation.

Select Build | Symmetry | Primitive Cell from the menu bar.

2. To optimize the structure of BN cubic

It is not necessary to perform geometry optimization before calculating elastic constants, so you can generate Cij data for experimentally observed structures. However, more consistent results are obtained if you perform full geometry optimization, including cell optimization, and then calculate the elastic constants for the structure corresponding to the theoretical ground state.

The accuracy of the elastic constants, especially of the shear constants, depends strongly on the quality of the SCF calculation and, in particular, on the quality of the Brillouin zone sampling and the degree of convergence of wavefunctions. Therefore, you should use the Fine setting for SCF tolerance and k-point sampling and a Fine derived FFT grid.

Now you will set up the geometry optimization.

Click the CASTEP button CASTEP on the Modules toolbar and select Calculation from the dropdown list or choose Modules | CASTEP | Calculation from the menu bar.

This opens the CASTEP Calculation dialog.

CASTEP Calculation dialog, Setup tab
CASTEP Calculation dialog, Setup tab

On the Setup tab, set the Task to Geometry Optimization, the Quality to Fine, and the Functional to GGA and PW91.

Click the More... button to open the CASTEP Geometry Optimization dialog. Check the Optimize cell checkbox and close the dialog.

Choose the Job Control tab on the CASTEP Calculation dialog and select the Gateway on which you wish to run the CASTEP job.

Click the Run button.

After optimization, the structure should have cell parameters of about a = b = c = 2.574 Å.

Right-click in the 3D Viewer and select Lattice Parameters from the shortcut menu.

The lattice parameters are shown. Now you can go on to calculate the elastic constants of the optimized structure.

3. To calculate the elastic constants of BN

Select the Setup tab on the CASTEP Calculation dialog. Select Elastic Constants from the Task dropdown list and click the More... button.

This opens the CASTEP Elastic Constants dialog.

CASTEP Elastic Constants dialog
CASTEP Elastic Constants dialog, Elastic Constants tab

Increase the Number of steps for each strain from 4 to 6 and close the dialog. Ensure that BN CASTEP GeomOpt/BN.xsd is the active document and click the Run button on the CASTEP Calculation dialog.

If BN CASTEP GeomOpt/BN.xsd is made active before the CASTEP Elastic Constants dialog is opened, the Strain pattern grid on this dialog will contain values.

CASTEP results for the Elastic Constants task are returned as a set of .castep output files. Each of them represents a geometry optimization run with a fixed cell, for a given strain pattern and strain amplitude. The naming convention for these files is:


where m is the current strain pattern and n is the current strain amplitude for the given pattern.

CASTEP can use these results to analyze the calculated stress tensors for each of these runs and generate a file with information about elastic properties.

Click the CASTEP button CASTEP on the Modules toolbar and select Analysis from the dropdown list or choose Modules | CASTEP | Analysis from the menu bar.

Select the Elastic constants option. The results file from the Elastic Constants job for BN should be displayed automatically in the Results file selector. Click the Calculate button.

A new text document, BN Elastic Constants.txt, is created in the results folder.

The information in this document includes a summary of the input strains and calculated stresses, results of linear fitting for each strain pattern (including quality of the fit), the correspondence between calculated stresses and elastic constants for a given symmetry, a table of elastic constants (Cij) and elastic compliances (Sij). The derived properties, such as bulk modulus and its inverse, compressibility, Young modulus, and Poisson ratios for three directions, and the Lame constants that are needed for modeling the material as an isotropic medium, are also reported in this document.

3. Description of the elastic constants file

The results you obtain may vary slightly from those shown because of minor differences in the structure of the starting model.

Two strain patterns are required for this lattice type. For each strain pattern, there is a summary of calculated stresses as extracted from the respective .castep files:

  Elastic constants from Materials Studio: CASTEP

         Summary of the calculated stresses

  Strain pattern:      1

  Current amplitude: 1
  Transformed stress tensor (GPa) :
      -12.255858       -0.000000       -0.000000
       -0.000000      -14.094398        1.282706
       -0.000000        1.282706      -14.094398

  Current amplitude: 2
  Transformed stress tensor (GPa) :
      -13.154269       -0.000000       -0.000000
       -0.000000      -14.247919        0.764497
       -0.000000        0.764497      -14.247919

  Current amplitude: 3
  Transformed stress tensor (GPa) :
      -14.042415       -0.000000       -0.000000
       -0.000000      -14.414859        0.246872
       -0.000000        0.246872      -14.414859

  Current amplitude: 4
  Transformed stress tensor (GPa) :
      -14.910298       -0.000000       -0.000000
       -0.000000      -14.542428       -0.253306
       -0.000000       -0.253306      -14.542428

  Current amplitude: 5
  Transformed stress tensor (GPa) :
      -15.792607       -0.000000       -0.000000
       -0.000000      -14.699735       -0.769431
       -0.000000       -0.769431      -14.699735

  Current amplitude: 6
  Transformed stress tensor (GPa) :
      -16.686454       -0.000000       -0.000000
       -0.000000      -14.849504       -1.266391
       -0.000000       -1.266391      -14.849504

Any information about the connection between components of the stress, strain, and elastic constants tensors is provided. At this stage, each elastic constant is represented by a single compact index rather than by a pair of ij indices. The correspondence between the compact notation and the conventional indexing is provided later in the file:

Stress corresponds to elastic coefficients (compact notation):
   1  7  7  4  0  0

 as induced by the strain components:
   1  1  1  4  0  0

A linear fit of the stress-strain relationship for each component of the stress is given in the following format:

     Stress    Cij       value of          value of
     index    index       stress            strain
       1        1       -12.255858         -0.003000
       1        1       -13.154269         -0.001800
       1        1       -14.042415         -0.000600
       1        1       -14.910298          0.000600
       1        1       -15.792607          0.001800
       1        1       -16.686454          0.003000
C (gradient)     :     736.568500
Error on C       :       1.743124
Correlation coeff:       0.999989
Stress intercept :     -14.473650

       2        7       -14.094398         -0.003000
       2        7       -14.247919         -0.001800
       2        7       -14.414859         -0.000600
       2        7       -14.542428          0.000600
       2        7       -14.699735          0.001800
       2        7       -14.849504          0.003000
C (gradient)     :     125.203500
Error on C       :       1.761830
Correlation coeff:       0.999604
Stress intercept :     -14.474807

       3        7       -14.094398         -0.003000
       3        7       -14.247919         -0.001800
       3        7       -14.414859         -0.000600
       3        7       -14.542428          0.000600
       3        7       -14.699735          0.001800
       3        7       -14.849504          0.003000
C (gradient)     :     125.203500
Error on C       :       1.761830
Correlation coeff:       0.999604
Stress intercept :     -14.474807

       4        4         1.282706         -0.003000
       4        4         0.764497         -0.001800
       4        4         0.246872         -0.000600
       4        4        -0.253306          0.000600
       4        4        -0.769431          0.001800
       4        4        -1.266391          0.003000
C (gradient)     :     424.939214
Error on C       :       1.471105
Correlation coeff:       0.999976
Stress intercept :       0.000824

The gradient provides the value of the elastic constant (or a linear combination of elastic constants); the quality of the fit, indicated by the correlation coefficient, provides the statistical uncertainty of that value. The stress intercept value is not used in further analysis, it is simply an indication of how far the converged ground state was from the initial structure.

The results for all the strain patterns are then summarized:

  Summary of elastic constants

 id  i  j       Cij (GPa)
 1   1  1     736.56850 +/-   1.743
 4   4  4     424.93921 +/-   1.471
 7   1  2     125.20350 +/-   1.246

The errors are only provided when more than two values for the strain amplitude are used, since there is no statistical uncertainty associated with fitting a straight line to only two points.

Elastic constants are then presented in a conventional 6 × 6 tensor form, followed by a similar 6 x 6 representation of the compliances:

      Elastic Stiffness Constants Cij (GPa)

   736.56850   125.20350   125.20350     0.00000     0.00000     0.00000
   125.20350   736.56850   125.20350     0.00000     0.00000     0.00000
   125.20350   125.20350   736.56850     0.00000     0.00000     0.00000
     0.00000     0.00000     0.00000   424.93921     0.00000     0.00000
     0.00000     0.00000     0.00000     0.00000   424.93921     0.00000
     0.00000     0.00000     0.00000     0.00000     0.00000   424.93921

      Elastic Compliance Constants Sij (1/GPa)

   0.0014282  -0.0002075  -0.0002075   0.0000000   0.0000000   0.0000000
  -0.0002075   0.0014282  -0.0002075   0.0000000   0.0000000   0.0000000
  -0.0002075  -0.0002075   0.0014282   0.0000000   0.0000000   0.0000000
   0.0000000   0.0000000   0.0000000   0.0023533   0.0000000   0.0000000
   0.0000000   0.0000000   0.0000000   0.0000000   0.0023533   0.0000000
   0.0000000   0.0000000   0.0000000   0.0000000   0.0000000   0.0023533

The final part of the file contains the derived properties:

Bulk modulus    =   328.99183 +/-  1.014 (GPa)

Compressibility =     0.00304 (1/GPa)

  Axis   Young Modulus<poisson Ratios
   X      700.18784       Exy=  0.1453 Exz=  0.1453
   Y      700.18784       Eyx=  0.1453 Eyz=  0.1453
   Z      700.18784       Ezx=  0.1453 Ezy=  0.1453

      Elastic constants for polycrystalline material (GPa)

                                 Voigt       Reuss<hill
Bulk modulus              :     328.99183   328.99183   328.99183
Shear modulus (Lame Mu)   :     377.23653   367.57761   372.40707
Lame lambda               :      77.50081    83.94009    80.72045

Universal anisotropy index:       0.13121

Calculated elastic properties of crystals are significantly more sensitive to the accuracy of the electronic structure calculation than, for example, calculated lattice parameters and atomic coordinates. It is always necessary to check the convergence of the calculated properties with respect to the following parameters:

An additional consideration in elastic properties calculations is the choice of exchange-correlation potential. Modern potentials that are designed to reproduce solid state properties more accurately than traditional LDA or PBE functionals are PBESOL and Wu-Cohen; these are recommended for calculating elastic coefficients of solids.

This is the end of the tutorial.

Accelrys Materials Studio 8.0 Help: Wednesday, December 17, 2014
Legal Notices