Based on the resonating valence bond theory, the linear combinations of the main orbits occupied by liquid hydrogen electrons are selected as the basis sets to construct the Jastrow antisymmetrized geminal product. The resonating valence bond (RVB) wave function which takes into consideration electron correlation effects provides lower energy than the local density approximation (LDA) function. In order to improve the nodal accuracy of the variational trial wave function, the backflow correlations are suggested to be employed whenever rs1.75 or T 15000 K, the improved wave function has about 1 mHa/atom decrease in local energy with respect to the one without backflow effects at the VMC level, and has a lower variance simultaneity. After combining the coupled electron-ion Monte Carlo (CEIMC) method with the RVB wave function, the simulation results we have obtained are in good agreement with the experimental and other ab-initio ones; the deuterium principal Hugoniot curve passing through the error bars of various existing experiments conducted via different high-pressure technologies has a maximum compression of 4.48 at about 50.3GPa, but the phenomenon of apparent increase in compression ratio along the Hugoniot between 100120 GPa has not been found. The RVB wave function discussed in this paper when adopted the CEIMC method is not only quite suitable for the simulation of liquid hydrogen within a wide range of density and temperature (1.0 rs2.2, 2800 K T60000 K), but also can give some more applicable thermodynamic properties of hydrogen under shock loading.