A collection of 24 hexacoordinate Co(II) complexes was investigated by ab initio CASSCF + NEVPT2 + SOC calculations. In addition to the energies of spin–orbit multiplets (Kramers doublets, KD) their composition of the spins is also analyzed, along with the projection norm to the effective Hamiltonian. The latter served as the evaluation of the axial and rhombic zero-field splitting parameters and the g-tensor components. The fulfilment of spin-Hamiltonian (SH) formalism was assessed by critical indicators: the projection norm for the first Kramers doublet N(KD1) > 0.7, the lowest g-tensor component g1 > 1.9, the composition of KDs from the spin states |±1/2> and |±3/2> with the dominating percentage p > 70%, and the first transition energy at the NEVPT2 level 4Δ1. Just the latter quantity causes a possible divergence of the second-order perturbation theory and a failure of the spin Hamiltonian. The data set was enriched by the structural axiality Dstr and rhombicity Estr, respectively, evaluated from the metal–ligand distances Co-O, Co-N and Co-Cl corrected to the mean values. The magnetic data (temperature dependence of the molar magnetic susceptibility, and the field dependence of the magnetization per formula unit) were fitted simultaneously, either to the Griffith–Figgis model working with 12 spin–orbit kets, or the SH-zero field splitting model that utilizes only four (fictitious) spin functions. The calculated data were analyzed using statistical methods such as Cluster Analysis and the Principal Component Analysis.