Large value of elastic constant c33

phialberta · June 21, 2023, 2:49am

Dear Abinit community,
I want to determine the elastic and piezoelectric properties of Wurtzite-AlN . Both atomic positions and the lattice constants have been optimized before beginning DFPT calculations.
#AlN in hypothetical wurzite (hexagonal) structure
#Response function calculation for:

* rigid-atom elastic tensor

* rigid-atom piezoelectric tensor

* interatomic force constants at gamma

* Born effective charges

ndtset 3

Set 1 : Initial self-consistent run

kptopt1 1
tolvrs1 1.0d-18 #need excellent convergence of GS quantities for RF runs

Set 2 : Calculate the ddk wf’s - needed for piezoelectric tensor and

Born effective charges in dataset 3

getwfk2 -1
iscf2 -3 #this option is needed for ddk
kptopt2 2 #use time-reversal symmetry only for k points
nqpt2 1 #one wave vector will be specified
qpt2 0 0 0 #need to specify gamma point
rfelfd2 2 #set for ddk wf’s only
tolwfr2 1.0d-20 #only wf convergence can be monitored here

Set 3 : response-function calculations for all needed perturbations

getddk3 -1
getwfk3 -2
kptopt3 2 #use time-reversal symmetry only for k points
nqpt3 1
qpt3 0 0 0
rfphon3 1 #do atomic displacement perturbation
rfstrs3 3 #do strain perturbation
tolvrs3 1.0d-10 #need reasonable convergence of 1st-order quantities

#Common input data

acell COPY RELAXED RESULT FROM PREVIOUS CALCULATION

Here is a default value, for automatic testing : suppress it and fill with values from the previous run

acell 5.8963399059E+00 5.8963399059E+00 9.5183958753E+00 Bohr

rprim sqrt(0.75) 0.5 0.0 #Better to specify hexagonal primitive vectors
-sqrt(0.75) 0.5 0.0 #with high accuracy to be
0.0 0.0 1.0 #sure that the symmetry is recognized
#and preserved in the optimization process

#Definition of the atom types and atoms
ntypat 2
znucl 13 7
natom 4
typat 1 1 2 2

xred COPY RELAXED RESULT FROM PREVIOUS CALCULATION

Here is a set of default values, for automatic testing : suppress it and fill with values from the previous run

xred 3.3333333333E-01 6.6666666667E-01 0.0000000000E+00
6.6666666667E-01 3.3333333333E-01 5.0000000000E-01
3.3333333333E-01 6.6666666667E-01 3.8056006731E-01
6.6666666667E-01 3.3333333333E-01 8.8056006731E-01

#Gives the number of bands, explicitely (do not take the default)
nband 8 # For an insulator (if described correctly as an
# insulator by DFT), conduction bands should not
# be included in response-function calculations

#Definition of the plane wave basis set
ecut 20 # Maximum kinetic energy cutoff (Hartree)
ecutsm 0.5 # Smoothing energy needed for lattice paramete
# optimization. This will be retained for
# consistency throughout.

#Definition of the k-point grid
kptopt 1 # Use symmetry and treat only inequivalent points
ngkpt 4 4 4 # 4x4x4 Monkhorst-Pack grid
nshiftk 1 # Use one copy of grid only (default)
shiftk 0.0 0.0 0.5 # This choice of origin for the k point grid
# preserves the hexagonal symmetry of the grid,
# which would be broken by the default choice.

#Definition of the self-consistency procedure
diemac 9.0 # Model dielectric preconditioner
nstep 40 # Maxiumum number of SCF iterations

enforce calculation of forces at each SCF step

optforces 1
I checked all scf are convergenced . the stresses are relaxed. I also confirm that the convergence in term of WFK had achieved.
I get a large value of c33 elastic constant :
Elastic Tensor (clamped ion) (unit:10^2GP):

