We developed a molecular simulation method suitable for estimation of binding free energy, called the filling potential method, based on the concept of the Taboo search, which is a type of self-avoiding random walk consisting of a cycle of local-minimum searches and transition state searches. The filling potential method is an umbrella potential sampling method, and enables the ligand molecule to drift from its local minima automatically. In the case of the filling potential method, the umbrella potential is a combination of Gaussian-type repulsive potentials, which are located on the trajectory of the ligand. Without setting the reaction coordinates a priori, this method searches for and determines the suitable reaction coordinates by successive generation of umbrella potentials based on its trajectory analysis. The weighted histogram analysis for these trajectories gives the binding free energy of the ligand to the receptor protein. It was applied to the complex of thermolysin and its inhibitors in explicit water, and the free energy surfaces with the stable binding state and the energy barrier were examined. The calculated binding free energies agreed well with the experimental results.