NEB energy and single-point energy discrepancy

[miguel_san-miguel]

I have two adsorption intermediates of O2 on a SnO2 surface, each with distinct electronic structures. The more stable state has a magnetic moment (mag) of 2.00, while the other has a mag of 0.00, with an energy difference of approximately 0.4 eV. My goal is to determine the activation energy required to transition between these two states.

I conducted a NEB calculation. Initially, I used IBRION = 2 for the optimization of the images for about 20 steps (in my experience this helps to speed up the convergence), and later switched to IBRION = 1 until convergence was reached. Due to computational limitations, I was only able to use 5 images. The NEB calculation led to the identification of a transition state, confirmed by frequency calculations, with a barrier energy of approximately 0.036 eV and a magnetic moment of 0.00. This result seems to align with experimental data since the less stable form of O2 is not observed on the surface.

To ensure the accuracy of the result, I initiated a CI-NEB calculation using the previous NEB images. Surprisingly, within just 8 steps, it also produced a transition state with a mag of 0.00 and a barrier energy of around 0.039 eV, as indicated by the final step of the CI-NEB run.

However, when I performed a single-point calculation from scratch using the same final structure obtained from the CI-NEB run (ISTART = 0 and ICHARG = 2), the energy result significantly differed. The barrier energy was approximately 0.2 eV, and the magnetic moment was 1.80. I copied the CONTCAR from the CI-NEB to the POSCAR of the single-point run.

I'm seeking an explanation for this discrepancy and guidance on whether I should disregard the second single-point calculation from scratch in favor of the results from the first NEB calculation and the optimization energy from the CI-NEB run.

Your assistance in resolving this issue would be greatly appreciated.

[merzuk.kaltak]

Dear Miguel,
without any further input and output files I can only draw some general conclusions.
Magnetic systems often have a very shallow energy landscape. This means that there can be various different magnetic alignments with similar total energy.
If you initiate a specific MAGMOM in the INCAR to "nudge" vasp into a certain solution, for instance into an anti-ferromagnetic alignment.

Usually, the NEB method is used to study the transition path from one minimum to another.
My guess is that by using the WAVECAR and CHGCAR from the NEB images, vasp ends up into similar minima as found by the NEB calculation. The magnetic moments and energies are probably similar when starting from a preconverged WAVECAR file.

However, by setting ISTART=0 and ICHARG=2, you essentially tell vasp to ignore the pre-converged WAVECAR and CHGCAR and to start the calculation from scratch without using the knowledge of the NEB output. My guess is that you end up in an alternative magnetic solution.

It seems that the magnetic alignment is quite important to estimate the energy barrier in your case.

I cannot judge which barrier is "more" valid. Naively, I would consider the structures with the lowest energy more reliable. However, if you are interested in the barrier between to specific magnetic states, you should make sure that the first and the last NEB image corresponds to the magnetic solution you are studying and that the magnetic moments of the intermediate NEB images "interpolate" between the magnetic moments of the initial and final state.

[miguel_san-miguel]

Dear, Merzuk

Thank you for your prompt response. Unfortunately, due to the constraint of using only five images, the interpolation of magnetic properties between them appears rather abrupt. Specifically, between the first image (0, mag = 0) and the last image (6, mag = 2.0), the magnetic moment exhibits a sequence of values: 0, 2, 2, 2, 2. The first value represents the transition state (TS). This aligns with the expected but lacks the smoothness I would prefer.

I also attempted the dimer method, set up from the NEB run using the VTST script, and encountered a similar behavior. The energy of the last step in the optimized dimer corresponds to the energy barrier from the preceding NEB run. However, when performing a single-point calculation using the optimized dimer structure with ISTART = 0 and ICHARG = 2, again the resulting energy is ~0.2 higher. I did not specify the MAGMOM tag in any of my calculations, relying on the default settings, which seemed to reproduce the expected behavior for both the surface and adsorbates, consistent with experimental knowledge.

It appears that these results align with your earlier remarks about the initial magnetic setup and the shallow energy landscape. I am leaning towards placing more trust in the lower energy barrier, as it was derived conditionally from the other images in the NEB method and I think explains the fact that the least stable intermediate has never been observed at room temperature. I would expect that the really low barrier would grant its conversion to the stablest intermediate spontaneous with only the thermal energy.