Concentrated solutions of short blunt-ended DNA duplexes, as short as 6 base pairs, are known to order into the nematic liquid crystal phases. This self-assembly is due to the stacking interactions between duplex terminals that promotes their aggregation into polydisperse chains with a significant persistence length. Experiments show that liquid crystal phases form above a critical volume fraction depending on the duplex length. We introduce and investigate via numerical simulations, a coarse-grained model of DNA double-helical duplexes. Each duplex is represented as an hard quasi-cylinder whose bases are decorated with two identical reactive sites. The stacking interaction between terminal sites is modeled via a short-range square-well potential. We compare the numerical results with predictions based on a free energy functional and find satisfactory quantitative matching of the isotropic-nematic phase boundary and of the system structure. Comparison of numerical and theoretical results with experimental findings confirm that the DNA duplex self-assembly can be properly modeled via equilibrium polymerization of cylindrical particles. This insight enables us to estimate the stacking energy. © 2011 American Chemical Society.