dimer transition state optimization problem

[jasius]

Colleagues, I have a very peculiar situation with Dimer method (IBRION=44) which up till now I thought I mastered well. I am trying to dissociate O-OH bond on Cu/Pt. I proceed normal way, e.g. calculating initial guess frequencies, taking the negative vector due to the right bond vibration and appending it at the bottom of POSCAR. Routinely, I preoptimize these structures with some lose parameters, normally ENCUT=270 eV which suffices to arrive to even better guess. Then I use that structure and reoptimize it with ENCUT=400 for a quick convergence. Below is the INCAR and files also attached

Code: Select all

GGA       = RP ! R-Perdew-Burke-Ernzerhof
PREC      = NORMAL
ENCUT     =  270  #(I switch to 400 later after this one converges to a reasonable geometry)
NELM      =   55
NELMIN    =    4
NELMDL    =   -5
EDIFF     =    0.000001
EDIFFG    =   -0.01
ISIF      =    2 ! The cell shape and the volume is fixed
ISYM      =    0 ! Symmetry is OFF
LCORR     =    T
ISMEAR    =    0
SIGMA     =    0.05
LREAL     = Auto
ALGO     =  Fast

#-----------------------------------------------------------------------------------------
#dynamic
IBRION  =  44
#NFREE   =  1
#POTIM   =  0.015
NSW     =  2500

The last optimization step with ENCUT=270 returns very favorable low forces, only shown below for the 3 atoms that are involved in TS

Code: Select all

 POSITION                                       TOTAL-FORCE (eV/Angst)
 -----------------------------------------------------------------------------------
      5.07804      3.77335     10.34551         0.001183     -0.002604      0.003790
      3.56832      3.06916     10.84629         0.004773     -0.005111      0.004338
      3.92548      2.41079     11.47076        -0.003853      0.006871     -0.002578
...................
    total drift:                                0.047841     -0.017459     -0.164061

Also, the frequency calculation of that 270 eV preoptimized structure with ENCUT=400 already returns only one negative frequency although I understand that those are calculated using different ENCUT.

Then I take that last coordinate file with its vectors and resubmit for a final ENCUT=400 TS iDM optimization and the very first step forces become unreasonable and optimization quickly "explodes" (mind you starting at the same proptimized geometry). These are the step 1 forces with ENCUT=400 from the same structure that was optimized with ENCUT=270:

Code: Select all

 POSITION                                       TOTAL-FORCE (eV/Angst)
 -----------------------------------------------------------------------------------
      5.07804      3.77335     10.34551         0.545761      0.252850     -0.125457
      3.56832      3.06916     10.84629        -0.644273     -0.108701      0.126654
      3.92548      2.41079     11.47076         0.089353     -0.154830      0.152156
..................
    total drift:                               -0.006608      0.050827      0.077685

Optimizing from ENCUT=400 (without preoptimization), switching to Accurate setting or LREAL=.F. always results in the same strange behavior: convergence at ENCUT=270 and divergence with 400 eV.
I was wondering if I could get some leads. Notably, this TS exists on Pt(111) alone and I have located using exactly the same methods.

[manuel_engel1]

Hi Jonas,

Thanks for reaching out to us. This is indeed strange and I will start looking into it. However, since this is a rather involved calculation, it might take some time to figure out. In the meantime, I would like to ask a few questions:

Optimizing from ENCUT=400 (without preoptimization)

Did you run the entire procedure of frequency analysis and IDM with ENCUT=400?
In other words, did you also take the phonon eigenvectors from a calculation with ENCUT=400?

Do you by any chance have a fast and simple test case where the same behavior occurs?

Best,
Manuel

[manuel_engel1]

Small update: We were able to reproduce the behavior and are now looking into potential causes and solutions.

[jasius]

Thank you. I normally run ENCUT 270 frequencies since that is not crucial and saves time. I am glad you reproduced it, I will wait for your response

[manuel_engel1]

We are much closer to solving the issue after running multiple simulations based on your system. However, we are still exploring some possibilities before we can give a definitive answer. Could I ask you for a bit more information?

• Do you have a structure (POSCAR) of the molecule adsorbed on the surface at 270eV, but preferably at a higher cutoff (ENCUT=400) if you have it? From the structure file you shared with us (the TS), we see the same as you do: the molecule quickly dissociates on the surface at 400eV.
• Could you please also share the structure from which you started the TS search at ENCUT=270, i.e., at the start of your pre-optimization run?
• What is the barrier height along the reaction path at ENCUT=270, i.e., how much higher is the TS in energy compared to the stable adsorption configuration?

[jasius]

Attached is the initial guess for 270 eV optimization - as you can see it is very far away from the final result of preoptimization. Note, direct 400 eV optimization of this initial structure diverges, so 270 eV is needed as it is more forgiving and leads to the right structure.
I suppose your question #1 is related to #3, so I am running IRC search (IBRION=40) with 400 eV from my 270 eV converged TS structure to show you both endpoints