The effects of different protonation states of the hydroxamic acid (HA) inhibitors against the class I histone deacetylase enzymes (HDACs) have been studied using the state of the art quantum polarized ligand docking (QPLD) and molecular mechanics-generalized Born surface area (MM-GBSA) approaches. The binding modes of the inhibitors and their inter-molecular interactions with class I HDACs, in response to the protonation states of the inhibitors, are explored. Our results indicate that the different protonation states of the inhibitors exhibit significant differences in their interactions with the catalytic zinc metal ion and the other active site residues in the HDAC enzymes, which in turn affect the 'Histidine-Aspartate' charge relay mechanism. The QPLD calculations show that the protonated states of the inhibitors display higher scores in all the class I HDACs in this study, while the deprotonated forms present lower scores. The molecular electrostatic potentials and the other physico-chemical descriptors support the results. The MM-GBSA approach employed in the present work has been able to accurately calculate the relative binding free energies of the neutral and the protonated HA inhibitors; those were close to the experimental values. However, the MM-GBSA approach breaks down while calculating the binding free energies of the deprotonated inhibitors, which resulted in unrealistic values. Large energetic differences were found in the polar electrostatic solvation energy terms and the coulombic contributions in the deprotonated inhibitors. Thus improvements in the present solvation models and force fields become inevitable for the inclusions of charged states of inhibitors in computational drug discovery.