We evaluate nucleon four-point functions in the framework of lattice QCD in order to extract the first Mellin moment of double parton distributions (DPDs) in the unpolarized proton. In this first study, we employ an n(f) = 2+1 ensemble with pseudoscalar masses of m(pi) = 355 MeV and m(K) = 441 MeV. The results are converted to the scale mu = 2 GeV. Our calculation includes all Wick contractions, and for almost all of them a good statistical signal is obtained. We analyze the dependence of the DPD Mellin moments on the quark flavor and the quark polarization. Furthermore, the validity of frequently used factorization assumptions is investigated.