I read in a paper entitled" http://www.tandfonline.com/doi/abs/10.1080/08927022.2011.561430#.UjAD1TZmWSo". However i couldn't understand from this paper as i am new to MD and Autodock
In order to do it right and computationally cost efficient, is to combine your md trajectory with a mm-pbsa or mm-gbsa method.
To be more precise you can write the binding energy as
Delta G= Gab-Ga-Gb
Where Gab, Ga, Gb are the absolute binding energies of comlex, ligand and receptor, respectively.
It can be written, furthermore, as
Delta G=Delta G (gas) - Delta Delta G (solvation)
Where the first term is the change in absolute binding energy due to complex formation in gas phase, which includes the entropic contribution, and the second term is the change in the solvation energy, which inludes the change in energy for creation of a cavity of size of your molecule in solvent (and so the reorganization entropy of solvent), plus the electrostatic term calculated using either pb or gb method.
My experience is that the only contribution that shouls be treated carefully is the configational entropy, which includes the rotational, translational, conformational and vibrational terms. The first two terms are treated separately, and require in additional two more simulations of the ligand and receptor separately in solvent. The other two terms may be treated either separately or together. In the last case a quasi harmonic approximation is used. When treated separately, you have to consider different possible conformations due to bond rotations of ligand and the vibration term using normal mode analysis.