I personally changed the spin multiplicity for Fe2+ to 5, since I am assuming a 3d6 configuration for my system.

You are right, the Fe2+ should have spin multiplicity of 5. I will change "chgmulti(26,2)=7" to "chgmulti(26,2)=5" in the next release, though this change will not evidently affect the result.

As you suggested, I will add a warning when unusual atomic charged state is involving.

A few comments that might help other users:

*) .wfn files generated wtih ORCA
MultiWFN suggests not to use .wfn file generated with ORCA due to non-standardized format. This problem can be circumvented by generating a molden input file from ORCA and convert this with MultiWFN to .wfn, since MultiWFN can handle the non-standardized ORCA molden input file.

*) Charges higher than +-3
In order to include more extreme charges (+- 4) the population.f90 has to be changed in a couple of places within the genatmradfile subroutine. It worked for -4, but not for -5. (e.g. for metal carbides this gives Error: atmrad/C@.rad was not prepared! for charges < -4 on carbon). However, the iterative Hirshfeld is not applicable to these problems anyways, because the free anions – and thus proatoms – are not stable. Maybe it would help to include a warning about the limit of the method for non-organic molecules?

*) Transition metals
It rightly states in the population.f90 file that the spin multiplicity for transition metal ions are determined based on chemical intuition. Users should check this in particular the metals in their system. I personally changed the spin multiplicity for Fe2+ to 5, since I am assuming a 3d6 configuration for my system. But this only matters if you let Gaussian generate the densities.

Requirement of possible charged states can be found in "subroutine genatmradfile" in population.f90 of Multiwfn source code package. For example, in this subroutine you can find

!C,Si,Ge,Sn,Pb
chgmulti(6,0)=3
chgmulti(6,1)=2
chgmulti(6,2)=1
chgmulti(6,-1)=4
chgmulti(6,-2)=3
chgmulti(14,:)=chgmulti(6,:)
chgmulti(32,:)=chgmulti(6,:)
chgmulti(50,:)=chgmulti(6,:)
chgmulti(82,:)=chgmulti(6,:)

That means, for carbon, you should provide "C-2.wfn","C-1.wfn","C_0.wfn","C+1.wfn","C+2.wfn" in the "atmrad" folder, and the charge and spin multiplicity for calculating these atomic .wfn files should be (-2,3), (-1,4), (0,3), (1,2), (2,1), respectively. Similary, for Si, Ge, Sn, Pb, if they are involved in the current system you should also provide corresponding .wfn files.

If any needed atomic .wfn file cannot be found in the "atmrad" folder, Multiwfn will try to invoke Gaussian to calculate them; while if all needed .wfn files can be directly found, then Multiwfn will directly start H-I calculation without invoking Gaussian.

BTW: Please do not use Multiwfn 3.6 for H-I calculation, there is a known severe bug for large system. Please update to Multiwfn 3.7(dev).

I am wondering whether the iterative Hirshfeld atomic charges can be generated without access to Gaussian. Is it somehow possible to generate the proatom densities prior to the calculation?

I am using MutiWFN v3.6.

Thanks for the help!

Best,
nico