Metal (4f)−ligand (Cl 3p) bonding in LnCl63- (Ln = Ce to Yb) complexes has been studied on the basis of 4f→4f and Cl,3p→4f charge-transfer spectra and on the analysis of these spectra within the valence bond configuration interaction model to show that mixing of Cl 3p into the Ln 4f ligand field orbitals does not exceed 1%. Contrary to this, Kohn−Sham formalism of density functional theory using currently available approximations to the exchange-correlation functional tends to strongly overestimate 4f−3p covalency, yielding, for YbCl63-, a much larger mixing of Cl 3p→4f charge transfer into the f13 ionic ground-state wave function. Thus, ligand field density functional theory, which was recently developed and applied with success to complexes of 3d metals in our group, yields anomalously large ligand field splittings for Ln, the discrepancy with experiment increasing from left to the right of the Ln 4f series. It is shown that eliminating artificial ligand-to-metal charge transfer in Kohn−Sham calculations by a procedure described in this work leads to energies of 4f−4f transitions in good agreement with experiment. We recall an earlier concept of Ballhausen and Dahl which describes ligand field in terms of a pseudopotential and give a thorough analysis of the contributions to the ligand field from electrostatics (crystal field) and exchange (Pauli) repulsion. The close relation of the present results with those obtained using the first-principles based and electron density dependent effective embedding potential is pointed out along with implications for applications to other systems.