We present here a set of scalar-relativistic norm-conserving 4f-in-core pseudopotentials, together with complementary valence-shell Gaussian basis sets, for the lanthanide (Ln) series (Ce-Lu). The Goedecker, Teter, and Hutter (GTH) formalism is adopted with the generalized gradient approximation (GGA) for the exchange-correlation Perdew-Burke-Ernzerhof (PBE) functional. The 4f-in-core pseudopotentials are built through attributing 4f-subconfiguration 4fn (n = 1-14) for Ln (Ln = Ce-Lu) into the atomic core region, making it possible to circumvent the difficulty of the description of the open 4fn valence shell. A wide variety of computational benchmarks and tests have been carried out on lanthanide systems including Ln3+-containing molecular complexes, aqueous solutions, and bulk solids to validate the accuracy, reliability, and efficiency of the optimized 4f-in-core GTH pseudopotentials and basis sets. The 4f-in-core GTH pseudopotentials successfully replicate the main features of lanthanide structural chemistry and reaction energetics, particularly for nonredox reactions. The chemical bonding features and solvation shells, hydrolysis energetics, acidity constants, and solid-state properties of selected lanthanide systems are also discussed in detail by utilizing these new 4f-in-core GTH pseudopotentials. This work bridges the idea of keeping highly localized 4f electrons in the atomic core and efficient pseudopotential formalism of GTH, thus providing a highly efficient approach for studying lanthanide chemistry in multi-scale modeling of constituent-wise and structurally complicated systems, including electronic structures of the condensed phase and first-principles molecular dynamics simulations.