NEB and geometry relaxation with Electric field

[santanumahapatra]

Hello, I want to test the effect of electric field on NEB barrier between 2 monolayer structures of MoS2. For that, first I am doing a geometry relaxation of the 2 structures in the presence of Efield. Following is my INCAR file for geometry relaxation.

Code: Select all

SYSTEM  = 2H-MoS2
ISTART=0
ENCUT=600
ALGO=Fast
IBRION=2
ISIF=3
NSW=800
EDIFF=1.0e-05
ISYM=0
LREAL=.FALSE.
ISMEAR=0
SIGMA=0.1
PREC=Med
ADDGRID=.TRUE.
NWRITE=1
LCHARG=.FALSE.
LWAVE=.FALSE.
EFIELD_PEAD = 0.0 0.0 0.8
LCALCEPS=.FALSE.

and below INCAR file for NEB:

Code: Select all

SYSTEM = Ge136

#PREC = Medium
PREC = Normal
ISTART=0
ISYM = 0
#IOPT = 3
IBRION = 3
#IBRION = 0
POTIM = 0.1
#POTIM = 0
NSW = 100
ISIF = 3
#NELMIN = 5
#NELM = 400
#NELMIN = 4
EDIFF = 1.0e-05
EDIFFG = -0.01
#ENCUT = 400.0
ENCUT = 600
#ICHARG = 2
LWAVE = .False.
LCHARG = .False.
IALGO = 48
ALGO = Fast
LREAL = .FALSE.
EFIELD_PEAD = 0.0 0.0 0.8
LCALCEPS=.FALSE.
ISMEAR = 0
SIGMA = 0.1
IMAGES= 8
SPRING= -5.0
ICHAIN= 0
LCLIMB = .TRUE.

I am not getting any energy changes when I relax the structures in the presence of Efield as compared to their energies when relaxed without any Efield, which should not be the case. Is this the correct settings for calculations in the presence of Electric field? I am trying to apply a Efield of 0.8V/nm, however, EFIELD_PEAD tag has units of ev/Ang. Do I also need to multiply 0.8 with 1.6E-19 in the INCAR file?

[ahampel]

Dear santanumahapatra,

thank you for asking your question on the official VASP forum.

The tag EFIELD_PEAD is using linear response to calculate the effect of an applied electric field. However, this will not affect calculated forces in VASP. Only the properties described here: tutorials/latest/response/part1/#response-e01 or wiki/index.php/LCALCEPS here are calculated under the influence of the electric field. Therefore, the result of your ionic relaxation must be the same, and performing ionic relaxation using this tag is not possible.

However, it is possible in layered systems (slabs, interfaces, mono-layers) to apply an electric field wiki/index.php/EFIELD and compensate the induced charged cell via dipole corrections. See also wiki/index.php/Electrostatic_corrections . I think your system MoS2 might have enough vacuum out of plane to try this. This means you could try to remove all EFIELD_PEAD related tags and instead use something like this:

Code: Select all

EFIELD = 0.8   
IDIPOL = 3      
LDIPOL = .TRUE. 

I only have limited experience with this, but maybe someone else on the forum tried this and can share his or her experience.

Please attach your POSCAR, KPOINTS, and POTCAR file the next time as well. This will me allow to try the calculation directly and verify it.

Best regards,
Alex H.

[santanumahapatra]

Please find the relevant input files.

[santanumahapatra]

Can you also specify if I need to factor in the value of 1.6E-19 in the INCAR file for unit conversion?

Thanks

[ahampel]

According to wiki/index.php/EFIELD the electric force field is applied in units eV/Å . If you need to apply a conversion to 1V/Å times 1e , yes you have to multiply by the charge of an electron, i.e. 1.60217663 × 10-19 . See also:
https://w.vasp.at/forum/viewtopic.php?t ... 265b439c3a for a discussion.

Best,
Alex H.

[santanumahapatra]

The relaxation is not converging. Following is my INCAR file, remaining input files are same.

Code: Select all

SYSTEM  = 2H-MoS2
NELM=200
ISTART=0
ENCUT=600
ALGO=Normal
IBRION=2
ISIF=3
NSW=100
EDIFF=1.0e-05
ISYM=0
LREAL=.FALSE.
ISMEAR=0
SIGMA=0.1
#IALGO=48
PREC=Accurate
ADDGRID=.TRUE.
NWRITE=1
LCHARG=.FALSE.
LWAVE=.FALSE.
LDIPOL = .TRUE.
IDIPOL = 3   # for Z direction
EFIELD = 2   # eV/A**

[ahampel]

the relaxation was not the problem. Your electronic scf calculation was not even converging, correct? I tried now to run your example. I think that the dipole that is placed in the cell by VASP was put at an awkward position. After specifying manually its position to be in between two mono-layers:

Code: Select all

DIPOL = 0.0 0.0 0.5

the calculation seems to run. I am not an expert but the calculation is converging and the relaxation goes to a different structure. I started with EFIELD=0.5 first and then increased it. There is no guarantee here for correctness, as I am not an expert.

Please check https://www.vasp.at/wiki/index.php/DIPOL for more information.

Best,
Alex H.

[santanumahapatra]

Thank you for your advice.
As far as NEB is concerned, do I simply add the electric field tags and run a new NEB as I do? Or do I take each structure in the NEB path for 0 Electric field and run self consistent calculations on each of them by adding the new electric field tags? The reason why I am asking is because under the presence of an electric field, the new NEB path can be different than what we obtain without any bias. In ideal case, I would like to keep the NEB pathway same, just check the effect of bias on the highest point (peak) of the barrier.

