Do you know how to compute binding energy using Gaussian? I have optimized binding pocket with ligand using ONIOM methodology implemented in Gaussian16. For the high layer B3LYP/6-31+G(d,p) was used, while low layer with PM6 method. Then, I recomputed the system without ligand and hence used only PM6 for the aminoacids. I though the way would be to subtract the energies of the system with and without ligand but it actually makes no sense for ONIOM total energy equals -1029.031167au and PM6 energy is -1.261231au.

Thanks in advance !

  • 2
    $\begingroup$ Take o look at the discussion here. $\endgroup$
    – Camps
    Nov 1, 2022 at 13:32
  • $\begingroup$ @Camps Thanks, I will read that $\endgroup$
    – farmaceut
    Nov 1, 2022 at 17:29

1 Answer 1


You have to also compute the energy of the free ligand at the DFT level, and subtract that too. A binding process is the reaction of a protein with a free ligand to form the protein-ligand complex, rather than a protein gaining a ligand from nowhere. So the energy of the free ligand must enter the equation of the binding energy.

Some comparatively minor but still important points:

  1. The basis set 6-31+G(d,p) is unreliable for single point energies, although it may be enough for geometry optimizations. The single point energies must be computed with a basis set that is at least a triple-zeta basis set, such as 6-311+G(d,p) (or even better, def2-TZVP)
  2. Don't forget to add the dispersion correction
  3. As you are almost certainly interested in the binding process in water, as opposed to in vacuum, you should add a solvation model
  4. Make sure to fully sample the binding conformations (for example under the PM6 level) and compute their Boltzmann average. Using a single hand-drawn binding pose is almost always unacceptable.
  • $\begingroup$ Thanks for the answer @wzkchem5 ! Just before reading it, the enlightenment came to me mind that the energy contribution from the ligand is indeed skipped! Regarding the remaining: (1) There's a great number of amino-acids the ligand interacts with, hence I restrained myself to optimization with the given level of theory. Anyway, the idea of recomputing SP with triple-dzeta might be wise... (2) It is - GD3 (3) PCM water is included too (4) How could I do that? I have performed docking in Vina and picked pose comparable to the similar ligand... however, sampling would be truly a wiser choice $\endgroup$
    – farmaceut
    Nov 1, 2022 at 17:35

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .