Inconsistency in softening parameter for soft_coulomb species

The documentation states that the softening parameter a of the potential v(r)=−zval​/(a2+r2) is defined in the input file, but actually in the code this parameter is taken as a**2. In line 1032 of species_pot.F90: vl(ip) = vl(ip) -species%get_zval()/sqrt(r2+species%get_softening())

To fix it, the documentation should stay the same but in the code the value species%get_softening() should be squared.