5.3966224 1.2721801 2.4291300 0.0000056 0.0000554 0.0006669
1.2721801 5.3977411 2.4292820 0.0000192 0.0000184 0.0004970
2.4292622 2.4291530 13.2992211 0.0000026 0.0000139 0.0000556
0.0000007 0.0001086 0.0000183 3.0872582 0.0000238 0.0000000
0.0000000 0.0000000 0.0000000 0.0000238 3.0871503 0.0000008
0.0002840 0.0000908 0.0000237 0.0000000 0.0000008 2.0617493

Elastic Tensor (relaxed ion) (unit:10^2GP):
(at fixed electric field boundary condition)

4.7008113 1.5952688 2.8114591 0.0000055 0.0000553 0.0006680
1.5952689 4.7017499 2.8116392 0.0000193 0.0000185 0.0004956
2.8115913 2.8115101 12.5149861 0.0000026 0.0000139 0.0000551
0.0000006 0.0001088 0.0000183 2.9364613 0.0000239 -0.0000000
-0.0000001 0.0000001 0.0000000 0.0000239 2.9363546 0.0000008
0.0002851 0.0000894 0.0000232 -0.0000000 0.0000008 1.5522459

1329 GPa for clamped ion and 1251 GPa for the relaxed ion, while the experimental value is about 350 GPa.
Can someone tells me what wrong and should I take the value of clamped or relaxed ion?
Many thnaks

ebousquet · June 21, 2023, 7:21am

Dear phialberta,
What is the Stress tensor/Pressure at the end of the Step 1 “initial self-consistent run”?
Also, can you show us the input file you used for the relaxation? How much did you put for dilatmx and strfact flags?
Best wishes,
Eric

phialberta · June 21, 2023, 2:43pm

Hi Eric,
Thank you for your reply. here is the input
ionmov1 2 # Use BFGS algorithm for structural optimization
ntime1 5 # Maximum number of optimization steps
tolmxf1 1.0e-6 # Optimization is converged when maximum force
# (Hartree/Bohr) is less than this maximum
natfix1 2 # Fix the position of two symmetry-equivalent atoms
# in doing the structural optimization
iatfix1 1 2 # Choose atoms 1 and 2 as the fixed atoms (see discussion)

Set 2 : Lattice parameter relaxation (including re-optimization of

internal coordinates)

dilatmx2 1.05 # Maximum scaling allowed for lattice parameters
getxred2 -1 # Start with relaxed coordinates from dataset 1
getwfk2 -1 # Start with wave functions from dataset 1
ionmov2 2 # Use BFGS algorithm
ntime2 14 # Maximum number of optimization steps
optcell2 2 # Fully optimize unit cell geometry, keeping symmetry
tolmxf2 1.0e-6 # Convergence limit for forces as above
strfact2 100 # Test convergence of stresses (Hartree/bohr^3) by
# multiplying by this factor and applying force
# convergence test
natfix2 2
iatfix2 1 2

#Common input data

#Starting approximation for the unit cell
acell 5.91 5.91 9.49 #this is a guess, with the c/a
#ratio based on ideal tetrahedral
#bond angles

rprim 0.866025403784439 0.5 0.0 #hexagonal primitive vectors must be
-0.866025403784439 0.5 0.0 #specified with high accuracy to be
0.0 0.0 1.0 #sure that the symmetry is recognized
#and preserved in the optimization
#process

#Definition of the atom types and atoms
ntypat 2
znucl 13 7
natom 4
typat 1 1 2 2

#Starting approximation for atomic positions in REDUCED coordinates
#based on ideal tetrahedral bond angles
xred 1/3 2/3 0.0
2/3 1/3 0.5
1/3 2/3 0.375
2/3 1/3 0.875

#Gives the number of bands, explicitely (do not take the default)
nband 8 # For an insulator (if described correctly as an
# insulator by DFT), conduction bands should not
# be included in response-function calculations

#Definition of the plane wave basis set
ecut 20 # Maximum kinetic energy cutoff (Hartree)
ecutsm 0.5 # Smoothing energy needed for lattice parameter
# optimization. This will be retained for
# consistency throughout.

