I tried SCF convergence algorithms like SCF=QC, but got the same result. Although, you are right that I need the unstable wavefunction, which is evident from the fact that when I tried the calculation using scf=(qc,conver=2), it gave me a spin pop at Fe ~ 1.3 which is close to the Fe(III), but that is not much use to me.

Thanks

By the way, you may also try to use different SCF convergence algorithms, such as SCF=QC. According to my experience, SCF=QC has a higher probability than default algorithm to converge to a unstable wavefunction, perhaps the "unstable" wavefunction is what you actually need.

2.)  I prepared a Gaussian input using cdft.fchk, with wavefunction written in the input using Multiwfn.

In both cases guess=(read,only,save) gave the Fe(III) (as expected)  but everytime  it went to Fe(IV) during SCF. So, the problem still persists.

I have a blog article about CDFT: http://www.shanxitv.org/271 (written in Chinese)

I am trying to obtain Fe(III) electromer of iron-porphyrin system published in this paper:

Reference:  J. Phys. Chem. Lett. 2011, 2, 2229–2235; DOI:  10.1021/jz2007534

I am using Gaussian16.

Figure 2a, I am able to obtain 4I(IV) state, successfully and the spin densities match with the reported values.  But I could not obtain the Fe(III) state 4I(III) in gas phase, given in Figure 2 from the above-mentioned paper.

The problem is that while solving the SCF, every time it falls to Fe(IV) state.

Is there a way to generate the desired guess for Fe(III) state from Fe(IV) fchk/chk files using Multiwfn?  That guess can then be used in gaussian calculation.

Thanks in advance.