Thanks

[ahampel]

As I said before I am not an expert on the matter but this seems to be feasible. Did you check literature for the correct procedure:

https://www.researchgate.net/figure/NEB ... _378715437 i.e. https://onlinelibrary.wiley.com/doi/abs ... .202301074

or https://www.sciencedirect.com/science/a ... 802200253X

Maybe these two have some insights that you find helpful. I thought the idea is that you have two structures and you would like to understand how this reaction is influenced when applying an electric force field. I would probably first try to run NEB without efield and then add the efield only for the NEB, i.e. keep the pathway the same. However, I think it could happen that you run into trouble since the electric field could change the band structure too much.

Hope that helps.

• Alex H.

[santanumahapatra]

I am trying to find the effect of EFIELD on NEB between two known structures of MoTe2 (Hex and T-dash). However, it seems external electric field is hardly changing the barrier height, which should not be correct. I am attaching the relevant input files for different applied electric fields, can you please advice what can be wrong in the calculations?

Thanks

[ahampel]

Do you have the suspicion that something in VASP is not working correctly. Does changing the field strength change something or not? If the result is numerically the same I am happy to look into it, as this points to a problem. Otherwise this seems to be more of a research question, and I am unfortunately not the right person to help here. Maybe another user in the forum can assist?

Best regards,
Alex

[santanumahapatra]

If I want to calculate the NEB between 2 structures of different unit cell size, for eg, the HEX version of MoTe2 has 1 Mo and 2 Te in its unitcell,whereas the T-dash version has 2 Mo and 4 Te in its unit cell, how do I match the 2 structures appropriately? I am creating a supercell of POSCAR1 in x direction as 2X1 to match the lattices of both the strcutures, but I believe, I am making a mistake in this process, which is messing up the NEB.

Can you please advise?

Thanks

POSCAR1:
"mo1 te2"
1.0
3.5473923683 0.0000000000 0.0000000000
-1.7736961842 3.0721319082 0.0000000000
0.0000000000 0.0000000000 18.6102008820
Mo Te
1 2
Cartesian
0.000000000 0.000000000 9.305100441
1.773696184 1.024044037 11.110199928
1.773696184 1.024044037 7.500000000

POSCAR2:
Mo Te
1.0000000000000000
7.2446165085000001 0.0000000000000000 0.0000000000000000
-1.6429158851000001 3.0327857418000002 0.0000000000000000
0.0000000000000000 0.0000000000000000 18.7289733886999983
Mo Te
2 4
Cartesian
0.3830421862607776 0.2076164477063406 9.2713813716718807
1.5963522189902730 2.8251726635186185 9.4574880524968368
1.8370588990892405 0.9930125464258893 11.4270439186893391
0.1422670458099391 2.0396604545969450 7.3019046728498926
5.4119138684010126 0.9690348510399284 10.8459491636392116
3.8121609682510118 2.0638608929317126 7.8832058961026581

[ahampel]

Hi,

I will take a look at your last question, but let me quickly stress that I made a mistake in earlier posts regarding the position of DIPOL. DIPOL should denote the maximum of your density, i.e. where the monolayer is located: DIPOL = 0.0 0.0 0.5 . VASP definition is such that it takes this value add integrates to DIPOL + lz*0.5 to determine the dipole correction. Sry for that. But this gives quite different results even though both converge. However, with this the warning that VASP might give in the end regarding a field free region should disappear. Maybe this helps you already.

[ahampel]

Did you have a look at https://vasp.at/tutorials/latest/transi ... tes/part1/ ? I am not sure I understand correctly what you did so far. Can you provide me with the exact set of POSCAR files that you used for the NEB calculations? The initial and final structure need to have the same number of atoms and matching species. It is also best to start with the the unit cell, even though with ISIF=3 you might be able to also change the cell shape in a NEB calculation. If you provide me with all the input files I can attempt to run you calculation.

Best,
Alex

[santanumahapatra]

Thank you for your reply. I have given the 2 POSCAR files below:

POSCAR1:
Mo Te
1.000000000000000
5.3759672164800003 0.0000000000000000 0.0000000000000000
-0.8335742256200001 4.4548567822800003 0.0000000000000000
0.2504815557600000 -1.3381739828399999 18.5990361925800016
Mo Te
2 4
Direct
0.1883941867275493 0.6194728383173062 0.4999966569049662
0.6902974963180881 0.1123888951434459 0.5000246834665349
0.5756972766347251 0.6402997363059992 0.5904640454215395
0.5749501610152195 0.5852057637448329 0.4093078273280262
0.0754435710775931 0.1439912251012796 0.5910609826754648
0.0752198698268164 0.0880503809871494 0.4091464228034732

POSCAR2:
Mo Te
1.0
7.0947799683000001 0.0000000000000000 0.0000000000000000
-1.7736954010999999 3.0721296071999999 0.0000000000000000
0.0000000000000000 0.0000000000000000 18.6102008820000009
Mo Te
2 4
Direct
0.0000000000000000 1.0000000000000000 0.5000000000000000
0.5000000000000000 1.0000000000000000 0.5000000000000000
0.3333300000000000 0.3333300000000000 0.5969900000000000
0.3333300000000000 0.3333300000000000 0.4030100000000000
0.8333300000000000 0.3333300000000000 0.5969900000000000
0.8333300000000000 0.3333300000000000 0.4030100000000000

Basically, I wish to test the effect of an external electric field on the barrier height, but unfortunately, the Efield does not seem to affect much. I believe, the issue is because of the drastic difference in the lattice constants of the 2 endpoints.

Thanks