Issue with calculating adsorption energy

[yygao7]

Dear Shankar,

I attempted to calculate the adsorption energy using both VASP and JDFTx separately. The calculation process is as follows:

Fads = F( NO ) – F( ) – F( NO ) = -435.2356455 – (-409.2763019) – (-25.94538336) = -0.37987314 Eads = E( NO ) – E( ) – E( NO ) = -626.00511 – (-613.86243) – (-12.306743) = 0.164063 Where F and E represent the energies calculated by JDFTx and VASP in a vacuum, respectively.

The structure computed by JDFTx was optimized using VASP. It can be observed that there is a significant difference in the results obtained from the two software packages. Additionally, from the perspective of the adsorption configuration, the adsorption energy should be positive.

I have attached the relevant documentation for your reference. I wonder if there are any issues in my JDFTx input file. I am looking forward to your suggestions.

Thanks for your time!

Best, Cathy jdftx-gas.zip jdftx-gasonslab.zip jdftx-slab.zip

[ximik69]

Hi Cathy, I looked your files through and have several suggestions. First of all one can compare the results of different programs only if calculating under the same conditions. There is no info about parameters in vasp. There one should pay attention to exclusion or compensation of interactions between layers, incl. forming dipoles and compensation of them. What I have seen in jdftx input files - NO is optimized but slabs with it or w/o it are not. I didn't check whether distances d(N-O) are exactly the same for isolated and adsorbed molecule. Everything should either be optimized or not. Mixing opt with noopt can give epic differences.

If I would calculate adsorbtion, I'd add coulomb-truncation slab 0 0 1 and coulomb-truncation embed 0 0 0.5 to slabs.

Cheers.

[yygao7]

Thank you for your kind advice!

As you suggested, I added the coulomb trunaction commands to the input file, and recalculated the energy of the optimized structure. However, when I calculating the energy of the slab with gas NO, the calculation crashed with the following error:

"---------- Setting up coulomb interaction ---------- Setting up double-sized grid for truncated Coulomb potentials: R = [ 35.742 -17.8721 0.104129 ] [ 0 30.949 0.11749 ] [ 0 0 75.6497 ] unit cell volume = 83682.1 G = [ 0.175793 0.101515 -0.000399633 ] [ 0 0.203017 -0.000315303 ] [ 0 -0 0.0830563 ] Chosen fftbox size, S = [ 144 144 320 ] Integer grid location selected as the embedding center: Grid: [ 0 0 80 ] Lattice: [ 0 0 0.5 ] Cartesian: [ 0.0260322 0.0293726 18.9124 ] Constructing Wigner-Seitz cell: 14 faces (6 quadrilaterals, 8 hexagons) Range-separation parameter for embedded mesh potentials due to point charges: 0.628518 bohrs. Lattice direction 001 is not perpendicular to the other two basis vectors. End date and time: Fri Dec 22 12:36:16 2023 (Duration: 0-0:00:00.55) Failed."

However, it didn’t occur when calculating the energy of the slab.

Would you please give me some hint where my computations might have gone wrong?

Thank you for your time!

Best, Cathy gasonslab-coulomb.zip

[ximik69]

It seems you have z-axis almost but not exactly perpendicular to other axes.

If it is isolated slab, which is the case, you can adjust by hand z axis in your NO-on-slab system. Good question whether it slightly changes geometry of slab+NO or not. Or how to change the unit cell without changing geometry of slab and adsorbate. At least angles are very close to 90 thus I think you can neglect changes. Just export model to vesta or diamond, change unit cell angles to 90 and then insert data into jdftx input file.

Hope that helps.

Best, Igor.

[shankar1729]

Hi Cathy,

Igor's suggestions on adding coulomb truncation would definitely help improve the accuracy of your simulation. However, likely your energy difference is emerging from the gas-phase NO calculation in JDFTx. You don't have any spin settings set, so you are calculating an unpolarized configuration with a higher energy than the ground state with an unpaired spin.

Make sure you set spintype z-spin and elec-initial-magnetization 1 yes to constrain to one unpaired spin for the gas-phase NO. (You could also leave the magnetization unconstrained and use fermi smearing, but I would recommend constraining the spin state for molecules when you know them explicitly.)

Best, Shankar

[yygao7]

Dear Igor,

Thanks for your advice!

I will adjust the z axis to ensure that the third direction is exactly orthogonal to the first two and recalculate again.

Thanks again for helping me.

Best, Cathy

[yygao7]

Dear Shankar,

Thank you for pointing out my problem.

I've encountered another issue now. When I do fixed potential calculations, the ElecMinimize can’t converge with the following error:

“ElecMinimize: Step failed along negative gradient direction. ElecMinimize: Probably at roundoff error limit. (Stopping) Setting wave functions to eigenvectors of Hamiltonian”

But the calculation in vacuum and solution converged well.

The slab contains iron. I suspect it might be due to the absence of elec-initial-magnetization settings in my calculations, or could there be other issues in my input file? I am looking forward to your advice.

Thanks in advance!

Best, Cathy constantpotential.zip

[shankar1729]

Hi Cathy,

Your smearing is likely too small for the chosen supercell and k-point sampling, leading to a bumpy energy landscape for optimizing the occupations. Use a larger smearing ~ 0.01 Eh, possibly with the Cold or MP1 options if you are concerned about the electronic entropy being unphysical.

