For the first time, the Fick diffusion coefficient matrix of a quaternary liquid mixture is sampled consistently by means of molecular dynamics simulation. The required phenomenological diffusion coefficient and thermodynamic factor matrices the mixture water + methanol + ethanol + 2-propanol are determined following the Green-Kubo formalism and Kirkwood-Buff theory. Further, a system size correction methodology for the Fick diffusion coefficient of multicomponent mixtures is proposed. Ten compositions are studied under ambient conditions and validated by analyzing the ternary limits of the quaternary Fick diffusion matrix. Because of complex intermolecular interactions due to the presence of hydrogen bonding, the elements of the Fick diffusion coefficient matrix exhibit a significant composition dependence. The magnitude of several cross coefficients indicates important coupling effects mainly affecting the diffusive flux of water. These effects are explained in the light of the structural information given by the radial distribution functions of the mixture. This work that solely rests on molecular dynamics simulation techniques to predict the Fick diffusion coefficient matrix of quaternary mixtures is expected to be a significant step forward for the understanding of multicomponent diffusion.