Constant potential calculations for large molecules

[oupengfei1989]

Dear jDFTx developers

We are trying to calculate the adsorption energy of a large molecule on a metal surface using jDFTx, and we have already successfully reproduced similar results for vacuum and solvation calculations, if compared to our DFT calculations. However, when we add the applied potential, even though under 0 V vs SHE, we find a large discrepancy between the ones without the applied potentials.

We'd like to double-check multiple details with you, to confirm the calculations are perfectly right.

1. If the applied potential in experiments is -1.3 V vs SHE, so we should set mu = -(-1.3+4.44)/27.12 = -0.11578 V in the input file, right?
2. The formula that we used to calculate the adsorption energy is deltaGcpm = G(Adenine/Cu(100)Q3)-G(Cu100Q1)-G(Adenine*Q2)-(Q1+Q2-Q3)mu, in which the Q1, Q2, and Q3 are the net charges determined by constant potential calculations, e.g., if the Adenine gains 0.5 electrons compared to its state in the vacuum, we thought the Q3 = -0.5.
3. And G(Adenine/Cu(100)Q2), G(Cu100Q1), and G(Adenine*Q2) are the solvation free energies for combined Adenine/Cu(100) (after Adenine adsorbed on Cu(100) surface), and pristine Cu(100) and Adenine molecules in the vacuum, respectively. All the input structures are directly from the converged DFT calculations and we have not enabled the ionic optimizations in our jDFTx calculations, only the electronic SCF.
4. What is the reference energy for Adenine, should it be exactly the same energy under different conditions, this is what we have learned from the ideal gas of H2, however, we thought should be different for Adenine, so we directly took the solvation free energy from each calculation.

We really appreciate your help in resolving these issues and look forward to your reply.

[shankar1729]

To answer your questions:

1. Yes, that is the correct conversion, but note that the mu input to the code is technically -0.11578 Hartrees, not Volts.
2. No, if G are the grand free energies you read from the log or Ecomponents files, the formula is just the difference of the G's. Don't add/subtract extra mu N terms. The whole point of grand-canonical ensemble is that you can easily compute free energy differences between systems with different electron components. (Also, note that the G is not a solvation free energy, as that would be an energy difference between solvated and vacuum cases.)
3. Sure, you can do fixed geometry calculations and only do the grand-canonical SCF/elec-minimize without ionic optimization.
4. For a molecule, you could do fixed-potential calculations, as you seem to have done, or you could just do a neutral solvated calculation to get Helmholtz energy F and then compute G = F - mu N for any mu you care about. (N is the nElectrons of that neutral calculation and mu is the absolute Hartree input you gave to target-mu.)

Finally, note that the adsorption energy at fixed potential of 0V need not be close to the neutral case. This is only likely if the PZC of the surface/adsorbed config is close to 0V. Otherwise, there is nothing special about 0V, and could be a configuration which is quite far from neutral.

Hope that makes sense!

Best, Shankar

[oupengfei1989]

Hi Shankar

Thank you very much for the response, super helpful!

1. Sorry for the typo, I intend to say Hartrees, not the Volts.
2. Just for clarification, the adsorption free energy is defined as the difference of Gs between the post and prior adsorption, right? If we are talking about the electrochemical reaction that involves a H+ or any other charged species, this equation is then used to calculate the reaction energy: deltaGcpm = G(Q3)-G(Q1)-G(Q2)-(Q1+Q2-Q3)mu, right?
3. Would you expect a large difference if we directly use the structures determined from DFT calculations? Since we do not include any effects of solvation and applied potential, the adsorption behavior might be totally different. Will you recommend us to do ionic optimization if we have enough time and resources?
4. My understanding for this point is that you can directly use the G retrieved from the out files under certain conditions, just make sure these numbers are under the same conditions; or you can do a neutral solvated calculation to get Helmholtz energy F and then compute G = F - mu N for any mu you care about. However, I also find these two numbers are not the same, and I suppose they should be the same.

Thank you very much again!

Best wishes

Pengfei

[shankar1729]

Hi Pengfei,

The reaction energy should just be deltaG = G3 - G1 - G2. The electron number terms (your Q mu products) are already included in G.

And yes, if you can, self-consistent structures with solvation and potential would be preferable. In some cases, these can change the adsorption configuration. But, applied potential calculations at fixed geometry will still be better than CHE alone, so what you are doing is at least a good starting point.

On the final point, G = F - mu N, with F and N taken from the neutral calculation should work approximately only for a molecule, not for a surface. Essentially, this equation is a linear extrapolation of the grand free energy from the neutral state, which is the essence of the CHE method. So the deviations of the self-consistent G from the above case for the metal surface and metal surfaces with adsorbates is the main effect that the fixed potential calculations are capturing for you beyond the CHE.

Best, Shankar

[oupengfei1989]

Hi Shankar

Thank you so much for all these detailed explanations.

After the ionic minimization, we found the G approaches the one using the N taken from the neutral calculation, which confirms your point, and I guess it gives us confidence about our results as well.

And it seems the adsorption energy after the ionic optimization reaches a more reasonable value.

Best wishes

Pengfei

[jiangyimin123456]

Hi Shankar

I have some doubts regarding constant potential calculations in JDFTx and I'm looking for accurate information. I kindly request your guidance and expertise to help clarify my queries.

1. How to calculate the reaction free energy deltaGcpm under a constant potential using JDFTx. For example: Reaction: a+b→c+d, deltaGcpm = Gc + Gd - Ga - Gb where Ga, Gb, Gc, and Gd are the energies obtained by fixing the same potential. Is my understanding correct?
2. Setting a constant potential E (parameter mu) E_SHE = -(mu-mu_SHE), There are two ways to set mu_SHE a. calibrating mu_SHE by comparig theoretical calculations of the PZC of various metal surfaces against experimental measurements. For the CANDLE solvation model, mu_SHE = -4.66 eV (from the CANDLE paper). b. mu_SHE from Neutral.out ( this is an example from the "Fixed-potential calculations" section on the JDFTx official website)
   3. E_SHE is relative to the Standard Hydrogen Electrode (SHE). If I want to correspond it to the Reversible Hydrogen Electrode (RHE) referenced in experiments, E_RHE = -(mu - mu_SHE) - 0.059*pH, is my calculation correct?

Best wishes

yimin

[shankar1729]

Hi Yimin,

1. That's right.
2. a is correct. The mu you get from Neutral.out is the potential of zero charge (PZC) of the specific system, not SHE. Note that the tutorial you reference plots V - PZC on the axis, not V vs SHE etc.
3. Yes, that looks correct. The SHE to RHE conversion is what you would normally have to do, that is not changed by the GC-DFT formalism in any way.

Best, Shankar