Phosphoric acid (H3PO4) and water (H2O) are well-known inorganic compounds, which have caught scientists' attention for their proton transfer mechanisms. Liquid H3PO4 has the highest intrinsic proton conductivity of all materials and H2O is the most important substance in nature and environment. Therefore, molecular studies of these compounds and their mixtures provide important insights into their structures and dynamics. We present a molecular dynamics study of hydrogen bonding in liquid H3PO4-H2O systems at different concentrations using GROMACS software package. A number of different generic force fields such as OPLS, GAFF, GROMOS 54, GROMOS 87, GROMOS 96 are tested and benchmarked for suitability in terms of agreement between various static and dynamical simulated and experimental properties: liquid densities, radial distribution functions (RDF), diffusion constants [1,2]. The most promising results for the description of hydrogen bonding have been achieved using the GAFF force field: density is determined with 3.15% relative error and the covalent bond peak (1.00Å) and hydrogen bond peak (1.58Å) in $g_{HX(X=O, P)}(r)$ close to experimental values of 0.98Å and 1.58Å, respectively. SPC/EF model is used for modeling water. A number of systems having 1:0, 3:1, 1:1, 1:3, 1:6 and 0:1 H3PO4-H2O ratios have been constructed and compared in terms of RDF, diffusivity and its activation energy. The RDF analysis shows that first and second solvation shells in 3:1, 1:1 and 1:3 H3PO4-H2O mixtures are almost identical to those of pure H3PO4 rather than H2O, with maximums slowly but gradually shifting from those of pure H3PO4 towards the ones' accordingly seen in H2O. A special analysis tool written in C++ for the hydrogen bonding in terms of bond lengths and angles for different hydrogen bonded species in these systems is developed and used for a more sophisticated analysis.
