Phosphorylation is a major post-translational modification, which is involved in the regulation of the dynamics and function of Intrinsically Disordered Proteins (IDPs). We recently characterized a phenomenon, which we termed n-Phosphate collaborations (nP-collabs), where bulk cations form stable bridges between several phosphoresidues in all-atom molecular dynamic simulations. nP-collabs were found to be sensitive to the combination of forcefields and cation types. Here, we attempt to assess the physical relevance of these nP-collabs by evaluating the strength of the cation/phosphate interaction through osmotic coefficient ( {varphi}) calculations on the model 2Na+HPO42- and 2K+HPO42- salts, using different classical forcefields for phosphorylations. All forcefields were found to overestimate the strength of the interaction to various degrees. We thus designed new parameters for CHARMM36m and AmberFF99SB-ILDN using the Electronic Continuum Correction (ECC) approach, which provide remarkable agreement for {varphi} values for both cation types and over a range of concentrations. We provide a preliminary test of these ECC parameters for phosphorylations by simulating the sevenfold-phosphorylated rhodopsin peptide 7PP and comparing secondary chemical shifts to experimental data. Conformational ensembles resulting from the ECC-derived phosphorylated forcefields display both qualitative and quantitative improvements with regards to full-charge forcefields. We thus conclude that long-lasting nP-collabs are artifacts for classical forcefields born from the lack of explicit polarization, and propose a possible computational strategy for the extensive parameterization of phosphorylations. The presence of long-lived nP-collabs in simulations produced using classical forcefields is therefore a serious concern for the accurate modelling of multiphosphorylated peptides and IDPs, which are at the center of research questions regarding neurodegenerative diseases such as Alzheimer\'s or Parkinson\'s.