UT theoretical chemistry code forum

Hi,

I have a question regarding an NEB calculation I have been carrying out...

My system is a PdCu(111) alloy and I am using NEB to find the transition state for dissociative hydrogen adsorption via Pd.

Initially I minimised the energy of the two extrema, with the initial state being a physisorbed hydrogen molecule effectively on the atop Pd site (albeit slightly more displaced in the z direction than when bound) and the final being one hydrogen atom in an fcc and one in a hcp hollow site (H-Pd-H angle of 180 deg).

The problem I am having is two-fold:

1) the energy of the final-1 and final-2 images are lower than the final image...
Naturally I thought this was due to not minimising the final image correctly. However, I then took the final-1/2 images, relaxed them and they returned to the same energy as the final image! For what reason might this happen? Do the NEB forces artificially create this behaviour? If I proceed to use the ELAST tag will this help or should this situation never occur in the first place?

I am confident there is a transition state here, as I took the maximum energy image and used the dimer method to find the TS which, after vibrational analysis seems sound enough.

2) Similarly with the First image when looking at a pure Cu(111) surface (not with the PdCu(111) surface), this is not lower than, but is much much closer to the transition barrier than expected and compared to previous reports of this mechanism. Could this be due to a poor choice of initial image? The image is where H-H is place ~2 angstroms above the surface, would you recommend something else? I had thought of using an image that has hydrogen in the gas phase and approaching the surface, although as this is a barrierless transition, would it affect the NEB calculation?

Upon visualisation of both of these scenarios, it is clear to see there is significant z-displacement by the metal atoms such that as the system moves from the initial to final state, the compression between layers in the slab increases. This is particularly notable in the PdCu(111) case as Pd sinks a bit into the surface causing a compression between layers. There are 5 layers of 4x4 metal atoms with only the bottom layer fixed in place.

I hope I've provided enough information, if not please ask!

My input/output files are attached.

Thanks for your time,
I would really appreciate any help/guidance you can provide.

Matt

Issues like this are usually related to small problems with the setup of the calculation. (you don't need to consider the elast tag).

Often, people will minimize the endpoints with different computational settings as the band. You must make sure that you use the same number of k-points and plane wave cutoff. Also, check that you have some frozen atoms which are in constant positions for every image along the band, including the endpoints. With 5 layers, it is typical to free the bottom two layers.

I can't tell the exact problem from your description, but check for simple geometric problems and that you endpoint INCAR settings are the same as for the NEB (except for the NEB settings).

If you can't get to the bottom of it, you can post a calculation and I'll take a look.

Thanks for your reply.

The only discrepancies I found between the input for the NEB and the endpoints was the use of the IDIPOL = 3 dipole correction. I had not used this before in the NEB, so now I have included it and it has made no difference...

Would you be able to take a quick look at my files?

I see my previous attempt to upload them failed, can I dropbox or email them to you?

Best wishes,
Matt

I've fixed the server so that posts/attachments up to 256MB should be allowed. Dropbox would also be fine.

Here's the initial and final state calculations along with the NEB. There's a few post processed figures/visualisations for materials studio in there too...

Thanks again.

Matt

I would start by rerunning your final state with vasp.5.3.3 instead of 5.2.12, so it is consistent with the Initial state and the NEB calculation. Something is different about the final state calculation; the bulk is somewhat expanded compared to the initial state.

I also noticed that you have atoms at the bottom of your slab which are not frozen, and other atoms close to the surface which are (due to the rotated periodic geometry that you choose). You might want to make a better choice of those constraints at the bottom of the slab.