Geometry optimization and the electronic structure calculations of Pu^Z+ complexes (Z = 3-6) in water solution have been performed, within the framework of the DMol³ and Relativistic Discrete-Variational (RDV) methods. For the simulation of Pu^Z+ molecular environment in aqueous solution we used 22 and 32 water molecules randomly distributed around cation. To model the effect of bulk solvent environment we used COSMO (Conductor-like Screening Model) potential for water (ε = 78.54). The obtained results showed that this approach allows the modeling of water dissociation and the formation of hydrolysis products. Our previously suggested scheme for the calculation of interaction energies between selected fragments of multi-molecular systems provides the quantitative estimation of the interaction strengths between plutonium in various oxidation states and each ligand in the first and second coordination shells in water solution.