We construct a model for double parton distribution functions (dPDFs) based on the notion of self-similarity, pursued earlier for small x physics at HERA. The most general form of dPDFs contains total thirteen parameters to be fitted from data of proton-proton collision at LHC. It is shown that the constructed dPDF does not factorize into two single PDFs in conformity with QCD expectation, and it satisfies the condition that at the kinematic boundary x 1 + x 2 = 1 (where x 1 and x 2 are the longitudinal fractional momenta of two partons), the dPDF vanishes.