The role played by water in biological processes spans over multiple scales. At the macromolecule scale, water plays an essential role in determining their structure, stability and function . It is also now widely acknowledged that water molecules have a significant effect in ligand binding, if not sometimes the main driven force in biomolecular recognition . When a small molecule binds to its target, it also have to compete with water molecules present in the binding site. And depending on the context, water molecules will be either retained or displaced, favoring or opposing ligand binding . The displacement of tightly bound water molecules leads to a gain in entropy and the overall displacement process is favorable when the ligand compensates for the loss of the water-protein interaction enthalpy, but it does not seem to be always the case . This has some important implications in lead-optimization and how functional groups are affecting water positions . But in some cases, water molecules are so tightly bound and conserved across family members, that they might be considered as part of the binding site. In addition to their individual contributions, they also need to be considered as a network. Water molecules present in binding sites are forming networks through hydrogen bonds, and the removal of one water molecule or more have consequences on the position and free energy of the remaining ones in a non-additive way . Thus, it is important to consider water molecules as a whole, and not as discrete and isolated entities, in which each individual water plays a role but depends on its interactions with the others in highly cooperative interactions. In the context of ligand binding, this has a direct effect on its binding, when reorganization or replacement of water networks can be often observed, contributing to the enthalpy/entropy . Mutations can also disrupt water networks and affect ligand binding indirectly, via concerted changes in the water network . Over the past decades, a lot of efforts were dedicated to the development of methods for identifying and characterizing water molecules in ligand binding site. This is reflected by the number of tools available, and still increasing, for such purpose. Those methods can be classified into four different groups: (1) empirical and knowledgebased like Consolv , WaterScore , AcquaAlta and Wa-terDock , (2) Statistical and molecular mechanics with GRID , SZMAP and 3D-RISM in combination with placevent or GAsol , (3) Simulation-based methods, either based on Molecular Dynamics (MD) simulations such as the post-processing methods IFST (used in WaterMap , STOW and GIST ), WATCLUST , WATsite , SPAM or based on Monte-Carlo (MC) simulations , like JAWS , grand-canonical Monte-Carlo (GCMC) simulations and Double Decoupling or by combining MC steps with MD (MC/MD) simulations , and very recently, (4) deep neural network approaches made their first appearances, promising instantaneous hydration free energies from static structures . For more detailed information about the various methods, you can refer to the following recent reviews on the subject: Bodnarchuk , Biedermannová and Schneider 52 , Spyrakis et. al. , Graves et. al. , Nittinger et. al. and Hu et. al. . More rapid methods perform acceptable placement predic-tions, but they trade accuracy for speed, therefore their energy estimates are usually only qualitative . Indeed, the accurate determination of individual water molecules position and properties requires more expensive calculations to sample multiple configurations, in part also due to the necessity of accounting for the water-water interactions essential to reproduce the correct water network. For that, methods based on MD simulations using explicit solvent models are considered as a de facto gold standard because they provide both positions and orientations sampling that can be used to perform quantitative estimate on thermodynamic properties, but are significantly more expensive than empirical methods. Moreover, conventional MD simulations have some significant limitations, such as the inability of sampling accurately water molecule positions in cavities that are not easily accessible and therefore have limited exchange with the bulk solvent (i.e., buried or enclosed pockets). This can lead to either the formation of vacuum or abnormal low density regions which lead to artifacts . Different strategies and protocols have been developed to circumvent this issue , at the cost of increased computational cost. One of the principal motivations for the development of the previously presented methods is the characterization of water molecules to guide drug design efforts, and the improvement of molecular modeling predictions. Computationally intense free energy calculation methods like Double Decoupling , Free Energy Perturbations (FEP) or Thermodynamic Integration (TI) , are used when precision and detailed thermodynamics are a priority for the study. They have been successfully applied for numerous sets of protein-ligand complexes , and such positive outcome is often linked to a proper treatment of buried water molecules after alchemical transformations, which is at the origin of large hysteresis observed during calculations . However, those methods are not yet suitable for high throughput applications such as docking. In docking applications, water molecules are very often neglected for complexity reasons and the lack of a standard protocol for their treatment. In practice, water molecules are commonly treated implicitly, with the possibility for users to chose special protocols when an explicit treatment is necessary. Efforts have been made to improve the modeling of desolvation effect during docking simulations with the goal of improving accuracy of large virtual screening libraries and reducing the number of false positives. While docking software are using different strategies, studies show in average an improvement of the docking pose accuracy and energy ranking . Some methods are using explicit water molecules during docking . Other approaches were explored by incorporating more accurate hydration information from MD simulations during docking . or as a post-processing step for re-ranking ligands . The choice of these protocols is dependent on the number of ligands to be docked, structural data availability, time constraints, and computational resources accessibility. Inclusion of more accurate and computational expensive methods into the docking pipeline is one of the main issue of a wider adoption of such methods. In particular, the need to run MD simulations prior to docking presents computational and technical challenges that can discourage common users to incorporate hydration information into their docking protocols, and limits its general applicability in a high-throughput fashion. Having an accurate description of the solvent-related energy components of the target structure during docking is essential to increase the accuracy of results, and reduce false positives. The motivation for development of WaterKit was the exploration of an alternative physics-based method to rapidly sample water configurations and generate high quality energy models of individual molecules. Such method needs to provide comparable accuracy to MD-based methods, providing full thermodynamic profiling of multiple hydration shells with much smaller computational needs. The protocol described here has been designed to be compatible and easily integrated with the AutoDock docking protocol .