#Definition of the k-point grid
ngkpt 4 4 4 # 4x4x4 Monkhorst-Pack grid
nshiftk 1 # Use one copy of grid only (default)
shiftk 0.0 0.0 0.5 # This choice of origin for the k point grid
# preserves the hexagonal symmetry of the grid,
# which would be broken by the default choice.

#Definition of the self-consistency procedure
diemac 9.0 # Model dielectric preconditioner
nstep 40 # Maxiumum number of SCF iterations
tolvrs 1.0d-18 # Strict tolerance on (squared) residual of the
# SCF potential needed for accurate forces and
# stresses in the structural optimization, and
# accurate wave functions in the RF calculations

enforce calculation of forces at each SCF step

optforces 1

For the Stress tensor at the end of Step 1 cycle 1 . I get this :
At SCF step 14 vres2 = 5.08E-19 < tolvrs= 1.00E-18 =>converged.

Cartesian components of stress tensor (hartree/bohr^3)
sigma(1 1)= 6.53252159E-05 sigma(3 2)= 0.00000000E+00
sigma(2 2)= 6.53252159E-05 sigma(3 1)= 0.00000000E+00
sigma(3 3)= -1.89501998E-04 sigma(2 1)= 0.00000000E+00

— !ResultsGS
iteration_state: {dtset: 1, itime: 1, icycle: 1, }
comment : Summary of ground state results
lattice_vectors:

[ 5.1182101, 2.9550000, 0.0000000, ]
[ -5.1182101, 2.9550000, 0.0000000, ]
[ 0.0000000, 0.0000000, 9.4900000, ]
lattice_lengths: [ 5.91000, 5.91000, 9.49000, ]
lattice_angles: [ 90.000, 90.000, 120.000, ] # degrees, (23, 13, 12)
lattice_volume: 2.8705942E+02
convergence: {deltae: 3.553E-15, res2: 5.084E-19, residm: 5.046E-19, diffor: 6.772E-11, }
etotal : -2.54844551E+01
entropy : 0.00000000E+00
fermie : 1.86263777E-01
cartesian_stress_tensor: # hartree/bohr^3
[ 6.53252159E-05, 0.00000000E+00, 0.00000000E+00, ]
[ 0.00000000E+00, 6.53252159E-05, 0.00000000E+00, ]
[ 0.00000000E+00, 0.00000000E+00, -1.89501998E-04, ]
pressure_GPa: 5.7716E-01
xred :
[ 3.3333E-01, 6.6667E-01, 0.0000E+00, Al]
[ 6.6667E-01, 3.3333E-01, 5.0000E-01, Al]
[ 3.3333E-01, 6.6667E-01, 3.7500E-01, N]
[ 6.6667E-01, 3.3333E-01, 8.7500E-01, N]
cartesian_forces: # hartree/bohr
[ -1.72389559E-23, 2.98587474E-24, -8.68313987E-03, ]
[ -1.03433735E-23, 2.98587474E-24, -8.68313987E-03, ]
[ 1.37911647E-23, -2.98587474E-24, 8.68313987E-03, ]
[ 1.37911647E-23, -2.98587474E-24, 8.68313987E-03, ]

Regards,

phialberta · June 21, 2023, 8:58pm

Hi Eric,
Thank you for your reply. here is the input file:

#Structural optimization run

ndtset 2 # There are 2 datasets in this calculation

Set 1 : Internal coordinate optimization

ionmov1 2 # Use BFGS algorithm for structural optimization
ntime1 5 # Maximum number of optimization steps
tolmxf1 1.0e-6 # Optimization is converged when maximum force
# (Hartree/Bohr) is less than this maximum
natfix1 2 # Fix the position of two symmetry-equivalent atoms
# in doing the structural optimization
iatfix1 1 2 # Choose atoms 1 and 2 as the fixed atoms (see discussion)