Obtaining Piezoelectric Stress Tensor for a Strained Monolayer with GGA+SOC by DFPT

[Nishat Tasnim Hiramony]

Dear users and developers, (I submitted the same post just a while ago but did not zip all the required files and hence I am resubmitting it)

I have been trying to obtain the piezoelectric stress tensor for monolayer CaAlGaSe4 at 6% tensile strain as reported using DFPT with GGA+SOC in the following paper:
https://doi.org/10.1039/D1TC01165K

My goal is to benchmark my method by trying to reproduce the results reported in the aforementioned paper, since I am also working on a strained monolayer that has potential for piezoelectric applications.

I have followed the following steps for obtaining the piezoelectric stress tensor for CaAlGaSe4 at 6% tensile strain:

1. I optimized the unstrained structure. I optimized the structure in all three dimensions, and did not restrict the change of cell parameter along the z-axis using any patch during compilation of VASP.
2. I applied 6% tensile strain to the structure and then relaxed it using fixed-cell structural optimization.
3. I formed an orthorhombic supercell by transforming the structure obtained in the previous step (as also done in the ref. paper).
4. I tried to obtain the piezoelectric stress tensor for the structure obtained in step-3 with GGA+SOC by DFPT. The POSCAR, INCAR and KPOINTS files used for this step and the corresponding OSZICAR, OUTCAR and stdout files are provided in the .zip file attached.

Piezoelectric IO for CaAlGaSe4.zip

The k-points, EDIFF, ENCUT have been set according to those used in the reference paper.

After following these steps, I obtained e11 = 3.3757e-10 C/m (Reference: 3.071e-10 C/m) and e31 = 0.1015e-10 C/m (Reference: -0.066e-10 C/m). It is worth mentioning that the "linear response to electric field" calculation for directions 1 and 3 was stuck at a certain low accuracy till the very end (You could see them if I was able to attach the OSZICAR, log.txt), which comes off to me as a warning that something is not right.

My results are far off from the results reported and hence I tweaked some parameters and re-ran my calculations but I have not been able to obtain reliable results so far. As such, I am at a loss as to what I have been doing wrong.

Under the circumstances, I would be extremely grateful if someone could kindly help me out by pointing out what I should change in my method or in the input tags so as to obtain the correct results. I would very much appreciate any kind help.

[pedro_melo]

Dear Nishat Tasnim Hiramony,

My apologies, I tried to move your previous post, but it ended up in the wrong topic. I paste it here for safekeeping.

Regarding your results, are the values of e11 and e31 that you posted coming from the calculation with the reference parameters, or from another calculation of yours in which you changes some values?

How do the structural parameters you obtained compare with those from the reference?

Kind regards,
Pedro Melo


