Large energy differences across successive geometry optimization steps in SIESTA

I am trying to calculate adsorption energies of acetophenone on a 6x6x3 nickel surface with about 12 Angstrom of vacuum above it using SIESTA 5.

I have arranged the molecule so that the ring plane is parallel to the surface plane, and the atom of the molecule closest to the top is about 3.5 Angstrom away.

I have been running several tests, varying parameters and geometric constraints, always keeping the nickel fixed, but I always obtain results that do not seem energetically reasonable to me because I notice too much energy difference between steps.

The energy values I observe for the first steps are

Target entalphy Free energy
-498938.1929 -499060.7102
-498940.3767 -499068.3973
-498941.1644 -499073.3934
-498941.2172 -499073.0439
-498942.9969 -499072.1696
-498942.3203 -499064.1844
-498943.1381 -499070.8697
-498944.7343 -499073.6971
-498944.7837 -499073.5593
-498944.9385 -499073.9446
-498946.0347 -499075.1088

The values according to the output are in eV units. Since these are weak interactions, the energy changes between the different steps seem exaggerated to me.

I have inspected the intermediate geometries and they seems to be reasonable, it is hard to tell the difference by simple inspection. The moving distances are reasonably small and I started the process from an optimized acetophenone.

The data that I believe are most relevant from the simulation are

XC.Functional GGA
XC.Authors PBE
DFTD3 true
PAO.BasisSize DZP
Mesh.Cutoff 175 Ry
Spin non-polarized
%block LatticeVectors   #In Angstroms
14.906518   0.0      0.0
7.4532590  12.909423 0.0
0.0         0.0     30.0855605
%endblock LatticeVectors
MD.VariableCell FALSE


Given the size of the SYSTEM (see the lattice vectors), I understand that using only 1 k-point should give at least reasonable results. Is this correct?

To choose the Mesh.Cutoff I analyzed the convergence both when calculating the surface and the molecule, and I was moderately generous. The value is quite larger than the one suggested in the pseudopotentials.

What I am missing here?

• I run some simulation using Ni, and the lower value of mesh cutoff I work with is 500Ry. I think that 175Ry is too low. Also, are you optimizing the system?
– Camps
Dec 6, 2023 at 17:11
• @Camps Thank you! Maybe I considered converged values that are not trully converged. I have not a clear criteria. Do you think that this may be the problem? In such case, could you provide some advice in an answer? Dec 6, 2023 at 17:39
• I am optimizing the geometry (relaxing) the system with the Niquel atom fixed. Dec 6, 2023 at 17:40

In order to calculate the adsorption energies of cetophenone on a nickel surface, you need to do three calculations.

Considering the same cell parameters:

• Run a calculation with only the nickel surface. Keep the optimized energy as $$E_{Ni}$$.
• Run the calculations with only the cetophenone molecule. Keep the optimized energy as $$E_{cet}$$.
• Run the calculation with the "complex": nickel surface + cetophenone. Keep the optimized energy as $$E_{Ni+cet}$$.
• Calculate the binding or adsorption as follows: $$E_{bind}=E_{Ni+cet}-E_{cet}-E_{Ni}$$.

Notes:

• In order to get more precise values, you need to run the calculations using the BSSE correction.
• Run the calculations in the other two possible conformations: with the cetophenone perpendicular to the nickel surface with the O-CH down and up.
• For this latter conformation, I imagine you will need a bigger cell, keep in mind this if you want to compare all the results.

Also, take a look at this question/answers.