In this work, we develop a DNN prediction model for the packing density of small organic molecules in an amorphous bulk phase and conduct a hyperscreening of 1.5 million candidate compounds. Our interest in this target property originates from our ongoing in silico discovery and design efforts for polymers with high refractive index (RI) to be used in optic and optoelec-tronic applications . We previously established an RI modeling protocol based on the Lorentz-Lorenz equation and parametrized with the polarizability and number density . For the number density, we introduced a hybrid physics-based/data-derived prediction model using the van der Waals volume computed via the Slonimskii method and the packing coefficient from a support vector regression machine learning model . An alternative and commonly employed route to computing the (number) density, is the use of molecular dynamics (MD) simulations, which we recently started exploring in our study of high-RI polymers. However, as these MD calculations are computationally expensive and technically challenging, they are not particularly well suited for the large-scale assessment of compounds in the course of high-throughput screening studies. To bypass this problem, we develop a DNN surrogate model for the MD density predictions. It allows us to rapidly and accurately compute the density values of the 1.5 million molecules of a virtual screening library we create for proof of concept. For this, we perform MD simulations on a subset of 100,000 compounds, use the results to train our DNN model, and subsequently employ it to compute the packing density of the remaining 1.4 million molecules. We mine the density results to identify patterns that lead to desirable outcomes (i.e., different density regimes). We also evaluate the learning curve for the density prediction to assess the dependence of training set size and model accuracy.