"""
Obtaining Piezoelectric Stress Tensor for a Strained Monolayer via DFPT
#3 Unread post by Nishat Tasnim Hiramony » Mon Jul 17, 2023 7:36 am

Dear users and developers,

I have been trying to obtain the piezoelectric stress tensor for monolayer CaAlGaSe4 at 6% tensile strain as reported using DFPT with GGA+SOC in the following paper:
https://doi.org/10.1039/D1TC01165K

My goal is to benchmark my method by trying to reproduce the results reported in the aforementioned paper, since I am also working on a strained monolayer that has potential for piezoelectric applications.

I have followed the following steps for obtaining the piezoelectric stress tensor for CaAlGaSe4 at 6% tensile strain:

1. I optimized the unstrained structure. I optimized the structure in all three dimensions, and did not restrict the change of cell parameter along the z-axis using any patch during compilation of VASP.
2. I applied 6% tensile strain to the structure and then relaxed it using fixed-cell structural optimization.
3. I formed an orthorhombic supercell by transforming the structure obtained in the previous step (as also done in the ref. paper).
4. I tried to obtain the piezoelectric stress tensor for the structure obtained in step-3 with GGA+SOC by DFPT. The POSCAR, INCAR and KPOINTS files used for this step are copied and pasted below. And I cannot seem to attach the OSZICAR, OUTCAR and log.txt files. The k-points, EDIFF, ENCUT have been set according to those used in the reference paper.

INCAR

! initialization
System = b2-CaAlGaSe4
ISTART = 0
ICHARG = 2

! electronic optimization
PREC = Accurate
! ALGO = Fast
ENCUT = 500
ISMEAR = 0; SIGMA = 0.05
! IALGO = 38
LREAL = .FALSE.
! ADDGRID = .TRUE.
LWAVE = .FALSE.
LCHARG = .FALSE.
LEPSILON = .TRUE.
IBRION = 8
EDIFF = 1E-08
NELM = 150
! LMAXMIX = 4
! LASPH = .T.
! AMIX = 0.1
! BMIX = 0.01
ISYM = 0

! ionic relaxation
LSORBIT = .TRUE.
LNONCOLLINEAR = .TRUE.
GGA_COMPAT = .FALSE.
MAGMOM = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

POSCAR

b2-CaAlGaSe4
1.0
7.3766465187 0.0000000000 0.0000000000
0.0000022684 4.2589178085 0.0000000000
0.0000000000 0.0000000000 27.6449050903
Se Ca Ga Al
8 2 2 2
Direct
0.833332300000000 0.499956608000000 0.324383303500000
0.333332300000000 0.999956608000000 0.324383303500000
0.166667804000000 0.499956250000000 0.675616696500000
0.666667819000000 0.999956250000000 0.675616696500000
0.166666761000000 0.499958038000000 0.444889411500000
0.666666746000000 0.999958038000000 0.444889411500000
0.833333135000000 0.499957740000000 0.558187082500000
0.333333135000000 0.999957740000000 0.558187082500000
0 0.999958158000000 0.501147374500000
0.500000000000000 0.499958158000000 0.501147374500000
0.166665420000000 0.499956667000000 0.359913736500000
0.666665435000000 0.999956667000000 0.359913736500000
0.833334565000000 0.499956429000000 0.641773537500000
0.333334565000000 0.999956429000000 0.641773537500000

KPOINTS

k-points
0
M
5 10 1
0 0 0

After following these steps, I obtained e11 = 3.3757e-10 C/m (Reference: 3.071e-10 C/m) and e31 = 0.1015e-10 C/m (Reference: -0.066e-10 C/m). It is worth mentioning that the "linear response to electric field" calculation for directions 1 and 3 was stuck at a certain low accuracy till the very end (You could see them if I was able to attach the OSZICAR, log.txt), which comes off to me as a warning that something is not right.

My results are far off from the results reported and hence I tweaked some parameters and re-ran my calculations but I have not been able to obtain reliable results so far. As such, I am at a loss as to what I have been doing wrong.

Under the circumstances, I would be extremely grateful if someone could kindly help me out by pointing out what I should change in my method or in the input tags so as to obtain the correct results. I would very much appreciate any kind help.
"""

[Nishat Tasnim Hiramony]

Dear Pedro,

Thanks for your reply. I have obtained the posted values of e11 and e31 using the input files that I have provided in the .zip file. The k-points, EDIFF and ENCUT used are same as those in the reference paper. The paper used GGA+SOC by DFPT for obtaining e_ij and so did I.

The reference paper focuses on two different materials and it did not mention the structural parameters for one of them, which is CaAlGaSe4, the material I am currently trying to benchmark my results with. And hence, I could not compare the structural parameters. However, I also obtained the elastic constants for the same material, and the obtained C11 and C12 match exactly with the reference paper. Hence, I am assuming that my structure is okay.

Thanks and regards.

[Nishat Tasnim Hiramony]

Dear Pedro,

My apologies, I just noticed that they did mention the lattice constant for the hexagonal unit cell, which is 4.02 Å, which is the same as the one that I obtained. Please note that the POSCAR file I provided here is for the orthorhombic supercell (used for obtaining the piezoelectric tensor), not the hexagonal unit cell.

Thanks and regards.

[Nishat Tasnim Hiramony]

Dear Pedro,

I would also like to mention that the posted values of e11 and e31 have been calculated from the OUTCAR in the following way:

Total Piezoelectric Tensor (in C/m^2) = PIEZOELECTRIC TENSOR (including local field effects) for field in x, y, z (C/m^2) + PIEZOELECTRIC TENSOR IONIC CONTR for field in x, y, z (C/m^2)

=> Total Piezoelectric Tensor (in C/m) = Total Piezoelectric Tensor (in C/m^2) * 27.6449050903e-10 (this is the "c" latt. parameter)

I have come to this method after benchmarking my results against another paper by the same author. However, since that material was not strained and the INCAR file for that material was already provided in the paper, I am trying to re-benchmark my method.

Thanks and regards.

[Nishat Tasnim Hiramony]

An update: the negative sign of e31 can be produced by reversing the z direction. Otherwise, the values agree to a reasonable degree in my opinion. So, I'm guessing the method is okay.