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.
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 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.
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.
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.
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 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
=============================================== 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 (GPa) 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.