Mechanical / Structural Properties

The Mechanical / structural properties module provides cross-code elastic property workflows for VASP and Quantum ESPRESSO.

It is available from:

VQkit
 -> Post-processing
 -> [7] Application-oriented analysis
 -> [3] Mechanical / structural properties

The current menu contains:

[1] Prepare elastic calculations
    Choose energy-strain or stress-strain; VASP / Quantum ESPRESSO.

[2] Analyze completed elastic calculations
    Auto-detect method and obtain C_ij, S_ij and mechanical properties.

[3] Directional elasticity: high-accuracy 3D plots
    200 x 200 3D surfaces with ELATE-compatible Powell optimization
    for shear modulus and Poisson ratio.

The module supports two independent elastic-property routes:

  • energy-strain fitting;

  • stress-strain fitting.

The underlying physics layer is code-independent. VQkit applies the same Voigt-strain description to either VASP or Quantum ESPRESSO structures.

Elastic tensor

The elastic stiffness tensor is written in Voigt notation as:

\[\boldsymbol{\sigma} = \mathbf{C}\boldsymbol{\epsilon}\]

where:

  • \(\boldsymbol{\sigma}\) is the six-component stress vector;

  • \(\boldsymbol{\epsilon}\) is the six-component engineering-strain vector;

  • \(\mathbf{C}\) is the \(6\times6\) stiffness tensor.

The compliance tensor is:

\[\mathbf{S} = \mathbf{C}^{-1}\]

VQkit reports:

C_ij in GPa
S_ij in GPa^-1

Voigt strain convention

The six strain components are interpreted as:

epsilon_1 = epsilon_xx
epsilon_2 = epsilon_yy
epsilon_3 = epsilon_zz
epsilon_4 = gamma_yz
epsilon_5 = gamma_xz
epsilon_6 = gamma_xy

The strain tensor is therefore:

\[\begin{split}\boldsymbol{\epsilon} = \begin{bmatrix} \epsilon_1 & \epsilon_6/2 & \epsilon_5/2 \\ \epsilon_6/2 & \epsilon_2 & \epsilon_4/2 \\ \epsilon_5/2 & \epsilon_4/2 & \epsilon_3 \end{bmatrix}\end{split}\]

The strained lattice is generated from:

