Binding energy with ONIOM

Question

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 !

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.