Elastic constant

import mdapy as mp
mp.__version__

'1.0.3'

MD to calculate elastic constant

al = mp.build_crystal('Al', 'fcc', 4.05)
calc = mp.NEP("../../tests/input_files/UNEP-v1.txt")

elas = mp.get_elastic_constant(al, calc)

elas.print()

Elastic tensor (GPa):
  121.27     54.01     54.01      0.00      0.00      0.00
   54.01    121.27     54.01      0.00      0.00      0.00
   54.01     54.01    121.27      0.00      0.00      0.00
    0.00      0.00      0.00     42.89      0.00      0.00
    0.00      0.00      0.00      0.00     42.89      0.00
    0.00      0.00      0.00      0.00      0.00     42.89

elas.print_vrh()

Voigt-Reuss-Hill averages:
  Property                              Voigt     Reuss      Hill  Unit
  ------------------------------------------------------------------
  Bulk modulus K                        76.43     76.43     76.43  GPa
  Shear modulus G                       39.19     38.64     38.91  GPa
  Young's modulus E                                         99.80  GPa
  Poisson's ratio nu                                       0.2824  -
  Voigt = upper bound  |  Reuss = lower bound  |  Hill = best estimate
  Voigt-Reuss gap (K+G): 0.55 GPa  (larger = more anisotropic)

elas.voigt.round(2)

array([[121.27,  54.01,  54.01,   0.  ,   0.  ,   0.  ],
       [ 54.01, 121.27,  54.01,   0.  ,   0.  ,   0.  ],
       [ 54.01,  54.01, 121.27,   0.  ,   0.  ,   0.  ],
       [  0.  ,   0.  ,   0.  ,  42.89,   0.  ,   0.  ],
       [  0.  ,   0.  ,   0.  ,   0.  ,  42.89,   0.  ],
       [  0.  ,   0.  ,   0.  ,   0.  ,   0.  ,  42.89]])

DFT to calculate elastic constant

• Perform full DFT optimization of the structure, including both atomic positions and cell parameters.

• Generate a set of deformed structures based on the optimized configuration.

• Compute the stress for each deformed structure using DFT (with energy minimization).

from mdapy.elastic import DeformedStructureSet, ElasticTensor, strain_from_deformation

# Load the DFT-optimized structure.
# Both cell and energy should be fully relaxed beforehand.
system = mp.System('opt.POSCAR')

# Generate a set of deformed structures.
# Note: mdapy does not account for symmetry, so even a simple FCC structure
# will produce 24 deformed configurations by default.
dfm_ss = DeformedStructureSet(system)

# Save each deformed structure.
# These structures can then be used for subsequent DFT calculations to obtain stress.
# Here, only energy minimization is considered as an example.
for i, j in enumerate(dfm_ss.deformed_systems):
    j.write_poscar(f'{i}.POSCAR')

# After completing the DFT calculations, collect the stress for each deformed structure.
strain_list, stress_list = [], []
for defo in dfm_ss.deformations:
    stress_list.append(stress)  # Stress corresponding to each deformed structure
    strain_list.append(strain_from_deformation(defo))

# Fit the elastic tensor using the collected strain–stress data.
elas = ElasticTensor.from_independent_strains(
    strain_list,
    stress_list,
    eq_stress=equi_stress  # Stress of the initial optimized (equilibrium) structure
)