Regarding change in spin densities derived from CDFT guess

Hello all,
I have performed an optimization using CDFT with desired spin densities (input attached). But when I am removing the cdft from input meaning doing a single point without cdft but using the guess from cdft+opt job (SCF_GUESS = read), spin densities changes back to what I had before doing CDFT OPT.
However, doing a SCF_GUESS = read and skip_scfman = true will keep the spin densities same but this will not serve the purpose as I need to get the converged SCF and I need SCF energy.
I also know that doing a single point while keeping the CDFT and guess=read will give the desired spin densities because the electronic structure is constrained.
Now, Is there a way to do a single point calculation while solving the SCF and using the electronic structure from CDFT OPT (SCF_GUESS = read) while keeping the spin densities unchanged? This calculation should not be a CDFT.

Input for CDFT OPT is as follows:

$comment
Fe(III) optimization with cdft 
$end

$molecule
0 4
Fe    0.000000    0.000000    0.000000
 O   -0.439492    0.039527    1.765463
 H   -1.410810    0.105902    1.798017
 C   -3.306617   -0.738118    3.302109
 C   -3.499722    0.418966    2.390203
 C   -3.159447    1.802757    2.809928
 H   -3.271417   -1.686101    2.754503
 H   -4.182390    0.293200    1.550596
 H   -3.152625    2.495463    1.962082
 H   -2.171958    1.846551    3.286249
 H   -3.876924    2.202427    3.545776
 H   -4.122510   -0.828471    4.038171
 H   -2.377258   -0.647756    3.876194
 C    0.466130   -1.698851   -2.777456
 S    0.500943    0.021444   -2.234090
 H    1.219399   -2.295820   -2.254824
 H    0.678462   -1.718608   -3.848982
 H   -0.518663   -2.144885   -2.607145
 N   -1.589303   -1.185995   -0.370126
 N    1.144834   -1.608840    0.295720
 N   -1.159449    1.618758   -0.347101
 C   -3.197933    0.502898   -1.117279
 C   -0.585073   -3.342400    0.208058
 C    3.312949   -0.506656    0.569151
 C    0.440530    3.364282    0.268570
 C   -2.811279   -0.804766   -0.858236
 C    0.718550   -2.907495    0.393725
 C    2.866221    0.805149    0.543328
 C   -0.808752    2.925107   -0.146762
 C   -3.666423   -1.950424   -0.977562
 C    1.837954   -3.768682    0.649795
 C    3.691935    1.956718    0.778585
 C   -1.901848    3.784795   -0.504183
 C   -2.952251   -3.025303   -0.541354
 C    2.944280   -2.977172    0.692841
 C    2.872270    3.041147    0.747620
 C   -2.906675    2.981432   -0.949595
 C   -1.649963   -2.541487   -0.182049
 C    2.498884   -1.627500    0.491783
 C    1.549419    2.555986    0.466164
 C   -2.437917    1.628942   -0.839023
 N    1.568684    1.193141    0.352294
 H   -4.201541    0.659758   -1.501468
 H   -0.774534   -4.406361    0.314032
 H    0.578260    4.435020    0.384964
 H    4.373886   -0.670619    0.731564
 H    1.770333   -4.842158    0.770704
 H    3.973469   -3.263726    0.866198
 H   -3.261038   -4.060891   -0.479261
 H   -4.686714   -1.919255   -1.337833
 H   -1.881876    4.864603   -0.432350
 H   -3.887064    3.263196   -1.311614
 H    3.125537    4.084166    0.886914
 H    4.758598    1.920216    0.958005
$end

$rem
Jobtype			opt
EXCHANGE		M06
Basis			cc-pVDZ-full
MAX_SCF_CYCLES		1000
GEOM_OPT_MAX_CYCLES	300
Symmetry		false
SYM_IGNORE		true
cdft			true
cdft_postdiis		true
cdft_prediis		true
cdft_thresh		5
cdft_becke_pop		true
$end

$cdft
1.0
1 1 1 s
0.0
1 2 3 s
1.0
1 4 13 s
1.0
1 14 54 s
$end

I am not sure that I understand, but I think what you are describing is that the electronic structure collapses to a different state (with different spin density) when normal, unconstrained SCF is used. This is expected because CDFT implements a constraint in order to prevent you from obtaining the true ground state. Fully-relaxed solution need not (and generally does not) satisfy the constraint(s) so electronic structure will change.

The CDFT calculation should give you an energy on the final optimization cycle, what is wrong with that energy?