I was performing MD simulation for protein-ligand complex following GROMACS tutorial.
I compared the energy of two ligands (S1 and S2) to a protein, using GROMACS LJ-SR (Lennard-Jones Short Range) and Coul-SR (Coulombic Short Range) option.
Below is the plot.
As you can see the LJ-SR is complete opposite of Coul-SR. Are they expected to be the opposite? Which is more appropriate choice to represent real binding energy?
Obi-Wan voice: "These are not the energies you're looking for."
If I understand your goal, you're trying to compare two ligands and determine which is better at binding to a given protein (and, ideally, quantify "how much better" one is compared to the other). If that is the case, then the quantity you're interested in is the relative free energy of binding, $\Delta\Delta G_\text{bind}$. The tutorial you're using explicitly states (regarding the energies you plot above):
It is important to note that this quantity is NOT a free energy or a binding energy. In fact, most force fields are not parametrized in such a way that this quantity is actually physically meaningful.
The energies you're looking at here only provide a small part of the total picture for the free energy. Among other things, entropic concerns can dwarf the pure energetics. The quantities you're looking at here may be useful for other things: for example, you might be able to characterize different binding poses based on these energies, and certainly bound vs. unbound would give significantly different values for these.
The simplest way to get binding free energies from MD simulations is probably to use MM-PBSA (or a related method). Here is a review article from a few years ago on those methods. These require some additional calculations, but a tool for doing MM-PBSA with Gromacs exists: https://rashmikumari.github.io/g_mmpbsa/.
A more exact way to calculate the relative free energy of binding would be alchemical free energy simulations, which combine MD simulations from non-physical ensembles to estimate the difference in free energy between different states. Here is a "living review" (intended to be regularly updated) of best practices for that kind of simulation; I recommend that as a starting point.
Another approach is to calculate the dissociation rate ($k_\text{off}$), using rare events methods, and use the assumption that $k_\text{on}$ is typically near-constant for reversible binders. In my experience, this tends to be a more challenging calculation than the alchemical approach, but it may be more forgiving if your initial binding pose is not perfect (depending on specifics of the algorithm used).
Just complementing the @dwhswenson answer, when studying ligand/protein binding process, there is no best energy to look for.
In the binding process there are several mechanism that will influence it. Knowing the influence of such and each mechanism, you can design your ligand accordingly to increase the binding.
One of the advantage of using force fields is that the energy can be "partitioned" in several terms each one representing one interaction mechanism.
For example, the GLIDE score function (a roughly free binding energy) is partitioned in the following terms:
Note: GLIDE is a docking software.
The CaFE plugin calculates the binding energy for a ligand/protein complex using the MM/PBSA method and returns the contribution of several terms (Elec, Vdw, PB, SA, Gas, Sol, Pol, Npol).
