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 1.0 

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.