Variational methods are used for targeting specific correlation effects by tailoring the configuration space. Independent sets of correlation orbitals, embedded in partitioned correlation functions (PCFs), are produced from multiconfiguration Hartree-Fock (MCHF) and DiracHartree-Fock (MCDHF) calculations. These non-orthogonal functions span configuration state function (CSF) spaces that are coupled to each other by solving the associated generalized eigenvalue problem. The Hamiltonian and overlap matrix elements are evaluated using the biorthonormal orbital transformations and efficient counter-transformations of the configuration interaction eigenvectors [1]. This method was successfully applied for describing the total energy of the ground state of beryllium [2]. Using this approach, we demonstrated the fast energy convergence in comparison with the conventional SD-MCHF method optimizing a single set of orthonormal one-electron orbitals for the complete configuration space. In the present work, we investigate the Partitioned Correlation Function Interaction (PCFI) approach for the two lowest states of neutral lithium, i.e. 1s 2 2s 2 S and 1s 2 2p 2 P o . For both states, we evaluate the total energy, as well as the expectation values of the specific mass shift operator, the hyperfine structure parameters and the transition probabilities using different models for tailoring the configuration space. We quantify the “constraint effect” due to the use of fixed PCF eigenvector compositions and illustrate the possibility of a progressive deconstraint, up to the non-orthogonal configuration interaction limit case. The PCFI estimation of the position of the quartet system relative to the ground state of B I will also be presented. The PCFI method leads to an impressive improvement in the convergence pattern of all the spectroscopic properties. As such, Li I, Be I and B I constitute perfect benchmarks for the PCFI method. For larger systems, it becomes hopeless to saturate a single common set of orthonormal orbitals and the PCFI method is a promising approach for getting high quality correlated wave functions. The present study constitutes a major step in the current developments of both atsp2K and grasp2K packages that adopt the biorthonormal treatment for estimating energies, isotope shifts, hyperfine structures and transition probabilities.