Firstly, do the geometry optimization for the doped system and calculate the charge density of the optimized configuration, i.e., \rho (dopant+host);
Secondly, calculate the charge density of a configuration without dopant, meanwhile, the atomic positions in this configuration are same to those in the first step. The obtained charge density is denoted as \rho (host)
Thirdly, calculate the charge density of a configuration containing only the dopant, here the atomic positions of dopant are same to those in the first step. The obtained charge density is denoted as \rho (dopant).
Finally, the charge difference can be estimated as \rho (dopant+host) - \rho (host) - \rho (dopant).