inversion symmetry of slab doesn't break with applied EFIELD

[sophie_weber]

Hello,

I am trying to perform calculations simulating an applied electric field perpendicular to an antiferromagnetic vacuum-terminated metallic slab (RuO2) which has inversion symmetry using the sawtooth potential (EFIELD) along with a dipole correction (IDIPOL=3, LDIPOL=True) along the direction of vacuum. At first I thought that my results were reasonable, but then I noticed that when I compare the converged OUTCAR file, in particular the atom-projected magnetization, with a finite electric field of 0.005 eV/angstrom (OUTCAR_p005E), the magnitude of the magnetic moments on the Ru, while different from the calculation with no EFIELD (OUTCAR_0E), are still exactly inversion-symmetric. I would expect that since an electric field breaks inversion symmetry, the magnetic moments of inversion-related sublattices on the top and bottom of the slab should no longer be identical.

Attached are the input files, as well as the OUTCAR files with and without electric fields. If someone could provide insight into what I'm doing wrong (or if the calculations seem correct, which the inversion breaking of the electric field doesn't seem to show up in the projected magnetic moments), I would greatly appreciate it.
Thank you in advance!

[henrique_miranda]

I took a quick look at the OUTCAR files, but could not immediately see that the inversion symmetry is not broken.
But you know what to look for so I would suggest that you try two things:

Add the NWRITE=3 INCAR tag to your input file
This will increase the verbosity of the OUTCAR file according to
https://www.vasp.at/wiki/index.php/NWRITE
but then you can look for the symmetry operations that VASP will use.
You should look that adding the EFIELD, IDIPOL = 3 and LDIPOL = .TRUE. INCAR tags will lead to less symmetries being used.
If that is not the case please report back in this thread.

Turn off symmetry by setting ISYM=0 or ISYM=-1

Let me know if any of these suggestions solves your problem.

[sophie_weber]

Hello,

Thank you very much for the prompt reply. The reason I say the inversion symmetry is not broken in the presence of the electric field is that the magnetic moments (they are printed at the very end of the OUTCAR) on the Ru atoms on the top and bottom surface of the slab, which are related by inversion symmetry in the absence of electric field, still have identical magnitudes. They should have different magnitudes to the presence of the polar vector of the electric field which breaks the inversion symmetry. I tried turning off symmetries but this did not resolve the issue.

I also today found something quite concerning which I think may point to a possible error in the EFIELD implementation in the >6 versions of VASP (or at least, I am not using the correct INCAR flags).
If I run a completely identical calculation with an old version of VASP, version 5.4.4 specifically, for the calculations without an electric field the 6.2 and 5.4 versions give basically identical results, but for the 0.005 eV/angstrom electric field, the magnetic moments do exactly as I expect, i.e. the magnitudes of the moments related by inversion symmetry in the structure are no longer equivalent due to the presence of the electric field. Now, the top and bottom surface Ru moments are quite different.

I know there were some bugs in magnetic symmetry in older versions of VASP that were fixed, so I'm very surprised that the newest versions of VASP would introduce a symmetry error that did not exist in the older version, but that does indeed seem to be the case based on what ought to happen physically.

Please let me know if you have insight into this, this is quite important to resolve.

I attach the files with the OUTCARs for both vasp versions, with and without electric fields, to compare.

Thanks a lot again for your help,
Sophie

[sophie_weber]

Sorry, but one more update; I see now that if I turn symmetries off in version VASP 5.4, the inversion symmetry of the magnetic moments in the OUTCAR is broken but only by an extremely tiny amount, and only for one sublattice (note the differences between sublattice 3 and sublattice 13). This is in start contrast to when the symmetry is turned on, and the change is quite substantial, which is opposite to what I'd expect you'd think that with symmetries off the changes would be more substantial if anything!) For the same value of electric field, with VASP version 6.2, the moments stay exactly the same in the bottom and top half of the slabs. So now I'm even more confused about which combination of VASP version and symmetry on/off is giving the correct physics.

