A new realization of the constant-pH molecular dynamics simulation method is proposed. Molecular dynamics simulation is performed for a protein in the most probable ionization microstate of the current conformation taking into account the potential of mean force of protein molecule in the water-proton bath in equilibrium titration conditions (MD-pH-ET). It is shown that: 1) the optimal one is the simulation of the protein in the most probable ionization for a given conformation, taking into account the correction ionization potential of mean force, which results in the librium ensemble of ionization states; 2) new method MD-pH-ET allow one to carry out an optimization of protein structure and the total free energy of a protein in the aqueous solution at constant pH, and to calculate the pH-dependent properties. Method MD-pH-ET possesses the unique features: 1) it uses the most precise and computational-effective realization of calculation of the electrostatic energy of a protein in water solution, the model of continuous dielectric media with Poisson equation and the generalized Born method with "ideal" Born atomic radii; 2) it uses the same model of the potential energy surface in the ionization-conformational phase space, both for calculating the potential energy of the protein and atomic forces and for determining the most probable ionization states; 3) it calculates the total free energy of the protein in the aqueous solution in proton reservoir under the conditions of equilibrium titration. The workability of the new method MD-pH-ET is demonstrated for the molecule of protein BPTI.