Transition state calculation for wealy bound molecules

[microomran]

Dear Dr.Graeme,

I am T.S. for different elementary reaction on Ni (111) .While most of calculation are going well and confirmed with one imaginary vibration,those of weakly bound molecules (namely CH4 and CO2) are not going well .

CH4:
After 150 iteration using QM:
0 0.000000 -200.038300 0.000000
1 0.076095 -200.030600 0.007700
2 0.218902 -197.373200 2.665100
3 0.262813 -195.285600 4.752700
4 0.157157 -195.424700 4.613600
5 0.217225 -195.084300 4.954000
6 0.252924 -197.042300 2.996000
7 0.412879 -197.557800 2.480500
8 0.176803 -197.794100 2.244200
9 0.009294 -197.983800 2.054500

Then switching to dimer (as below 0.5 eV) I got some positive curvature that persists after about 400 iteration:
Step Force Torque Energy Curvature Angle
1 0.99496 43.43601 -194.76306 22.99635 61.11114
1 0.99496 41.38738 -194.76306 -24.38334 46.97009
1 0.99496 8.24616 -194.76306 -18.51177 10.10026
1 0.99496 6.46426 -194.76306 -20.71544 3.24645
2 0.78132 7.54465 -195.38218 -3.71035 12.77097
2 0.78132 24.48226 -195.38218 -46.18370 36.89431
2 0.78132 13.16304 -195.38218 -32.24095 18.45483
2 0.78132 11.05970 -195.38218 10.83804 8.76179
....
20 1.07599 0.58094 -195.22214 0.05091 1.44860
21 1.16995 0.91533 -195.19354 0.07023 0.78856
22 1.23400 0.49295 -195.16451 0.07168 1.38721
23 1.32429 0.81208 -195.13735 0.05768 0.69044
24 1.37741 0.42034 -195.11009 0.02901 1.24705
25 1.45780 0.71129 -195.08484 0.03989 0.58221
26 1.50184 0.35152 -195.05947 0.06210 1.26089
27 1.57108 0.63423 -195.03676 0.06528 0.51063

Similarly with CO2 ,I have even got high forces when starting calculation that may not converge (after 50 iteration QD) :
0 0.000000 -199.171500 0.000000
1 16.472382 -197.844800 1.326700
2 0.809860 -198.928400 0.243100
3 2.732337 -194.419200 4.752300
4 6.463525 -193.745700 5.425800
5 7.247404 -194.691200 4.480300
6 7.349327 -195.932400 3.239100
7 7.528407 -196.360200 2.811300
8 7.612267 -196.832200 2.339300
9 0.009683 -197.728900 1.442600

For both cases I tried to repeat initial co-adsorption calculation but no significant change.I attached the INCAR files for both cNEB and dimer.
My current approach is to reduce forces <0.5 eV using QM and then switch to dimer.It worked with all other species.
Any advice to adjust those calculations?

Thank you in advance

[graeme]

Before getting into the dimer calculations, the NEB calculations do not looks right. I see jumps of over 2 eV between images in the first calculation and oscillations in the energy and high forces in the second. Check these paths to see if they are reasonable. I would not be surprised if you are looking at molecular dissociation in the gas phase. I don't have the data to look at the pathways, but there is some problem with the paths even before you are starting the dimer calculations.

[microomran]

Hi Dr.Graeme

Yes,as you expected it is gas dissociation.Specifically ,I am facing the problem with H2,CH4 and CO2 which are physically adsorbed with week adsorption energies(-0.02 to -0.06 eV).The two main problems -as you mentioned- are either I got is either jump between image (+2eV) or high forces after my 1st 50 iterations using QM/cNEB.I tried to repeat co-adsorption in case of CO2 dissociation 2 times and start by QM again but forces still high.

The issue that I have used the same INCAR setting /methodology with all other molecules and it worked well ,confirmed by only 1 imaginary vibration frequency.Should I repeat Coadsorption again and redo cNEB? ; Can I reduce the number of images that I usually use (8 images) for the purpose of troubleshooting?

I attached that the data for your kind advice.

Thanks

[graeme]

The problem is your initial paths. You have mechanisms in which molecules are being broken in the gas phase and then a fragment adsorbs on the surface (in the case of CO2 and CH4) or just left in gas phase (H2). If you are trying to understand dissociative adsorption, you need to find stable final states for the products on the surface. To give some examples, you can look at publications #6 and #24 on http://theory.cm.utexas.edu/henkelman/pubs/ .

[microomran]

Dear Dr.Graeme,

