Thank you for the comment. 

The model I use for Phaser is coming from Sculptor, so it does not have heteroatoms.

I introduce heteroatom (Fe) only after Phaser. I use Phaser-EP to find it, eventhough it is obvious from the density. Then I merge Fe and Phaser.pdb in Coot after which I feed it to ReadySet. So the pdb that I give to Refine either have Fe with occupancy of 1, <1 or >1 (I tried both). I use PDB editor to add charge, but the outcome is still the same. After refinement Fe occupancy is either 1 or very close.