\[\mathbf{A}' = \mathbf{A} \left(\mathbf{I}+\boldsymbol{\epsilon}\right)^T\]

where \(\mathbf{A}\) is the original lattice matrix.

Default strain grid

The current default strain values are:

-0.015
-0.010
-0.005
 0.000
 0.005
 0.010
 0.015

Users may enter custom strain values.

The strain set must contain:

  • at least three points;

  • zero strain.

[1] Prepare elastic calculations

Choose:

[1] Prepare elastic calculations

VQkit first asks for the elastic method:

[1] Energy-strain
    Fits curvature of total energy versus strain

[2] Stress-strain
    Fits sigma_i = sum_j C_ij epsilon_j

Then select:

[1] VASP
[2] Quantum ESPRESSO

Source structures

For VASP, VQkit searches for:

POSCAR
CONTCAR

For Quantum ESPRESSO, VQkit searches for:

*.in

The selected source structure must contain explicit lattice vectors.

Metric lattice-family detection

VQkit calculates:

  • a, b, c;

  • alpha, beta, gamma;

  • cell volume.

It then performs a metric-lattice classification.

Possible detected families include:

cubic
tetragonal
orthorhombic
hexagonal
trigonal
monoclinic
triclinic

This is not a space-group determination.

The detector uses only lattice lengths and angles.

For example:

  • approximately equal a=b=c and 90-degree angles -> cubic;

  • a=b != c with 90-degree angles -> tetragonal;

  • three unequal lengths with 90-degree angles -> orthorhombic;

  • a=b with 120/60-degree gamma -> hexagonal.

General elastic mode

VQkit provides a:

General / asymmetric

elastic mode.

This mode fits the complete symmetric:

\[6\times6\]

stiffness tensor with:

\[21\]

independent elastic constants.

No crystal-symmetry equalities are imposed.

This is the recommended treatment for:

  • doped systems;

  • vacancies;

  • defects;

  • substituted structures;

  • reconstructed cells;

  • chemically asymmetric cells;

  • disordered structures;

  • systems whose atomic symmetry is lower than their cell shape suggests.

This is important because a metrically cubic supercell is not necessarily elastically cubic after doping or defect formation.

Symmetry-reduced mode

VQkit also supports:

Metric-symmetry reduced

analysis for pristine, standardized high-symmetry crystals.

The nominal independent elastic-constant counts are:

cubic        : 3
hexagonal    : 5
tetragonal   : 6
orthorhombic : 9
trigonal     : 6
monoclinic   : 13
triclinic    : 21
general      : 21

The user may also choose a crystal family manually.

A manual or metric-reduced choice should only be used when the assumed elastic symmetry is scientifically justified.

Energy-strain method

For the energy-strain route, VQkit generates strain modes according to the chosen elastic family.

For cubic materials, the reduced set includes patterns associated with:

C11
C11_C12
C44

For hexagonal materials, representative modes include:

C11
C33
C44
C11_C12
C11_C33_C13

For tetragonal materials:

C11
C33
C44
C66
C11_C12
C11_C33_C13

For orthorhombic materials, nine independent mode combinations are generated.

For general low-symmetry analysis, VQkit generates the full 21-component basis.

Energy-strain physics

For small strain:

\[E(\epsilon) = E_0 + V_0 \frac{1}{2} \boldsymbol{\epsilon}^T \mathbf{C} \boldsymbol{\epsilon}\]

The numerical energy data are fitted quadratically.

The resulting curvature is converted to GPa using:

\[1\ \mathrm{eV/\AA^3} = 160.21766208\ \mathrm{GPa}\]

VASP preparation

For VASP, VQkit creates a calculation directory tree and writes strained POSCAR files.

It also creates an elastic INCAR profile with high-precision settings including the current overrides:

EDIFF  = 1E-08
IBRION = 2
NSW    = 100
ISIF   = 2
EDIFFG = -0.01
LREAL  = .FALSE.
PREC   = Accurate
LWAVE  = .FALSE.
LCHARG = .FALSE.
ADDGRID = .TRUE.

Where available, VQkit copies:

KPOINTS
POTCAR
job.sh

from the source directory into every generated calculation folder.

Quantum ESPRESSO preparation

For QE, VQkit inherits the numerical and pseudopotential settings from the selected pw.x input.

It updates the input to use:

calculation = 'relax'

and writes the strained cell explicitly using:

CELL_PARAMETERS angstrom

with atomic positions written as:

ATOMIC_POSITIONS crystal

For stress-strain calculations, QE stress output is enabled with:

tstress = .true.

The user should review inherited QE settings before running the calculations.

Output directories

The default energy-strain workflow directory is:

VQkit_elastic

The default stress-strain workflow directory is:

VQkit_elastic_stress

Inside the workflow, each strain mode contains subdirectories for the strain values.

A typical structure is:

VQkit_elastic/
├── C11/
│   ├── strain_m0p0150/
│   ├── strain_m0p0100/
│   ├── strain_m0p0050/
│   ├── strain_z0p0000/
│   ├── strain_p0p0050/
│   ├── strain_p0p0100/
│   └── strain_p0p0150/
└── ...

Manifest file

Every prepared elastic workflow contains:

VQkit_ELASTIC_MANIFEST.json

The manifest records information such as:

  • method;

  • code;

  • source structure;

  • detected metric lattice family;

  • chosen analysis family;

  • symmetry treatment;

  • nominal independent constants;

  • reference volume;

  • strain values;

  • strain modes;

  • calculation directories.

This manifest allows VQkit to auto-detect the correct analysis method later.

Stress-strain method

The stress-strain workflow uses six independent applied engineering-strain directions:

eps1_xx
eps2_yy
eps3_zz
eps4_yz
eps5_xz
eps6_xy

For each applied strain direction, VQkit later fits the full stress response.

This follows generalized Hooke’s law:

\[\sigma_i = \sum_j C_{ij}\epsilon_j\]

For a general tensor, six applied strain modes are sufficient because each calculation provides all six stress components.

This can require substantially fewer calculations than a 21-mode general energy-strain workflow.

Stress-strain accuracy

Stress-derived elastic constants are sensitive to:

  • plane-wave cutoff;

  • k-point density;

  • Pulay stress;

  • SCF convergence;

  • residual stress;

  • ionic relaxation quality.

VQkit explicitly warns users to verify stress convergence.

[2] Analyze completed elastic calculations

Choose:

[2] Analyze completed elastic calculations

VQkit asks for the workflow directory, for example:

VQkit_elastic

The method is read from:

VQkit_ELASTIC_MANIFEST.json

and analysis is automatically routed to:

  • energy-strain;

  • stress-strain.

Energy-strain fitting

VQkit reads the completed total energies from each strain directory.

The quadratic energy fit is stored with:

  • quadratic coefficient;

  • linear coefficient;

  • constant coefficient;

  • fitted elastic quantity in GPa;

  • \(R^2\);

  • number of data points.

The fit table is written to:

ENERGY_STRAIN_FITS.dat

Stress-strain fitting

For the stress-strain route, each stress component is fitted linearly against the corresponding applied strain.

The fit relation is:

\[\sigma_i = C_{ij}\epsilon_j + b_i\]

For each fit, VQkit records:

  • strain mode;

  • applied component;

  • response component;

  • slope in GPa;

  • intercept in GPa;

  • \(R^2\);

  • number of points.

The fit table is written to:

STRESS_STRAIN_FITS.dat

Raw stress stiffness tensor

For stress-strain analysis, the direct fitted tensor is written as:

RAW_STRESS_STIFFNESS_TENSOR.dat

Because numerical fitting can produce small reciprocity differences:

\[C_{ij} \ne C_{ji}\]

VQkit reports the maximum raw reciprocity mismatch.

The final reported stiffness tensor is reciprocity-symmetrized before optional crystal-family projection.

Elastic stiffness tensor

The final stiffness tensor is written to:

ELASTIC_TENSOR.dat

in GPa.

It has the form:

\[\begin{split}\mathbf{C} = \begin{bmatrix} C_{11}&C_{12}&C_{13}&C_{14}&C_{15}&C_{16}\\ C_{12}&C_{22}&C_{23}&C_{24}&C_{25}&C_{26}\\ C_{13}&C_{23}&C_{33}&C_{34}&C_{35}&C_{36}\\ C_{14}&C_{24}&C_{34}&C_{44}&C_{45}&C_{46}\\ C_{15}&C_{25}&C_{35}&C_{45}&C_{55}&C_{56}\\ C_{16}&C_{26}&C_{36}&C_{46}&C_{56}&C_{66} \end{bmatrix}\end{split}\]

Compliance tensor

The compliance tensor is:

\[\mathbf{S} = \mathbf{C}^{-1}\]

and is written to:

COMPLIANCE_TENSOR.dat

with units:

GPa^-1

Voigt-Reuss-Hill properties

VQkit calculates polycrystalline elastic properties using Voigt, Reuss, and Hill averaging.

The output contains:

B_V
B_R
B_H

for bulk modulus and:

G_V
G_R
G_H

for shear modulus.

The Hill averages are:

\[B_H = \frac{B_V+B_R}{2}\]
\[G_H = \frac{G_V+G_R}{2}\]

Young’s modulus

The Hill Young’s modulus is:

\[E_H = \frac{9B_HG_H} {3B_H+G_H}\]

Poisson ratio

The Hill Poisson ratio is:

\[\nu_H = \frac{3B_H-2G_H} {2(3B_H+G_H)}\]

Pugh ratio

VQkit reports:

\[\frac{B_H}{G_H}\]

as the Pugh ratio.

This quantity is commonly used as an empirical ductility/brittleness descriptor, but it should not be treated as an absolute classification for every material.

Elastic anisotropy

The mechanical-property summary reports:

  • universal anisotropy index \(A_U\);

  • shear anisotropy \(A_G\);

  • bulk anisotropy \(A_B\).

A fully isotropic solid approaches zero anisotropy for the corresponding measures.

Mechanical stability

VQkit evaluates the eigenvalues of the symmetric stiffness tensor.

A mechanically stable positive-definite tensor requires:

\[\lambda_i > 0\]

for all stiffness-matrix eigenvalues.

The report therefore gives:

STABLE

or:

UNSTABLE

based on positive definiteness.

For cubic systems, VQkit also reports the corresponding cubic stability criteria when applicable.

Mechanical output summary

The main analysis writes:

ELASTIC_TENSOR.dat
COMPLIANCE_TENSOR.dat
MECHANICAL_PROPERTIES.txt

and one of:

ENERGY_STRAIN_FITS.dat

or:

STRESS_STRAIN_FITS.dat

For stress-strain analysis, it also writes:

RAW_STRESS_STIFFNESS_TENSOR.dat

The MECHANICAL_PROPERTIES.txt summary records:

  • method;

  • code;

  • detected metric family;

  • elastic-analysis family;

  • symmetry treatment;

  • fitted component count;

  • reference volume;

  • Voigt/Reuss/Hill moduli;

  • Young’s modulus;

  • Poisson ratio;

  • Pugh ratio;

  • anisotropy measures;

  • stiffness eigenvalues;

  • stability result;

  • missing-calculation count where applicable.

[3] Directional elasticity: high-accuracy 3D plots

Choose:

[3] Directional elasticity: high-accuracy 3D plots

This workflow analyzes elastic anisotropy directly from a completed \(6\times6\) stiffness tensor.

Input sources

VQkit offers:

[1] Read elastic tensor from file
[2] Read completed VQkit elastic workflow
[3] Enter 6 x 6 elastic tensor manually

For file input, conventional filenames include:

ELASTIC_TENSOR
ELASTIC_TENSOR.dat
ELASTIC_TENSOR.in

The input tensor must be symmetric within numerical tolerance.

If the maximum:

\[|C_{ij}-C_{ji}|\]

exceeds approximately:

1e-3 GPa

the file is rejected.

Positive-definite requirement

Directional analysis requires a mechanically stable positive-definite stiffness tensor.

VQkit checks the stiffness eigenvalues before plotting.

If the minimum eigenvalue is non-positive, directional analysis is stopped.

Young’s modulus

For loading direction \(\mathbf{a}\), the directional Young’s modulus is:

\[E(\mathbf{a}) = \left[ S_{ijkl} a_i a_j a_k a_l \right]^{-1}\]

where \(S_{ijkl}\) is the fourth-order compliance tensor.

Linear compressibility

Directional linear compressibility is evaluated from:

\[\beta(\mathbf{a}) = S_{ijkk}a_i a_j\]

The 3D visualization uses the same compliance tensor.

Shear modulus

For loading direction \(\mathbf{a}\) and a perpendicular direction \(\mathbf{b}\):

\[G(\mathbf{a},\mathbf{b}) = \left[ 4S_{ijkl} a_i b_j a_k b_l \right]^{-1}\]

For each loading direction, VQkit searches over perpendicular orientations to obtain:

  • minimum shear modulus;

  • maximum shear modulus.

Poisson ratio

The directional Poisson ratio is:

\[\nu(\mathbf{a},\mathbf{b}) = - \frac{ S_{ijkl}b_i b_j a_k a_l }{ S_{ijkl}a_i a_j a_k a_l }\]

VQkit obtains directional:

  • minimum Poisson ratio;

  • maximum Poisson ratio.

A negative minimum Poisson ratio indicates an auxetic direction.

High-resolution grid

The default 3D surface resolution is:

200 x 200

in spherical angular space.

Young’s modulus and linear compressibility are evaluated directly on this grid.

ELATE-compatible optimization

For shear modulus and Poisson-ratio surfaces, VQkit uses an ELATE-compatible Powell optimization over the transverse orientation angle.

This reproduces the directional-envelope convention used by ELATE-style elastic visualization.

The code also contains an analytic transverse-plane eigensolver used for cross-checking and directional-extrema searches.

Automatic directional extrema

Before plotting, VQkit searches automatically for:

E_min
E_max
beta_min
beta_max
G_min
G_max
nu_min
nu_max

For properties involving a transverse vector, the report also includes the corresponding:

a direction
b direction

and verifies their perpendicularity.

Directional output directory

The default output directory is:

VQkit_directional_elasticity

VQkit writes:

VQkit_DIRECTIONAL_ELASTICITY_summary.txt

and high-resolution 3D surfaces for:

Young's modulus
Linear compressibility
Minimum shear modulus
Maximum shear modulus
Minimum Poisson ratio
Maximum Poisson ratio

The Young’s-modulus output stem is:

VQkit_Young_3D

The linear-compressibility output stem is:

VQkit_Linear_Compressibility_3D

The shear outputs include:

VQkit_Shear_Minimum_3D
VQkit_Shear_Maximum_3D

and the Poisson outputs follow the corresponding minimum/maximum naming convention.

Each surface is saved as:

.png
.pdf
.csv

The CSV files contain Cartesian surface coordinates:

x
y
z

Directional summary

The summary records:

  • minimum and maximum Young’s modulus;

  • minimum and maximum linear compressibility;

  • minimum and maximum shear modulus;

  • minimum and maximum Poisson ratio;

  • associated directional vectors;

  • whether an auxetic direction is detected.

For example:

Auxetic direction detected: YES

is written if:

\[\nu_{\min} < 0\]

Energy-strain versus stress-strain

Energy-strain

Advantages:

  • based on total-energy curvature;

  • robust when highly converged total energies are available.

Considerations:

  • general 21-component analysis can require many calculations;

  • quadratic fits must be well behaved;

  • strain range must remain in the elastic regime.

Stress-strain

Advantages:

  • only six strain directions are required for a fully general tensor;

  • each calculation provides all six stress responses.

Considerations:

  • stress is especially sensitive to numerical convergence;

  • Pulay stress can affect the tensor;

  • residual stress must be controlled;

  • accurate stress output is essential.

Scientific cautions

Small-strain regime

The elastic tensor is a linear-response quantity.

Strains should be small enough that the material remains in the harmonic elastic regime.

Too-large strains can introduce anharmonic effects.

Too-small strains can make the signal comparable to numerical noise.

The default ±1.5% range is a practical starting point, not a universal guarantee.

Internal relaxation

The current generated VASP/QE elastic workflows allow ionic relaxation at fixed strained cell.

This corresponds to relaxed-ion elastic response.

If clamped-ion elastic constants are required, the computational setup must be modified accordingly.

Symmetry

Metric symmetry is not equivalent to atomic or chemical symmetry.

For doped, defected, or reconstructed systems, use the general 21-component mode unless stronger symmetry is independently established.

Convergence

Elastic constants can be more sensitive than ordinary total energies.

Converge:

  • plane-wave cutoff;

  • k-point mesh;

  • electronic convergence;

  • force threshold;

  • stress accuracy.

Stability

A positive-definite stiffness matrix is a necessary condition for the directional elastic analysis implemented here.

An unstable tensor should not be used to generate physically interpreted directional moduli.

Summary

The VQkit Mechanical / structural properties module provides:

Preparation
-> VASP / QE
-> energy-strain or stress-strain
-> general or symmetry-reduced tensor

Analysis
-> C_ij
-> S_ij
-> Voigt/Reuss/Hill moduli
-> Young's modulus
-> Poisson ratio
-> Pugh ratio
-> anisotropy indices
-> mechanical stability

Directional elasticity
-> E
-> linear compressibility
-> G_min / G_max
-> nu_min / nu_max
-> 200 x 200 high-resolution 3D surfaces

The workflow is designed to support both high-symmetry crystals and fully asymmetric doped/defective structures without silently imposing unjustified elastic symmetry.