Piezoelectric tensor

Dear community, i’m trying to perform a d33 calculation of BaTiO3 but the anaddb code failed and give this message :
Calculation of the internal-strain tensor

instrflag=1, so extract the internal strain constant from the 2DTE

src_file: anaddb.F90
src_line: 840
mpi_rank: 0
message: |
DDB file must contain both uniaxial and shear strain for piezoelectric, Check your calculations

abinit_abort: decision taken to exit. Check above messages for more info

this is my input file :

ndtset 3
paral_rf2 1
paral_rf3 1

#Definition of the elementary cell
rprim 1.0 0.0 0.0
0.0 1.0 0.0
0.0 0.0 1.0
acell 7.6152197110E+00 7.6152197110E+00 7.6152197110E+00 # relaxation value

#Definition of the atoms

xred 0.000000 0.000000 0.000000
0.500000 0.500000 0.500000
0.000000 0.500000 0.500000
0.500000 0.000000 0.500000
0.500000 0.500000 0.000000

nband 32 # total valence electron/2 + 10
natom 5
ntypat 3
znucl 56 22 8
typat 1 2 3 3 3

#Definition of the SCF procedure
nstep 200
fftalg 401
diemac 4.0

#Definition of the plane wave basis set
ecut 32
pawecutdg 90
ecutsm 0.5
dilatmx 1.05

ngkpt 6 6 6
nshiftk 1
shiftk 0.5 0.5 0.5

#Ground state calculation
kptopt1 1 # Automatic generation of k points, taking
# into account the symmetry
toldfe1 1.0d-12 # SCF stopping criterion

#Response Function calculation : d/dk
rfelfd2 2 # Activate the calculation of the d/dk perturbation
rfdir2 1 1 1 # Need to consider the perturbation in the x-direction only
# This is rather specific, due to the high symmetry of the AlAs crystal
# In general, just use rfdir 1 1 1
# In the present version of ABINIT (v4.6), symmetry cannot be used
# to reduce the number of ddk perturbations
nqpt2 1
qpt2 0.0 0.0 0.0 # This is a calculation at the Gamma point
getwfk2 -1 # Uses as input the output wf of the previous dataset
kptopt2 2 # Automatic generation of k points,
# using only the time-reversal symmetry to decrease
# the size of the k point set.
iscf2 -3
tolwfr2 1.0d-22 # Must use tolwfr for non-self-consistent calculations
# Here, the value of tolwfr is very low.

#Response Function calculation : electric field perturbation and phonons
rfphon3 1 # Activate the calculation of the atomic dispacement perturbations
rfatpol3 1 5 # All the atoms will be displaced
rfelfd3 3 # Activate the calculation of the electric field perturbation
rfdir3 1 1 1 # All directions are selected. However, symmetries will be used to decrease
# the number of perturbations, so only the x electric field is needed
# (and this explains why in the second dataset, rfdir was set to 1 0 0).
nqpt3 1
qpt3 0.0 0.0 0.0 # This is a calculation at the Gamma point
getwfk3 -2 # Uses as input wfs the output wfs of the dataset 1
getddk3 -1 # Uses as input ddk wfs the output of the dataset 2
kptopt3 2 # Automatic generation of k points,
# using only the time-reversal symmetry to decrease
# the size of the k point set.
toldfe3 1.0d-12
pawxcdev3 0

add to conserve old < 6.7.2 behavior for calculating forces at each SCF step

optforces 1

is there any suggestions ?

Dear Ayman,
If I’m not mistaken, I do not see any perturbation for the strain in your input ( see DFPT - abinit), right? If so, it is normal that the piezo cannot be calculated.
Best wishes,

I noticed that mistake and I added rfstrs = 3 but I get a negative value of the d33 which should be positive as the experimental value !!
Best regards

OK, did you check that all of your calculations (ground state = deep relaxation of the forces and stress, DFTP scf or ecut and kpoints) are well converged?

Yeah, i did all of that, maybe the experimental value is negative !