Thanks for your reply and share those two similar examples.One main idea-as you said- I got is that my final states are not stable on surface ((although they are converged).I will try fix that to see how it will work.

Thanks

[microomran]

Dear Prof.Graeme,

Updating our previous discussion,I have tried to ensure final configuration for these gas dissociation reaction.

H2:
While it worked for H2 dissociation H2-->H*+H* with 0.16eV in agreement with literature expectation ,I could not confirm the transition state by either dynamic matrix or VASP IBRION=5.I tried to use all non-frozen atom in case of IBRION and change setting(for DISPLACECAR and INCAR) in case of dynamic matrix but always 3 imaginary images persists.

CH4:
CH4-->CH3+H
After 200 steps with QM I got
0 0.000000 -200.036100 0.000000
1 0.019338 -199.853800 0.182300
2 0.047769 -199.766900 0.269200
3 0.110265 -199.533300 0.502800
4 0.148254 -198.479200 1.556900
5 0.073031 -198.744400 1.291700
6 0.061109 -198.789900 1.246200
7 0.009320 -199.254300 0.781800

As forces become <0.5 eV I switched to dimer but it is not converging after 1200 step.Forces are fluctuating and curvture are not consistently negative .
I am thinking to use optimizer (Can I use after I run neb2dim the vfin.pl results folder or should I repeat first 200 steps cNEB again?)

What is your kind advice?

Thanks you

[graeme]

I don't see the NEB calculations here. Start with a 'standard' NEB calculation using about 5-8 images. These are relatively long paths because the initial state has a gas phase molecule. Make sure that the climbing image is fairly close to the transition state with a force on the order of 0.1 eV/Ang. Then you can switch to the dimer calculation, or continue the NEB calculation until the force is below 0.01 eV/Ang. If you are too far away from the saddle, the dimer will roam the landscape looking for a saddle rather than converging directly to one that is nearby.

[microomran]

Thank you very much for your answer and patience.

I attached the initial cNEB calculation here.results are shown in Results subfolder.

I came across
viewtopic.php?f=2&t=3252&p=9530&hilit=d ... uble#p9530

it seems he faced the same problem too.

As per that previous case, I continue using cNEB IOPT=3.Is that O.K?

[graeme]

Yes, exactly. Try to get the force of the 04 image lower in the NEB calculations - below 0.1 eV/Ang. If you have the resources, you can then continue running the NEB calculation using a more aggressive optimizer, such as IOPT=1 (LBFGS), or switch to the dimer method.

The good news is that your path does now look reasonable.

[microomran]

Thanks for clarification and good news

I tried to reduce the residual forces to less than 0.1 eV/A as below (attached)
0 0.000000 -200.036100 0.000000
1 0.009956 -198.167800 1.868300
2 0.041363 -198.085300 1.950800
3 0.095117 -197.881700 2.154400
4 0.089032 -196.932200 3.103900
5 0.091099 -197.228300 2.807800
6 0.093513 -197.288400 2.747700
7 0.009320 -199.254300 0.781800

However,the saddle point barrier jumped to 3.1 which I do not know why it happened.

Now,I will switch to IOPT=1 (LBFGS) and see how it will go..

[microomran]

Please find attached calculations

[graeme]

The increase in the barrier is because you are using different settings in the calculations of the endpoints as for the NEB. If you reconverge the endpoints with the lower energy cutoff (and any other differences to the NEB calculation) then your barriers will be reasonable. The calculation appears to be converging nicely.

[microomran]

Dear Prof.Graeme,

After reach below 0.1 force ,I switched to IOPT=1(LBFGS) as discussed ,After 450 steps the forces seems are not going down or very very slow (attached) .Should I consider FIRE or try DIMER better?

Thanks for your help
Omran

[graeme]

If you look at your path, you can see the the CH4 molecule is dissociating over a different surface atom than you were anticipating in your final state. If you minimize the structure in state 04 you will find a different final state. The path connecting that state and the final state you choose will involve diffusion of the CH3 group.

To resolve this path, you can break it up into the two steps. After you find the intermediate minimum, separately calculate the dissociative adsorption step to a minimized 04 structure. Then, optionally, calculate the CH3 diffusion step to your chosen final state.

Use IBRION or IOPT 3 with a sufficiently small POTIM for the initial paths. Then, you can continue with a larger POTIM when the forces drop below about 0.2 eV/Ang. If you choose to switch to IBRION or IOPT = 1, make sure that you have a sufficiently small POTIM or INVCURV. The calculation that you post is unstable, as you can see by the spikes in the force.

[microomran]

Dear Prof.Dr.Graeme,

Thank you very much for your suggestion,

(a)I have put your latest suggestion about intermediate convergence into action and waiting the results..just to clarify (please correct me if I am wrong):
1-I will minimize image 04 and consider it as my final state instead of my earlier suggested final state (Coadsorption)
2-Do cNEB calculation between CH4 (IS) and newly converged final state ->get T.S and energy barrier from that path.

However,I did get the purpose of second part "Then, optionally, calculate the CH3 diffusion step to your chosen final state." and how it change the result (either T.S. or barrier)..Can you elaborate?

(b)As you suggested earlier,after I brought the forces to below 0.1 eV/Ang using QM ,I switched to Dimer(instead of going for cNEB again).
The dimer calculation went O.K until convergence.However ,I got a strange result with T.S. barrier -5.15 eV ! (Yes,-ve) .The result was strange (not logical),so I repeated calculation again but I got the same results(attached). Can you help me to understand what happened?

Thank you for your suggestions and help

Regards,
Omran