Unrelated, your chosen number of processes (64) is not efficient for a calculation with nStates = 10. Please see the Getting started page and numerous other threads here on how to optimally use the hybrid MPI + threads parallelization.

Best, Shankar

[yygao7]

Dear Shankar,

Thanks for your reply, and it helps a lot!

I changed the smearingWidth to 0.01 Eh and it converged successfully. But I'm puzzled as to why, when I replace the Fe element with other elements such as Co, Ni, Cu, the convergence can be achieved using the same smearing width as before. The supercell and k-point sampling are the same.

Additionally, thank you for addressing the issues I encountered with parallelization. I have now utilized the -n and -c options in accordance with the instructions on the 'Getting Started' page to enhance the speed of my calculations.

Best, Cathy

[shankar1729]

Glad to hear that: the lack of convergence is difficult to predict. The low smearing can create a problematic energy landscape, but it depends on the electronic structure near the Fermi level. Additionally, even in a problematic energy landscape, you may not always encounter the convergence difficulty.

Best, Shankar

[yygao7]

Dear Shankar,

Thank you for all your responses; they are very helpful in guiding me through the use of JDFTx.

Have a nice day!

Best, Cathy

[yygao7]

Dear Shankar,

I attempted to calculate the adsorption energy of NO under different potential using JDFTx. For example, at -1 V and 1 V respectively, the calculation process is as follows: ΔΦ ads（-1 V） = Φ( NO ) – Φ( ) – F( NO –μ𝑁 ) = [-902.64 – (-878.1127723) – (-25.95870932 + 11 0.1345) ] 27.2 = -1.312850102 ΔΦ ads（1 V） = Φ( NO ) – Φ( ) – F( NO –μ𝑁 ) = [-876.6860627 – (-852.9383807) – (-25.95870932 + 11 0.208) ] 27.2 = -2.09450377 The results show that the higher the potential, the more negative the adsorption energy. This result is very strange. I wonder if there are any issues in my calculation. I am looking forward to your suggestions. Thanks for your time! Best, Cathy

[shankar1729]

Hi Cathy,

Could you expand on why you state that this result is strange? Couldn't the nominally negative NO adsorbate bind more strongly to a more positively charged electrode at a higher potential?

Best, Shankar

[yygao7]

Dear Shankar,

I also calculated the adsorption energy of H under different potentials, the calculation process is as follows: ΔΦ ads（-1 V） = Φ( H ) – Φ( ) – [ F(H2)/2 –μ𝑁 ] = [-878.57 – (-878.1127723) – (-1.165788444/2 + 0.1345) ] 27.2 = -0.26 eV ΔΦ ads（1 V） = Φ( H ) – Φ( ) – [ F(H2)/2 –μ𝑁 ] = [-853.356976 – (-852.9383807) – (-1.165788444/2 + 0.208) ] 27.2 = -1.19 eV

The change of H adsorption energy with potential shows the same trend as NO. In theory, the adsorption energy of a hydrogen atom with a positive charge should be more negative at a lower potential. I'm not sure whether these results are reasonable or not. I have attached the relevant documentation for your reference. slab.zip slab+H.zip

Could you please point out any mistakes in my calculations or understanding. I‘m looking forward to your advice.

Thanks for your time! Best, Cathy

[yygao7]

Dear Shankar,

In addition, I would like to confirm whether eq(1) or eq(2) is correct when calculating the change of free energy during NO hydrogenation, as I found that the formulas in different literatures are not uniform. ΔΦ ads（1 V） = Φ( NOH ) – Φ(NO) – [F(H2)/2 –μ𝑁] = [-877.0534895 – (-876.6860627) – (-1.165788444/2 + 0.208) ] 27.2 = -2.09 (eq1) ΔΦ ads（1 V） = Φ( NOH ) – Φ(NO) – [F(H2)/2 –μ𝑁] + eφ = [-877.0534895 – (-876.6860627) – (-1.165788444/2 + 0.208) ] 27.2 + 1 = -1.09 (eq2)

Thanks again for your time! I am looking forward to your reply.

Best, Cathy

[shankar1729]

Hi Cathy,

I think equation 1 vs 2 correspond to adsorption of H* vs H+; it's the reference state that is changing between them. Perhaps, you are counting the energy to transfer the e- in one of the cases as well that leads to the issues you are finding?

Best, Shankar

[yygao7]

Dear Shankar,

Thank you for your reply. That is to say, if I want to calculate the adsorption of H+ in aqueous solution, I should use eq(2), right?

Best, Cathy

[shankar1729]

Essentially, you want Φ( NOH ) – Φ(NO) - Φ(H+). Note that you don't need to balance electrons since you are in the grand canonical ensemble.

Now, Φ(H+) = Φ(H2)/2 at SHE by definition, and because H+ has no electrons, Φ(H+) is independent of electrode potential (at the same pH). Hence, Φ(H+) = F(H2)/2 –μ(SHE)𝑁 (with 𝑁=1) at all potentials. Finally, since μ = μ(SHE) - eφ by definition of φ, you end up with eq. 2 written above.

[yygao7]

Dear Shankar,

Thanks for your reply！

I seem to understand… So if I want to calculate the adsorption of NO by Φ( NO ) – Φ( ) – Φ(NO), then Φ(NO) = F(NO) –μ𝑁 (𝑁 = 11), with μ = μ(SHE) – eφ. This is different from calculating H+, am I right?

Best, Cathy

[shankar1729]

Yes, that's right. And this would correspond to NO gas at standard temperature and pressure as the reference state for NO.