The sign is positive if I remember well, to be checked, how much did you get for the amplitude?
I also see that you have a dilatmx of 1.05, which is large, check that the ground state WF calculation has a stress as small as at the end of the relaxation.

Sorry I did not check that it is a physics problem that you have here, you try to compute the piezo response of the cubic phase of BaTiO3 which is not piezoelectric. You have to do the calculations in one of the 3 ferroelectric phases, ideally in R3m phase to have no remaining unstable mode in the structure but the other phases can be calculated knowing the remaining instabilities in some directions.

Hi Eric,
I too got a negative value of d33 while calculating the piezoelectric properties of AlN in superlattice (two cells along z) . I add that I obtained a value close to the experimental when I calculated it taking a single cell. Regarding your answer , I have checked if there are any negative frequencies; I got this
Phonon frequencies in cm-1 :

  • 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00
  • 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00
  • 0.000000E+00 0.000000E+00 3.505561E+02 3.842994E+02 3.842994E+02
  • 3.938523E+02 3.938523E+02 3.938526E+02 3.938526E+02 4.009857E+02
  • 4.009863E+02 5.211343E+02 5.211343E+02 5.359462E+02
    chkph3 : WARNING -
    The dynamical matrix was incomplete : phonon frequencies may be wrong …

Best regards,

Dear phialberta,
Your phonon frequencies show that your phonon calculation is incomplete as stated by the warning. I expect it is because the rfatpol does not include all of the atoms N, i.e. it should be rfatpol = 1 N, can you check this? Of course any further calculation using phonon perturbation will be incorrect…
Best wishes,

Dear Eric,
Thank you for your help. Effectively , I add rfatpol 1 N , the matrix was completed and i get the correct value of d (piezoelectric constant). I have another question about fixing atoms during atoms relaxation GS calculation; following the elastic tutorial of wurtzite AlN, I fixed the Al positions to be at reduced c axis coordinates. Is it necessary to fix them in the next superlattice calculations? or just relax all the atoms in all directions?


Dear phialberta,
Good that your problem is fixed now :slightly_smiling_face:
For the atom fixed, you can relax all the atoms I do not remember this atom fix in the tuto, thank you for pointing toward it.
Best wishes,

1 Like

Hello again,
I come back to this question because I have just discovered that the wave functions have not converged despite the fact that I tried to decrease the tolwfr from 1.0d-20 to 1.0d-18. I get this ;
scprqt: WARNING -
** nstep= 200 was not enough SCF cycles to converge;**
** maximum residual= 2.656E+00 exceeds tolwfr= 1.000E-18**

Have you any suggestions about this problem?

Best regards

Dear phialberta,
For insulators, the residual of the WF might looks bad if you have empty bands in your calculations. Often, the last unocciped bands are not converged but they do not contribute to the energy so that we can disregard them. You have to check how many occupied bands you have in your system and either you put nband to the number of unoccupied bands in your DFPT calculation (occopt=1 only) or you can also use the flag nbdbuf and put it to 2 or more (up to the number of unoccupied bands, occopt=1 or metallic cases), see the description of this input variable, it has been designed for that problem.
This “problem” is often evidenced when the energy and potential residuals are converged but not the WF one.
Hopefully this is the problem you have, otherwise this means that the SCF has to be tuned with diemix, etc kind of convergence flags.
Best wishes,

Hi Eric,
Thank you for your response. I have just decrease the nband taking into account only the occupied one and set nbdbuf to 4 , the wf had converged but I got a negative piezoelectric constant d ?
Is it the problem with nbdbuf 4 ?


Dear @phialberta,
The nbdbuf is useful if you have bands than the occupied ones, but if you set nband to the number of occupied ones, then you should not use nbdbuf as it will disregard the last occupied bands then, which you want to take into account in the residual of the SCF.
Regarding the negative d, couldn’t it be negative, are the elastic constants negative too?
You have to check if the relaxation is well done with the pressure and check if your calculation is converged, see this post:

Best wishes,

Hi Eric,
Effectively, I recalculated it with nband equals to the number of occupied ones without using nbdbuf and I get the correct value of d.
Many thanks,
Best regards