Transition state calculation for wealy bound molecules

Vasp transition state theory tools

Moderator: moderators

microomran
Posts: 28
Joined: Fri Aug 03, 2018 7:06 am

Transition state calculation for wealy bound molecules

Post by 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
Attachments
DIMCAR.txt
(3.79 KiB) Downloaded 855 times
graeme
Site Admin
Posts: 2256
Joined: Tue Apr 26, 2005 4:25 am
Contact:

Re: Transition state calculation for wealy bound molecules

Post by 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
Posts: 28
Joined: Fri Aug 03, 2018 7:06 am

Re: Transition state calculation for wealy bound molecules

Post by 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
Attachments
H_H_discn.zip
(76.86 MiB) Downloaded 1154 times
CO_3rdtrial.zip
CO2
(5.61 MiB) Downloaded 875 times
CH4_trial.zip
(27.46 MiB) Downloaded 823 times
graeme
Site Admin
Posts: 2256
Joined: Tue Apr 26, 2005 4:25 am
Contact:

Re: Transition state calculation for wealy bound molecules

Post by 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
Posts: 28
Joined: Fri Aug 03, 2018 7:06 am

Re: Transition state calculation for wealy bound molecules

Post by 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
Posts: 28
Joined: Fri Aug 03, 2018 7:06 am

Re: Transition state calculation for weakly bound molecules

Post by 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
Attachments
H_H_Compressed.tar.gz
(761.75 MiB) Downloaded 1126 times
CH4_Converge.tar.gz
(13.24 KiB) Downloaded 824 times
graeme
Site Admin
Posts: 2256
Joined: Tue Apr 26, 2005 4:25 am
Contact:

Re: Transition state calculation for wealy bound molecules

Post by 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
Posts: 28
Joined: Fri Aug 03, 2018 7:06 am

Re: Transition state calculation for weakly bound molecules

Post by 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?
Attachments
cNEB_CH4.tar.gz
(56.16 MiB) Downloaded 1120 times
graeme
Site Admin
Posts: 2256
Joined: Tue Apr 26, 2005 4:25 am
Contact:

Re: Transition state calculation for wealy bound molecules

Post by 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
Posts: 28
Joined: Fri Aug 03, 2018 7:06 am

Re: Transition state calculation for weakly bound molecules

Post by 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
Posts: 28
Joined: Fri Aug 03, 2018 7:06 am

Re: Transition state calculation for wealy bound molecules

Post by microomran »

Please find attached calculations
Attachments
001_NEB02.tar.gz
(33.53 MiB) Downloaded 1085 times
graeme
Site Admin
Posts: 2256
Joined: Tue Apr 26, 2005 4:25 am
Contact:

Re: Transition state calculation for wealy bound molecules

Post by 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
Posts: 28
Joined: Fri Aug 03, 2018 7:06 am

Re: Transition state calculation for weakly bound molecules

Post by 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
Attachments
NEB04_image4.png
NEB04_image4.png (105.19 KiB) Viewed 27386 times
NEB04.tar.gz
(45.09 MiB) Downloaded 1056 times
NEB03.tar.gz
(23.61 KiB) Downloaded 750 times
graeme
Site Admin
Posts: 2256
Joined: Tue Apr 26, 2005 4:25 am
Contact:

Re: Transition state calculation for wealy bound molecules

Post by 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
Posts: 28
Joined: Fri Aug 03, 2018 7:06 am

Re: Transition state calculation for weakly bound molecules

Post by 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
Attachments
CH4.tar.gz
(106.03 MiB) Downloaded 1009 times
Post Reply