I attach the two OUTCAR files with symmetries off (ISYM=0)

thanks,
Sophie

[henrique_miranda]

Thanks for your quick feedback.
I played around with your input files of my side a little bit:
1. I did not see any difference in the symmetry operations that are detected depending on EFIELD and IDIPOL. I think that when we apply an external electric field by modifying the potential, we are breaking some symmetries of the system. This is something we might need to fix in VASP. In the meantime, you can either turn off the symmetries completely with ISYM=-1 or slightly break the initial magmom such that you obtain the same effect. The advantage of the second option is that you still keep some symmetries, so your calculation is faster.
2. It seems that you are running into some convergence issues in this calculation. There might be many reasons for it, but I would suggest you address these convergence problems before drawing any conclusions about different results between different versions of VASP. The point is that if the calculation is slowly converging then the code might stop at a point where the energy convergence criteria is satisfied, but you are not converged. And the different versions of VASP change slightly the path you take to get to convergence.
3. You might want to consider using a different lattice in your POSCAR which you can generate using the following script https://gist.github.com/henriquemiranda ... af6484f890 (see attached POSCAR).
This will make it clearer (at least to me) in which direction you should have k-point sampling (no need for k-points in the z direction).
Also which direction is z which is important for IDIPOL=3 [EDIT: this statement is wrong as pointed out by sophie_weber in the post below, IDIPOL refers to the lattice vector and not Cartesian direction]
4. You have a metallic system(?) then it might be better to use ISMEAR=1 instead of ISMEAR=0 but please do your own testing.

As a rule of thumb, I would say that the calculation should converge below the default 60 SCF steps, at least that would be ideal.
Can you check if with some of the modifications above you get a better SCF convergence?

[henrique_miranda]

If this new POSCAR with a new corresponding KPOINTS file + ISMEAR=1 instead of ISMEAR=0 does not help with convergence, then you might consider trying some of the suggestions in this page:
https://www.vasp.at/wiki/index.php/Trou ... onvergence

[sophie_weber]

Hello,

Thanks very much for the helpful responses.

I am pretty confident that the fairly long time to scf convergence is not the issue here; I've worked with a lot of large, metallic systems, and it is not so uncommon that in these compounds it takes up to a few hundred steps to converge. But I've done quite a few tests across different settings, different versions of VASP, etc., and in these systems (including this one) the outputs for any calculations without electric field are essentially identical. So It think that's a good indicator its' properly converged.

I agree however that it's good to eliminate variables so I'm running tests now with an insulating system, FeF2, which is isostructural, and takes less time to converge (my main concern is that I know that the relative changes in magnetisation for this system will be smaller), so in this sense It may be harder to debug the EFIELD issue. In any case they are running now, and I'll send you the results as soon as they finish.

A couple other comment's questions.
1. To be safe I'll rotate the cell so the surface normal of the slab is along z (this is a (110) slab hence why the lattice vectors looked the way they did), but I'm quite sure the k-mesh and IDIPOL directions should be fine, because this is a cell with all right angles, and I'm using direct rather than cartesian coordinates. Also, I think the vasp wiki specifics explicitly that IDIPOL directions are determined by lattice vectors, not cartesian coordinates.
https://www.vasp.at/wiki/index.php/IDIPOL
2. " I did not see any difference in the symmetry operations that are detected depending on EFIELD and IDIPOL. I think that when we apply an external electric field by modifying the potential, we are breaking some symmetries of the system. This is something we might need to fix in VASP. "
Is it not the case that the symmetry routine in VASP only looks at the symmetry of the POSCAR? if that's the case, then IDIPOL and EFIELD settings should not matter for the symmetry operations that are detected. But if indeed the symmetry printed out should take into account the symmetry reduction of the electric field, then that's a huge problem...if the EFIELD VASP implementation doesn't break inversion symmetry how does it possibly mimic any sort of realistic electric field properties? I'll do tests with symmetries turned off, but if the implementation of the electric field itself fails to incorporate inversion breaking then I don't see how this implementation as is can be used.
3. Are you aware of any changes to the EFIELD implementation between VASP 5.4.4 and VASP6.2? The results with symmetries turned off seem somewhat similar (though not what I would expect) but the results with symmetries turned on are completely different.

Thanks again, I'll update you as soon as I have more results

[henrique_miranda]

1. You are correct of course, I apologize for my confusion (I added a note in my post above so other users don't get confused by my comment).
2. The symmetry routines look at the lattice, positions and velocities (from POSCAR) but also MAGMOM (which is in the INCAR). These symmetry operations are then written to the OUTCAR and are used to reduce the number of k-points and to symmetrize the charge density, for example. It is a problem that the applied electric field is not reducing the amount of symmetries in VASP, and we will fix this in the code. The implementation can be used, but you need to either turn off all the symmetries (ISYM=0) or break the symmetries in another way (for example, by setting slightly different value for MAGMOM for atoms that are related by symmetry). It is of course likely that you will get wrong results when VASP is applying symmetries that should not exist in the system because they were broken by the external electric field so for now I recommend setting NWRITE=3 and inspecting the symmetry operations manually.
3. No. I looked at the output files you sent me and I could not find a difference between VASP 5.4.4 and VASP 6.x.x. I will check if I reproduce your results and if so, what is the reason for this difference. My first suspicion is the one I mentioned above (slow convergence + different path to convergence due to small changes in the code between VASP 5.4.4 and VASP 6.x.x, but it might be something else completely).

[sophie_weber]

Thanks a lot for the clarifications, I appreciate it.

I'm sorry to be annoying, but I just want to make sure I'm understanding correctly; "The symmetry routines look at the lattice, positions and velocities (from POSCAR) but also MAGMOM (which is in the INCAR). These symmetry operations are then written to the OUTCAR and are used to reduce the number of k-points and to symmetrize the charge density, for example. It is a problem that the applied electric field is not reducing the amount of symmetries in VASP, and we will fix this in the code. "
But if the symmetry assumed by VASP takes into account the MAGMOM and the POSCAR, my question still holds; In my INCAR files I always initialise the magnetic moments to be inversion symmetric/identical at top and bottom slabs because the point/idea is that when the electric field is included in the calculation, I want to see how much it changes the moments, which will be inversion symmetric in the absence of the electric field. Moreover, note that in the calculation with 0.005 eV/Angstrom used in VASP 5.4.4, with symmetries left on, yields an OUTCAR which has very different magnetic moments for inversion-symmetric sites. It is only in VASP 6.2 that the OUTCAR keeps inversion-symmetric magnetic moments.

This makes me suspicious that simply turning symmetries off will solve the trick, because it seems like there is something else problematic in the implementation of the EFIELD; it gives inversion-asymmetric results for both ISYM=0 and symmetries on for VASP 5.4.4, but with very different magnitudes. In VASP 6.2, it only gives inversion-asymmetric results with symmetries turned off.
Relatedly, what is the EFIELD doing if not breaking the inversion symmetry? Isn't it just a sawtooth potential?

If you could let me know your results trying to reproduce these that would be great. I'll check with this insulating isostructural system in the meantime to try and eliminate the convergence issue possibility and update you asap they should be done tomorrow.

Thanks again for your help!

[henrique_miranda]

I would also expect that since VASP is enforcing the symmetries (while it shouldn't when you apply an external electric field) the final result would also be symmetric, but as you reported it is not.
It is strange that it is when you turn off the symmetries in VASP 5.4.4 that the result becomes symmetric again.
It is more likely is that the correct result is the one you obtain when the symmetries are turned off: the electric field might be coupling only lightly to the magnetization.

That being said, I need to find out what is the reason for this change in behavior between VASP 5.4.4 and VASP 6.
For that, I tried running the INCAR files that you shared, but these calculations take too long to run to be practical to track down what causes the different behavior between vasp 5.4.4 and vasp 6
Are you able to reproduce this same issue on a smaller system? That would make bisecting the code changes between the different versions much easier.
Maybe one could simply remove a few atoms and vacuum from the POSCAR you shared before.
Would that make sense?

[sophie_weber]

Apologies, the calculations are taking me longer than expected. I tried using a smaller slab and it is converging quicker now, and in fact, at least with the symmetries turned off, VASP 5.4 VASP 6.2 seem to produce the same results.
However, the changes in magnetization seem like that may be random, rather than systematic, in terms of the strength of the electric field, which is physically not reasonable.

These should be done over the weekend and I'll send you the output Monday.

Thanks a lot for your time and help,
Sophie

[sophie_weber]

Hello,

Apologies for the delay. I ran tests on a slab structure (FeF2 versus RuO2, but it is magnetically and crystallographically isostructural, and converges easier than RuO2) which is half the thickness of the previous tests I showed you. With this four-layer slab I was able to converge calculations with symmetries turned off, and when I do this, indeed I get the same output whether I use version VASP 5.4 or version 6. However, if you look back at the EFIELD_vaspversion_compare.zip file I sent you, with symmetries on there is clearly a huge difference in output, with inversion symmetry being maintained in v 6.2 in the magnetic moments in the OUTCAR, but not in version 5.4, with an applied field. So I think obviously there was some change in the implementation of EFIELD between these two versions, even if keeping symmetries on is not the correct procedure.

Separately, even with symmetries turn off and converged calculations with the smaller slab (the attached outputs are with VASP 5.4.4, but as I said in the case 6.2 yields the same results), the implementation of EFIELD seems problematic. the changes in magnetic moments seem completely random with the strength of the electric field (see attached plot, where I've just plotted the total magnetic moment of the two Fe sublattices in the top layer of the slab, as a function of EFIELD value). At the very least, the changes in magnetization should be monotonic with the strength of the applied field, if the calculations are converged (which as far as I can tell, they seem to be). So this makes me think that (1) the calculations are not converged, event though they appear to be, or (2) the EFIELD implementation, even with symmetries turned off, does not reliably mimic an applied electric field.

If you have any insight or thoughts based on the output files (or if it's possible to do your own tests with the smaller slab) I'd greatly appreciate it.

Thanks,
Sophie

[henrique_miranda]

Sorry for my delay in answering.
I tried running your calculations locally, but I still could not reach a definitive conclusion or have a suggestion.
The data you show about the magnetic moments oscillating wildly with changing electric field seems to indicate that there is a problem reaching convergence or the parameters of the calculation are not precise enough to resolve these moments clearly.
But I cannot make a defenitive statement about this.

What the EFIELD implementation does is rather simple, you can verify it yourself: it adds a contribution to the local potential to simulate the electric field.
Because we use periodic boundary conditions, at some point it needs to have a step in the potential which it tries to put in the vaccum region (there is some logic in VASP to find where that is).
You can inspect this by looking at the LOCPOT file or visualizing it or reading from the vaspout.h5 file if you have vasp.6.4.3 by using the WRT_POTENTIAL tag
https://www.vasp.at/wiki/index.php/WRT_POTENTIAL
you can try the attached python file (it's not very carefully written, but maybe it helps).
(I uploaded a new version of the python script)

I know this is not very helpful.
I will try these calculations a bit more and let you know if I find something that can help.

[henrique_miranda]

The problem of convergence comes from the electrostatic corrections

Code: Select all

 EFIELD=0.05
 IDIPOL=3
 LDIPOL=.TRUE.
 DIPOL=0.5 0.5 0.228045

One trick that you could try to greatly improve convergence is to run the calculations in two steps:
1. with these tags commented out (probably will converge in less than 60 SCF iterations)
2. and then another one with the tags uncommented, starting from the WAVECAR of the previous run (probably converged in less than 20 interations).
With this trick, you can also generate the plot of dipole vs electric field much more efficiently.

Let me know if this helps.

PS: you might use this same trick to see if you still end up with different results when using vasp 5 or 6 (always with ISYM=-1).

[sophie_weber]

Thanks a lot; I will try this and let you know if it helps!

Best,
Sophie